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