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
111 if(msgFlag)
112 {
113 cout<<"TrackExt: ******************* Start a event *******************"<<endl;
114 cout<<"run= "<<runNumber<<"; event= "<<eventNumber<<endl;
115 }
116
117
118 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
119 if(!mucDigiCol) {
120 log << MSG::FATAL << "Could not find MUC digi" << endreq;
121 return( StatusCode::SUCCESS);
122 }
123
124
127
128 bool setTrk=false;
129
130 int parID;
131 if(myParticleName=="e") parID=0;
132 else if(myParticleName=="mu") parID=1;
133 else if(myParticleName=="pi") parID=2;
134 else if(myParticleName=="kaon") parID=3;
135 else if(myParticleName=="proton"||myParticleName=="anti_proton") parID=4;
136
137 if(myInputTrk=="Mdc")
138 {
139 SmartDataPtr<RecMdcTrackCol> aMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
140 if(!aMdcTrackCol)
141 {
142 log << MSG::WARNING << "Can't find RecMdcTrackCol in TDS!" << endreq;
143 return( StatusCode::SUCCESS);
144 }
146
147 }
148 else if(myInputTrk=="Kal")
149 {
150 SmartDataPtr<RecMdcKalTrackCol> aMdcKalTrackCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
151 if(!aMdcKalTrackCol)
152 {
153 log << MSG::WARNING << "Can't find RecMdcKalTrackCol in TDS!" << endreq;
154 return( StatusCode::SUCCESS);
155 }
157 }
158 else
159 {
160 log << MSG::WARNING << "Wrong type of inputTrk:" << myInputTrk << endreq;
161 return( StatusCode::SUCCESS);
162 }
163
165 if(setTrk)
166 {
168 {
169
171
172 for(int i=0; i<5; i++)
173 {
175 {
177
184 double tofInMdc = aExtMdcTrack.
GetTrkTof();
185
186 if(msgFlag)
187 {
188 cout<<"Start From:"<<position.x()<<' '<<position.y()<<' '
189 <<position.z()<<endl;
191 cout<<
"Start Error matrix:"<<
error<<endl;
192 cout<<"Path before start:"<< pathInMDC << endl;
193 }
194
195 G4String aParticleName(
parName[i]);
197 if(!aParticleName.contains("proton"))
198 {
199 if(charge>0) aParticleName += "+";
200 else aParticleName += "-";
201 }
202 else
203 {
204 if(charge>0) aParticleName = "proton";
205 else aParticleName = "anti_proton";
206 }
207
208 if(msgFlag)
209 {
210 cout<<"Charge: "<<charge<<endl;
211 cout<<"Particle: "<<aParticleName<<endl;
212 }
213
217
218 extSteppingAction->
Reset();
221 bool m_trackstatus = false;
222 int trk_startpart = 0;
223 while(!m_trackstatus)
224 {
225
226 trk_startpart++;
227 if(trk_startpart>20)
228 {cout<<"-------has modified more than 20 times---------"<<endl;break;}
229 if(myExtTrack->
Set(position,
momentum,error,aParticleName,pathInMDC,tofInMdc))
230 {
233 m_trackstatus = extSteppingAction->
TrackStop();
234 }
235 else
236 m_trackstatus = true;
237 }
238 }
239
240 }
241
243
244 if(msgFlag) cout<<"will add aExtTrack!"<<endl;
245 if(aExtTrackCol)
246 {
247 if(aExtTrack) aExtTrackCol->add(aExtTrack);
248 else if(msgFlag) cout<<"No aExtTrack!"<<endl;
249 }
250 else
251 {
252 if(msgFlag) cout<<"No aExtTrackCol!"<<endl;
253 }
254 if(msgFlag) cout<<"add a aExtTrack!"<<endl;
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
286 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
287
288
289 DataObject *extTrackCol;
290 eventSvc()->findObject("/Event/Recon/RecExtTrackCol",extTrackCol);
291 if(extTrackCol != NULL) {
292 dataManSvc->clearSubTree("/Event/Recon/RecExtTrackCol");
293 eventSvc()->unregisterObject("/Event/Recon/RecExtTrackCol");
294 }
295
296
297
298 StatusCode sc = eventSvc()->registerObject("/Event/Recon/RecExtTrackCol", aExtTrackCol);
299 if(sc!=StatusCode::SUCCESS) {
300 log << MSG::FATAL << "Could not register RecExtTrackCol in TDS!" << endreq;
301 return( StatusCode::FAILURE);
302 }
303
304
305
306 SmartDataPtr<RecExtTrackCol> aExtTrkCol(eventSvc(),"/Event/Recon/RecExtTrackCol");
307 if (!aExtTrkCol) {
308 log << MSG::FATAL << "Can't find RecExtTrackCol in TDS!" << endreq;
309 return( StatusCode::FAILURE);
310 }
311
312 RecExtTrackCol::iterator iterOfExtTrk;
313 int j=1;
314
315 for(iterOfExtTrk = aExtTrkCol->begin();iterOfExtTrk!=aExtTrkCol->end();iterOfExtTrk++)
316 {
317 if(myResultFlag)
318 {
319 for(int i=0; i<5; i++)
320 {
321
322 cout<<"##########track"<<j<<": "<<"("<<i<<")"<<endl;
323 cout<<"******TOF1:******"<<endl;
324 cout<<"VolumeName: "<<(*iterOfExtTrk)->tof1VolumeName(i)<<"\t"
325 <<"VolumeNumber: "<<(*iterOfExtTrk)->tof1VolumeNumber(i)<<"\t"<<endl
326 <<"Position: "<<(*iterOfExtTrk)->tof1Position(i)<<"\t"
327 <<"Momentum: "<<(*iterOfExtTrk)->tof1Momentum(i)<<"\t"<<endl
328 <<"Error matrix: "<<(*iterOfExtTrk)->tof1ErrorMatrix(i)
329 <<"Error z: "<<(*iterOfExtTrk)->tof1PosSigmaAlongZ(i)<<"\t"
330 <<"Error Tz: "<<(*iterOfExtTrk)->tof1PosSigmaAlongT(i)<<"\t"
331 <<"Error x: "<<(*iterOfExtTrk)->tof1PosSigmaAlongX(i)<<"\t"
332 <<"Error y: "<<(*iterOfExtTrk)->tof1PosSigmaAlongY(i)<<endl
333 <<"Tof: "<<(*iterOfExtTrk)->tof1(i)<<"\t"
334 <<"PathOF: "<<(*iterOfExtTrk)->tof1Path(i)
335 <<endl;
336 cout<<"******TOF2:******"<<endl;
337 cout<<"VolumeName: "<<(*iterOfExtTrk)->tof2VolumeName(i)<<"\t"
338 <<"VolumeNumber: "<<(*iterOfExtTrk)->tof2VolumeNumber(i)<<"\t"<<endl
339 <<"Position: "<<(*iterOfExtTrk)->tof2Position(i)<<"\t"
340 <<"Momentum: "<<(*iterOfExtTrk)->tof2Momentum(i)<<"\t"<<endl
341 <<"Error matrix: "<<(*iterOfExtTrk)->tof2ErrorMatrix(i)
342 <<"Error z: "<<(*iterOfExtTrk)->tof2PosSigmaAlongZ(i)<<"\t"
343 <<"Error Tz: "<<(*iterOfExtTrk)->tof2PosSigmaAlongT(i)<<"\t"
344 <<"Error x: "<<(*iterOfExtTrk)->tof2PosSigmaAlongX(i)<<"\t"
345 <<"Error y: "<<(*iterOfExtTrk)->tof2PosSigmaAlongY(i)<<endl
346 <<"Tof: "<<(*iterOfExtTrk)->tof2(i)<<"\t"
347 <<"PathOF: "<<(*iterOfExtTrk)->tof2Path(i)
348 <<endl;
349
350
351 cout<<"******EMC:******"<<endl
352 <<"VolumeName: "<<(*iterOfExtTrk)->emcVolumeName(i)<<"\t"
353 <<"VolumeNumber: "<<(*iterOfExtTrk)->emcVolumeNumber(i)<<"\t"<<endl
354 <<"Position: "<<(*iterOfExtTrk)->emcPosition(i)<<"\t"
355 <<"Momentum: "<<(*iterOfExtTrk)->emcMomentum(i)<<"\t"<<endl
356 <<"Error matrix: "<<(*iterOfExtTrk)->emcErrorMatrix(i)
357 <<"Error theta: "<<(*iterOfExtTrk)->emcPosSigmaAlongTheta(i)<<"\t"
358 <<"Error phi: "<<(*iterOfExtTrk)->emcPosSigmaAlongPhi(i)<<"\t"
359 <<"EMC path: "<<(*iterOfExtTrk)->emcPath(i)
360 <<endl;
361
362
363 cout<<"******MUC:******"<<endl
364 <<"VolumeName: "<<(*iterOfExtTrk)->mucVolumeName(i)<<"\t"
365 <<"VolumeNumber: "<<(*iterOfExtTrk)->mucVolumeNumber(i)<<endl
366 <<"Position: "<<(*iterOfExtTrk)->mucPosition(i)<<"\t"
367 <<"Momentum: "<<(*iterOfExtTrk)->mucMomentum(i)<<"\t"<<endl
368 <<"Error matrix: "<<(*iterOfExtTrk)->mucErrorMatrix(i)
369 <<"Error z: "<<(*iterOfExtTrk)->mucPosSigmaAlongZ(i)<<"\t"
370 <<"Error Tz: "<<(*iterOfExtTrk)->mucPosSigmaAlongT(i)<<"\t"
371 <<"Error x: "<<(*iterOfExtTrk)->mucPosSigmaAlongX(i)<<"\t"
372 <<"Error y: "<<(*iterOfExtTrk)->mucPosSigmaAlongY(i)
373 <<endl;
374
375 cout<<"*******MUC KALMANFILTER***********"<<endl;
376 cout<<"Chisq is "<<(*iterOfExtTrk)->MucKalchi2(i)<<endl;
377 cout<<"Nfit is "<<(*iterOfExtTrk)->MucKaldof(i)<<endl;
378 cout<<"chiL "<<(*iterOfExtTrk)->MucKalchi2()<<endl;
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400 }
401 }
402 j++;
403
404 }
405
406 if(msgFlag) cout<<"****************** End a event! ****************"<<endl<<endl;
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
533
534return StatusCode::SUCCESS;
535}
**********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 Set_which_tof_version(int version)
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)