//*-- Author : Stefano Spataro (spataro@lns.infn.it) // *************************************************************************** //_HADES_CLASS_DESCRIPTION ////////////////////////////////////////////////////////////////////////////// // HPidTofRec // In case of experiments where there exists no start hit and the TOF of the // particles has to be recalculated. This class will reconstruct it, // and the new value (plus the new beta and new mass2), will be written inside HPidTrackData // To mark those objects as changed one can ask to HPidTrackData::getTofRecFlag(algorythm). // The result should be 1 or 2 depending on the method which has been used, -1 if the // reconstructor failed (such as in cases when there are only positive charged tracks and no // leptons), 0 (default value) if the reconstructor did not run. // If the start time reconstruction was done starting from a TOFINO track with multiplicity>1 // then the Flag will be 10 for lepton ID, 20 for pi- ID- // Anyhow the original tof can be retrieved direct HPidHitData::getTof(). // // To setup the constructor call HPidTofRec(Text_t *name,Text_t *title, const Char_t *select, const Char_t *track), // where "select" contains the key words for configuration (just "," separated list) // nospline : switch off spline // nohires : switch off high resolution kickplane // nolowres : switch off low resolution kickplane // nork : switch off Runge Kutta tracking // noreft : switch off Reference Trajectories tracking // norich : switch off rich mode // ntuple : switch on create,fill and write ntuple for debugging // chiin : switch on the chi2>-1 cut for inner MDC segments // chiout : switch on the chi2>-1 cut for outer MDC segments // common : use spline as reconstructor for all the trackings // skip : skip event in DST if start is not reconstructed // sim : switch on simulation mode // nopidparam : do not get rich parameters from HPidAlgStandCutsPar // startrk : catStartHit will be filled from RK // startspline: catStartHit will be filled from SPLINE // // if common mode is switched on, "track" is the chosen tracking algorythm for global tof reconstruction: // spline : spline // lowres : low resolution kickplane // hires : high resolution kickplane // rk : Runge Kutta tracking // reft : Reference Trajectories tracking // // defaults (empty string) are: // spline ON // low resolution kickplane ON // high resolution kickplane ON // RungeKutta ON // Reference Trajectories ON // fRichMode ON // ntuple OFF // chi cut (inner MDC) mode OFF // chi cut (outer MDC) mode OFF // skip event OFF // common OFF // simulation OFF // pidparam ON // startrk ON // startspline OFF // // Usage: // // HPidTofRec *pidtofrec = new HPidTofRec("tofreconstructor","tofreconstructor",option,track); // pidtofrec->setQualityRichParams(200.,4.,2.8,5); // HTaskSet * pidtofrectask = new HTaskSet("",""); // pidtofrectask->connect(pidtofrec); // masterTaskSet->add(pidtofrec); // ////////////////////////////////////////////////////////////////////////////// // *************************************************************************** #include "hpidtofrec.h" #include "hades.h" #include "hruntimedb.h" #include "hevent.h" #include "hcategory.h" #include "hmatrixcategory.h" #include "hspectrometer.h" #include "hdetector.h" #include "TMath.h" #include "hpidtrackcandsim.h" #include "hpidgeanttrackset.h" #include "hpidtrackdata.h" #include "hpidhitdata.h" #include "hpidphysicsconstants.h" #include "hpidalgstandcutspar.h" #include "hstarthit.h" #include "hlocation.h" #include "hlinearcategory.h" #include "hstartdef.h" #include "hgeantdef.h" #include "showertofinodef.h" #include "tofinodef.h" #include "hshowerhittof.h" #include "htofinocalsim.h" #include ClassImp(HPidTofRec) HPidTofRec::HPidTofRec() { // Default constructor richParam = NULL; pidalg = -1111; fRichQparams[0] = 0.; fRichQparams[1] = 0.; fRichQparams[2] = 1000.; fRichQparams[3] = 0.; nCommonFlag = -1; nStartAlg = ALG_SPLINE; for (Short_t nSs = 0; nSs < 6; nSs ++) { nMdcSeg[nSs] = 0; for (Short_t nMm = 0; nMm < 4; nMm ++) { if (gHades->getSetup()->getDetector("Mdc")->getModule(nSs,nMm)) { nMdcSeg[nSs] ++; } } } hit = 0; richcorr = 0; bRichMode = kTRUE; bMinusMode = kTRUE; bFillNTuple = kFALSE; bChiCutInner = kFALSE; bChiCutOuter = kFALSE; bSkip = kFALSE; bIsCommon = kFALSE; bIsSimulation = kFALSE; bPidParam = kTRUE; bSetStartHit = kFALSE; bReset = kFALSE; for (Int_t nAlg = 0; nAlg < MAX_MOM_ALGS; nAlg ++) { // switch on all the tracking algorythms bTrackAlg[nAlg] = kTRUE; } iterTrackCand = 0; } HPidTofRec::HPidTofRec(const Text_t *name,const Text_t *title, const Char_t *select, const Char_t *track) :HReconstructor(name,title) { // Constructor TString sel = select; sel.ToLower(); TString tracking = track; tracking.ToLower(); bRichMode = kTRUE; bMinusMode = kTRUE; bFillNTuple = kFALSE; bChiCutInner = kFALSE; bChiCutOuter = kFALSE; bSkip = kFALSE; bIsCommon = kFALSE; bIsSimulation = kFALSE; bPidParam = kTRUE; bSetStartHit = kFALSE; bReset = kFALSE; for (Int_t nAlg = 0; nAlg < MAX_MOM_ALGS; nAlg ++) { // switch on all the tracking algorythms bTrackAlg[nAlg] = kTRUE; } richParam = NULL; pidalg = -1111; fRichQparams[0] = 0.; fRichQparams[1] = 0.; fRichQparams[2] = 1000.; fRichQparams[3] = 0.; nCommonFlag = -1; nStartAlg = ALG_SPLINE; if (sel.Contains("nospline")) bTrackAlg[ALG_SPLINE] = kFALSE; if (sel.Contains("nohires")) bTrackAlg[ALG_KICK123] = kFALSE; if (sel.Contains("nolowres")) bTrackAlg[ALG_KICK] = kFALSE; if (sel.Contains("nork")) bTrackAlg[ALG_RUNGEKUTTA] = kFALSE; if (sel.Contains("noreft")) bTrackAlg[ALG_REFT] = kFALSE; if (sel.Contains("norich")) bRichMode = kFALSE; if (sel.Contains("ntuple")) bFillNTuple = kTRUE; if (sel.Contains("chiin")) bChiCutInner = kTRUE; if (sel.Contains("chiout")) bChiCutOuter = kTRUE; if (sel.Contains("skip")) bSkip = kTRUE; if (sel.Contains("sim")) bIsSimulation = kTRUE; if (sel.Contains("nopidparam")) bPidParam = kFALSE; if (sel.Contains("startrk")) { nStartAlg = ALG_RUNGEKUTTA; bSetStartHit = kTRUE; } if (sel.Contains("startspline")) { nStartAlg = ALG_SPLINE; bSetStartHit = kTRUE; } if (sel.Contains("common")) { if (tracking.Contains("spline")) nCommonFlag = ALG_SPLINE; if (tracking.Contains("lowres")) nCommonFlag = ALG_KICK; if (tracking.Contains("hires")) nCommonFlag = ALG_KICK123; if (tracking.Contains("rk")) nCommonFlag = ALG_RUNGEKUTTA; if (tracking.Contains("reft")) nCommonFlag = ALG_REFT; if (nCommonFlag!=-1) { bIsCommon = kTRUE; } else Warning("init","HPidTofRec: No common tracking setted -> Switching Common Mode OFF"); } for (Short_t nSs = 0; nSs < 6; nSs ++) { nMdcSeg[nSs] = 0; for (Short_t nMm = 0; nMm < 4; nMm ++) { if (gHades->getSetup()->getDetector("Mdc")->getModule(nSs,nMm)) nMdcSeg[nSs] ++; } } hit = 0; richcorr = 0; sDir = "./"; sFile = ""; iterTrackCand = 0; } HPidTofRec::~HPidTofRec(void){ // Destructor. if( iterTrackCand ) { delete iterTrackCand; iterTrackCand = 0; } } Bool_t HPidTofRec::init(void){ Info("init()", "PID Tof Reconstructor default initialization *** "); if(gHades->getCurrentEvent()->getCategory(catGeantKine) && gHades->getEmbeddingMode() > 0) { isEmbedding = kTRUE; } else { isEmbedding = kFALSE; } pShowerHitTofCat = 0; pTofinoCalCat = 0; if(isEmbedding){ pShowerHitTofCat = (HCategory*)(gHades->getCurrentEvent()->getCategory(catShowerHitTofTrack)); if(!pShowerHitTofCat) Warning("init","No HShowerHitTof Input -> Switching usage of this category OFF"); pTofinoCalCat = (HMatrixCategory*)(gHades->getCurrentEvent()->getCategory(catTofinoCal)); if(!pTofinoCalCat) Warning("init","No HTofinoCal Input -> Switching usage of this category OFF"); } pTrackCand = gHades->getCurrentEvent()->getCategory(catPidTrackCand); if (!pTrackCand) { Error("init","HPidTofRec: No HPIDTrackCand Input -> Switching HPidTofRec OFF"); } else { iterTrackCand= (HIterator *)pTrackCand->MakeIterator(); } if(bSetStartHit) { pStartHitCat = (HLinearCategory *)(gHades->getCurrentEvent()->getCategory(catStartHit)); if (!pStartHitCat) { pStartHitCat = new HLinearCategory("HStartHit"); if (!pStartHitCat) { Error("init()","Could not create catStartHit!"); return kFALSE; } else { gHades->getCurrentEvent()->addCategory(catStartHit,pStartHitCat,"Start"); } } } cout<<"Spline: " << bTrackAlg[ALG_SPLINE] <<" Kick lowres: "<< bTrackAlg[ALG_KICK] <<" Kick hires: " << bTrackAlg[ALG_KICK123] <<" RungeKutta: " << bTrackAlg[ALG_RUNGEKUTTA] <<" RefTraj: " << bTrackAlg[ALG_REFT] <getRuntimeDb() ->getContainer(PIDALGSTANDCUTSPAR_NAME)) == NULL) { Error("HPidTofRec::init()", "Cannot get parameters: %s", PIDALGSTANDCUTSPAR_NAME); return kFALSE; } if (bTrackAlg[ALG_RUNGEKUTTA]) { pidalg = ALG_RUNGEKUTTA; } else if (bTrackAlg[ALG_SPLINE]) { pidalg = ALG_SPLINE; } else if (bTrackAlg[ALG_REFT]) { pidalg = ALG_REFT; } else if (bTrackAlg[ALG_KICK123]) { pidalg = ALG_KICK123; } else if (bTrackAlg[ALG_KICK]) { pidalg = ALG_KICK; } if (-1111 == pidalg) { Error("HPidTofRec::init()","no valid momentum algorithm activated"); return kFALSE; } richParam->setContext(pidalg); } else { // ! bPidParam cout<<"Rich Quality parameters: RingPatMat "<< fRichQparams[0] <<" RingPadNr " << fRichQparams[1] <<" Centroid " << fRichQparams[2] <<" RingAverageCharge " << fRichQparams[3] < return kFALSE (analysis can not be done in such case) //and print error message if (bPidParam) { return richParam->checkContext(pidalg); } else { return kTRUE; } } Int_t HPidTofRec::execute(void){ // Execute function //------------------------------------------------------------------------- // create the ntuples in the first event loop if (bIsFirst && bFillNTuple) { TString dstFile = gHades->getOutputFile()->GetName(); if (sFile == "") { sFile = sDir + "/pidtofrec_" + dstFile(dstFile.Last('/') + 1,dstFile.Last('.') - dstFile.Last('/') - 1) + ".root"; }else { sFile = sDir + "/" + sFile; } cout << " *** HPidTofRec: Ntuple output file:\t" << sFile << "\t ***" << endl; bIsFirst = kFALSE; } //------------------------------------------------------------------------- //------------------------------------------------------------------------- // create helper arrays and init them Bool_t bIsRichRec [2][MAX_MOM_ALGS]; Bool_t bIsMinusRec[2][MAX_MOM_ALGS]; Float_t fStartTime [2][MAX_MOM_ALGS]; Float_t fRecTof [2][MAX_MOM_ALGS]; Short_t nRecTofMult[2][MAX_MOM_ALGS]; Float_t fRecStart [2][MAX_MOM_ALGS]; Short_t nIndex [2][MAX_MOM_ALGS]; Short_t nRecSys [2][MAX_MOM_ALGS]; Short_t nRecPol [2][MAX_MOM_ALGS]; Short_t nRecSec [2][MAX_MOM_ALGS]; Short_t nRecIndex [2][MAX_MOM_ALGS]; Bool_t bIsRec [2]; //------------------------------------------------------------------------- //------------------------------------------------------------------------- // find a good particle to recalulate the start time. // In embedding mode this loop runs twice, once for real // and once for tracks with sim TOF measurements. The // start time correction will be done for sim /real tracks // separately. This is not fully consistent, but the only way // to mimic somehow the resolution effect for the simulated // tracks. Bool_t bIsRealTrack = kTRUE; for(Int_t n = 0; n < 2; n ++) { bIsRec[n] = kFALSE; for(Int_t nAlg = 0; nAlg < MAX_MOM_ALGS; nAlg ++) { bIsRichRec [n][nAlg] = kFALSE; bIsMinusRec[n][nAlg] = kFALSE; nIndex [n][nAlg] = 0; fStartTime [n][nAlg] = 0.; nRecSys [n][nAlg] = -10; nRecPol [n][nAlg] = 0; nRecSec [n][nAlg] = -1; fRecTof [n][nAlg] = 0.; nRecTofMult[n][nAlg] = -1; fRecStart [n][nAlg] = 0.; nRecIndex [n][nAlg] = -1; } } HPidTrackCandSim *pCand = 0; Int_t nMax = 1; if(isEmbedding) nMax = 2; // loop n==0 -> real tracks , n == 1 -> sim tracks for(Int_t n = 0; n < nMax; n ++ ) { iterTrackCand->Reset(); while ((pCand = (HPidTrackCandSim *)iterTrackCand->Next()) != 0) // begin of PidTrackCand iterator { //------------------------------------------------------------------------- // in embedding mode we have to take care if we // use real or embedded tracks for the recalibration // of the start time. if(isEmbedding) { bIsRealTrack = isRealTrack(pCand); if(n == 0) { // real tracks if(!bIsRealTrack){ continue; } } else { // sim tracks if( bIsRealTrack){ continue; } } } //------------------------------------------------------------------------- for (Int_t nAlg = 0; nAlg < MAX_MOM_ALGS; nAlg ++) // loop over tracking algorythms { if(bReset){ //------------------------------------------------------------------------- // reset to old values pCand->getTrackData()->nTofRecFlag [nAlg] = 0; pCand->getTrackData()->fTofRecTof [nAlg] = pCand->getHitData()->getTof(); pCand->getTrackData()->fTofRecBeta [nAlg] = pCand->getTrackData()->fCorrectedBeta[nAlg]; pCand->getTrackData()->fTofRecMassSquared[nAlg] = pCand->getTrackData()->fMassSquared [nAlg]; //------------------------------------------------------------------------- } if (bTrackAlg[nAlg]) // selection of algorythm activated for reconstruction { nIndex[n][nAlg] ++; pCand->getTrackData()->nTofRecFlag[nAlg] = -1; // reset HPidTofRec flag if ((pCand->getTrackData()->bIsAccepted[nAlg]) // check on isXXXAccepted() && !((bChiCutInner) && // check on chi for inner mdc segment (pCand->getHitData()->fInnerMdcChiSquare == -1) ) && !((bChiCutOuter) && // check on chi for outer mdc segment (nAlg != ALG_KICK) && (pCand->getHitData()->fOuterMdcChiSquare == -1) ) && !((nAlg == ALG_SPLINE) && // skip spline tracks with Qspline==-1 (pCand->getTrackData()->fSplineChiSquare == -1) ) && ( pCand->getHitData()->iSystem != -1) // skip tracks without META correlation /* // commented line, to use TofinoMult>1 && !((pCand->getHitData()->iSystem == 0) && // skip system==0 && TofinoMult!=1 (pCand->getHitData()->iTofinoMult != 1) ) */ ) { Short_t tofMult =-1; if (pCand->getHitData()->iSystem==0) tofMult = pCand->getHitData()->getTofinoMult(); else tofMult = 1; if ((bRichMode) && (pCand->getHitData()->hasRingCorrelation[nAlg])) { // check of RICH correlation if (bFillNTuple) // fill RICH debug histos { Float_t richcorr_ntuple[17] = { pCand->getHitData() ->fMdcTheta, pCand->getHitData() ->fRichTheta, pCand->getHitData() ->fMdcPhi, pCand->getHitData() ->fRichPhi, pCand->getHitData() ->nSector, pCand->getHitData() ->iSystem, pCand->getTrackData()->nPolarity[nAlg], pCand->getTrackData()->fMomenta[nAlg], nAlg, pCand->getHitData() ->fInnerMdcChiSquare, pCand->getHitData() ->fOuterMdcChiSquare, pCand->getHitData() ->fMdcZCoord, pCand->getHitData() ->fMdcRCoord, pCand->getHitData() ->nRingPatMat, pCand->getHitData() ->nRingPadNr, pCand->getHitData() ->fRingCentroid, pCand->getHitData() ->nRingAmplitude }; richcorr->Fill(richcorr_ntuple); } // enf of RICH debug histos Bool_t bGoodRing = kFALSE; if (bPidParam) { Int_t mdc_sec = pCand->getHitData()->nSector; if (// positrons ( (pCand->getHitData()->nRingPatMat > richParam->getValue(2, 0, mdc_sec, 0)) && (pCand->getHitData()->fRingCentroid < richParam->getValue(2, 0, mdc_sec, 1)) && (pCand->getHitData()->nRingPadNr >= richParam->getValue(2, 0, mdc_sec, 2)) && (pCand->getHitData()->nRingAmplitude / pCand->getHitData()->nRingPadNr >= richParam->getValue(2, 0, mdc_sec, 3)) ) #ifdef THIS_IS_BEAUTIFUL_BUT_OBSOLETE // positrons and electrons have same thr || // electrons ( (pCand->getHitData()->nRingPatMat > richParam->getValue(3, 0, mdc_sec, 0)) && (pCand->getHitData()->fRingCentroid < richParam->getValue(3, 0, mdc_sec, 1)) && (pCand->getHitData()->nRingPadNr >= richParam->getValue(3, 0, mdc_sec, 2)) && (pCand->getHitData()->nRingAmplitude / pCand->getHitData()->nRingPadNr >= richParam->getValue(3, 0, mdc_sec, 3)) ) #endif // THIS_IS_BEAUTIFUL_BUT_OBSOLETE ) { bGoodRing = kTRUE; } } else { // ! PidParam if ( (pCand->getHitData()->nRingPatMat > fRichQparams[0]) && // quality parameters (pCand->getHitData()->nRingPadNr >= fRichQparams[1]) && (pCand->getHitData()->fRingCentroid < fRichQparams[2]) && ((pCand->getHitData()->nRingAmplitude / pCand->getHitData()->nRingPadNr) >= fRichQparams[3]) ) { bGoodRing = kTRUE; } } // ! PidParam if (bGoodRing) { if ( (!bIsRichRec[n][nAlg]) // if there are no other rich correlated tracks || ( bIsRichRec[n][nAlg] && // or if new one has lower pad multiplicity (tofMult < nRecTofMult[n][nAlg] ) // priority on tracks with lower pad multiplicity ) || ( bIsRichRec[n][nAlg] && // priority on sectors with more MDC chambers (tofMult == nRecTofMult[n][nAlg] ) && (nMdcSeg[pCand->getHitData()->nSector] > nMdcSeg[nRecSec[n][nAlg]]) ) || ( bIsRichRec[n][nAlg] && // or if new one TOF, old one TOFINO (tofMult == nRecTofMult[n][nAlg] ) && (nMdcSeg[pCand->getHitData()->nSector] == nMdcSeg[nRecSec[n][nAlg]]) && (pCand->getHitData()->iSystem == 1 && nRecSys[n][nAlg] == 0 ) ) || ( bIsRichRec[n][nAlg] && // or if new one pol==1 and old pol ==1 (tofMult == nRecTofMult[n][nAlg] ) && (nMdcSeg[pCand->getHitData()->nSector] == nMdcSeg[nRecSec[n][nAlg]]) && (pCand->getHitData()->iSystem == nRecSys[n][nAlg] ) && (pCand->getTrackData()->nPolarity[nAlg] == -1 && nRecPol[n][nAlg] == 1 ) ) || ( bIsRichRec[n][nAlg] && // faster particle (tofMult == nRecTofMult[n][nAlg] ) && (nMdcSeg[pCand->getHitData()->nSector] == nMdcSeg[nRecSec[n][nAlg]]) && (pCand->getHitData()->iSystem == nRecSys[n][nAlg]) && (pCand->getTrackData()->nPolarity[nAlg] == nRecPol[n][nAlg]) && (pCand->getHitData()->getTof() < fRecTof[n][nAlg]) ) ) { bIsRec [n] = kTRUE; bIsRichRec [n][nAlg] = kTRUE; nRecIndex [n][nAlg] = nIndex[n][nAlg]; fRecTof [n][nAlg] = pCand->getHitData()->getTof(); nRecSys [n][nAlg] = pCand->getHitData()->iSystem; nRecTofMult[n][nAlg] = tofMult; nRecPol [n][nAlg] = pCand->getTrackData()->nPolarity[nAlg]; nRecSec [n][nAlg] = pCand->getHitData()->nSector; fRecStart [n][nAlg] = pCand->getHitData()->getTof() * pCand->getTrackData()->fCorrectedBeta[nAlg] * TMath::C() * 1e-6 * sqrt(pCand->getTrackData()->fMomenta[nAlg] * pCand->getTrackData()->fMomenta[nAlg] + HPidPhysicsConstants::mass(3) * HPidPhysicsConstants::mass(3)) / (TMath::C() * 1e-6 * pCand->getTrackData()->fMomenta[nAlg]) - pCand->getHitData()->getTof(); } } // end of lepton recognition } // end of rich loop if (bMinusMode && (!bIsRichRec[n][nAlg]) && (pCand->getTrackData()->nPolarity[nAlg] == -1)) { if ( ( !bIsMinusRec[n][nAlg]) // if there are no other negative charged tracks || ( bIsRichRec[n][nAlg] && // or if new one has lower pad multiplicity (tofMult < nRecTofMult[n][nAlg] ) // priority on tracks with lower pad multiplicity ) || ( bIsMinusRec[n][nAlg] && // priority on sectors with more MDC chambers (tofMult == nRecTofMult[n][nAlg] ) && (nMdcSeg[pCand->getHitData()->nSector] > nMdcSeg[nRecSec[n][nAlg]])) || ( bIsMinusRec[n][nAlg] && // or if new one TOF, old one TOFINO (tofMult == nRecTofMult[n][nAlg] ) && (nMdcSeg[pCand->getHitData()->nSector] == nMdcSeg[nRecSec[n][nAlg]]) && (pCand->getHitData()->iSystem == 1 && nRecSys[n][nAlg] == 0) ) || (bIsMinusRec[n][nAlg] && // faster particle (tofMult == nRecTofMult[n][nAlg] ) && (nMdcSeg[pCand->getHitData()->nSector] == nMdcSeg[nRecSec[n][nAlg]]) && (pCand->getHitData()->iSystem == nRecSys[n][nAlg]) && (pCand->getHitData()->getTof() < fRecTof[n][nAlg])) ) { bIsRec [n] = kTRUE; bIsMinusRec[n][nAlg] = kTRUE; nRecIndex [n][nAlg] = nIndex[n][nAlg]; fRecTof [n][nAlg] = pCand->getHitData()->getTof(); nRecSys [n][nAlg] = pCand->getHitData()->iSystem; nRecTofMult[n][nAlg] = tofMult; nRecPol [n][nAlg] = pCand->getTrackData()->nPolarity[nAlg]; nRecSec [n][nAlg] = pCand->getHitData()->nSector; fRecStart [n][nAlg] = pCand->getHitData()->getTof() * pCand->getTrackData()->fCorrectedBeta[nAlg] * TMath::C() * 1e-6 * sqrt(pCand->getTrackData()->fMomenta[nAlg] * pCand->getTrackData()->fMomenta[nAlg] + HPidPhysicsConstants::mass(8) * HPidPhysicsConstants::mass(8)) / (TMath::C() * 1e-6 * pCand->getTrackData()->fMomenta[nAlg]) - pCand->getHitData()->getTof(); } } // end of pi- recognition } // end of "quality" check of the track } // end of loop over different tracking } } // end of HPidTrackCand iterator } //------------------------------------------------------------------------- //------------------------------------------------------------------------- Short_t nRecAlg = 0; // loop n==0 -> real tracks , n == 1 -> sim tracks for(Int_t n = 0; n < nMax; n ++) { // time/beta/mass2 recalculation for (Int_t nAlg = 0; nAlg < MAX_MOM_ALGS; nAlg ++) { nIndex[n][nAlg] = 0; } iterTrackCand->Reset(); while ((pCand = (HPidTrackCandSim *)iterTrackCand->Next()) != 0) // begin of PidTrackCand iterator { //------------------------------------------------------------------------- // in embedding mode we have to take care if we // use real or embedded tracks for the recalibration // of the start time. if(isEmbedding) { bIsRealTrack = isRealTrack(pCand); if(n == 0) { // real tracks if(!bIsRealTrack){ continue; } } else { // sim tracks if( bIsRealTrack){ continue; } } } //------------------------------------------------------------------------- for (Int_t nAlg = 0; nAlg < MAX_MOM_ALGS; nAlg ++) // loop over tracking algorythms { if (bTrackAlg[nAlg]) // selection of algorythm activated for reconstruction, plus isXXXAccepted { nIndex[n][nAlg] ++; if ((pCand->getTrackData()->bIsAccepted[nAlg]) && // check on isXXXAccepted() (pCand->getHitData()->iSystem != -1) && (bIsRec[n]) ) { if (bIsCommon) nRecAlg = nCommonFlag; else nRecAlg = nAlg; if (nRecTofMult[n][nRecAlg] ==1) { if ( bIsRichRec[n][nRecAlg]) pCand->getTrackData()->nTofRecFlag[nAlg] = 1; if ((!bIsRichRec[n][nRecAlg]) && (bIsMinusRec[n][nRecAlg])) pCand->getTrackData()->nTofRecFlag[nAlg] = 2; } else if (nRecTofMult[n][nRecAlg] >1) { if ( bIsRichRec[n][nRecAlg]) pCand->getTrackData()->nTofRecFlag[nAlg] = 10; if ((!bIsRichRec[n][nRecAlg]) && (bIsMinusRec[n][nRecAlg])) pCand->getTrackData()->nTofRecFlag[nAlg] = 20; } if (pCand->getTrackData()->nTofRecFlag[nAlg] > 0) // if the start time was reconstructed { pCand->getTrackData()->fTofRecTof[nAlg] = pCand->getHitData()->getTof() + fRecStart[n][nRecAlg]; pCand->getTrackData()->fTofRecBeta[nAlg] = pCand->getHitData()->getTof() * pCand->getTrackData()->fCorrectedBeta[nAlg] / pCand->getTrackData()->fTofRecTof[nAlg]; pCand->getTrackData()->fTofRecMassSquared[nAlg] = pCand->getTrackData()->fMomenta[nAlg] * pCand->getTrackData()->fMomenta[nAlg] * (1. / pCand->getTrackData()->fTofRecBeta[nAlg] / pCand->getTrackData()->fTofRecBeta[nAlg] - 1); } if (bFillNTuple) // fill debug ntuples { if (bIsSimulation) { Float_t hit_ntuple[29] = { pCand->getTrackData()->fMomenta[nAlg], pCand->getTrackData()->nPolarity[nAlg], pCand->getHitData()->nSector, pCand->getHitData()->iSystem, pCand->getHitData()->fMdcTheta, pCand->getHitData()->fMdcPhi, pCand->getTrackData()->nTofRecFlag[nAlg], nAlg, pCand->getHitData()->getTof(), pCand->getTrackData()->fTofRecTof[nAlg], pCand->getTrackData()->fCorrectedBeta[nAlg], pCand->getTrackData()->fTofRecBeta[nAlg], pCand->getTrackData()->fMassSquared[nAlg], pCand->getTrackData()->fTofRecMassSquared[nAlg], nRecSys[n][nAlg], nRecPol[n][nAlg], nRecSec[n][nAlg], pCand->getHitData()->fInnerMdcChiSquare, pCand->getHitData()->fOuterMdcChiSquare, pCand->getHitData()->fMdcZCoord, pCand->getHitData()->fMdcRCoord, pCand->getHitData()->iTofinoMult, (Float_t)(nRecIndex[n][nAlg] == nIndex[n][nAlg]), pCand->getHitData()->hasRingCorrelation[nAlg], pCand->getTrackData()->fSplineChiSquare, pCand->getGeantTrackSet()->getCorrelationFlag(), pCand->getGeantTrackSet()->getTotalMomentum(), pCand->getGeantTrackSet()->getGeantPID(), pCand->getGeantTrackSet()->getGeantParentID() }; hit->Fill(hit_ntuple); } // end of ntuple filling sim else { Float_t hit_ntuple[25] = { pCand->getTrackData()->fMomenta[nAlg], pCand->getTrackData()->nPolarity[nAlg], pCand->getHitData()->nSector, pCand->getHitData()->iSystem, pCand->getHitData()->fMdcTheta, pCand->getHitData()->fMdcPhi, pCand->getTrackData()->nTofRecFlag[nAlg], nAlg, pCand->getHitData()->getTof(), pCand->getTrackData()->fTofRecTof[nAlg], pCand->getTrackData()->fCorrectedBeta[nAlg], pCand->getTrackData()->fTofRecBeta[nAlg], pCand->getTrackData()->fMassSquared[nAlg], pCand->getTrackData()->fTofRecMassSquared[nAlg], nRecSys[n][nAlg], nRecPol[n][nAlg], nRecSec[n][nAlg], pCand->getHitData()->fInnerMdcChiSquare, pCand->getHitData()->fOuterMdcChiSquare, pCand->getHitData()->fMdcZCoord, pCand->getHitData()->fMdcRCoord, pCand->getHitData()->iTofinoMult, (Float_t)(nRecIndex[n][nAlg] == nIndex[n][nAlg]), pCand->getHitData()->hasRingCorrelation[nAlg], pCand->getTrackData()->fSplineChiSquare }; hit->Fill(hit_ntuple); } // end of ntuple filling exp } // end of ntuple filling } // cut on isXXXAccepted() sys!=-1 and tof recalcutation } // end of algorythm loop } }// end of loop on HPidTrackCand } //------------------------------------------------------------------------- //------------------------------------------------------------------------- // skip non recalibrated events if wanted if ( (!bIsRec[0] && !bIsRec[1]) && bSkip) { return kSkipEvent; } //------------------------------------------------------------------------- //------------------------------------------------------------------------- // fill the recalibrated start time if wanted if(bSetStartHit) { HLocation loc; loc.set(2,0,0); HStartHit* pHit = (HStartHit*)pStartHitCat->getObject(loc); if(!pHit){ pHit = (HStartHit*)pStartHitCat->getNewSlot(loc); } if(pHit){ pHit = new (pHit) HStartHit; pHit->setModule(0); pHit->setStrip(0); if(bIsRec[0]){ pHit->setTime(-fRecStart[0][nStartAlg]); } else { pHit->setTime(0); } pHit->setStartMult(1); pHit->setFlag(kTRUE); } else { Error("execute()","NULL pointer for HStartHit slot retrieved!"); } } //------------------------------------------------------------------------- return 0; } Bool_t HPidTofRec::finalize(void){ if(bFillNTuple) { TFile *r = TFile::Open(sFile,"RECREATE"); hit ->Write(); richcorr->Write(); r->Save(); r->Close(); r->Delete(); hit ->Delete(); richcorr->Delete(); } return kTRUE; } Bool_t HPidTofRec::isRealTrack(HPidTrackCand* pCand) { // checks by the TOF measurement of the // HPidTrackCand object is related to a real track // or embedded track. Return kTRUE if the track // is real. //------------------------------------------------------------------------- if(isEmbedding) { HPidTrackCandSim* pCandSim = (HPidTrackCandSim*) pCand; HPidGeantTrackSet* pGeantTrackSet = pCandSim->getGeantTrackSet(); if(pCand->getHitData()->iSystem == 1) { // TOF hit for(UInt_t i = 0; i < pGeantTrackSet->getNTOFTracks() ; i ++){ if(pGeantTrackSet->getTOFTrack(i) > 0) { return kFALSE; } } } if(pCand->getHitData()->iSystem == 0) { //------------------------------------------------------------------------- // for the TOFINO we have to do some effort // because the TOFINO track number is not used by // HPidTrackCand if(pShowerHitTofCat && pTofinoCalCat) { // TOFINO hit HShowerHitTof* pShowerhittof = (HShowerHitTof*)pShowerHitTofCat->getObject(pCand->getHitData()->getIndShower()); if(pShowerhittof) { HLocation loc; loc.set(3,pShowerhittof->getSector(),0,pShowerhittof->getTofinoCell()); HTofinoCalSim* ptofinocalsim = 0; if(loc[1] >= 0) { ptofinocalsim = (HTofinoCalSim*) pTofinoCalCat->getObject(loc); if(ptofinocalsim) { if(ptofinocalsim->getTrack(0) > 0) { return kFALSE; } } else { Error("execute()","NULL pointer received for TofinoCal object!"); } } } else { Error("execute()","NULL pointer received for ShowerHitTof object with index = %i!",pCand->getHitData()->getIndShower()); } } } //------------------------------------------------------------------------- } //------------------------------------------------------------------------- return kTRUE; }