96{
97
98 MsgStream log(
msgSvc(), name());
99 log << MSG::INFO << "execute()" << endreq;
100 int eventNumber, runNumber;
101
102 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
103 if (!eventHeader)
104 {
105 log << MSG::FATAL << "Could not find Event Header" << endreq;
106 return( StatusCode::FAILURE);
107 }
108 runNumber = eventHeader->runNumber();
109 eventNumber = eventHeader->eventNumber();
110 if(runNumber==13219&&eventNumber==3031231) return( StatusCode::SUCCESS);
111
112 if(msgFlag)
113 {
114 cout<<"TrackExt: ******************* Start a event *******************"<<endl;
115 cout<<"run= "<<runNumber<<"; event= "<<eventNumber<<endl;
116 }
117
118
119 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
120 if(!mucDigiCol) {
121 log << MSG::FATAL << "Could not find MUC digi" << endreq;
122 return( StatusCode::SUCCESS);
123 }
124
125
128
129 bool setTrk=false;
130
131 int parID;
132 if(myParticleName=="e") parID=0;
133 else if(myParticleName=="mu") parID=1;
134 else if(myParticleName=="pi") parID=2;
135 else if(myParticleName=="kaon") parID=3;
136 else if(myParticleName=="proton"||myParticleName=="anti_proton") parID=4;
137
138 if(myInputTrk=="Mdc")
139 {
140 SmartDataPtr<RecMdcTrackCol> aMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
141 if(!aMdcTrackCol)
142 {
143 log << MSG::WARNING << "Can't find RecMdcTrackCol in TDS!" << endreq;
144 return( StatusCode::SUCCESS);
145 }
147
148 }
149 else if(myInputTrk=="Kal")
150 {
151 SmartDataPtr<RecMdcKalTrackCol> aMdcKalTrackCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
152 if(!aMdcKalTrackCol)
153 {
154 log << MSG::WARNING << "Can't find RecMdcKalTrackCol in TDS!" << endreq;
155 return( StatusCode::SUCCESS);
156 }
158 }
159 else
160 {
161 log << MSG::WARNING << "Wrong type of inputTrk:" << myInputTrk << endreq;
162 return( StatusCode::SUCCESS);
163 }
164
166 if(setTrk)
167 {
169 {
170
172
173 for(int i=0; i<5; i++)
174 {
176 {
178
185 double tofInMdc = aExtMdcTrack.
GetTrkTof();
186
187 if(msgFlag)
188 {
189 cout<<"Start From:"<<position.x()<<' '<<position.y()<<' '
190 <<position.z()<<endl;
192 cout<<
"Start Error matrix:"<<
error<<endl;
193 cout<<"Path before start:"<< pathInMDC << endl;
194 }
195
196 G4String aParticleName(
parName[i]);
198 if(!aParticleName.contains("proton"))
199 {
200 if(charge>0) aParticleName += "+";
201 else aParticleName += "-";
202 }
203 else
204 {
205 if(charge>0) aParticleName = "proton";
206 else aParticleName = "anti_proton";
207 }
208
209 if(msgFlag)
210 {
211 cout<<"Charge: "<<charge<<endl;
212 cout<<"Particle: "<<aParticleName<<endl;
213 }
214
217 extSteppingAction->
Reset();
220 bool m_trackstatus = false;
221 int trk_startpart = 0;
222 while(!m_trackstatus)
223 {
224
225 trk_startpart++;
226 if(trk_startpart>20)
227 {cout<<"-------has modified more than 20 times---------"<<endl;break;}
228 if(myExtTrack->
Set(position,
momentum,error,aParticleName,pathInMDC,tofInMdc))
229 {
232 m_trackstatus = extSteppingAction->
TrackStop();
233 }
234 else
235 m_trackstatus = true;
236 }
237 }
238
239 }
240
242
243 if(msgFlag) cout<<"will add aExtTrack!"<<endl;
244 if(aExtTrackCol)
245 {
246 if(aExtTrack) aExtTrackCol->add(aExtTrack);
247 else if(msgFlag) cout<<"No aExtTrack!"<<endl;
248 }
249 else
250 {
251 if(msgFlag) cout<<"No aExtTrackCol!"<<endl;
252 }
253 if(msgFlag) cout<<"add a aExtTrack!"<<endl;
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272 }
273 }
274
275
276
277
278
279
280
281
282
283
284
285 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
286
287
288 DataObject *extTrackCol;
289 eventSvc()->findObject("/Event/Recon/RecExtTrackCol",extTrackCol);
290 if(extTrackCol != NULL) {
291 dataManSvc->clearSubTree("/Event/Recon/RecExtTrackCol");
292 eventSvc()->unregisterObject("/Event/Recon/RecExtTrackCol");
293 }
294
295
296
297 StatusCode sc = eventSvc()->registerObject("/Event/Recon/RecExtTrackCol", aExtTrackCol);
298 if(sc!=StatusCode::SUCCESS) {
299 log << MSG::FATAL << "Could not register RecExtTrackCol in TDS!" << endreq;
300 return( StatusCode::FAILURE);
301 }
302
303
304
305 SmartDataPtr<RecExtTrackCol> aExtTrkCol(eventSvc(),"/Event/Recon/RecExtTrackCol");
306 if (!aExtTrkCol) {
307 log << MSG::FATAL << "Can't find RecExtTrackCol in TDS!" << endreq;
308 return( StatusCode::FAILURE);
309 }
310
311 RecExtTrackCol::iterator iterOfExtTrk;
312 int j=1;
313
314 for(iterOfExtTrk = aExtTrkCol->begin();iterOfExtTrk!=aExtTrkCol->end();iterOfExtTrk++)
315 {
316 if(myResultFlag)
317 {
318 for(int i=0; i<5; i++)
319 {
320
321 cout<<"##########track"<<j<<": "<<"("<<i<<")"<<endl;
322 cout<<"******TOF1:******"<<endl;
323 cout<<"VolumeName: "<<(*iterOfExtTrk)->tof1VolumeName(i)<<"\t"
324 <<"VolumeNumber: "<<(*iterOfExtTrk)->tof1VolumeNumber(i)<<"\t"<<endl
325 <<"Position: "<<(*iterOfExtTrk)->tof1Position(i)<<"\t"
326 <<"Momentum: "<<(*iterOfExtTrk)->tof1Momentum(i)<<"\t"<<endl
327 <<"Error matrix: "<<(*iterOfExtTrk)->tof1ErrorMatrix(i)
328 <<"Error z: "<<(*iterOfExtTrk)->tof1PosSigmaAlongZ(i)<<"\t"
329 <<"Error Tz: "<<(*iterOfExtTrk)->tof1PosSigmaAlongT(i)<<"\t"
330 <<"Error x: "<<(*iterOfExtTrk)->tof1PosSigmaAlongX(i)<<"\t"
331 <<"Error y: "<<(*iterOfExtTrk)->tof1PosSigmaAlongY(i)<<endl
332 <<"Tof: "<<(*iterOfExtTrk)->tof1(i)<<"\t"
333 <<"PathOF: "<<(*iterOfExtTrk)->tof1Path(i)
334 <<endl;
335 cout<<"******TOF2:******"<<endl;
336 cout<<"VolumeName: "<<(*iterOfExtTrk)->tof2VolumeName(i)<<"\t"
337 <<"VolumeNumber: "<<(*iterOfExtTrk)->tof2VolumeNumber(i)<<"\t"<<endl
338 <<"Position: "<<(*iterOfExtTrk)->tof2Position(i)<<"\t"
339 <<"Momentum: "<<(*iterOfExtTrk)->tof2Momentum(i)<<"\t"<<endl
340 <<"Error matrix: "<<(*iterOfExtTrk)->tof2ErrorMatrix(i)
341 <<"Error z: "<<(*iterOfExtTrk)->tof2PosSigmaAlongZ(i)<<"\t"
342 <<"Error Tz: "<<(*iterOfExtTrk)->tof2PosSigmaAlongT(i)<<"\t"
343 <<"Error x: "<<(*iterOfExtTrk)->tof2PosSigmaAlongX(i)<<"\t"
344 <<"Error y: "<<(*iterOfExtTrk)->tof2PosSigmaAlongY(i)<<endl
345 <<"Tof: "<<(*iterOfExtTrk)->tof2(i)<<"\t"
346 <<"PathOF: "<<(*iterOfExtTrk)->tof2Path(i)
347 <<endl;
348
349
350 cout<<"******EMC:******"<<endl
351 <<"VolumeName: "<<(*iterOfExtTrk)->emcVolumeName(i)<<"\t"
352 <<"VolumeNumber: "<<(*iterOfExtTrk)->emcVolumeNumber(i)<<"\t"<<endl
353 <<"Position: "<<(*iterOfExtTrk)->emcPosition(i)<<"\t"
354 <<"Momentum: "<<(*iterOfExtTrk)->emcMomentum(i)<<"\t"<<endl
355 <<"Error matrix: "<<(*iterOfExtTrk)->emcErrorMatrix(i)
356 <<"Error theta: "<<(*iterOfExtTrk)->emcPosSigmaAlongTheta(i)<<"\t"
357 <<"Error phi: "<<(*iterOfExtTrk)->emcPosSigmaAlongPhi(i)<<"\t"
358 <<"EMC path: "<<(*iterOfExtTrk)->emcPath(i)
359 <<endl;
360
361
362 cout<<"******MUC:******"<<endl
363 <<"VolumeName: "<<(*iterOfExtTrk)->mucVolumeName(i)<<"\t"
364 <<"VolumeNumber: "<<(*iterOfExtTrk)->mucVolumeNumber(i)<<endl
365 <<"Position: "<<(*iterOfExtTrk)->mucPosition(i)<<"\t"
366 <<"Momentum: "<<(*iterOfExtTrk)->mucMomentum(i)<<"\t"<<endl
367 <<"Error matrix: "<<(*iterOfExtTrk)->mucErrorMatrix(i)
368 <<"Error z: "<<(*iterOfExtTrk)->mucPosSigmaAlongZ(i)<<"\t"
369 <<"Error Tz: "<<(*iterOfExtTrk)->mucPosSigmaAlongT(i)<<"\t"
370 <<"Error x: "<<(*iterOfExtTrk)->mucPosSigmaAlongX(i)<<"\t"
371 <<"Error y: "<<(*iterOfExtTrk)->mucPosSigmaAlongY(i)
372 <<endl;
373
374 cout<<"*******MUC KALMANFILTER***********"<<endl;
375 cout<<"Chisq is "<<(*iterOfExtTrk)->MucKalchi2(i)<<endl;
376 cout<<"Nfit is "<<(*iterOfExtTrk)->MucKaldof(i)<<endl;
377 cout<<"chiL "<<(*iterOfExtTrk)->MucKalchi2()<<endl;
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399 }
400 }
401 j++;
402
403 }
404
405 if(msgFlag) cout<<"****************** End a event! ****************"<<endl<<endl;
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533return StatusCode::SUCCESS;
534}
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
ObjectVector< RecExtTrack > RecExtTrackCol
void SetTrackId(int trackId)
void SetParType(int aParType=2)
double GetParticleCharge() const
double GetTrackLength() const
bool SetMdcKalTrkCol(RecMdcKalTrackCol *aPointer)
bool SetMdcRecTrkCol(RecMdcTrackCol *aPointer)
const Hep3Vector GetPosition() const
const HepSymMatrix GetErrorMatrix() const
const Hep3Vector GetMomentum() const
void SetMsgFlag(bool aFlag)
void SetExtTrackPointer(RecExtTrack *aExtTrack)
void InfmodMuc(Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
void SetMucDigiColPointer(MucDigiCol *rawdigicol)
void TrackExtrapotation()
ExtSteppingAction * GetStepAction()
bool Set(const Hep3Vector &xv3, const Hep3Vector &pv3, const HepSymMatrix &err, const std::string &particleName, const double pathInMDC, const double tofInMdc)