69{
70 MsgStream log(
msgSvc(), name());
71 log << MSG::INFO << "in execute()" << endreq;
72
73 int numOfChargedTrks = 0;
74 int numOfKalTrks = 0;
75 int numOfDedx = 0;
76 int numOfExt = 0;
77 int numOfTof = 0;
78 int numOfShower = 0;
79 int numOfMatchedShower = 0;
80 int numOfMucTrks = 0;
81
82
83 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
84 if (!eventHeader)
85 {
86 log << MSG::FATAL << "Could not find Event Header" << endreq;
87 return( StatusCode::FAILURE);
88 }
89 log << MSG::INFO << ": retrieved event: " << eventHeader->eventNumber() << " run: " << eventHeader->runNumber() << endreq;
90
91
93 RecMdcTrackCol::iterator beginOfMdcTrkCol = RecMdcTrackCol::iterator(NULL);
94 if(!recMdcTrackCol)
95 {
96 log<<MSG::INFO<<"Could not find RecMdcTrackCol!"<<endreq;
97 }
98 else
99 {
100 beginOfMdcTrkCol = recMdcTrackCol->begin();
101 numOfChargedTrks = recMdcTrackCol->size();
102 }
103
105 RecMdcKalTrackCol::iterator beginOfMdcKalTrackCol = RecMdcKalTrackCol::iterator(NULL);
106 if(!mdcKalTrackCol)
107 {
108 log<<MSG::INFO<<"Could not find MdcKalTrackCol!"<<endreq;
109 }
110 else
111 {
112 numOfKalTrks = mdcKalTrackCol->size();
113 beginOfMdcKalTrackCol = mdcKalTrackCol->begin();
114 }
115
116
118 RecMdcDedxCol::iterator beginOfDedxCol = RecMdcDedxCol::iterator(NULL);
119 if(!mdcDedxCol)
120 {
121 log<<MSG::INFO<<"Could not find MdcDedxCol!"<<endreq;
122 }
123 else {
124 numOfDedx = mdcDedxCol->size();
125 beginOfDedxCol = mdcDedxCol->begin();
126 }
127
128
130 RecExtTrackCol::iterator beginOfExtTrackCol = RecExtTrackCol::iterator(NULL);
131 if(!extTrackCol)
132 {
133 log<<MSG::INFO<<"Could not find RecExtTrackCol!"<<endreq;
134 }
135 else
136 {
137 beginOfExtTrackCol = extTrackCol->begin();
138 numOfExt = extTrackCol->size();
139 }
140
142 RecTofTrackCol::iterator beginOfTofTrackCol = RecTofTrackCol::iterator(NULL);
143 if(!TofTrackCol)
144 {
145 log<<MSG::INFO<<"Could not find RecTofTrackCol!"<<endreq;
146 }
147 else
148 {
149 beginOfTofTrackCol = TofTrackCol->begin();
150 numOfTof = TofTrackCol->size();
151 }
152
154 RecEmcShowerCol::iterator beginOfEmcTrackCol = RecEmcShowerCol::iterator(NULL);
155 if(!emcRecShowerCol)
156 {
157 log<<MSG::INFO<<"Could not find RecEmcShowerCol!"<<endreq;
158 }
159 else
160 {
161 beginOfEmcTrackCol = emcRecShowerCol->begin();
162 numOfShower = emcRecShowerCol->size();
163 }
164
166 RecMucTrackCol::iterator beginOfMucTrackCol = RecMucTrackCol::iterator(NULL);
167 if(!mucTrackCol)
168 {
169 log<<MSG::INFO<<"Could not find RecMucTrackCol!"<<endreq;
170 }
171 else
172 {
173 beginOfMucTrackCol = mucTrackCol->begin();
174 numOfMucTrks = mucTrackCol->size();
175 }
176
177
178
179 int MaxHit = numOfChargedTrks+numOfShower+numOfTof;
180
181 vector<bool> dEdxMatched;
182 vector<bool> kalTrkMatched;
183 vector<bool> extMatched;
184 vector<bool> TofMatched;
185 vector<bool> emcMatched;
186 vector<bool> mucMatched;
187
188 for(int i=0;i<MaxHit;i++)
189 {
190 dEdxMatched.push_back(false);
191 kalTrkMatched.push_back(false);
192 extMatched.push_back(false);
193 TofMatched.push_back(false);
194 emcMatched.push_back(false);
195 mucMatched.push_back(false);
196 }
197
198
199
200
202 if(numOfChargedTrks+numOfShower)
203 {
204 for(int i=0;i<numOfChargedTrks;i++)
205 {
207 int mdcTrkID = (*(beginOfMdcTrkCol+i))->trackId();
209 aNewEvtRecTrack->
setMdcTrack(*(beginOfMdcTrkCol+i));
210
211 for(int j=0;j<numOfKalTrks;j++)
212 if( !kalTrkMatched[j] && mdcTrkID==(*(beginOfMdcKalTrackCol+j))->trackId() )
213 {
215
216 kalTrkMatched[j]=true;
217 }
218
219 for(int j=0;j<numOfDedx;j++)
220 if( !dEdxMatched[j] && mdcTrkID==(*(beginOfDedxCol+j))->trackId() )
221 {
222 aNewEvtRecTrack->
setMdcDedx(*(beginOfDedxCol+j));
223
224 dEdxMatched[j]=true;
225 }
226
227 for(int j=0;j<numOfExt;j++)
228 if( !extMatched[j] && mdcTrkID==(*(beginOfExtTrackCol+j))->trackId() )
229 {
230 aNewEvtRecTrack->
setExtTrack(*(beginOfExtTrackCol+j));
231
232 extMatched[j]=true;
233 }
234 for(int j=0;j<numOfTof;j++)
235 if( !TofMatched[j] && mdcTrkID==(*(beginOfTofTrackCol+j))->trackID() )
236 {
237 aNewEvtRecTrack->
addTofTrack(*(beginOfTofTrackCol+j));
238
239 TofMatched[j]=true;
240 }
241 for(int j=0;j<numOfMucTrks;j++)
242 if( !mucMatched[j] && mdcTrkID==(*(beginOfMucTrackCol+j))->trackId() )
243 {
244 aNewEvtRecTrack->
setMucTrack(*(beginOfMucTrackCol+j));
245
246 mucMatched[j]=true;
247 }
248
249 if(m_Output==1){
250 myExted = -1.;
251 myMatched = -1.;
252 myEnergyMatched = -1.;
253 }
255 {
256 if(m_Output==1)myExted = 1.;
257
258 Hep3Vector extPos = (aNewEvtRecTrack->
extTrack())->emcPosition();
259 double extTheta = extPos.theta();
260 double extPhi = extPos.phi();
261
262 double dAngleMin = 9999.;
263 int indexOfEmcNearest = -1;
264
265 for(int j=0;j<numOfShower;j++)
266 {
267 if(!emcMatched[j])
268 {
269 HepPoint3D emcPot = (*(beginOfEmcTrackCol+j))->position();
270 Hep3Vector emcPos(emcPot.x(),emcPot.y(),emcPot.z());
271 double dAngle = extPos.angle(emcPos);
272 if(dAngle<dAngleMin)
273 {
274 dAngleMin = dAngle;
275 indexOfEmcNearest = j;
276 }
277 }
278 }
279 if(indexOfEmcNearest>=0)
280 {
281 HepPoint3D emcPot = (*(beginOfEmcTrackCol+indexOfEmcNearest))->position();
282 double theta = emcPot.theta();
283 double phi = emcPot.phi();
284 double deltaTheta = (extTheta-theta);
286 if(myActiveCutFlag)
287 {
288 double tanLamda = (*(beginOfMdcTrkCol+i))->helix(4);
289 double kappa = (*(beginOfMdcTrkCol+i))->helix(2);
290 double p = sqrt(1+tanLamda*tanLamda)/fabs(kappa);
291 double aThetaCut = thetaCut(p,myEmcThetaNSigmaCut,myParticleId);
292 double aPhiCUt = phiCut(p,myEmcPhiNSigmaCut,myParticleId);
293 if(fabs(deltaTheta)<aThetaCut && fabs(deltaPhi)<aPhiCUt)
294 {
295 aNewEvtRecTrack->
setEmcShower(*(beginOfEmcTrackCol+indexOfEmcNearest));
296
297 emcMatched[indexOfEmcNearest]=true;
298 numOfMatchedShower++;
299 if(m_Output==1)myMatched = 1.;
300 }
301 }
302 else if(fabs(deltaTheta)<myEmcThetaCut && fabs(deltaPhi)<myEmcPhiCut)
303 {
304 aNewEvtRecTrack->
setEmcShower(*(beginOfEmcTrackCol+indexOfEmcNearest));
305
306 emcMatched[indexOfEmcNearest]=true;
307 numOfMatchedShower++;
308 if(m_Output==1)myMatched = 1.;
309 }
310 }
311 }
312 aNewEvtRecTrackCol->push_back(aNewEvtRecTrack);
313 if(m_Output==1){
314 myEnergyMatched = 0.;
315 myNTuple1->write();
316 }
317 }
318
319 int id=0;
320 for(int i=0;i<numOfShower;i++)
321 {
322 if(!emcMatched[i])
323 {
325 aNewEvtRecTrack->
setTrackId(
id+numOfChargedTrks);
327
328
329
330 HepPoint3D emcPos((*(beginOfEmcTrackCol+i))->position());
331 double emcPhi = emcPos.phi();
332 emcPhi = emcPhi < 0 ? emcPhi+CLHEP::twopi : emcPhi;
333
334 int module=9999, nphi1=9999, nphi2=9999, nphi1_mrpc=9999;
335 module = (*(beginOfEmcTrackCol+i))->module();
336 if(module==1) {
337 nphi1 = (int)(emcPhi*88./CLHEP::twopi);
338 nphi2 = (int)(emcPhi*88./CLHEP::twopi+0.5);
339 nphi2 += 88;
340 } else {
341 nphi1 = (int)(emcPhi*48./CLHEP::twopi+0.5);
342
343
344
345 double fi_endtof_mrpc = atan2(emcPos.y(),emcPos.x())-3.141592654/2.;
346 if(fi_endtof_mrpc<0) fi_endtof_mrpc += 2*3.141592654;
347 nphi1_mrpc=((int)(fi_endtof_mrpc/(3.141593/18)))+1;
348
349 if(nphi1_mrpc%2==0){nphi1_mrpc=nphi1_mrpc/2;}
350 else {nphi1_mrpc=(nphi1_mrpc+1)/2;}
351
352 if(emcPos.z()>0) {(nphi1_mrpc -= 19); nphi1_mrpc=nphi1_mrpc*(-1);}
353
354
355 }
356
357 int dPhiMin = 9999;
358 int partid_dPhiMin = 7;
359 bool mrpcmatch = false;
360
361 int indexOfTofNearest = -1;
362 for(int j=0;j<numOfTof;j++) {
363 if( !TofMatched[j] ) {
364
365 Identifier id((*(beginOfTofTrackCol+j))->tofID());
367
368 if(!is_mrpc)
369 {
373 im += layer * 88;
374
375 if(module!=barrel_ec) continue;
376
377 int dPhi;
378 if(layer == 0) {
379 dPhi =
abs(im - nphi1);
380 } else {
381 dPhi =
abs(im - nphi2);
382 }
383 if(module==1) {
384 dPhi = dPhi>=44 ? 88-dPhi : dPhi;
385 } else {
386 dPhi = dPhi>=24 ? 48-dPhi : dPhi;
387 }
388
389 if(dPhi < dPhiMin) {
390 dPhiMin = dPhi;
391 indexOfTofNearest = j;
392 }
393 }
394 else
395 {
398 int module_mrpc = (int)im/25;
399
400 if(module==1) continue;
401 if(module==0 && (barrel_ec==5 || barrel_ec==6)) continue;
402 if(module==2 && (barrel_ec==3 || barrel_ec==4)) continue;
403
404 int dphi;
405 dphi =
abs(module_mrpc-nphi1_mrpc);
406
407 if(dphi < dPhiMin)
408 {
409 dPhiMin=dphi;
410 partid_dPhiMin=barrel_ec;
411 indexOfTofNearest = j;
412 mrpcmatch = true;
413 }
414 }
415
416 }
417 }
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433 if(indexOfTofNearest>=0) {
434 HepPoint3D tofPos(0,0.867,(*(beginOfTofTrackCol+indexOfTofNearest))->zrhit());
435 double matWindow = 0;
436 double eShower = (*(beginOfEmcTrackCol+i))->e5x5();
437 if(eShower>0) {
438 matWindow = 0.01277+0.004268/sqrt(eShower);
439 }
440 }
441
442
443 if(indexOfTofNearest>=0 && dPhiMin<=1) {
444
445 double eShower = (*(beginOfEmcTrackCol+i))->e5x5();
446 double matWindow = 0;
447 if(eShower>0) {
448 matWindow = 0.01277+0.004268/sqrt(eShower);
449 }
450 matWindow = 0.2;
451
452 HepPoint3D tofPos(0,0.867,(*(beginOfTofTrackCol+indexOfTofNearest))->zrhit());
453 if(fabs(tofPos.theta()-emcPos.theta())<matWindow || module!=1) {
454 aNewEvtRecTrack->
addTofTrack(*(beginOfTofTrackCol+indexOfTofNearest));
455 TofMatched[indexOfTofNearest]=true;
456
457
458 }
459 }
460
461 aNewEvtRecTrackCol->push_back(aNewEvtRecTrack);
462 id++;
463 }
464 }
465 }
466
467
468 if(myMsgFlag) cout<<"Charged tracks: "<<numOfChargedTrks<<" Ext track: "<<numOfExt <<" Shower:"<<numOfShower<<" Matched shower:"<<numOfMatchedShower<<endl;
469
470
471 DataObject *aDataObject1;
473 if ( aDataObject1 == NULL ) {
476 if ( sc.isFailure() ) {
477 log << MSG::FATAL << "Could not register EvtRecObject in TDS!" << endreq;
478 return StatusCode::FAILURE;
479 }
480 }
481
482 DataObject* aDataObject2;
484 if ( aDataObject2 == NULL ) {
486 aRecEvent->
setTotalTracks(numOfChargedTrks+numOfShower-numOfMatchedShower);
490 if ( sc.isFailure() ) {
491 log << MSG::FATAL << "Could not register EvtRecEvent in TDS!" << endreq;
492 return( StatusCode::FAILURE);
493 }
494 }
495 else {
497 aRecEvent->
setTotalTracks(numOfChargedTrks+numOfShower-numOfMatchedShower);
500 }
501
502 DataObject* aDataObject3;
504 if ( aDataObject3 != NULL ) {
505
506 SmartIF<IDataManagerSvc> dataMgrSvc(eventSvc());
509 }
510
512 if ( sc.isFailure() ) {
513 log << MSG::FATAL << "Could not register EvtRecTrackCol in TDS!" << endreq;
514 return( StatusCode::FAILURE);
515 }
516
517
519 if(!evtRecTrackCol)
520 {
521 log<<MSG::FATAL<<"Could not find RecTrackListCol in TDS!"<<endreq;
522 }
523
524 EvtRecTrackCol::iterator itOfRecTrkListCol = evtRecTrackCol->begin();
525 if(myMsgFlag)
526 for(int i=0;itOfRecTrkListCol!=evtRecTrackCol->end();itOfRecTrkListCol++,i++)
527 {
528 cout<<"******************** Track "<<i<<": *********************"<<endl;
529 cout<<"Track ID: "<<(*itOfRecTrkListCol)->trackId()<<endl;
530 if((*itOfRecTrkListCol)->isMdcTrackValid())
531 {
532 RecMdcTrack* mdcTrk = (*itOfRecTrkListCol)->mdcTrack();
533 double kappa = mdcTrk->
helix(2);
534 double tanl = mdcTrk->
helix(4);
535 cout<<"Mdc kappa, tanl: "<<kappa<<", "<<tanl<<endl;
536 }
537 if((*itOfRecTrkListCol)->isExtTrackValid())
538 {
539 RecExtTrack* extTrack = (*itOfRecTrkListCol)->extTrack();
541 cout<<"Ext Emc pos:"<<emcPos.x()<<","<<emcPos.y()<<","<<emcPos.z()<<endl;
542 }
543 if((*itOfRecTrkListCol)->isEmcShowerValid())
544 {
545 RecEmcShower* emcTrack = (*itOfRecTrkListCol)->emcShower();
547 double x = emcPos.x();
548 double y = emcPos.y();
549 double z = emcPos.z();
550 cout<<
"Emc rec pos:"<<
x<<
","<<y<<
","<<z<<endl;
551 }
552 }
553
554 return StatusCode::SUCCESS;
555}
double abs(const EvtComplex &c)
ObjectVector< EvtRecTrack > EvtRecTrackCol
HepPoint3D position() const
const Hep3Vector emcPosition() const
const HepVector helix() const
......
void setTotalTracks(const int tottks)
void setTotalNeutral(const int nneu)
void setTotalCharged(const int nchrg)
void setMucTrack(const RecMucTrack *trk)
void setMdcTrack(const RecMdcTrack *trk)
void setTrackId(const int trkId)
void setMdcKalTrack(const RecMdcKalTrack *trk)
void setMdcDedx(const RecMdcDedx *trk)
void setEmcShower(const RecEmcShower *shower)
void addTofTrack(const SmartRef< RecTofTrack > trk)
void setExtTrack(const RecExtTrack *trk)
static int phi_module(const Identifier &id)
static int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
static bool is_mymrpc(const Identifier &id)
static int layer(const Identifier &id)
_EXTERN_ std::string Event
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecTrackCol
_EXTERN_ std::string RecExtTrackCol
_EXTERN_ std::string RecMdcDedxCol
_EXTERN_ std::string RecTofTrackCol
_EXTERN_ std::string RecMdcTrackCol
_EXTERN_ std::string RecMdcKalTrackCol
_EXTERN_ std::string RecMucTrackCol
_EXTERN_ std::string RecEmcShowerCol