X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=CORRFW%2FAliCFTrackCutPid.cxx;h=39607355edfe9d165f750c7294e87c819b7c3a67;hb=fc0882f371b96f1a712eaa5f3475a2ae575fb254;hp=14ae79d5ae480c593b901669b1b3266c8283a581;hpb=107a31009f303a729cd6513c1931d8f64e3a0e54;p=u%2Fmrichter%2FAliRoot.git diff --git a/CORRFW/AliCFTrackCutPid.cxx b/CORRFW/AliCFTrackCutPid.cxx index 14ae79d5ae4..39607355edf 100644 --- a/CORRFW/AliCFTrackCutPid.cxx +++ b/CORRFW/AliCFTrackCutPid.cxx @@ -56,6 +56,7 @@ // + #include "AliCFTrackCutPid.h" #include "AliLog.h" #include @@ -73,6 +74,7 @@ AliCFTrackCutPid::AliCFTrackCutPid() : fMinDiffProbability(0.001), fgParticleType(10), fgIsComb(kTRUE), + fgIsAOD(kFALSE), fCheckResponse(kFALSE), fCheckSelection(kTRUE), fIsPpriors(kFALSE), @@ -82,7 +84,8 @@ AliCFTrackCutPid::AliCFTrackCutPid() : fNbins(101), fDetRestr(-1), fiPartRestr(-1), - fDetProbRestr(1) + fDetProbRestr(1), + fProbThreshold(0.) { // //Default constructor @@ -107,6 +110,7 @@ AliCFTrackCutPid::AliCFTrackCutPid(const Char_t* name, const Char_t* title) : fMinDiffProbability(0.001), fgParticleType(10), fgIsComb(kTRUE), + fgIsAOD(kFALSE), fCheckResponse(kFALSE), fCheckSelection(kTRUE), fIsPpriors(kFALSE), @@ -116,7 +120,8 @@ AliCFTrackCutPid::AliCFTrackCutPid(const Char_t* name, const Char_t* title) : fNbins(101), fDetRestr(-1), fiPartRestr(-1), - fDetProbRestr(1) + fDetProbRestr(1), + fProbThreshold(0.) { // //Constructor @@ -141,6 +146,7 @@ AliCFTrackCutPid::AliCFTrackCutPid(const AliCFTrackCutPid& c) : fMinDiffProbability(c.fMinDiffProbability), fgParticleType(c.fgParticleType), fgIsComb(c.fgIsComb), + fgIsAOD(c.fgIsAOD), fCheckResponse(c.fCheckResponse), fCheckSelection(c.fCheckSelection), fIsPpriors(c.fIsPpriors), @@ -150,7 +156,8 @@ AliCFTrackCutPid::AliCFTrackCutPid(const AliCFTrackCutPid& c) : fNbins(c.fNbins), fDetRestr(c.fDetRestr), fiPartRestr(c.fiPartRestr), - fDetProbRestr(c.fDetProbRestr) + fDetProbRestr(c.fDetProbRestr), + fProbThreshold(c.fProbThreshold) { // //Copy constructor @@ -180,9 +187,10 @@ AliCFTrackCutPid& AliCFTrackCutPid::operator=(const AliCFTrackCutPid& c) AliCFCutBase::operator=(c) ; this->fCut=c.fCut; this->fMinDiffResponse=c.fMinDiffResponse; - this->fMinDiffProbability=c.fMinDiffProbability; + this->fMinDiffProbability=c.fMinDiffProbability; this->fgParticleType=c.fgParticleType; this->fgIsComb=c.fgIsComb; + this->fgIsAOD=c.fgIsAOD; this->fCheckResponse=c.fCheckResponse; this->fCheckSelection=c.fCheckSelection; this->fIsPpriors=c.fIsPpriors; @@ -193,6 +201,7 @@ AliCFTrackCutPid& AliCFTrackCutPid::operator=(const AliCFTrackCutPid& c) this->fDetRestr=c.fDetRestr; this->fiPartRestr=c.fiPartRestr; this->fDetProbRestr=c.fDetProbRestr; + this->fProbThreshold=c.fProbThreshold; for(Int_t i=0; i< kNdets ; i++ ) { this->fDets[i]=c.fDets[i]; @@ -394,8 +403,9 @@ Int_t AliCFTrackCutPid::GetID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][Al } }else{ Double_t calcprob[5]; - CombPID(status,pid,calcprob); + CombPID(status,pid,calcprob); iPart = Identify(calcprob); + } @@ -415,6 +425,21 @@ Int_t AliCFTrackCutPid::GetID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][Al return iPart; } +//_________________________________________________________________________________ +Int_t AliCFTrackCutPid::GetAODID(AliAODTrack *aodtrack) const +{ +// +// Identifies the AOD Track using the combined pid responses +// + + Double_t combpid[AliPID::kSPECIES]; + for(Int_t i=0; i< AliPID::kSPECIES; i++) { + combpid[i]= aodtrack->PID()[i]; + if(!fhCombResp[i]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called",i)); + else fhCombResp[i]->Fill(combpid[i]); + } + return Identify(combpid); +} //__________________________________ Bool_t AliCFTrackCutPid::Check(const Double_t *p, Int_t iPsel, Double_t minDiff) const { @@ -443,25 +468,31 @@ Int_t AliCFTrackCutPid::Identify(Double_t pid[AliPID::kSPECIES]) const // The identification is actually performed here with possible // checks on the det responses and/or probabilities // + Int_t iPart = -1; - + + AliDebug(2,Form("calc response bef: %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4])); + AliDebug(2,Form("priors : %f %f %f %f %f",fPriors[0],fPriors[1],fPriors[2],fPriors[3],fPriors[4])); + AliPID getpid(pid,kTRUE); + getpid.SetPriors(fPriors); - if(fgIsComb) { - Double_t priors[5]={0.2,0.2,0.2,0.2,0.2}; - getpid.SetPriors(priors); + Double_t probability[AliPID::kSPECIES]={0.,0.,0.,0.,0.}; + for(Int_t iP=0; iPFill(probability[iP]); } - else getpid.SetPriors(fPriors); - - - AliPID::EParticleType sel = getpid.GetMostProbable(); - Double_t probability[AliPID::kSPECIES]; - for(Int_t iP=0; iPfCut) iPart= (Int_t)sel; - AliDebug(2,Form("calc response : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4])); - AliDebug(2,Form("probabilities : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4])); + + if (fProbThreshold > 0.) { + if (probability[fgParticleType] >= fProbThreshold) iPart=fgParticleType; + } + else { + AliPID::EParticleType sel = getpid.GetMostProbable(); + if(getpid.GetProbability(sel,fPriors)>fCut) iPart= (Int_t)sel; + AliDebug(2,Form("probabilities : %f %f %f %f %f",probability[0],probability[1],probability[2],probability[3],probability[4])); + } if(fCheckResponse && !Check(pid,iPart, fMinDiffResponse)) iPart=kCheckResp; if(fCheckSelection && !Check(probability,iPart,fMinDiffProbability)) iPart=kCheckProb; @@ -480,13 +511,8 @@ Int_t AliCFTrackCutPid::IdentifyQA(const Double_t pid[AliPID::kSPECIES], Int_t i AliDebug(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4])); AliPID getpid(pid,kTRUE); - - if(fgIsComb) { - Double_t priors[5]={0.2,0.2,0.2,0.2,0.2}; - getpid.SetPriors(priors); - } - else getpid.SetPriors(fPriors); - + getpid.SetPriors(fPriors); + AliPID::EParticleType sel = getpid.GetMostProbable(); Double_t probability[AliPID::kSPECIES]; for(Int_t iP=0; iPClassName()); - if (className.CompareTo("AliESDtrack") != 0) { - AliError("obj must point to a AliESDtrack "); - return kFALSE ; - } - - AliESDtrack *esdTrack = (AliESDtrack *)track; + if (className.CompareTo("AliESDtrack") == 0) { + AliESDtrack *esdTrack = dynamic_cast(track); ULong_t status[kNdets+1]={0,0,0,0,0,0}; Double_t pid[kNdets+1][AliPID::kSPECIES]; TrackInfo(esdTrack,status,pid); if(fIsPpriors) SetPPriors(esdTrack); - if(GetID(status,pid)==fgParticleType) return kTRUE; - else return kFALSE; + if(GetID(status,pid)==fgParticleType) sel = kTRUE; + } + + if (className.CompareTo("AliAODTrack") == 0) { + AliAODTrack *aodtrack = dynamic_cast(track); + if(GetAODID(aodtrack) == fgParticleType) sel = kTRUE; + } + + return sel; + } //__________________________________ void AliCFTrackCutPid::CombPID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][AliPID::kSPECIES],Double_t *combpid) const @@ -535,16 +566,15 @@ void AliCFTrackCutPid::CombPID(ULong_t status[kNdets+1],Double_t pid[kNdets+1][ // Bool_t isdet=kFALSE; - Double_t sum=0.; Double_t prod[AliPID::kSPECIES]={1.,1.,1.,1.,1.}; - Double_t comb[AliPID::kSPECIES]={0.,0.,0.,0.,0.}; - Double_t priors[AliPID::kSPECIES]={0.2,0.2,0.2,0.2,0.2}; ULong_t andstatus =0; if(fIsDetAND) { andstatus = StatusForAND(status); AliDebug(1,Form("AND combination %lu",andstatus)); } + + //Products of single detector responses for(Int_t j=0; j trk status %lu and status %lu -> trk-ANDdetector status combination %lu",status[kNdets],andstatus,status[kNdets]&andstatus)); AliDebug(1,Form("In det %i -> particle %i response is %f",i,j,pid[i][j])); } - } - else { + } else { prod[j]*=pid[i][j]; isdet=kTRUE; AliDebug(2,Form("In det %i -> particle %i response is %f",i,j,pid[i][j])); + + if(fIsQAOn){ + if(!fhResp[i][j]) {AliDebug(1,Form("no pointer to the histo fhResp%i%i, check if pidcut->Init() was called",i,j));} + else fhResp[i][j]->Fill(pid[i][j]); + + if(!fhProb[i][j]) {AliDebug(1,Form("no pointer to the histo fhProb%i%i, check if pidcut->Init() was called",i,j));} + else { + AliPID detprob(pid[i],kTRUE); + detprob.SetPriors(fPriors); + fhProb[i][j]->Fill(detprob.GetProbability((AliPID::EParticleType)j)); + } + } } }//combined mode }//loop on dets }//loop on species - - if(fIsQAOn) { - for(Int_t iqa =0; iqa < kNdets; iqa++){ - if(!fDets[iqa]) continue; - AliPID normresp(pid[iqa]); - normresp.SetPriors(priors); - for(Int_t ip =0; ip< AliPID::kSPECIES; ip++){ - if(!fhResp[iqa][ip]) {AliDebug(1,Form("no pointer to the histo fhResp%i%i, check if pidcut->Init() was called",iqa,ip));} - else fhResp[iqa][ip]->Fill(normresp.GetProbability((AliPID::EParticleType)ip)); - }//loop on part - }//loop on dets - }//if qa - - - Double_t add = 0; for(Int_t isumm=0; isumm<5; isumm++) add+=prod[isumm]; - if(add>0) AliDebug(1,Form("calculated comb response: %f %f %f %f %f",prod[0]/add,prod[1]/add,prod[2]/add,prod[3]/add,prod[4]/add)); - AliDebug(1,Form("the ESDpid response: %f %f %f %f %f",pid[kNdets][0],pid[kNdets][1],pid[kNdets][2],pid[kNdets][3],pid[kNdets][4])); + + //no detectors found, then go to ESD pid... if(!isdet) { AliWarning("\n !! No detector found for the combined pid --> responses are from the ESDpid !! \n"); Double_t sumesdpid=0; - for(Int_t nn=0; nnInit() was called",nn)); - else fhCombResp[nn]->Fill(combpid[nn]); - }//fIsQAOn - }//loop on species - if(sumesdpid > 0) for(Int_t ih = 0; ih < AliPID::kSPECIES; ih++) { - if(!fhCombResp[ih]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called \n",ih)); - else fhCombProb[ih]->Fill(fPriors[ih]*combpid[ih]/sumesdpid); - } - else AliDebug(1,"priors or ESDpid are zero, please check them"); - }// end no det - - else{ - for(Int_t k=0; kInit() was called \n",k)); - else fhCombResp[k]->Fill(prod[k]); - } - AliDebug(1,Form("species: %i priors %f and combined prod %f",k,fPriors[k],prod[k])); - comb[k]=fPriors[k]*prod[k]; - AliDebug(1,Form("combined probability %i %f",k,comb[k])); - sum+=comb[k]; - } + for(Int_t nn=0; nnInit() was called",k)); + else fhCombResp[k]->Fill(combpid[k]); + } + }//loop on species + } + return; + } - if(sum == 0) { - AliDebug(1,"Check the detector responses or the priors, their combined products are zero"); - return; - } - for(Int_t n=0; nInit() was called",n)); - else fhCombProb[n]->Fill(combpid[n]); - } - } - }//end else + Double_t add = 0; for(Int_t isumm=0; isumm<5; isumm++) add+=prod[isumm]; + if(add>0) { + for(Int_t ip =0; ip < AliPID::kSPECIES; ip++) { + combpid[ip] = prod[ip]/add; + if(fIsQAOn) { + if(!fhCombResp[ip]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() was called",ip)); + else fhCombResp[ip]->Fill(combpid[ip]); + + } + } + AliDebug(1,Form("calculated comb response: %f %f %f %f %f",combpid[0],combpid[1],combpid[2],combpid[3],combpid[4])); + } else { + AliDebug(1,"single detector responses are inconsistent, please check them...."); + return; + } + AliDebug(1,Form("the ESDpid response: %f %f %f %f %f",pid[kNdets][0],pid[kNdets][1],pid[kNdets][2],pid[kNdets][3],pid[kNdets][4])); + } //__________________________________________ // @@ -652,29 +672,58 @@ void AliCFTrackCutPid::DefineHistograms() // //QA histo booking // - char *detect[kNdets]={"ITS","TPC","TRD","TOF","HMPID"}; - char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"}; - - for(Int_t iDet =0; iDet< kNdets; iDet++) - { - if(!fDets[iDet]) continue; - for(Int_t iP =0; iP < AliPID::kSPECIES; iP++){ - fhResp[iDet][iP] = new TH1F(Form("rDet%iPart%i",iDet,iP),Form("%s %s response ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax); - fhProb[iDet][iP] = new TH1F(Form("pDet%iPart%i",iDet,iP),Form("%s %s probability ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax); - } - } - -if(fgIsComb) - { - for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++) + if(fgIsAOD){ + const char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"}; + + for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++) { - fhCombResp[iPart] = new TH1F(Form("rCombPart%i",iPart),Form(" %s combined response ",partic[iPart]),fNbins,fXmin,fXmax); - Printf(Form( "rCombPart%i is booked!!",iPart)); - fhCombProb[iPart] = new TH1F(Form("pCombPart%i",iPart),Form("%s combined probability ",partic[iPart]),fNbins,fXmin,fXmax); - Printf(Form( "rCombProb%i is booked!!",iPart)); - } - } + fhCombResp[iPart] = new TH1F(Form("%s_rCombPart%i",GetName(),iPart),Form(" %s combined response (AODTrack) ",partic[iPart]),fNbins,fXmin,fXmax); + fhCombResp[iPart]->SetXTitle(Form(" %s combined response ",partic[iPart])); + fhCombResp[iPart]->SetYTitle("entries"); + AliDebug(1,Form( "%s is booked!!",fhCombResp[iPart]->GetName())); + fhCombProb[iPart] = new TH1F(Form("%s_pCombPart%i",GetName(),iPart),Form("%s combined probability (AODTrack) ",partic[iPart]),fNbins,fXmin,fXmax); + fhCombProb[iPart]->SetXTitle(Form(" %s combined probability ",partic[iPart])); + fhCombProb[iPart]->SetYTitle("entries"); + AliDebug(1,Form( "%s is booked!!",fhCombProb[iPart]->GetName())); + } + } + + + else { + const char *detect[kNdets]={"ITS","TPC","TRD","TOF","HMPID"}; + const char *partic[AliPID::kSPECIES]={"electron","muon","pion","kaon","proton"}; + + for(Int_t iDet =0; iDet< kNdets; iDet++) + { + if(!fDets[iDet]) continue; + for(Int_t iP =0; iP < AliPID::kSPECIES; iP++){ + fhResp[iDet][iP] = new TH1F(Form("%s_rDet%iPart%i",GetName(),iDet,iP),Form("%s response for %s ",detect[iDet],partic[iP]),fNbins,fXmin,fXmax); + fhResp[iDet][iP]->SetXTitle(Form(" %s response ",partic[iP])); + fhResp[iDet][iP]->SetYTitle("entries"); + fhProb[iDet][iP] = new TH1F(Form("%s_pDet%iPart%i",GetName(),iDet,iP),Form("%s calculated probability for %s",detect[iDet],partic[iP]),fNbins,fXmin,fXmax); + fhProb[iDet][iP]->SetXTitle(Form(" %s probability ",partic[iP])); + fhProb[iDet][iP]->SetYTitle("entries"); + } + } + + + if(fgIsComb) + { + for(Int_t iPart =0; iPart < AliPID::kSPECIES; iPart++) + { + fhCombResp[iPart] = new TH1F(Form("%s_rCombPart%i",GetName(),iPart),Form(" %s combined response ",partic[iPart]),fNbins,fXmin,fXmax); + fhCombResp[iPart]->SetXTitle(Form(" %s response ",partic[iPart])); + fhCombResp[iPart]->SetYTitle("entries"); + AliDebug(1,Form( "%s is booked!!",fhCombResp[iPart]->GetName())); + fhCombProb[iPart] = new TH1F(Form("%s_pCombPart%i",GetName(),iPart),Form("%s combined probability ",partic[iPart]),fNbins,fXmin,fXmax); + fhCombProb[iPart]->SetXTitle(Form(" %s response ",partic[iPart])); + fhCombProb[iPart]->SetYTitle("entries"); + AliDebug(1,Form( "%s is booked!!",fhCombProb[iPart]->GetName())); + } + } + } + } //___________________________________________________ @@ -686,7 +735,7 @@ void AliCFTrackCutPid::AddQAHistograms(TList *qalist) if(!fIsQAOn) return; DefineHistograms(); - if(fgIsComb){ + if(fgIsComb || fgIsAOD){ for(Int_t iPart =0; iPartAdd(fhCombResp[iPart]); qalist->Add(fhCombProb[iPart]); @@ -702,3 +751,4 @@ void AliCFTrackCutPid::AddQAHistograms(TList *qalist) } } +