//_HADES_CLASS_DESCRIPTION //////////////////////////////////////////////////////////////////////////// //*-- AUTHOR : G. Kornakov // // HParticleT0Reco // // Reconstructor to recalulate the T0 form the reconstructed particle // to increase the time resolution. By default T0 is calulated form all particles // flagged as kIsUsed. The user has to take care that enough used particles // are avalilable. As standard hadrons and leptons are flagged. As alternative // all candidates with fitted inner/outer segment,RKchi2<1000,MetaHit and beta!=-1 // can be used when setUseFlaggedCandidates(kFALSE) is used. It is assumed // that metaQA normalization , momentum correction for systematics and path length // correction is applied before (HParticleCandFiller). // // Particles with id 2, 3 ( e+, e-) (without ues of RICH/SHOWER) // 5, 6 (mu+,mu-) // 8, 9 (pi+,pi-) // 11,12 (K+ , K-) // 14 (p) // 45 (d) // 46 (t) // 49 (he3) // // are identified by beta vs momentum and beta vs MDC dedx simultaniously. Each // Track is tested for all hyposisis. The PID with smallest deviation from // theoretical curves is chosen. beta vs mom is weighted higher of beta vs dedx. // The procedure is performed in 2 steps. In the first iteration a rough temporay ID is // performed to retrieve T0. After correction of TO in the second iteration the PIDs // are asigned. Circular buffers for the mean and sigmas of the different START // strips are filled to improve resolution and to take care of variation with time. // Particle ID takes into account the errors of time measurement. //-------------------------------------------------------------------------- // USAGE: // // HParticleT0Reco t0Reco("apr12"); // t0Reco.setUseFlaggedCandidates(kTRUE); // // // 2 possible ways : // a. As Reconstructor add it to standard task list // b. With HLoop : // // before eventLoop : t0Reco.init(); // inside eventloop : t0Reco.execute(); // // //////////////////////////////////////////////////////////////////////////// #include "hparticlet0reco.h" #include "hades.h" #include "hevent.h" #include "hphysicsconstants.h" #include "hhistmap.h" #include "hhistconverter.h" #include "hparticletool.h" #include "hcategory.h" #include "hparticlecand.h" #include "hparticlecandsim.h" #include "hstart2hit.h" #include "hstart2cal.h" #include "hstartdef.h" ClassImp(HParticleT0Reco) HParticleT0Reco::HParticleT0Reco(const Text_t* name,const Text_t* title, const Text_t* beamtime) : HReconstructor(name,title) { run = beamtime; run.ToLower(); c = 299.792458; t0rpc = 0.; t0tof = 0.; fUseFlagged =kTRUE; fisSimulation=kFALSE; fisEmbedding =kFALSE; } HParticleT0Reco::~HParticleT0Reco() { } Int_t HParticleT0Reco::execute(){ static Int_t checkSim =0; static Int_t checkEmb =0; if(checkSim == 0 || checkEmb == 0){ if(candCat) { HParticleCand* pCand; for(Int_t j=0;jgetEntries();j++) { pCand = (HParticleCand*)candCat->getObject(j); if(!pCand) continue; HParticleCandSim* csim = dynamic_cast(pCand); if(csim){ checkSim = 1; fisSimulation = kTRUE; Int_t tr =-1; Bool_t found=kFALSE; for(Int_t i=0; i<2; i++){ tr = csim->getGeantTrackInnerMdc(i); if(tr==gHades->getEmbeddingRealTrackId()) {found = kTRUE; break;} tr = csim->getGeantTrackOuterMdc(i); if(tr==gHades->getEmbeddingRealTrackId()) {found = kTRUE; break;} } if(found) { fisEmbedding = kTRUE; checkEmb =1; break; } } } } } fill(); return 0; } Bool_t HParticleT0Reco::init() { //--------------------------------------------------------------------------- candCat = gHades->getCurrentEvent()->getCategory(catParticleCand); catStart = gHades->getCurrentEvent()->getCategory(catStart2Hit); catStartcal = gHades->getCurrentEvent()->getCategory(catStart2Cal); fEvHeader = (HEventHeader*)gHades->getCurrentEvent()->getHeader(); if (!candCat) { Error("init","No ParticleCand Category"); return kFALSE; } if (!catStart ) { Error("init","No Start2Hit Category"); return kFALSE; } /* if (!catStartcal ) { Error("init","No Start2Cal Category"); return kFALSE; } */ //--------------------------------------------------------------------------- //--------------------------------------------------------------------------- // build buffers list temp(1,0.); for(Int_t i=0;i > temp2; for(Int_t j=0;j setDefaultPar(run); } else eLoss = HEnergyLossCorrPar::getObject() ; //--------------------------------------------------------------------------- return kTRUE; } void HParticleT0Reco::fillMETA(Float_t t0, Float_t err,Int_t sys) { buff.t0 [sys].push_back(t0); buff.t0err[sys].push_back(err); } void HParticleT0Reco::fillBuffer() { for(Int_t system = 0; system < 2; system++){ fillMETA(getT0(system),getT0Err(system),system); // If the buffer is already larger than 1000, empty the first! if(buff.t0 [system].size() > 200)buff.t0 [system].pop_front(); if(buff.t0err[system].size() > 200)buff.t0err[system].pop_front(); } Float_t startHitTime; HStart2Hit* sH=0; if((sH = (HStart2Hit *) catStart->getObject(0)) != NULL){ Int_t flag = sH->getCorrFlag(); startHitTime = sH->getTime(); if(catStartcal && flag >= 0) { HStart2Cal* sHc=0; for(Int_t i = 0; i < catStartcal->getEntries(); i ++) { if((sHc = (HStart2Cal *) catStartcal->getObject(i)) != NULL){ // find the corresponding start hits to the used Int_t m = sHc->getModule(); Int_t st = sHc->getStrip(); if( m!=0 && m!=1 ) continue; // skip VETO startvalues& startval = start[m][st]; for(Int_t j = 1;j <= sHc->getMultiplicity() && j <= sHc->getMaxMultiplicity(); j ++) { if( fabs(sHc->getTime(j)-startHitTime) < 0.3) { if(m == 0 || m == 1) { buff.meanqstrip[m][st].push_back(sHc->getWidth(j)); if(m == 0) { Float_t t0r = (-1.*(startHitTime - sHc->getTime(j)) + startval.h->Interpolate(sHc->getWidth(j) - startval.meanq) ); buff.t0forrpcstrip[m][st].push_back(getT0(0) - t0r); } if(m == 1) { Float_t t0r = (+1.*(startHitTime - sHc->getTime(j)) + startval.h->Interpolate(sHc->getWidth(j) - startval.meanq) ); buff.t0forrpcstrip[m][st].push_back(getT0(0) - t0r); } buff.t0forrpcstriperr [m][st].push_back( getT0Err(0) ); if(buff.t0forrpcstrip [m][st].size()>50 )buff.t0forrpcstrip [m][st].pop_front(); if(buff.t0forrpcstriperr[m][st].size()>50 )buff.t0forrpcstriperr[m][st].pop_front(); if(buff.meanqstrip [m][st].size()>100)buff.meanqstrip [m][st].pop_front(); } } } // endloop hits cal } // startcal } // endloop start cal } // startcal + corrflag>=0 } // starthit } void HParticleT0Reco::fillCorrections() { t0tof = 0.; t0rpc = 0.; t0startcont = 0.; t0startconterr = 0.; // local vars Float_t errtmp = 0.; Float_t err[2] = {0,0}; Float_t t0 [2] = {0,0}; for(Int_t sys = 0; sys < 2; sys++){ if(buff.t0[sys].size()>2) { list::iterator ite = buff.t0err[sys].begin(); for(list::iterator it = buff.t0[sys].begin();it != buff.t0[sys].end(); it++) { errtmp = *ite; t0 [sys] += (*it)/(errtmp*errtmp); err[sys] += 1./(errtmp*errtmp); ite++; } if(err[sys]) t0[sys]/=err[sys]; else t0[sys] = 0.; } else t0[sys] = 0.; } t0tof = t0[1]; t0rpc = t0[0]; for(Int_t m=0;m::iterator it = buff.meanqstrip[mod][strip].begin();it!=buff.meanqstrip[mod][strip].end();it++) { if((*it)!=0) { qmean += (*it); n++; } } if(n>10) { qmean/=n; } else qmean=0.; } else qmean=0.; return qmean; } Float_t HParticleT0Reco::getOnlineStripCorrection(Int_t mod, Int_t strip) { Float_t tcorr = 0; Int_t n = 0; Float_t errtmp = 0.; Float_t err = 0.; if(buff.t0forrpcstrip[mod][strip].size()) { list::iterator ite = buff.t0forrpcstriperr[mod][strip].begin(); for(list::iterator it = buff.t0forrpcstrip[mod][strip].begin();it!=buff.t0forrpcstrip[mod][strip].end();it++) { if((*it)!=0) { errtmp = (*ite); tcorr += (*it)/(errtmp*errtmp); err += 1./(errtmp*errtmp); n++; } ite++; } if(err) { tcorr/=err; } else tcorr=0.; } else tcorr=0.; return tcorr; } void HParticleT0Reco::fill() { // This function fills the particle buffers for TOF and RPC. // It suppose that the meta matching was already recalculated!!!!! // It uses tracks with some minimum quality criteria. // Also it does a mild cut on vertex : dist>(1./mom)*1200.+10. vector vparticleList; for(UInt_t i = 0;igetVertexReco().getPos(); if(candCat) { for(Int_t j=0;jgetEntries();j++) { pCand = (HParticleCand*)candCat->getObject(j); if(!pCand) continue; if(fUseFlagged) { if (!pCand->isFlagBit(Particle::kIsUsed)) continue; } else { if(pCand->getBeta()==-1 || pCand->getChi2() > 500 || !pCand->isFlagAND(4, Particle::kIsAcceptedHitInnerMDC, Particle::kIsAcceptedHitOuterMDC, Particle::kIsAcceptedHitMETA, Particle::kIsAcceptedRK) ) continue; } if(pCand->isFakeRejected()) continue; if(pCand->getChi2()>200.) continue; if(pCand->getMetaMatchQuality()>3.) continue; Float_t mom = pCand->getMomentum(); Float_t theta = pCand->getTheta(); Int_t sys = pCand->getSystemUsed(); if(sys<0) continue; HGeomVector base, dir; HParticleTool::calcSegVector(pCand->getZ(), pCand->getR(), pCand->getPhi(), theta, base, dir); Float_t dist = HParticleTool::calculateMinimumDistanceStraightToPoint(base,dir,vertex); if(dist>(1./mom)*1200.+ 10.) continue; //-------------------------------------------------------------- // All the possible hypothesis! // 5 sigma cut for(UInt_t i = 0;i6.) continue; Float_t dtime = val.deltaTime; Float_t sigma = val.sigmaTime; //-------------------------------------------------------------- tofvalues tofval; tofval.c = pCand; tofval.index = j; tofval.times = dtime; tofval.error = 1./sigma/sigma; // Push back values! if(pCand->getSelectedMeta()==Particle::kRpcClst) { eventRPC.push_back(tofval); } else { eventTOF.push_back(tofval); } } } fillCorrections(); fillBuffer(); calculateStartT0Cont(); // Improves start time measurement by using making the walk correction and by dynamic resolution. correctBeta(); // Just fix beta+mass in particleCand, measured time is not going to be modified, neither START. setPIDs(); // Function which sets the PIDs of the already corrected particles. } void HParticleT0Reco::setPIDs() { // Function which loops over candidates onces corrected by T0 and sets the best hypo just to the best traces. HParticleCand* pCand=0; vector vparticleList; for(UInt_t i = 0;igetEntries();j++) { pCand = (HParticleCand*)candCat->getObject(j); if(!pCand) continue; pCand->setPID(-1); if(fUseFlagged) { if (!pCand->isFlagBit(Particle::kIsUsed)) continue; } else { if(pCand->getBeta()==-1 || pCand->getChi2() > 1000 || !pCand->isFlagAND(4, Particle::kIsAcceptedHitInnerMDC, Particle::kIsAcceptedHitOuterMDC, Particle::kIsAcceptedHitMETA, Particle::kIsAcceptedRK) ) continue; } if( pCand->isFakeRejected() ) continue; if( pCand->getMetaMatchQuality() > 4.5 ) continue; Int_t sys = pCand->getSystemUsed(); if(sys<0) continue; //-------------------------------------------------------------- // All the possible hypothesis! for(UInt_t i = 0;i4.5) continue; //-------------------------------------------------------------- pCand->setPID(val.id); if(val.id>0) pCand->calc4vectorProperties(HPhysicsConstants::mass(val.id)); } } } void HParticleT0Reco::correctBeta() { // is a two step procedure; first recalculate all betas+mass, // then use them to modify the category values. Othwerwise, // the modified beta value would modify the output of the // loop if done in just one single step! vector betas; HParticleCand* pCand=0; if(candCat) { for(Int_t j=0;jgetEntries();j++){ pCand = (HParticleCand*)candCat->getObject(j); if(!pCand) continue; betas.push_back(getBeta(pCand)); } for(Int_t j=0;jgetEntries();j++){ pCand = (HParticleCand*)candCat->getObject(j); if(!pCand) continue; if(pCand->getBeta()!=-1 && pCand->isFlagAND(4, Particle::kIsAcceptedHitInnerMDC, Particle::kIsAcceptedHitOuterMDC, Particle::kIsAcceptedHitMETA, Particle::kIsAcceptedRK) ) { pCand->setBeta(betas[j]); Float_t mom = pCand->getMomentum(); Float_t mass2 = mom*mom/betas[j]/betas[j]*(1.-betas[j]*betas[j]); pCand->setMass2(mass2); } } } } void HParticleT0Reco::calculateStartT0Cont() { t0startcont = 0.; t0startconterr = 0.; Float_t startHitTime = 0.; // As it is used for META correction is 0 by deffinition! Float_t fstartSimWidth = -1000; if(catStart) { HStart2Hit* sH=0; if((sH = (HStart2Hit *) catStart->getObject(0)) != NULL){ Int_t flag = sH->getCorrFlag(); startHitTime = sH->getTime(); //------------------------------------------------------- // pure simulation fstartSimWidth = sH->getSimWidth(); if(fstartSimWidth<0) fstartSimWidth = 0.06; if(fisSimulation && !fisEmbedding){ t0startconterr += fstartSimWidth; t0startcont = 0; return; } //------------------------------------------------------- if(catStartcal && flag>=0) { Bool_t found = kFALSE; HStart2Cal* sHc=0; for(Int_t i=0;i< catStartcal->getEntries();i++) { if((sHc = (HStart2Cal *) catStartcal->getObject(i)) != NULL){ // find the corresponding start hits to the used Int_t m = sHc->getModule(); Int_t st = sHc->getStrip(); if( m!=0 && m!=1 ) continue; // skip VETO startvalues& startval = start[m][st]; for(Int_t j=1;j<=sHc->getMultiplicity()&&j<=sHc->getMaxMultiplicity();j++) { if( fabs(sHc->getTime(j)-startHitTime)<0.3) { found = kTRUE; // This start cal was used! Float_t errtmp = startval.he->Interpolate(sHc->getWidth(j)- startval.meanq)*1.2; if(errtmp) t0startconterr+= 1./errtmp/errtmp; else { t0startconterr += 1./0.2/0.2; errtmp = 0.2; } if(m==0) { t0startcont += (-1.*(startHitTime - sHc->getTime(j)) + startval.h->Interpolate(sHc->getWidth(j) - startval.meanq) + 2*startval.scor )*1./errtmp/errtmp; } if(m==1) { t0startcont += (+1.*(startHitTime - sHc->getTime(j)) + startval.h->Interpolate(sHc->getWidth(j) - startval.meanq) + 2*startval.scor )*1./errtmp/errtmp; } } // match arround starthit } // endloop hits startcal } // startcal } // endloop startcal } // startCal + corrflag >-1 } // startHit } // catStart } Float_t HParticleT0Reco::getBeta(HParticleCand* pCand) { Float_t t0 = 0; Float_t err = 0; Float_t t02 = 0; Float_t err2 = 0; if(eventTOF.size()+eventRPC.size() <1 || pCand->getBeta() == -1) return pCand->getBeta(); t0 += t0startcont; err += t0startconterr; // Read FLAGs in Hit // Read time in Hit for(UInt_t i=0;igetSelectedMeta()==Particle::kRpcClst ) t0 += val.times*val.error; else t0 += (val.times + (t0tof-t0rpc) )*val.error*0.8; err += val.error; } } for(UInt_t i=0;igetSelectedMeta()==Particle::kRpcClst ) t0 += (val.times - (t0tof-t0rpc) )*val.error*0.5; else t0 += (val.times)*val.error; err += val.error; } } // Remove outliers and introduce assymetry factor in order to get the right velocity // This is due to a non gaussian shape of momentum distribution. if(err!=0) t0 = t0/err; else t0 = 0.; t02 += t0startcont; err2 += t0startconterr; // Read FLAGs in Hit // Read time in Hit Float_t assymetryFact = 0; for(UInt_t i=0;igetSelectedMeta()==Particle::kRpcClst ) { if( fabs(t0-val.times)<6.*sqrt(1./val.error) ) { assymetryFact = (t0-val.times) < 0 ? 2 : 10; t02 += val.times* val.error * pow(1./(1. + fabs(t0-val.times)*sqrt(val.error)/assymetryFact), 2); err2 += val.error * pow(1./(1. + fabs(t0-val.times)*sqrt(val.error)/assymetryFact), 2) ; } else continue; } else { if( fabs(t0 - val.times - (t0tof-t0rpc))<6.*sqrt(1./val.error) ) { assymetryFact = (t0-val.times- (t0tof-t0rpc)) < 0 ? 2 : 10; t02 += (val.times + (t0tof-t0rpc) )*val.error* 0.8* pow(1./(1. + fabs(t0-val.times- (t0tof-t0rpc))*sqrt(val.error)/assymetryFact), 2); err2 += val.error *0.8* pow(1./(1. + fabs(t0-val.times- (t0tof-t0rpc))*sqrt(val.error)/assymetryFact), 2) ; } else continue; } } } for(UInt_t i=0;igetSelectedMeta()==Particle::kRpcClst ) { if( fabs(t0-val.times + (t0tof-t0rpc))<6.*sqrt(1./val.error) ) { assymetryFact = (t0-val.times + (t0tof-t0rpc)) < 0 ? 2 : 10; t02 += (val.times + (-t0tof+t0rpc) )*val.error*0.5*pow(1./(1. + fabs(t0-val.times + (t0tof-t0rpc))*sqrt(val.error)/assymetryFact), 2); err2 += val.error* 0.5*pow(1./(1. + fabs(t0-val.times + (t0tof-t0rpc))*sqrt(val.error)/assymetryFact), 2); } else continue; } else { if( fabs(t0-val.times)<6.*sqrt(1./val.error) ) { assymetryFact = (t0-val.times) < 0 ? 2 : 10; t02 += (val.times)*val.error * pow(1./(1. + fabs(t0-val.times)*sqrt(val.error)/assymetryFact), 2); err2 += val.error * pow(1./(1. + fabs(t0-val.times)*sqrt(val.error)/assymetryFact), 2); } else continue; } } } if(err!=0) t02 = t02/err2; else t02 = 0.; Float_t oldtof = pCand->getDistanceToMetaHit()/pCand->getBeta()/c; Float_t newBeta = pCand->getDistanceToMetaHit()/(oldtof - t02)/c; return newBeta; } Float_t HParticleT0Reco::getT0(Int_t sys) { Float_t t0 = 0; Float_t err = 0; Float_t t02 = 0; Float_t err2 = 0; Float_t assymetryFact = 0; vector& event = *eventMETA[sys]; if(event.size()<1) return t0; for(UInt_t i=0;i0.) t02 = t02/err2; else t02 = 0.; return t02; } Float_t HParticleT0Reco::getT0Err(Int_t sys) { Float_t err = 0; vector& event = *eventMETA[sys]; if(event.size()<1) return 100.; for(UInt_t i=0;iGet(Form("h_q_maxdt_mod%i_strip%i",i,j)); // he[i][j] = (TH1D*)fileWC->Get(Form("h_q_sdt_mod%i_strip%i" ,i,j)); // HHistConverter::fillArray(h [i][j],astartdt [i][j],Form("h_q_maxdt_mod%i_strip%i",i,j),10,10, kFALSE); // HHistConverter::fillArray(he[i][j],astartsdt[i][j],Form("h_q_sdt_mod%i_strip%i" ,i,j),10,10, kFALSE); // } // } // // // // // ofstream out; // out.open("myhists.txt"); // // for(Int_t i=0;i<2;i++) { // for(Int_t j=0;j<16;j++) { // HHistConverter::writeArray(out,h [i][j]->GetName(),astartdt [i][j],10); // } // } // for(Int_t i=0;i<2;i++) { // for(Int_t j=0;j<16;j++) { // HHistConverter::writeArray(out,he[i][j]->GetName(),astartsdt[i][j],10); // } // } // out.close(); // // // out.open("myhists_array.txt"); // // Int_t imax = 2; // Int_t jmax = 16; // Int_t kmax = astartdt[0][0].GetSize(); // // out<0)out<0&&k0)out<0&&kSetDirectory(0); start[i][j].he->SetDirectory(0); start[i][j].meanqall = tmp[i][j]; } } }