169 #include "TProfile.h"
174 #include "TGeoManager.h"
177 #include "StTrackNode.h"
178 #include "StContainers.h"
179 #include "StPrimaryVertex.h"
180 #include "StVertex.h"
181 #include "StMeasuredPoint.h"
182 #include "StDedxPidTraits.h"
183 #include "StTrackPidTraits.h"
184 #include "StTrackGeometry.h"
185 #include "StTrackDetectorInfo.h"
186 #include "StGlobalTrack.h"
187 #include "StParticleTypes.hh"
188 #include "StThreeVector.hh"
189 #include "StThreeVectorF.hh"
190 #include "StThreeVectorD.hh"
191 #include "StLorentzVectorD.hh"
192 #include "StGlobals.hh"
193 #include "StGetConfigValue.hh"
194 #include "StTimer.hh"
195 #include "StMemoryInfo.hh"
196 #include "StMessMgr.h"
198 #include "St_DataSet.h"
199 #include "St_DataSetIter.h"
200 #include "StEventTypes.h"
201 #include "StTpcDedxPidAlgorithm.h"
202 #include "StParticleDefinition.hh"
203 #include "StPhysicalHelix.hh"
204 #include "StHelixD.hh"
206 #include "PhysicalConstants.h"
207 #include "SystemOfUnits.h"
209 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
210 #include "StMuDSTMaker/COMMON/StMuEvent.h"
211 #include "StMuDSTMaker/COMMON/StMuDst.h"
212 #include "StMuDSTMaker/COMMON/StMuTrack.h"
213 #include "StMuDSTMaker/COMMON/StMuMtdCollection.h"
214 #include "StMuDSTMaker/COMMON/StMuMtdRawHit.h"
215 #include "StMuDSTMaker/COMMON/StMuMtdHeader.h"
216 #include "StMuDSTMaker/COMMON/StMuMtdHit.h"
217 #include "StMuDSTMaker/COMMON/StMuBTofHit.h"
218 #include "StMuDSTMaker/COMMON/StMuBTofPidTraits.h"
219 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
220 #include "StMuDSTMaker/COMMON/StMuMtdPidTraits.h"
222 #include "StEventMaker/StEventMaker.h"
223 #include "StAssociationMaker/StAssociationMaker.h"
224 #include "StMcEventMaker/StMcEventMaker.h"
225 #include "StDetectorDbMaker/St_MagFactorC.h"
226 #include "tables/St_vertexSeed_Table.h"
228 #include "StBTofPidTraits.h"
229 #include "StBTofCollection.h"
230 #include "StBTofHit.h"
231 #include "StBTofRawHit.h"
232 #include "StBTofHeader.h"
233 #include "StBTofUtil/tofPathLength.hh"
234 #include "StBTofUtil/StBTofGeometry.h"
235 #include "StBTofUtil/StBTofDaqMap.h"
236 #include "StBTofUtil/StBTofHitCollection.h"
237 #include "StMtdPidTraits.h"
239 #include "StMtdUtil/StMtdGeometry.h"
240 #include "StMtdMatchMaker.h"
249 doPrintMemoryInfo = kFALSE;
250 doPrintCpuInfo = kFALSE;
253 mMinFitPointsPerTrack=15;
262 memset(mVDrift,0,
sizeof(mVDrift));
265 index2Primary.clear();
268 mLockBField = kFALSE;
277 fZReso =
new TF1(
"fZReso",
"sqrt([0]/x/x+[1])",0,100);
278 fZReso->SetParameters(148.7,1.654);
279 fPhiReso =
new TF1(
"fPhiReso",
"sqrt([0]/x/x+[1])",0,100);
280 fPhiReso->SetParameters(9.514e-4,7.458e-6);
282 mEventCounterHisto = NULL;
283 mCellsMultInEvent = NULL;
284 mHitsMultInEvent = NULL;
285 mHitsPrimaryInEvent = NULL;
286 mHitsGlobalInEvent = NULL;
287 mHitsMultPerTrack = NULL;
288 mHitsPosition = NULL;
289 for(
int i=0; i<mNBacklegs; i++)
291 mDaqOccupancy[i] = NULL;
292 mDaqOccupancyProj[i] = NULL;
294 mHitCorrModule[i] = NULL;
295 mDeltaHitFinal[i] = NULL;
299 mTrackNFitPts = NULL;
300 mTrackdEdxvsp = NULL;
301 mNSigmaPivsPt = NULL;
302 mCellsPerEventMatch1 = NULL;
303 mHitsPerEventMatch1 = NULL;
304 mCellsPerTrackMatch1 = NULL;
305 mTracksPerCellMatch1 = NULL;
306 mDaqOccupancyMatch1 = NULL;
307 mDeltaHitMatch1 = NULL;
308 mCellsPerEventMatch2 = NULL;
309 mHitsPerEventMatch2 = NULL;
310 mCellsPerTrackMatch2 = NULL;
311 mTracksPerCellMatch2 = NULL;
312 mDaqOccupancyMatch2 = NULL;
313 mDeltaHitMatch2 = NULL;
314 mCellsPerEventMatch3 = NULL;
315 mHitsPerEventMatch3 = NULL;
316 mCellsPerTrackMatch3 = NULL;
317 mTracksPerCellMatch3 = NULL;
318 mDaqOccupancyMatch3 = NULL;
319 mDeltaHitMatch3 = NULL;
320 mCellsPrimaryPerEventMatch3 = NULL;
322 hTofPhivsProj = NULL;
325 hMtdPhivsProj = NULL;
326 hMtddPhivsBackleg = NULL;
327 hMtddZvsBackleg = NULL;
330 StMtdMatchMaker::~StMtdMatchMaker(){
335 void StMtdMatchMaker::Clear(
const char*){
342 for(
int i=0;i<mNAllTrays;i++){
343 for(
int j=0;j<mNStrips;j++){
347 if(mHisto) bookHistograms();
349 return StMaker::Init();
353 void StMtdMatchMaker::bookHistograms(){
355 mEventCounterHisto =
new TH1D(
"eventCounter",
"eventCounter",20,0,20);
356 AddHist(mEventCounterHisto);
358 mCellsMultInEvent =
new TH1D(
"cellsPerEvent",
"cellsPerEvent",1000,0,1000);
359 AddHist(mCellsMultInEvent);
361 mHitsMultInEvent =
new TH1D(
"hitsPerEvent",
"hitsPerEvent",1000,0,1000);
362 AddHist(mHitsMultInEvent);
364 mHitsPrimaryInEvent =
new TH1D(
"hitsPrimaryPerEvent",
"hitsPrimaryPerEvent",1000,0,1000);
365 AddHist(mHitsPrimaryInEvent);
367 mHitsMultPerTrack =
new TH1D(
"hitsPerTrack",
"hitsPerTrack",10,0,10);
368 AddHist(mHitsMultPerTrack);
370 mHitsPosition =
new TH2D(
"hitsPosition",
"hitsPositions",1000,-500.,500.,1000,-500.,500);
371 AddHist(mHitsPosition);
377 for(
int i=0;i<mNBacklegs;i++) {
379 sprintf(hisname,
"Occupancy_Backleg_%d",i+1);
380 mDaqOccupancy[i] =
new TH1D(hisname,
";fired cell(module*12+stripId)",60,0,60);
381 sprintf(hisname,
"OccupancyProj_Backleg_%d",i+1);
382 mDaqOccupancyProj[i] =
new TH1D(hisname,
";projected cell",60,0,60);
384 AddHist(mDaqOccupancy[i]);
385 AddHist(mDaqOccupancyProj[i]);
389 for(
int i=0;i<mNBacklegs;i++) {
391 sprintf(hisname,
"Corr_Backleg_%d",i+1);
392 mHitCorr[i] =
new TH2D(hisname,
";project stripId;fired stripId",60,0,60,60,0,60);
393 sprintf(hisname,
"Corr_Backleg_%d_module",i+1);
394 mHitCorrModule[i] =
new TH2D(hisname,
";project moduleId;fired moduleId",6,0,6,6,0,6);
395 AddHist(mHitCorr[i]);
396 AddHist(mHitCorrModule[i]);
400 for(
int i=0;i<mNBacklegs;i++) {
402 sprintf(hisname,
"LocalYZ_Backleg_%d",i+1);
403 mDeltaHitFinal[i] =
new TH2D(hisname,
";localY;localZ",300,-15.,15.,320,-80.,80);
404 AddHist(mDeltaHitFinal[i]);
407 mTrackPtEta =
new TH2D(
"trackPtEta",
";p_{T} (GeV/c);#eta",500,0.,10.,60,-1.5,1.5);
408 mTrackPtPhi =
new TH2D(
"trackPtPhi",
";p_{T} (GeV/c);#phi",500,0.,10.,120,0.,2.*M_PI);
409 mTrackNFitPts =
new TH1D(
"trackNFitPts",
";nHitsFit",50,0.,50.);
410 mTrackdEdxvsp =
new TH2D(
"trackdEdxvsp",
";p_{T} (GeV/c);dE/dx",500,0.,10.,1000,0.,10.);
411 mNSigmaPivsPt =
new TH2D(
"nSigmaPivsPt",
";p_{T} (GeV/c);n#sigma_{#pi}",500,0.,10.,1000,-10.,10.);
414 mCellsPerEventMatch1 =
new TH1D(
"cellsPerEventMatch1",
"cellPerEventMatch1;# hits matched with track(s) per event",10,0,10);
415 mHitsPerEventMatch1 =
new TH1D(
"hitsPerEventMatch1",
"hitsPerEventMatch1;# tracks matched with hits per event",10,0,10);
416 mCellsPerTrackMatch1 =
new TH1D(
"cellsPerTrackMatch1",
"cellsPerTrackMatch1;# tracks matched with same hit",10,0,10);
417 mTracksPerCellMatch1 =
new TH1D(
"tracksPerCellMatch1",
"tracksPerCellMatch1;# hits matched with same track",10,0,10);
418 mDaqOccupancyMatch1 =
new TH1D(
"daqPerCellMatch1",
"daqPerCellMatch1;stripId",60,0,60);
419 mDeltaHitMatch1 =
new TH2D(
"deltaHitMatch1",
"deltaHitMatch1;localY;localZ",300,-15,15,800,-80.,80);
422 mCellsPerEventMatch2 =
new TH1D(
"cellsPerEventMatch2",
"cellPerEventMatch2;# hits matched with track(s) per event",10,0,10);
423 mHitsPerEventMatch2 =
new TH1D(
"hitsPerEventMatch2",
"hitsPerEventMatch2;# tracks matched with hits per event",10,0,10);
424 mCellsPerTrackMatch2 =
new TH1D(
"cellsPerTrackMatch2",
"cellsPerTrackMatch2;# tracks matched with same hit",10,0,10);
425 mTracksPerCellMatch2 =
new TH1D(
"tracksPerCellMatch2",
"tracksPerCellMatch2;# hits matched with same track",10,0,10);
426 mDaqOccupancyMatch2 =
new TH1D(
"daqPerCellMatch2",
"daqPerCellMatch2;stripId",60,0,60);
427 mDeltaHitMatch2 =
new TH2D(
"deltaHitMatch2",
"deltaHitMatch2;localY;localZ",300,-15,15,800,-80.,80);
430 mCellsPerEventMatch3 =
new TH1D(
"cellsPerEventMatch3",
"cellPerEventMatch3;# hits matched with track(s) per event",10,0,10);
431 mHitsPerEventMatch3 =
new TH1D(
"hitsPerEventMatch3",
"hitsPerEventMatch3;# tracks matched with hits per event",10,0,10);
432 mCellsPerTrackMatch3 =
new TH1D(
"cellsPerTrackMatch3",
"cellsPerTrackMatch3;# tracks matched with same hit",10,0,10);
433 mTracksPerCellMatch3 =
new TH1D(
"tracksPerCellMatch3",
"tracksPerCellMatch3;# hits matched with same track",10,0,10);
434 mDaqOccupancyMatch3 =
new TH1D(
"daqPerCellMatch3",
"daqPerCellMatch3;stripId",60,0,60);
435 mDeltaHitMatch3 =
new TH2D(
"deltaHitMatch3",
"deltaHitMatch3;localY;localZ",300,-15,15,800,-80.,80);
437 mCellsPrimaryPerEventMatch3 =
new TH1D(
"cellsPrimaryPerEventMatch3",
"cellsPrimaryPerEventMatch3",10,0,10);
439 hphivsz =
new TH2F(
"hphivsz",
"hphivsz",500,0,TMath::Pi()*2,500,-500,500);
440 hTofPhivsProj=
new TH2F(
"hTofPhivsProj",
"hTofPhivsProj",100,0,TMath::Pi()*2,100,0,TMath::Pi()*2);
441 hTofZvsProj=
new TH2F(
"hTofZvsProj",
"hTofzvsProj",600,-300,300,600,-300,300);
442 hMtdPhivsProj=
new TH2F(
"hMtdPhivsProj",
"hMtdPhivsProj;projected #phi; fired #phi",360,-TMath::Pi(),TMath::Pi(),180,0,2.*TMath::Pi());
443 hMtddPhivsBackleg=
new TH2F(
"hMtddPhivsBackleg",
"hMtddPhivsBackleg;backleg; d#phi",30,0,30,1000,-2.*TMath::Pi(),2.*TMath::Pi());
444 hMtddZvsBackleg =
new TH2F(
"hMtddZvsBackleg",
"hMtddZvsBackleg;backleg; dz",30,0,30,400,-200,200);
445 hMtdZvsProj=
new TH2F(
"hMtdZvsProj",
"hMtdzvsProj;projected z(cm); fired z(cm)",600,-300,300,600,-300,300);
447 AddHist(mTrackPtEta);
448 AddHist(mTrackPtPhi);
449 AddHist(mTrackNFitPts);
450 AddHist(mTrackdEdxvsp);
451 AddHist(mNSigmaPivsPt);
453 AddHist(mCellsPerEventMatch1);
454 AddHist(mHitsPerEventMatch1);
455 AddHist(mCellsPerTrackMatch1);
456 AddHist(mTracksPerCellMatch1);
457 AddHist(mDaqOccupancyMatch1);
458 AddHist(mDeltaHitMatch1);
460 AddHist(mCellsPerEventMatch2);
461 AddHist(mHitsPerEventMatch2);
462 AddHist(mCellsPerTrackMatch2);
463 AddHist(mTracksPerCellMatch2);
464 AddHist(mDaqOccupancyMatch2);
465 AddHist(mDeltaHitMatch2);
467 AddHist(mCellsPerEventMatch3);
468 AddHist(mHitsPerEventMatch3);
469 AddHist(mCellsPerTrackMatch3);
470 AddHist(mTracksPerCellMatch3);
471 AddHist(mDaqOccupancyMatch3);
472 AddHist(mDeltaHitMatch3);
474 AddHist(mCellsPrimaryPerEventMatch3);
477 AddHist(hTofPhivsProj);
478 AddHist(hTofZvsProj);
479 AddHist(hMtdZvsProj);
480 AddHist(hMtdPhivsProj);
481 AddHist(hMtddPhivsBackleg);
482 AddHist(hMtddZvsBackleg);
491 mYear= (Int_t)(runnumber/1000000) + 1999;
492 LOG_INFO <<
"Run year = " << mYear << endm;
497 LOG_INFO <<
"Initializing MTD Geometry:" << endm;
500 LOG_INFO <<
"TGeoManager (" << gGeoManager->GetName() <<
"," << gGeoManager->GetTitle() <<
") exists" << endm;
504 GetDataBase(
"VmcGeometry");
507 int year = GetDateTime().GetYear();
508 if(mGeomTag.Length()==0)
510 mGeomTag = Form(
"y%da",year);
511 if(year<2012) mGeomTag =
"y2014a";
512 LOG_INFO <<
"Load default geometry " << mGeomTag.Data() <<
" for year " << year << endm;
516 LOG_INFO <<
"Load input geometry " << mGeomTag.Data() <<
" for year " << year << endm;
519 TString ts = Form(
"$STAR/StarVMC/Geometry/macros/loadStarGeometry.C(\"%s\",1)",mGeomTag.Data());
521 gROOT->Macro(ts.Data(),&ierr);
529 mMtdGeom =
new StMtdGeometry(
"mtdGeometry",
"mtdGeometry in StMtdMatchMaker",gGeoManager);
530 mMtdGeom->SetCosmicFlag(mCosmicFlag);
531 mMtdGeom->SetELossFlag(mELossFlag);
532 mMtdGeom->SetNExtraCells(mNExtraCells);
533 if(mLockBField) mMtdGeom->SetLockBField(mLockBField);
534 mMtdGeom->Init(
this);
535 LOG_INFO<<
" Created a new mtdGeometry ..."<<endm;
537 assert(gMtdGeometry);
538 mMtdGeom = gMtdGeometry;
547 TDataSet* dbDataSet = this->GetDataBase(
"Calibrations/rhic");
550 vertexSeed_st* vSeed = ((St_vertexSeed*) (dbDataSet->FindObject(
"vertexSeed")))->GetTable();
558 LOG_INFO <<
"StMtdMatchMaker -- No Database for beamline" << endm;
561 LOG_INFO <<
"BeamLine Constraint (StMtdMatchMaker): " << endm;
562 LOG_INFO <<
"x(z) = " << x0 <<
" + " << dxdz <<
" * z" << endm;
563 LOG_INFO <<
"y(z) = " << y0 <<
" + " << dydz <<
" * z" << endm;
566 double pt = 88889999;
567 double nxy=::sqrt(dxdz*dxdz + dydz*dydz);
569 LOG_WARN <<
"StMtdMatchMaker:: Beam line must be tilted!" << endm;
602 mEvent=(
StEvent *) GetInputDS(
"StEvent");
606 mMuDst = (
StMuDst *)GetInputDS(
"MuDst");
610 if(mHisto) mEventCounterHisto->Fill(0);
612 if(doPrintCpuInfo) timer.start();
613 if(doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
616 if(mMuDstIn) cleanUpMtdPidTraits();
622 mtdCellHitVector daqCellsHitVec;
623 idVector validModuleVec;
624 if(!readMtdHits(daqCellsHitVec,validModuleVec))
return kStOK;
627 LOG_INFO<<
" Sect.A: =============================="<<endm;
628 LOG_INFO <<
" total # of cells =" << daqCellsHitVec.size() << endm;
629 for(
size_t iv=0;iv<validModuleVec.size();iv++){
630 LOG_INFO <<
" module# "<< validModuleVec[iv] <<
" Valid! "<<endm;
632 LOG_INFO<<
" daqCellsHitVec:"<<endm;
633 mtdCellHitVectorIter daqIter = daqCellsHitVec.begin();
634 for(
unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
635 int daqChn = (daqIter->module-1)*12 + (daqIter->cell);
636 int daqBackleg = daqIter->backleg;
637 LOG_INFO<<
" daq backleg:"<<daqBackleg<<
" chn:"<<daqChn<<endm;
642 mCellsMultInEvent->Fill(daqCellsHitVec.size());
643 if(daqCellsHitVec.size()) mEventCounterHisto->Fill(6);
647 LOG_INFO <<
"CPU time after Step A - loading hits: "
648 << timer.elapsedTime() <<
" sec"<<endm;
653 mtdCellHitVector allCellsHitVec;
654 Int_t nPrimaryHits(0);
656 project2Mtd(daqCellsHitVec,allCellsHitVec,nPrimaryHits);
659 mHitsMultInEvent->Fill(allCellsHitVec.size());
660 mHitsPrimaryInEvent->Fill(nPrimaryHits);
661 if(allCellsHitVec.size()) mEventCounterHisto->Fill(7);
665 LOG_INFO<<
" Sect.B: =============================="<<endm;
666 LOG_INFO<<
" allCellsHitVec:"<<endm;
667 mtdCellHitVectorIter proIter = allCellsHitVec.begin();
668 for(
unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
669 int proChn = (proIter->module-1)*12 + (proIter->cell);
670 int proBackleg = proIter->backleg;
671 LOG_INFO<<
" proj backleg:"<<proBackleg<<
" chn:"<<proChn<<endm;
677 LOG_INFO <<
"CPU time after Step B - projection:"
678 << timer.elapsedTime() <<
"sec" <<endm;
684 mtdCellHitVector matchHitCellsVec;
685 matchMtdHits(daqCellsHitVec,allCellsHitVec,matchHitCellsVec);
688 LOG_INFO<<
" Sect.C: =============================="<<endm;
689 LOG_INFO<<
" matchCellsHitVec:"<<endm;
690 mtdCellHitVectorIter matIter = matchHitCellsVec.begin();
691 for(
unsigned int imat=0;imat<matchHitCellsVec.size();imat++, matIter++) {
692 int matChn = (matIter->module-1)*12 + (matIter->cell);
693 int matBackleg = matIter->backleg;
694 LOG_INFO<<
" match backleg:"<<matBackleg<<
" chn:"<<matChn<<endm;
699 if(Debug()){ LOG_INFO <<
"C:before/after:"<< allCellsHitVec.size()<<
"/"<<matchHitCellsVec.size()<<endm;}
700 if(mHisto&&matchHitCellsVec.size()) mEventCounterHisto->Fill(8);
703 LOG_INFO <<
"CPU time after Step C - matching:"
704 <<timer.elapsedTime()<<
"sec" <<endm;
710 mtdCellHitVector singleHitCellsVec;
711 mtdCellHitVector multiHitsCellsVec;
712 sortSingleAndMultiHits(matchHitCellsVec,singleHitCellsVec,multiHitsCellsVec);
715 LOG_INFO<<
" Sect.D: =============================="<<endm;
716 LOG_INFO<<
" singleCellsHitVec:"<<endm;
717 mtdCellHitVectorIter singleIter = singleHitCellsVec.begin();
718 for(
unsigned int isingle=0;isingle<singleHitCellsVec.size();isingle++, singleIter++) {
719 int singleChn = (singleIter->module-1)*12 + (singleIter->cell);
720 int singleBackleg = singleIter->backleg;
721 LOG_INFO<<
" single backleg:"<<singleBackleg<<
" chn:"<<singleChn<<endm;
722 LOG_INFO<<
" LeTimeWest:"<<singleIter->leadingEdgeTime.first<<
" LeTimeEast:"<<singleIter->leadingEdgeTime.second<<endm;
724 LOG_INFO<<
" multiCellsHitVec:"<<endm;
725 mtdCellHitVectorIter multiIter = multiHitsCellsVec.begin();
726 for(
unsigned int imulti=0;imulti<multiHitsCellsVec.size();imulti++, multiIter++) {
727 int multiChn = (multiIter->module-1)*12 + (multiIter->cell);
728 int multiBackleg = multiIter->backleg;
729 LOG_INFO<<
" multi backleg:"<<multiBackleg<<
" chn:"<<multiChn<<endm;
735 LOG_INFO <<
" CPU time after Step D - sorting:"
736 <<timer.elapsedTime()<<
" sec"<<endm;
742 mtdCellHitVector finalMatchedCellsVec;
743 finalMatchedMtdHits(singleHitCellsVec,finalMatchedCellsVec);
746 LOG_INFO<<
" Sect.E: =============================="<<endm;
747 LOG_INFO<<
" finalMatchedCellsHitVec:"<<endm;
748 mtdCellHitVectorIter finalIter = finalMatchedCellsVec.begin();
749 for(
unsigned int ifinal=0;ifinal<finalMatchedCellsVec.size();ifinal++, finalIter++) {
750 int finalChn = (finalIter->module-1)*12 + (finalIter->cell);
751 int finalBackleg = finalIter->backleg;
752 LOG_INFO<<
" final backleg:"<<finalBackleg<<
" chn:"<<finalChn<<endm;
756 if(Debug()){ LOG_INFO <<
"E: before/after:"<< singleHitCellsVec.size() <<
"/" <<finalMatchedCellsVec.size()<<endm;}
758 LOG_INFO <<
"CPU time after Step E - final matched:"
759 <<timer.elapsedTime() <<
" sec"<<endm;
766 if(finalMatchedCellsVec.size()) mEventCounterHisto->Fill(10);
767 mCellsPerEventMatch3->Fill(finalMatchedCellsVec.size());
769 Int_t nValidSingleHitCells(0),nValidSinglePrimHitCells(0);
770 fillPidTraits(finalMatchedCellsVec,nValidSingleHitCells,nValidSinglePrimHitCells);
773 LOG_INFO<<
" Sect.F: =============================="<<endm;
774 LOG_INFO<<
" nValidSingleHitCells:"<<nValidSingleHitCells<<endm;
776 LOG_INFO <<
"#(daq):" <<daqCellsHitVec.size()
777 <<
"#(proj):"<<allCellsHitVec.size()
778 <<
"#(prim proj):" <<nPrimaryHits
779 <<
"#(single valid):"<<nValidSingleHitCells <<endm;
780 LOG_INFO <<
"F: before/after"<< finalMatchedCellsVec.size()<<
"/"<<nValidSingleHitCells <<endm;
784 LOG_INFO <<
"CPU time after Step F - final :"
785 << timer.elapsedTime()<<
"sec"<<endm;
792 if(theMtd->hitsPresent()){
793 LOG_DEBUG <<
" MtdCollection: hit container present. "<<endm;
795 LOG_INFO <<
"# of hits in this event:" << theMtd->hitsPresent() <<endm;
796 for(Int_t i=0;i<theMtd->hitsPresent();i++){
799 LOG_INFO <<
"backleg:"<<aHit->backleg()
800 <<
" module:"<<aHit->backleg()
801 <<
" cell:"<<aHit->cell()
802 <<
" tof:"<<aHit->tof()<<endm;
810 if(theMtd->hitsPresent()){
811 LOG_DEBUG <<
" MtdCollection: hit container present. "<<endm;
813 StSPtrVecMtdHit& tmpCellMtdVec = theMtd->mtdHits();
814 LOG_INFO <<
"# of hits in this event:" << tmpCellMtdVec.size()<<endm;
815 for(
size_t i=0;i<tmpCellMtdVec.size();i++){
817 LOG_INFO <<(*p)<<endm;
824 if(doPrintMemoryInfo){
825 StMemoryInfo::instance()->snapshot();
826 StMemoryInfo::instance()->print();
830 LOG_INFO<<
"CPU time for StMtdMatchMaker::Make():"
831 <<timer.elapsedTime()<<
"sec"<<endm;
834 LOG_INFO<<
"StMtdMatchMaker -- bye--bye"<<endm;
842 void StMtdMatchMaker::cleanUpMtdPidTraits()
844 index2Primary.clear();
845 for(Int_t ii=0;ii<(Int_t)mMuDst->
array(muPrimary)->GetEntries();ii++)
848 if(!pTrack)
continue;
850 if(index2Global<0)
continue;
851 index2Primary[index2Global] = ii;
854 UInt_t Nnodes = mMuDst->numberOfGlobalTracks();
855 for(UInt_t iNode=0;iNode<Nnodes;iNode++)
858 if(!theTrack)
continue;
859 if(theTrack->index2MtdHit()<0)
continue;
863 theTrack->setMtdPidTraits(pidMtd);
864 theTrack->setIndex2MtdHit(-999);
867 map<Int_t, Int_t>::iterator it = index2Primary.find(iNode);
868 if(it!=index2Primary.end()){
874 thePrimaryTrack->setMtdPidTraits(pidMtd);
875 thePrimaryTrack->setIndex2MtdHit(-999);
883 Bool_t StMtdMatchMaker::readMtdHits(mtdCellHitVector& daqCellsHitVec,idVector& validModuleVec){
888 LOG_INFO <<
"No Mudst ... bye-bye" <<endm;
892 int nMtdHits = mMuDst->numberOfMTDHit();
893 if(nMtdHits<=0)
return kFALSE;
894 for(Int_t i=0;i<nMtdHits;i++){
897 if(aHit->backleg()<=0||aHit->backleg()>mNBacklegs)
continue;
900 aHit->setIndex2Primary(-1);
901 aHit->setIndex2Global(-1);
902 aHit->setAssociatedTrackKey(-1);
903 int backlegId = aHit->backleg();
904 int moduleId = aHit->module();
905 int cellId = aHit->cell();
906 if(Debug()) {LOG_INFO <<
"A: fired hit in " <<
"backleg=" << backlegId <<
" module="<<moduleId<<
" cell="<<cellId<<endm;}
907 StructCellHit aDaqCellHit;
908 aDaqCellHit.backleg = backlegId;
909 aDaqCellHit.module= moduleId;
910 aDaqCellHit.cell=cellId;
911 aDaqCellHit.tot=aHit->tot();
912 aDaqCellHit.leadingEdgeTime=aHit->leadingEdgeTime();
913 aDaqCellHit.index2MtdHit=i;
914 aDaqCellHit.idTruth = aHit->idTruth();
915 daqCellsHitVec.push_back(aDaqCellHit);
918 int id=backlegId*100+moduleId;
919 if(find(validModuleVec.begin(),validModuleVec.end(),id) == validModuleVec.end())
920 validModuleVec.push_back(
id);
921 int hisIndex = backlegId - 1;
923 mDaqOccupancy[hisIndex]->Fill((moduleId-1)*12+cellId);
927 if(!mEvent||!(mEvent->mtdCollection())||!(mEvent->mtdCollection()->hitsPresent())){
928 if(!mEvent){LOG_INFO <<
"no StEvent" <<endm;
return kFALSE;}
929 else if(!(mEvent->mtdCollection())) {
930 LOG_INFO <<
"no MTD Collection" <<endm;
933 else if(!(mEvent->mtdCollection()->hitsPresent())){LOG_INFO <<
"no MTD hits present" <<endm;
return kFALSE;}
945 StSPtrVecMtdHit& mtdHits= theMtd->mtdHits();
947 for(
size_t i=0;i<mtdHits.size();i++){
950 if(aHit->backleg()<=0||aHit->backleg()>mNBacklegs)
continue;
953 aHit->setAssociatedTrack(0);
955 int backlegId = aHit->backleg();
956 int moduleId = aHit->module();
957 int cellId = aHit->cell();
959 if(Debug()) {LOG_INFO <<
"A: fired hit in " <<
"backleg=" << backlegId <<
" module="<<moduleId<<
" cell="<<cellId<<endm;}
960 LOG_DEBUG <<
"A: fired hit in " <<
"backleg=" << backlegId <<
" module="<<moduleId<<
" cell="<<cellId<<
" leadingWest="<<aHit->leadingEdgeTime().first<<
" leadingEast="<<aHit->leadingEdgeTime().second<<endm;
961 StructCellHit aDaqCellHit;
962 aDaqCellHit.backleg = backlegId;
963 aDaqCellHit.module= moduleId;
964 aDaqCellHit.cell=cellId;
965 aDaqCellHit.tot=aHit->tot();
966 aDaqCellHit.leadingEdgeTime=aHit->leadingEdgeTime();
967 aDaqCellHit.index2MtdHit=i;
968 aDaqCellHit.idTruth = aHit->idTruth();
969 daqCellsHitVec.push_back(aDaqCellHit);
972 int id=backlegId*100+moduleId;
973 if(find(validModuleVec.begin(),validModuleVec.end(),id) == validModuleVec.end())
974 validModuleVec.push_back(
id);
976 int hisIndex = backlegId - 1;
978 mDaqOccupancy[hisIndex]->Fill((moduleId-1)*12+cellId);
987 void StMtdMatchMaker::project2Mtd(mtdCellHitVector daqCellsHitVec,mtdCellHitVector& allCellsHitVec,Int_t& nPrimaryHits){
994 Nnodes = mMuDst->numberOfGlobalTracks();
995 for(UInt_t iNode=0;iNode<Nnodes;iNode++){
997 if(!theTrack)
continue;
1000 bool isPrimary=kFALSE;
1001 Int_t pIndex = -999;
1002 map<Int_t, Int_t>::iterator it = index2Primary.find(iNode);
1003 if(it!=index2Primary.end())
1004 pIndex = it->second;
1005 if(pIndex>=0) isPrimary=kTRUE;
1010 int vtxIndex = ((
StMuTrack *)mMuDst->
array(muPrimary)->UncheckedAt(pIndex))->vertexIndex();
1014 if(vertex) pVtx = vertex->position();
1017 if(matchTrack2Mtd(daqCellsHitVec,theTrack->
outerHelix(),theTrack->
charge(),allCellsHitVec,iNode,pVtx)){
1019 if(isPrimary) nPrimaryHits++;
1024 Nnodes = mEvent->trackNodes().size();
1025 for(UInt_t iNode=0;iNode<Nnodes;iNode++){
1026 StSPtrVecTrackNode& nodes=mEvent->trackNodes();
1028 if(!theTrack)
continue;
1031 StSPtrVecTrackPidTraits& traits = theTrack->pidTraits();
1032 for (StSPtrVecTrackPidTraitsIterator it = traits.begin(); it != traits.end(); it++)
1034 if( (*it)->detector() == kMtdId )
1041 if (mEvent->primaryVertex()) pVtx = mEvent->primaryVertex()->position();
1042 else pVtx.set(0,0,0);
1043 bool isPrimary =kFALSE;
1048 if(pTrack->vertex()) pVtx = pTrack->vertex()->position();
1051 StSPtrVecTrackPidTraits& ptraits = pTrack->pidTraits();
1052 for (StSPtrVecTrackPidTraitsIterator it = ptraits.begin(); it != ptraits.end(); it++)
1054 if( (*it)->detector() == kMtdId )
1062 if(matchTrack2Mtd(daqCellsHitVec,theTrack->outerGeometry()->helix(),theTrack->geometry()->charge(),allCellsHitVec,iNode,pVtx)){
1064 if(isPrimary) nPrimaryHits++;
1074 bool StMtdMatchMaker::matchTrack2Mtd(mtdCellHitVector daqCellsHitVec,
const StPhysicalHelixD &helix, Int_t gq, mtdCellHitVector& allCellsHitVec,
unsigned int iNode,
StThreeVectorD pVtx){
1076 if(mMuDstIn) mField = mMuDst->
event()->runInfo().magneticField();
1077 else mField = mEvent->runInfo()->magneticField();
1083 float mFieldFromGeom = mMtdGeom->GetFieldZ(100,100,0.);
1084 if(fabs(mFieldFromGeom-mField/10.)>0.2){
1085 LOG_WARN<<
"Wrong magnetc field mField = "<<mField/10.<<
" mFieldFromGeom = "<<mFieldFromGeom<<
" check the magF input!!!"<<endm;
1092 mMtdGeom->HelixCrossCellIds(helix,idVec,pathVec,crossVec,tofVec);
1094 mMtdGeom->HelixCrossCellIds(helix,pVtx,idVec,pathVec,crossVec,tofVec);
1097 if((idVec.size()!=pathVec.size()) || idVec.size()!=crossVec.size() || idVec.size()!=tofVec.size()){
1099 LOG_INFO<<
"Inconsistent size idVec = "<<idVec.size()<<
" pathVec = "<<pathVec.size()<<
" crossVec = "<<crossVec.size()<<
" tofVec = "<<tofVec.size()<<endm;
1106 LOG_INFO<<
"dca:"<<dca.mag()<<endm;
1107 LOG_INFO<<
"idVec:"<<idVec.size()<<endm;
1108 LOG_INFO<<
"iNode:"<<iNode<<endm;
1109 if(idVec.size()>1) LOG_INFO<<
"MARK: one track cross two modules!"<<endm;
1111 for (UInt_t i = 0; i < idVec.size(); i++) {
1114 Int_t iCell = -9999;
1115 mMtdGeom->DecodeCellId(idVec[i],iBL,iMod,iCell);
1116 StMtdGeoModule *mMtdGeoModule = mMtdGeom->GetGeoModule(iBL,iMod);
1117 Double_t local[3]={0,0,0};
1118 Double_t global[3]={crossVec[i].x(),crossVec[i].y(),crossVec[i].z()};
1119 if(mMtdGeoModule) mMtdGeoModule->MasterToLocal(global,local);
1122 StructCellHit cellHit;
1123 cellHit.backleg = iBL;
1124 cellHit.module = iMod;
1125 cellHit.cell = iCell;
1126 cellHit.trackIdVec.push_back((Int_t)iNode);
1127 cellHit.hitPosition = crossVec[i];
1128 cellHit.zhit = ol.z();
1129 cellHit.yhit = ol.y();
1130 cellHit.theta = ol.theta();
1131 cellHit.pathLength = pathVec[i];
1132 cellHit.expTof2MTD = tofVec[i];
1133 allCellsHitVec.push_back(cellHit);
1137 if(idVec.size()>1) LOG_INFO<<
"iBL:iMod:iCell="<<iBL<<
" "<<iMod<<
" "<<iCell<<endm;
1144 void StMtdMatchMaker::matchMtdHits(mtdCellHitVector& daqCellsHitVec,mtdCellHitVector& allCellsHitVec,mtdCellHitVector& matchHitCellsVec){
1145 StructCellHit cellHit;
1146 mtdCellHitVectorIter daqIter = daqCellsHitVec.begin();
1148 for(
unsigned int idaq=0;idaq<daqCellsHitVec.size();idaq++, daqIter++) {
1149 mtdCellHitVectorIter proIter = allCellsHitVec.begin();
1150 StMtdGeoModule *geoModule = mMtdGeom->GetGeoModule(daqIter->backleg,daqIter->module);
1151 if(daqIter->backleg==9)
1155 LOG_WARN <<
"Geometry module not available for BL = " << daqIter->backleg <<
" Mod = " << daqIter->module << endm;
1161 for(
unsigned int ipro=0;ipro<allCellsHitVec.size();ipro++, proIter++) {
1163 int daqIndex = (daqIter->module-1)*12 + (daqIter->cell);
1164 int proIndex = (proIter->module-1)*12 + (proIter->cell);
1165 int hisIndex = daqIter->backleg - 1;
1168 double stripPhiCen = 0.;
1169 int channel = daqIter->cell;
1170 stripPhiCen = geoModule->GetCellPhiCenter(channel);
1172 double mLeTimeWest = daqIter->leadingEdgeTime.first;
1173 double mLeTimeEast = daqIter->leadingEdgeTime.second;
1175 double stripZCen = modCen.z() - (mLeTimeWest-mLeTimeEast)/2./gMtdCellDriftV*1000.;
1177 double daqphi = stripPhiCen;
1178 double daqz = stripZCen;
1180 hMtdZvsProj->Fill(proIter->hitPosition.z(),daqz);
1181 hMtdPhivsProj->Fill(proIter->hitPosition.phi(),daqphi);
1182 hMtddPhivsBackleg->Fill(hisIndex,proIter->hitPosition.phi()-daqphi);
1183 hMtddZvsBackleg->Fill(hisIndex,proIter->hitPosition.z()-daqz);
1187 if(daqIter->backleg==proIter->backleg) {
1189 if(hisIndex>=0&&hisIndex<30) {
1190 mHitCorr[hisIndex]->Fill(proIndex,daqIndex);
1191 mHitCorrModule[hisIndex]->Fill(proIter->module-1,daqIter->module-1);
1193 LOG_WARN <<
" weird tray # " << daqIter->module<< endm;
1209 Int_t ibackleg = daqIter->backleg;
1210 Int_t imodule = daqIter->module;
1211 Int_t icell = daqIter->cell;
1213 double zdaq = (daqIter->leadingEdgeTime.second-daqIter->leadingEdgeTime.first)/2./mVDrift[(ibackleg-1)*gMtdNModules+imodule-1][icell]*1e3;
1214 bool isMatch =
false;
1217 if(daqIter->backleg==proIter->backleg){
1218 if(daqIter->module==proIter->module) isMatch =
true;
1222 if(daqIter->backleg==proIter->backleg && daqIter->module!=proIter->module){
1224 if((daqIter->module-1>0)&&daqIter->module-1==proIter->module) isMatch =
true;
1226 if(abs(daqIter->module-proIter->module)<=1) isMatch =
true;
1228 if((daqIter->module+1<6)&&daqIter->module+1==proIter->module) isMatch =
true;
1235 cellHit.backleg = daqIter->backleg;
1236 cellHit.module = daqIter->module;
1237 cellHit.cell = daqIter->cell;
1238 cellHit.hitPosition = proIter->hitPosition;
1239 cellHit.trackIdVec = proIter->trackIdVec;
1240 cellHit.zhit = proIter->zhit;
1241 cellHit.yhit = proIter->yhit;
1242 cellHit.tot = daqIter->tot;
1243 cellHit.leadingEdgeTime = daqIter->leadingEdgeTime;
1244 cellHit.index2MtdHit = daqIter->index2MtdHit;
1245 cellHit.theta = proIter->theta;
1246 cellHit.pathLength = proIter->pathLength;
1247 cellHit.expTof2MTD = proIter->expTof2MTD;
1248 cellHit.idTruth = daqIter->idTruth;
1249 matchHitCellsVec.push_back(cellHit);
1258 void StMtdMatchMaker::sortSingleAndMultiHits(mtdCellHitVector& matchHitCellsVec,mtdCellHitVector& singleHitCellsVec,mtdCellHitVector& multiHitsCellsVec)
1260 StructCellHit cellHit;
1261 StructCellHit cellHit_candidate;
1262 mtdCellHitVector tempVec = matchHitCellsVec;
1263 mtdCellHitVector erasedVec = tempVec;
1264 mtdCellHitVector multiHitsCellsVec_temp;
1265 while (tempVec.size() != 0)
1268 idVector trackIdVec;
1270 mtdCellHitVectorIter tempIter=tempVec.begin();
1271 mtdCellHitVectorIter erasedIter=erasedVec.begin();
1272 while(erasedIter!= erasedVec.end())
1274 if(tempIter->backleg == erasedIter->backleg &&
1275 tempIter->module == erasedIter->module &&
1276 tempIter->cell == erasedIter->cell)
1283 cellHit.cell = tempIter->cell;
1284 cellHit.module = tempIter->module;
1285 cellHit.backleg = tempIter->backleg;
1286 cellHit.hitPosition = tempIter->hitPosition;
1287 cellHit.trackIdVec.push_back(tempIter->trackIdVec.back());
1288 cellHit.zhit = tempIter->zhit;
1289 cellHit.matchFlag = -999;
1290 cellHit.yhit = tempIter->yhit;
1291 cellHit.tot = tempIter->tot;
1292 cellHit.leadingEdgeTime = tempIter->leadingEdgeTime;
1293 cellHit.index2MtdHit = tempIter->index2MtdHit;
1294 cellHit.theta = tempIter->theta;
1295 cellHit.pathLength = tempIter->pathLength;
1296 cellHit.expTof2MTD = tempIter->expTof2MTD;
1297 cellHit.idTruth = tempIter->idTruth;
1300 LOG_INFO<<
"track Info: "<<cellHit.zhit<<
" "<<cellHit.yhit<<endm;
1301 LOG_INFO<<
"hit Info: "<<cellHit.backleg<<
" "<<cellHit.module<<
" "<<cellHit.cell<<endm;
1303 multiHitsCellsVec_temp.push_back(cellHit);
1306 cellHit_candidate.cell = erasedIter->cell;
1307 cellHit_candidate.module = erasedIter->module;
1308 cellHit_candidate.backleg = erasedIter->backleg;
1309 cellHit_candidate.hitPosition = erasedIter->hitPosition;
1310 cellHit_candidate.trackIdVec.push_back(erasedIter->trackIdVec.back());
1311 cellHit_candidate.zhit = erasedIter->zhit;
1312 cellHit_candidate.matchFlag = -999;
1313 cellHit_candidate.yhit = erasedIter->yhit;
1314 cellHit_candidate.tot = erasedIter->tot;
1315 cellHit_candidate.leadingEdgeTime = erasedIter->leadingEdgeTime;
1316 cellHit_candidate.index2MtdHit = erasedIter->index2MtdHit;
1317 cellHit_candidate.theta = erasedIter->theta;
1318 cellHit_candidate.pathLength = erasedIter->pathLength;
1319 cellHit_candidate.expTof2MTD = erasedIter->expTof2MTD;
1320 cellHit_candidate.idTruth = erasedIter->idTruth;
1321 multiHitsCellsVec_temp.push_back(cellHit_candidate);
1324 LOG_INFO<<
"track Info: "<<cellHit_candidate.zhit<<
" "<<cellHit_candidate.yhit<<endm;
1325 LOG_INFO<<
"hit Info: "<<cellHit_candidate.backleg<<
" "<<cellHit_candidate.module<<
" "<<cellHit_candidate.cell<<endm;
1330 trackIdVec.push_back(erasedIter->trackIdVec.back());
1331 erasedVec.erase(erasedIter);
1337 cellHit.cell = tempIter->cell;
1338 cellHit.module = tempIter->module;
1339 cellHit.backleg = tempIter->backleg;
1340 cellHit.hitPosition = tempIter->hitPosition;
1341 cellHit.trackIdVec = trackIdVec;
1342 cellHit.zhit = tempIter->zhit;
1343 cellHit.matchFlag = -999;
1344 cellHit.yhit = tempIter->yhit;
1345 cellHit.tot = tempIter->tot;
1346 cellHit.leadingEdgeTime = tempIter->leadingEdgeTime;
1347 cellHit.index2MtdHit = tempIter->index2MtdHit;
1348 cellHit.theta = tempIter->theta;
1349 cellHit.pathLength = tempIter->pathLength;
1350 cellHit.expTof2MTD = tempIter->expTof2MTD;
1351 cellHit.idTruth = tempIter->idTruth;
1356 cellHit.matchFlag = 1;
1357 singleHitCellsVec.push_back(cellHit);
1362 bool isMatchToSingleTrack = kFALSE;
1363 idVector tmpIdVec = trackIdVec;
1364 sort(tmpIdVec.begin(),tmpIdVec.end());
1365 tmpIdVec.erase(unique(tmpIdVec.begin(),tmpIdVec.end()),tmpIdVec.end());
1366 if(tmpIdVec.size()==1)
1368 isMatchToSingleTrack = kTRUE;
1371 LOG_INFO<<
"A hit is matched to one track that is projected to two modules!"<<endm;
1374 for(
int iTrack = 0; iTrack < nTracks; iTrack++)
1376 if(isMatchToSingleTrack) multiHitsCellsVec_temp[iTrack].matchFlag = 1;
1377 multiHitsCellsVec.push_back(multiHitsCellsVec_temp[iTrack]);
1382 LOG_WARN <<
"D: no tracks extrapolate to matched cell ... should not happen!" << endm;
1387 LOG_DEBUG <<
"D: backleg=" << cellHit.backleg <<
" imodule=" << cellHit.module <<
" icell=" << cellHit.cell <<
"\ttrackid:";
1388 idVectorIter ij=trackIdVec.begin();
1389 while (ij != trackIdVec.end()) { LOG_DEBUG<<
" " << *ij<<endm; ij++; }
1391 tempVec = erasedVec;
1392 multiHitsCellsVec_temp.clear();
1397 LOG_INFO<<
"multiHitsCellsVec size="<<multiHitsCellsVec.size()<<endm;
1402 mtdCellHitVector tempVec_2Trck = multiHitsCellsVec;
1403 mtdCellHitVector erasedVec_2Trck = tempVec_2Trck;
1404 while(tempVec_2Trck.size()!=0)
1406 StructCellHit Cellhit;
1409 vector<StThreeVectorD> vPosition;
1410 vector<Int_t> vchannel, vbackleg, vmodule, vcell;
1411 vector<Int_t> vflag;
1412 vector<Float_t> vzhit, vyhit;
1415 vector<Double_t> vtheta;
1416 vector<Double_t> vpathLength;
1417 vector<Double_t> vexpTof2MTD;
1418 vector<Int_t> vindex2MtdHit;
1419 vector<Int_t> vidTruth;
1421 mtdCellHitVectorIter temp_Iter=tempVec_2Trck.begin();
1422 mtdCellHitVectorIter erased_Iter=erasedVec_2Trck.begin();
1424 while(erased_Iter!= erasedVec_2Trck.end())
1426 if(temp_Iter->backleg == erased_Iter->backleg &&
1427 temp_Iter->module == erased_Iter->module &&
1428 temp_Iter->cell == erased_Iter->cell)
1431 vbackleg.push_back(erased_Iter->backleg);
1432 vmodule.push_back(erased_Iter->module);
1433 vcell.push_back(erased_Iter->cell);
1434 vflag.push_back(erased_Iter->matchFlag);
1435 vPosition.push_back(erased_Iter->hitPosition);
1436 vTrackId.push_back(erased_Iter->trackIdVec.back());
1437 vzhit.push_back(erased_Iter->zhit);
1438 vyhit.push_back(erased_Iter->yhit);
1439 vtot.push_back(erased_Iter->tot);
1440 vtdc.push_back(erased_Iter->leadingEdgeTime);
1441 vindex2MtdHit.push_back(erased_Iter->index2MtdHit);
1442 vtheta.push_back(erased_Iter->theta);
1443 vpathLength.push_back(erased_Iter->pathLength);
1444 vexpTof2MTD.push_back(erased_Iter->expTof2MTD);
1445 vidTruth.push_back(erased_Iter->idTruth);
1448 LOG_INFO<<
"ntracks ="<<ntracks<<endm;
1449 LOG_INFO<<
"tack INFO: "<<vzhit[ntracks-1]<<
" "<<vyhit[ntracks-1]<<endm;
1450 LOG_INFO<<
"hits INFO: "<<vbackleg[ntracks-1]<<
" "<<vmodule[ntracks-1]<<
" "<<vcell[ntracks-1]<<endm;
1452 erasedVec_2Trck.erase(erased_Iter);
1460 LOG_INFO<<
"What? one hit just match with one track in multimatch!!!!!!! "<<endm;
1467 Int_t thiscandidate = -99;
1468 Int_t thisMatchFlag = 0;
1469 Float_t maxR = 9999.;
1471 for(
int j=0; j<ntracks; j++)
1473 Int_t ibackleg = vbackleg[j];
1474 Int_t imodule = vmodule[j];
1475 Int_t icell = vcell[j];
1476 Int_t iidTruth = vidTruth[j];
1478 Float_t trkLocalY = vyhit[j];
1479 Float_t hitLocalY = -999.;
1480 if(mYear<=2016) hitLocalY = mMtdGeom->GetGeoModule(ibackleg,imodule)->GetCellLocalYCenter(icell,ibackleg,iidTruth);
1481 else hitLocalY = mMtdGeom->GetGeoModule(ibackleg,imodule)->GetCellLocalYCenter(icell);
1482 Float_t dy = fabs(trkLocalY-hitLocalY);
1484 Float_t trkGlobalZ = vPosition[j].z();
1485 Float_t hitGlobalZ =
getMtdHitGlobalZ(vtdc[j].first, vtdc[j].second, imodule);
1486 Float_t dz = fabs(trkGlobalZ-hitGlobalZ);
1489 Float_t dr = abs(imodule-projMod)*gMtdRadiusDiff;
1492 if(mCosmicFlag) rr = dy;
1493 else rr = sqrt(dy*dy+dz*dz+dr*dr);
1496 LOG_INFO<<
"hit information::"<<
" "<<ibackleg<<
" "<<imodule<<
" "<<icell<<endm;
1497 LOG_INFO<<
"track information::"<<
" "<<vTrackId[j] <<
" " << trkGlobalZ<<
" "<<trkLocalY<<endm;
1498 LOG_INFO<<
"distance::"<<
" "<<dz<<
" "<<dy<<
" "<<rr<<endm;
1505 thisMatchFlag = (vflag[j]==1? 1 : 7);
1509 if (thiscandidate>=0)
1511 Cellhit.backleg = vbackleg[thiscandidate];
1512 Cellhit.module = vmodule[thiscandidate];
1513 Cellhit.cell = vcell[thiscandidate];
1514 Cellhit.trackIdVec.push_back(vTrackId[thiscandidate]);
1515 Cellhit.hitPosition = vPosition[thiscandidate];
1516 Cellhit.matchFlag = thisMatchFlag;
1517 Cellhit.zhit = vzhit[thiscandidate];
1518 Cellhit.yhit = vyhit[thiscandidate];
1519 Cellhit.tot = vtot[thiscandidate];
1520 Cellhit.leadingEdgeTime = vtdc[thiscandidate];
1521 Cellhit.index2MtdHit = vindex2MtdHit[thiscandidate];
1522 Cellhit.theta = vtheta[thiscandidate];
1523 Cellhit.pathLength = vpathLength[thiscandidate];
1524 Cellhit.expTof2MTD = vexpTof2MTD[thiscandidate];
1525 Cellhit.idTruth = vidTruth[thiscandidate];
1526 singleHitCellsVec.push_back(Cellhit);
1529 LOG_INFO<<
"flag = "<<thisMatchFlag<<endm;
1530 LOG_INFO<<
"Pick track " << vTrackId[thiscandidate] <<
" with local z" << vzhit[thiscandidate] << endm;
1534 tempVec_2Trck = erasedVec_2Trck;
1540 void StMtdMatchMaker::finalMatchedMtdHits(mtdCellHitVector& singleHitCellsVec,mtdCellHitVector& finalMatchedCellsVec){
1541 mtdCellHitVector tempVec = singleHitCellsVec;
1542 mtdCellHitVector erasedVec = tempVec;
1544 while (tempVec.size() != 0)
1546 StructCellHit cellHit;
1549 vector<StThreeVectorD> vPosition;
1550 vector<Int_t> vchannel, vbackleg, vmodule, vcell;
1551 vector<Float_t> vzhit, vyhit;
1554 vector<Double_t> vtheta;
1555 vector<Double_t> vpathLength;
1556 vector<Double_t> vexpTof2MTD;
1557 vector<Int_t> vindex2MtdHit;
1558 vector<Int_t> vflag;
1559 vector<Int_t> vidTruth;
1560 mtdCellHitVectorIter tempIter=tempVec.begin();
1561 mtdCellHitVectorIter erasedIter=erasedVec.begin();
1562 while(erasedIter!= erasedVec.end())
1564 if(tempIter->trackIdVec.back() == erasedIter->trackIdVec.back())
1567 vbackleg.push_back(erasedIter->backleg);
1568 vmodule.push_back(erasedIter->module);
1569 vcell.push_back(erasedIter->cell);
1570 vPosition.push_back(erasedIter->hitPosition);
1571 vTrackId.push_back(erasedIter->trackIdVec.back());
1572 vzhit.push_back(erasedIter->zhit);
1573 vyhit.push_back(erasedIter->yhit);
1574 vtot.push_back(erasedIter->tot);
1575 vtdc.push_back(erasedIter->leadingEdgeTime);
1576 vindex2MtdHit.push_back(erasedIter->index2MtdHit);
1577 vtheta.push_back(erasedIter->theta);
1578 vflag.push_back(erasedIter->matchFlag);
1579 vpathLength.push_back(erasedIter->pathLength);
1580 vexpTof2MTD.push_back(erasedIter->expTof2MTD);
1581 vidTruth.push_back(erasedIter->idTruth);
1584 LOG_INFO<<
"flag 1 ::"<<
" "<<nCells<<
" "<<erasedIter->matchFlag<<endm;
1587 erasedVec.erase(erasedIter);
1596 cellHit.backleg = vbackleg[0];
1597 cellHit.module = vmodule[0];
1598 cellHit.cell = vcell[0];
1599 cellHit.trackIdVec.push_back(vTrackId[0]);
1600 cellHit.hitPosition = vPosition[0];
1601 cellHit.matchFlag = vflag[0];
1602 cellHit.zhit = vzhit[0];
1603 cellHit.yhit = vyhit[0];
1604 cellHit.tot = vtot[0];
1605 cellHit.leadingEdgeTime = vtdc[0];
1606 cellHit.index2MtdHit = vindex2MtdHit[0];
1607 cellHit.theta = vtheta[0];
1608 cellHit.pathLength = vpathLength[0];
1609 cellHit.expTof2MTD = vexpTof2MTD[0];
1610 cellHit.idTruth = vidTruth[0];
1611 finalMatchedCellsVec.push_back(cellHit);
1615 LOG_INFO<<
"flag 2::"<<vflag[0]<<endm;
1621 LOG_DEBUG <<
"E: ibackleg=" << cellHit.backleg <<
" imodule=" << cellHit.module <<
" icell=" << cellHit.cell <<
"\ttrackid:";
1622 idVectorIter ij=vTrackId.begin();
1623 while (ij != vTrackId.end()) { LOG_DEBUG <<
" " << *ij; ij++; }
1630 Int_t thiscandidate = -99;
1631 Int_t thisMatchFlag = 0;
1634 vector<Int_t> ttCandidates;
1635 for (Int_t i=0;i<nCells;i++)
1637 pair<Double_t,Double_t> tt = vtot[i];
1638 if(tt.first<40.) ttCandidates.push_back(i);
1640 if (ttCandidates.size()==1)
1644 LOG_INFO<<
"flag 3 ::"<<vflag[ttCandidates[0]]<<endm;
1646 thiscandidate = ttCandidates[0];
1647 thisMatchFlag = vflag[ttCandidates[0]] + 1;
1649 else if (ttCandidates.size()>1)
1652 vector<Int_t> ssCandidates;
1653 for(
size_t j=0;j<ttCandidates.size();j++)
1655 Int_t ibackleg = vbackleg[ttCandidates[j]];
1656 Int_t imodule = vmodule[ttCandidates[j]];
1657 Int_t icell = vcell[ttCandidates[j]];
1658 Int_t iidTruth = vidTruth[ttCandidates[j]];
1660 Float_t trkLocalY = vyhit[ttCandidates[j]];
1661 Float_t hitLocalY = -999.;
1662 if(mYear<=2016) hitLocalY = mMtdGeom->GetGeoModule(ibackleg,imodule)->GetCellLocalYCenter(icell,ibackleg,iidTruth);
1663 else hitLocalY = mMtdGeom->GetGeoModule(ibackleg,imodule)->GetCellLocalYCenter(icell);
1664 Float_t dy = fabs(trkLocalY-hitLocalY);
1666 Float_t trkGlobalZ = vPosition[ttCandidates[j]].z();
1667 Float_t hitGlobalZ =
getMtdHitGlobalZ(vtdc[ttCandidates[j]].first, vtdc[ttCandidates[j]].second, imodule);
1668 Float_t dz = fabs(trkGlobalZ-hitGlobalZ);
1670 Int_t projMod =
getProjModule(vzhit[ttCandidates[j]],vPosition[ttCandidates[j]].z());
1671 Float_t dr = abs(imodule-projMod)*gMtdRadiusDiff;
1674 if(mCosmicFlag) rr = dy;
1675 else rr = sqrt(dy*dy+dz*dz+dr*dr);
1679 ssCandidates.clear();
1680 ssCandidates.push_back(ttCandidates[j]);
1683 ssCandidates.push_back(ttCandidates[j]);
1685 if (ssCandidates.size()>=1)
1689 LOG_INFO<<
"flag 4::"<<ssCandidates.size()<<
" "<<vflag[ttCandidates[0]]<<endm;
1691 thiscandidate = ssCandidates[0];
1692 thisMatchFlag = vflag[ssCandidates[0]]+2;
1696 if (thiscandidate>=0)
1698 cellHit.backleg = vbackleg[thiscandidate];
1699 cellHit.module = vmodule[thiscandidate];
1700 cellHit.cell = vcell[thiscandidate];
1701 cellHit.trackIdVec.push_back(vTrackId[thiscandidate]);
1702 cellHit.hitPosition = vPosition[thiscandidate];
1703 cellHit.matchFlag = thisMatchFlag;
1704 cellHit.zhit = vzhit[thiscandidate];
1705 cellHit.yhit = vyhit[thiscandidate];
1706 cellHit.tot = vtot[thiscandidate];
1707 cellHit.leadingEdgeTime = vtdc[thiscandidate];
1708 cellHit.index2MtdHit = vindex2MtdHit[thiscandidate];
1709 cellHit.theta = vtheta[thiscandidate];
1710 cellHit.pathLength = vpathLength[thiscandidate];
1711 cellHit.expTof2MTD = vexpTof2MTD[thiscandidate];
1712 cellHit.idTruth = vidTruth[thiscandidate];
1714 finalMatchedCellsVec.push_back(cellHit);
1717 if(Debug()) { LOG_DEBUG <<
"E: ibackleg=" << cellHit.backleg <<
" imodule=" << cellHit.module <<
" icell=" << cellHit.cell <<
"\ttrackid:" << vTrackId[thiscandidate] << endm; }
1722 LOG_WARN <<
"E: no cells belong to this track ... should not happen!" << endm;
1725 tempVec = erasedVec;
1731 void StMtdMatchMaker::fillPidTraits(mtdCellHitVector& finalMatchedCellsVec,Int_t& nValidSingleHitCells,Int_t& nValidSinglePrimHitCells){
1733 for (
size_t ii=0; ii < finalMatchedCellsVec.size(); ii++){
1734 Int_t backleg = finalMatchedCellsVec[ii].backleg;
1735 Int_t module = finalMatchedCellsVec[ii].module;
1736 Int_t cell = finalMatchedCellsVec[ii].cell;
1738 if (finalMatchedCellsVec[ii].trackIdVec.size()!=1)
1739 LOG_WARN <<
"F: WHAT!?! mult.matched cell in single cell list " << backleg <<
" " << module <<
" " << cell << endm;
1741 Float_t trkLocalY = finalMatchedCellsVec[ii].yhit;
1742 Float_t trkLocalZ = finalMatchedCellsVec[ii].zhit;
1743 Float_t trkGlobalZ = finalMatchedCellsVec[ii].hitPosition.z();
1745 Int_t hitIdTruth = finalMatchedCellsVec[ii].idTruth;
1746 Float_t hitLocalY = -999.;
1747 if(mYear<=2016) hitLocalY = mMtdGeom->GetGeoModule(backleg,module)->GetCellLocalYCenter(cell,backleg,hitIdTruth);
1748 else hitLocalY = mMtdGeom->GetGeoModule(backleg,module)->GetCellLocalYCenter(cell);
1749 Float_t LeTimeWest = finalMatchedCellsVec[ii].leadingEdgeTime.first;
1750 Float_t LeTimeEast = finalMatchedCellsVec[ii].leadingEdgeTime.second;
1753 Float_t dy = trkLocalY - hitLocalY;
1754 Float_t dz = trkGlobalZ - hitGlobalZ;
1758 mTracksPerCellMatch3->Fill(finalMatchedCellsVec[ii].trackIdVec.size());
1759 mDeltaHitMatch3->Fill(dy, dz);
1760 int hisIndex = backleg-1;
1761 mDeltaHitFinal[hisIndex]->Fill(dy,dz);
1765 Int_t trackNode = finalMatchedCellsVec[ii].trackIdVec[0];
1768 LOG_DEBUG<<
"In StMuDst mode: mtd hit matched with track successfully : track nodeId:"<<finalMatchedCellsVec[ii].trackIdVec[0]<<
" mtd hitId:"<<finalMatchedCellsVec[ii].index2MtdHit<<endm;
1772 LOG_WARN <<
"Wrong global track!" << endm;
1775 StMuMtdHit *mtdHit = mMuDst->mtdHit(finalMatchedCellsVec[ii].index2MtdHit);
1776 if(mtdHit->backleg()!=backleg || mtdHit->module()!=module || mtdHit->cell()!=cell) {
1777 LOG_WARN <<
"Wrong hit in the MtdHitCollection!" << endm;
1780 nValidSingleHitCells++;
1782 mtdHit->setAssociatedTrackKey(gTrack->
id());
1783 mtdHit->setIndex2Global(trackNode);
1784 gTrack->setIndex2MtdHit(finalMatchedCellsVec[ii].index2MtdHit);
1787 pidMtd.setMatchFlag(finalMatchedCellsVec[ii].matchFlag);
1788 pidMtd.setYLocal(trkLocalY);
1789 pidMtd.setZLocal(trkLocalZ);
1790 pidMtd.setDeltaY(dy);
1791 pidMtd.setDeltaZ(dz);
1792 pidMtd.setThetaLocal(finalMatchedCellsVec[ii].theta);
1793 pidMtd.setPosition(finalMatchedCellsVec[ii].hitPosition);
1794 pidMtd.setPathLength(finalMatchedCellsVec[ii].pathLength);
1795 pidMtd.setExpTimeOfFlight(finalMatchedCellsVec[ii].expTof2MTD);
1796 gTrack->setMtdPidTraits(pidMtd);
1798 Int_t pIndex = -999;
1799 map<Int_t, Int_t>::iterator it = index2Primary.find(trackNode);
1800 if(it!=index2Primary.end()){
1801 pIndex = it->second;
1804 if(pTrack && pIndex>-1){
1805 mtdHit->setIndex2Primary(pIndex);
1806 pTrack->setIndex2MtdHit(finalMatchedCellsVec[ii].index2MtdHit);
1808 ppidMtd.setMatchFlag(finalMatchedCellsVec[ii].matchFlag);
1809 ppidMtd.setYLocal(trkLocalY);
1810 ppidMtd.setZLocal(trkLocalZ);
1811 ppidMtd.setDeltaY(dy);
1812 ppidMtd.setDeltaZ(dz);
1813 ppidMtd.setThetaLocal(finalMatchedCellsVec[ii].theta);
1814 ppidMtd.setPosition(finalMatchedCellsVec[ii].hitPosition);
1815 ppidMtd.setPathLength(finalMatchedCellsVec[ii].pathLength);
1816 ppidMtd.setExpTimeOfFlight(finalMatchedCellsVec[ii].expTof2MTD);
1817 pTrack->setMtdPidTraits(ppidMtd);
1822 LOG_INFO<<
"In StEvent mode: mtd hit matched with track successfully : track nodeId:"<<finalMatchedCellsVec[ii].trackIdVec[0]<<
" mtd hitId:"<<finalMatchedCellsVec[ii].index2MtdHit<<endm;
1825 StSPtrVecTrackNode &nodes = mEvent->trackNodes();
1828 LOG_WARN <<
"Wrong global track!" << endm;
1835 StSPtrVecMtdHit& mtdHits = theMtd->mtdHits();
1836 StMtdHit *mtdHit = mtdHits[finalMatchedCellsVec[ii].index2MtdHit];
1837 if(mtdHit->backleg()!=backleg || mtdHit->module()!=module || mtdHit->cell()!=cell) {
1838 LOG_WARN <<
"Wrong hit in the MtdHitCollection!" << endm;
1841 nValidSingleHitCells++;
1843 mtdHit->setAssociatedTrack(globalTrack);
1846 pidMtd->setMtdHit(mtdHit);
1847 pidMtd->setMatchFlag(finalMatchedCellsVec[ii].matchFlag);
1848 pidMtd->setYLocal(trkLocalY);
1849 pidMtd->setZLocal(trkLocalZ);
1850 pidMtd->setDeltaY(dy);
1851 pidMtd->setDeltaZ(dz);
1852 pidMtd->setThetaLocal(finalMatchedCellsVec[ii].theta);
1853 pidMtd->setPosition(finalMatchedCellsVec[ii].hitPosition);
1854 pidMtd->setPathLength(finalMatchedCellsVec[ii].pathLength);
1855 pidMtd->setExpTimeOfFlight(finalMatchedCellsVec[ii].expTof2MTD);
1856 globalTrack->addPidTraits(pidMtd);
1860 ppidMtd->setMtdHit(mtdHit);
1861 ppidMtd->setMatchFlag(finalMatchedCellsVec[ii].matchFlag);
1862 ppidMtd->setYLocal(trkLocalY);
1863 ppidMtd->setZLocal(trkLocalZ);
1864 ppidMtd->setDeltaY(dy);
1865 ppidMtd->setDeltaZ(dz);
1866 ppidMtd->setThetaLocal(finalMatchedCellsVec[ii].theta);
1867 ppidMtd->setPosition(finalMatchedCellsVec[ii].hitPosition);
1868 ppidMtd->setPathLength(finalMatchedCellsVec[ii].pathLength);
1869 ppidMtd->setExpTimeOfFlight(finalMatchedCellsVec[ii].expTof2MTD);
1870 primaryTrack->addPidTraits(ppidMtd);
1890 float geta = mtt.eta;
1891 int gnFtPts = mtt.nFtPts;
1892 int gnDedxPts = mtt.nDedxPts;
1894 if(geta<mMinEta||geta>
mMaxEta)
return kFALSE;
1895 if(gpt<
mMinPt)
return kFALSE;
1897 if (mtt.flag<=0 || mtt.flag>=1000)
return kFALSE;
1902 if (gnFtPts < mMinFitPointsPerTrack)
return kFALSE;
1906 if(mtt.nHitsPoss!=0) ratio = gnFtPts / (1.0*mtt.nHitsPoss);
1915 return (module-3)*gMtdCellLength + (leadingEastTime - leadingWestTime)/2./gMtdCellDriftV*1e3;
1922 for(
int i=0; i<5; i++)
1924 if(abs(global_z-(i-2)*gMtdCellLength-local_z)<5)
1935 MtdTrack::MtdTrack(
StTrack *stt){
1937 pt = -999.; eta = -999.; nFtPts = 0;
1938 nDedxPts = 0; flag = 0; nHitsPoss = 999;
1940 pt = stt->geometry()->momentum().perp();
1941 eta = stt->geometry()->momentum().pseudoRapidity();
1942 nFtPts = stt->fitTraits().numberOfFitPoints(kTpcId);
1946 if(PidAlgorithm.traits()){
1947 nDedxPts=PidAlgorithm.traits()->numberOfPoints();
1951 nHitsPoss = stt->numberOfPossiblePoints(kTpcId);
1957 pt = -999.; eta = -999.; nFtPts = 0;
1958 nDedxPts = 0; flag = 0; nHitsPoss = 999;
1961 eta = mut->
momentum().pseudoRapidity();
static StMuPrimaryVertex * primaryVertex()
return pointer to current primary vertex
static TObjArray * globalTracks()
returns pointer to the global tracks list
StThreeVectorF primaryVertexPosition(int vtx_id=-1) const
The StMuDst is supposed to be structured in 'physical events'. Therefore there is only 1 primary vert...
Int_t index2Global() const
Returns index of associated global track. If not in order can be set with StMuDst::fixTrackIndeces() ...
short id() const
Returns the track id(or key), is unique for a track node, i.e. global and primary tracks have the sam...
virtual Int_t Init()
initialize drifting velocity and histograms.
Float_t mMinFitPointsOverMax
minimum dE/dx fit points per track
UShort_t nHitsDedx() const
Return number of hits used for dEdx.
virtual void Clear(Option_t *option="")
User defined functions.
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
UShort_t nHitsFit() const
Return total number of hits used in fit.
bool validTrack(StTrack *track)
check track quality
short flag() const
Returns flag, (see StEvent manual for type information)
Short_t charge() const
Returns charge.
Float_t mMinEta
minimum ratio
StMtdMatchMaker(const char *name="MtdMatch")
Default constructor.
UShort_t nHitsPoss() const
Return number of possible hits on track.
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
const StThreeVectorF & momentum() const
Returns 3-momentum at dca to primary vertex.
Int_t getProjModule(Float_t local_z, Float_t global_z)
calculate module of projected position
Float_t mMinPt
maximum pseudorapidity
static TClonesArray * array(int type)
returns pointer to the n-th TClonesArray
Float_t mMaxEta
minimum pseudorapidity
virtual Int_t Make()
associate tracks with mtd hits
Int_t InitRun(int runnumber)
InitRun: initialize geometries (retrieve beam line constraint from database)
virtual Int_t Finish()
write QA histograms
StPhysicalHelixD outerHelix() const
Returns outer helix (last measured point)
Int_t mMindEdxFitPoints
minimum fit points per track
Float_t getMtdHitGlobalZ(Float_t leadingWestTime, Float_t leadingEastTime, Int_t module)
calculate global z of the MTD hits based on its timing information
Int_t FinishRun(int runnumber)
FinishRun: clean up BeamHelix (will be reinitialized at the next initRun)