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