00001
00002
00003
00004
00005
00006 #include "GaudiKernel/MsgStream.h"
00007 #include "GaudiKernel/AlgFactory.h"
00008 #include "GaudiKernel/ISvcLocator.h"
00009 #include "GaudiKernel/SmartDataPtr.h"
00010 #include "GaudiKernel/IDataProviderSvc.h"
00011 #include "GaudiKernel/IDataManagerSvc.h"
00012 #include "GaudiKernel/PropertyMgr.h"
00013 #include "EventModel/EventHeader.h"
00014
00015 #include "CLHEP/Matrix/SymMatrix.h"
00016 #include "CLHEP/Vector/ThreeVector.h"
00017
00018 #include "ExtEvent/RecExtTrack.h"
00019 #include "TrkExtAlg/TrkExtAlg.h"
00020 #include "TrkExtAlg/ExtSteppingAction.h"
00022
00023 #include "TrackUtil/Helix.h"
00024 #include "DetVerSvc/IDetVerSvc.h"
00026 #include "TrkExtAlg/ExtMdcTrack.h"
00027
00028 using namespace Event;
00029
00030 string parName[5]={"e","mu","pi","kaon","proton"};
00031
00032
00033 TrkExtAlg::TrkExtAlg(const string& name, ISvcLocator* pSvcLocator):Algorithm(name, pSvcLocator)
00034 {
00035 myParticleName = "pi";
00036 msgFlag = false;
00037 myResultFlag = false;
00038 myInputTrk= "Kal";
00039
00040 declareProperty("ParticleName",myParticleName);
00041 declareProperty("GeantGeomOptimization",myGeomOptimization=true);
00042 declareProperty("MessageFlag",msgFlag);
00043 declareProperty("ResultMessageFlag",myResultFlag);
00044 declareProperty("BFieldOn",myBFieldOn=true);
00045 declareProperty("InputTrk",myInputTrk);
00046 declareProperty("UseMucKalFilter",myUseMucKal=true);
00047 declareProperty("MucWindow",myMucWindow=6);
00048 }
00049
00050
00051 TrkExtAlg::~TrkExtAlg()
00052 {
00053 if(myExtTrack) delete myExtTrack;
00054 }
00055
00057 StatusCode TrkExtAlg::initialize()
00058 {
00059 MsgStream log(msgSvc(), name());
00060 log << MSG::INFO << "initialize()" << endreq;
00061
00062 IDetVerSvc* detVerSvc;
00063 StatusCode sc_det = service("DetVerSvc", detVerSvc);
00064 if( sc_det.isFailure() ) {
00065 log << MSG::FATAL << "can't retrieve DetVerSvc instance" << endreq;
00066 return sc_det;
00067 }
00068
00069 m_detVer = detVerSvc->phase();
00070 log << MSG::INFO << "** ~~~~~~ ** : retrieved DetectorStage = " << m_detVer << endreq;
00071
00072 myExtTrack = new Ext_track(msgFlag,myBFieldOn,myGeomOptimization,m_detVer,myUseMucKal,myMucWindow);
00073 myExtTrack->Initialization(msgFlag,myBFieldOn,myGeomOptimization,myUseMucKal,myMucWindow);
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100 return StatusCode::SUCCESS;
00101 }
00102
00103
00105 StatusCode TrkExtAlg::execute()
00106 {
00107
00108 MsgStream log(msgSvc(), name());
00109 log << MSG::INFO << "execute()" << endreq;
00110 int eventNumber, runNumber;
00111
00112 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00113 if (!eventHeader)
00114 {
00115 log << MSG::FATAL << "Could not find Event Header" << endreq;
00116 return( StatusCode::FAILURE);
00117 }
00118 runNumber = eventHeader->runNumber();
00119 eventNumber = eventHeader->eventNumber();
00120
00121 if(msgFlag)
00122 {
00123 cout<<"TrackExt: ******************* Start a event *******************"<<endl;
00124 cout<<"run= "<<runNumber<<"; event= "<<eventNumber<<endl;
00125 }
00126
00128 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
00129 if(!mucDigiCol) {
00130 log << MSG::FATAL << "Could not find MUC digi" << endreq;
00131 return( StatusCode::SUCCESS);
00132 }
00133
00135 ExtMdcTrack aExtMdcTrack;
00136 aExtMdcTrack.SetMsgFlag(msgFlag);
00137
00138 bool setTrk=false;
00139
00140 int parID;
00141 if(myParticleName=="e") parID=0;
00142 else if(myParticleName=="mu") parID=1;
00143 else if(myParticleName=="pi") parID=2;
00144 else if(myParticleName=="kaon") parID=3;
00145 else if(myParticleName=="proton"||myParticleName=="anti_proton") parID=4;
00146
00147 if(myInputTrk=="Mdc")
00148 {
00149 SmartDataPtr<RecMdcTrackCol> aMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00150 if(!aMdcTrackCol)
00151 {
00152 log << MSG::WARNING << "Can't find RecMdcTrackCol in TDS!" << endreq;
00153 return( StatusCode::SUCCESS);
00154 }
00155 setTrk=aExtMdcTrack.SetMdcRecTrkCol(aMdcTrackCol);
00156
00157 }
00158 else if(myInputTrk=="Kal")
00159 {
00160 SmartDataPtr<RecMdcKalTrackCol> aMdcKalTrackCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
00161 if(!aMdcKalTrackCol)
00162 {
00163 log << MSG::WARNING << "Can't find RecMdcKalTrackCol in TDS!" << endreq;
00164 return( StatusCode::SUCCESS);
00165 }
00166 setTrk=aExtMdcTrack.SetMdcKalTrkCol(aMdcKalTrackCol);
00167 }
00168 else
00169 {
00170 log << MSG::WARNING << "Wrong type of inputTrk:" << myInputTrk << endreq;
00171 return( StatusCode::SUCCESS);
00172 }
00173
00174 RecExtTrackCol *aExtTrackCol = new RecExtTrackCol;
00175 if(setTrk)
00176 {
00177 while(aExtMdcTrack.GetOneGoodTrk())
00178 {
00179
00180 RecExtTrack *aExtTrack = new RecExtTrack;
00181
00182 for(int i=0; i<5; i++)
00183 {
00184 if(aExtMdcTrack.ReadTrk(i))
00185 {
00186 aExtTrack->SetParType(i);
00187
00188 int trackID = aExtMdcTrack.GetTrackID();
00189 aExtTrack->SetTrackId(trackID);
00190 Hep3Vector position = aExtMdcTrack.GetPosition();
00191 Hep3Vector momentum = aExtMdcTrack.GetMomentum();
00192 HepSymMatrix error = aExtMdcTrack.GetErrorMatrix();
00193 double pathInMDC = aExtMdcTrack.GetTrackLength();
00194 double tofInMdc = aExtMdcTrack.GetTrkTof();
00195
00196 if(msgFlag)
00197 {
00198 cout<<"Start From:"<<position.x()<<' '<<position.y()<<' '
00199 <<position.z()<<endl;
00200 cout<<"Start Momentum:"<<momentum.x()<<' '<<momentum.y()<<' '<<momentum.z()<<endl;
00201 cout<<"Start Error matrix:"<<error<<endl;
00202 cout<<"Path before start:"<< pathInMDC << endl;
00203 }
00204
00205 G4String aParticleName(parName[i]);
00206 double charge = aExtMdcTrack.GetParticleCharge();
00207 if(!aParticleName.contains("proton"))
00208 {
00209 if(charge>0) aParticleName += "+";
00210 else aParticleName += "-";
00211 }
00212 else
00213 {
00214 if(charge>0) aParticleName = "proton";
00215 else aParticleName = "anti_proton";
00216 }
00217
00218 if(msgFlag)
00219 {
00220 cout<<"Charge: "<<charge<<endl;
00221 cout<<"Particle: "<<aParticleName<<endl;
00222 }
00223
00224 ExtSteppingAction *extSteppingAction;
00225 extSteppingAction = myExtTrack->GetStepAction();
00226 extSteppingAction->Set_which_tof_version(m_detVer);
00227
00228 extSteppingAction->Reset();
00229 extSteppingAction->SetMucDigiColPointer(mucDigiCol);
00230 extSteppingAction->SetExtTrackPointer(aExtTrack);
00231 bool m_trackstatus = false;
00232 int trk_startpart = 0;
00233 while(!m_trackstatus)
00234 {
00235
00236 trk_startpart++;
00237 if(trk_startpart>20)
00238 {cout<<"-------has modified more than 20 times---------"<<endl;break;}
00239 if(myExtTrack->Set(position,momentum,error,aParticleName,pathInMDC,tofInMdc))
00240 {
00241 myExtTrack->TrackExtrapotation();
00242 extSteppingAction->InfmodMuc(position,momentum,error);
00243 m_trackstatus = extSteppingAction->TrackStop();
00244 }
00245 else
00246 m_trackstatus = true;
00247 }
00248 }
00249
00250 }
00251
00252 aExtTrack->SetParType(parID);
00253
00254 if(msgFlag) cout<<"will add aExtTrack!"<<endl;
00255 if(aExtTrackCol)
00256 {
00257 if(aExtTrack) aExtTrackCol->add(aExtTrack);
00258 else if(msgFlag) cout<<"No aExtTrack!"<<endl;
00259 }
00260 else
00261 {
00262 if(msgFlag) cout<<"No aExtTrackCol!"<<endl;
00263 }
00264 if(msgFlag) cout<<"add a aExtTrack!"<<endl;
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283 }
00284 }
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00297
00298
00299 DataObject *extTrackCol;
00300 eventSvc()->findObject("/Event/Recon/RecExtTrackCol",extTrackCol);
00301 if(extTrackCol != NULL) {
00302 dataManSvc->clearSubTree("/Event/Recon/RecExtTrackCol");
00303 eventSvc()->unregisterObject("/Event/Recon/RecExtTrackCol");
00304 }
00305
00306
00307
00308 StatusCode sc = eventSvc()->registerObject("/Event/Recon/RecExtTrackCol", aExtTrackCol);
00309 if(sc!=StatusCode::SUCCESS) {
00310 log << MSG::FATAL << "Could not register RecExtTrackCol in TDS!" << endreq;
00311 return( StatusCode::FAILURE);
00312 }
00313
00314
00315
00316 SmartDataPtr<RecExtTrackCol> aExtTrkCol(eventSvc(),"/Event/Recon/RecExtTrackCol");
00317 if (!aExtTrkCol) {
00318 log << MSG::FATAL << "Can't find RecExtTrackCol in TDS!" << endreq;
00319 return( StatusCode::FAILURE);
00320 }
00321
00322 RecExtTrackCol::iterator iterOfExtTrk;
00323 int j=1;
00324
00325 for(iterOfExtTrk = aExtTrkCol->begin();iterOfExtTrk!=aExtTrkCol->end();iterOfExtTrk++)
00326 {
00327 if(myResultFlag)
00328 {
00329 for(int i=0; i<5; i++)
00330 {
00331
00332 cout<<"##########track"<<j<<": "<<"("<<i<<")"<<endl;
00333 cout<<"******TOF1:******"<<endl;
00334 cout<<"VolumeName: "<<(*iterOfExtTrk)->tof1VolumeName(i)<<"\t"
00335 <<"VolumeNumber: "<<(*iterOfExtTrk)->tof1VolumeNumber(i)<<"\t"<<endl
00336 <<"Position: "<<(*iterOfExtTrk)->tof1Position(i)<<"\t"
00337 <<"Momentum: "<<(*iterOfExtTrk)->tof1Momentum(i)<<"\t"<<endl
00338 <<"Error matrix: "<<(*iterOfExtTrk)->tof1ErrorMatrix(i)
00339 <<"Error z: "<<(*iterOfExtTrk)->tof1PosSigmaAlongZ(i)<<"\t"
00340 <<"Error Tz: "<<(*iterOfExtTrk)->tof1PosSigmaAlongT(i)<<"\t"
00341 <<"Error x: "<<(*iterOfExtTrk)->tof1PosSigmaAlongX(i)<<"\t"
00342 <<"Error y: "<<(*iterOfExtTrk)->tof1PosSigmaAlongY(i)<<endl
00343 <<"Tof: "<<(*iterOfExtTrk)->tof1(i)<<"\t"
00344 <<"PathOF: "<<(*iterOfExtTrk)->tof1Path(i)
00345 <<endl;
00346 cout<<"******TOF2:******"<<endl;
00347 cout<<"VolumeName: "<<(*iterOfExtTrk)->tof2VolumeName(i)<<"\t"
00348 <<"VolumeNumber: "<<(*iterOfExtTrk)->tof2VolumeNumber(i)<<"\t"<<endl
00349 <<"Position: "<<(*iterOfExtTrk)->tof2Position(i)<<"\t"
00350 <<"Momentum: "<<(*iterOfExtTrk)->tof2Momentum(i)<<"\t"<<endl
00351 <<"Error matrix: "<<(*iterOfExtTrk)->tof2ErrorMatrix(i)
00352 <<"Error z: "<<(*iterOfExtTrk)->tof2PosSigmaAlongZ(i)<<"\t"
00353 <<"Error Tz: "<<(*iterOfExtTrk)->tof2PosSigmaAlongT(i)<<"\t"
00354 <<"Error x: "<<(*iterOfExtTrk)->tof2PosSigmaAlongX(i)<<"\t"
00355 <<"Error y: "<<(*iterOfExtTrk)->tof2PosSigmaAlongY(i)<<endl
00356 <<"Tof: "<<(*iterOfExtTrk)->tof2(i)<<"\t"
00357 <<"PathOF: "<<(*iterOfExtTrk)->tof2Path(i)
00358 <<endl;
00359
00360
00361 cout<<"******EMC:******"<<endl
00362 <<"VolumeName: "<<(*iterOfExtTrk)->emcVolumeName(i)<<"\t"
00363 <<"VolumeNumber: "<<(*iterOfExtTrk)->emcVolumeNumber(i)<<"\t"<<endl
00364 <<"Position: "<<(*iterOfExtTrk)->emcPosition(i)<<"\t"
00365 <<"Momentum: "<<(*iterOfExtTrk)->emcMomentum(i)<<"\t"<<endl
00366 <<"Error matrix: "<<(*iterOfExtTrk)->emcErrorMatrix(i)
00367 <<"Error theta: "<<(*iterOfExtTrk)->emcPosSigmaAlongTheta(i)<<"\t"
00368 <<"Error phi: "<<(*iterOfExtTrk)->emcPosSigmaAlongPhi(i)<<"\t"
00369 <<"EMC path: "<<(*iterOfExtTrk)->emcPath(i)
00370 <<endl;
00371
00372
00373 cout<<"******MUC:******"<<endl
00374 <<"VolumeName: "<<(*iterOfExtTrk)->mucVolumeName(i)<<"\t"
00375 <<"VolumeNumber: "<<(*iterOfExtTrk)->mucVolumeNumber(i)<<endl
00376 <<"Position: "<<(*iterOfExtTrk)->mucPosition(i)<<"\t"
00377 <<"Momentum: "<<(*iterOfExtTrk)->mucMomentum(i)<<"\t"<<endl
00378 <<"Error matrix: "<<(*iterOfExtTrk)->mucErrorMatrix(i)
00379 <<"Error z: "<<(*iterOfExtTrk)->mucPosSigmaAlongZ(i)<<"\t"
00380 <<"Error Tz: "<<(*iterOfExtTrk)->mucPosSigmaAlongT(i)<<"\t"
00381 <<"Error x: "<<(*iterOfExtTrk)->mucPosSigmaAlongX(i)<<"\t"
00382 <<"Error y: "<<(*iterOfExtTrk)->mucPosSigmaAlongY(i)
00383 <<endl;
00384
00385 cout<<"*******MUC KALMANFILTER***********"<<endl;
00386 cout<<"Chisq is "<<(*iterOfExtTrk)->MucKalchi2(i)<<endl;
00387 cout<<"Nfit is "<<(*iterOfExtTrk)->MucKaldof(i)<<endl;
00388 cout<<"chiL "<<(*iterOfExtTrk)->MucKalchi2()<<endl;
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410 }
00411 }
00412 j++;
00413
00414 }
00415
00416 if(msgFlag) cout<<"****************** End a event! ****************"<<endl<<endl;
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544 return StatusCode::SUCCESS;
00545 }
00546
00548 StatusCode TrkExtAlg::finalize()
00549 {
00550 MsgStream log(msgSvc(), name());
00551 log << MSG::INFO << "finalize()" << endreq;
00552
00553
00554
00555
00556 return StatusCode::SUCCESS;
00557 }