//
+
#include "AliCFTrackCutPid.h"
#include "AliLog.h"
#include <TMath.h>
fMinDiffProbability(0.001),
fgParticleType(10),
fgIsComb(kTRUE),
+ fgIsAOD(kFALSE),
fCheckResponse(kFALSE),
fCheckSelection(kTRUE),
fIsPpriors(kFALSE),
fNbins(101),
fDetRestr(-1),
fiPartRestr(-1),
- fDetProbRestr(1)
+ fDetProbRestr(1),
+ fProbThreshold(0.)
{
//
//Default constructor
fMinDiffProbability(0.001),
fgParticleType(10),
fgIsComb(kTRUE),
+ fgIsAOD(kFALSE),
fCheckResponse(kFALSE),
fCheckSelection(kTRUE),
fIsPpriors(kFALSE),
fNbins(101),
fDetRestr(-1),
fiPartRestr(-1),
- fDetProbRestr(1)
+ fDetProbRestr(1),
+ fProbThreshold(0.)
{
//
//Constructor
fMinDiffProbability(c.fMinDiffProbability),
fgParticleType(c.fgParticleType),
fgIsComb(c.fgIsComb),
+ fgIsAOD(c.fgIsAOD),
fCheckResponse(c.fCheckResponse),
fCheckSelection(c.fCheckSelection),
fIsPpriors(c.fIsPpriors),
fNbins(c.fNbins),
fDetRestr(c.fDetRestr),
fiPartRestr(c.fiPartRestr),
- fDetProbRestr(c.fDetProbRestr)
+ fDetProbRestr(c.fDetProbRestr),
+ fProbThreshold(c.fProbThreshold)
{
//
//Copy constructor
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;
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];
}
}else{
Double_t calcprob[5];
- CombPID(status,pid,calcprob);
+ CombPID(status,pid,calcprob);
iPart = Identify(calcprob);
+
}
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
{
// 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; iP<AliPID::kSPECIES; iP++) {
+ probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
+ AliDebug(2,Form("prob %i %f",iP, probability[iP]));
+ if(fIsQAOn) fhCombProb[iP]->Fill(probability[iP]);
}
- else getpid.SetPriors(fPriors);
-
-
- AliPID::EParticleType sel = getpid.GetMostProbable();
- Double_t probability[AliPID::kSPECIES];
- for(Int_t iP=0; iP<AliPID::kSPECIES; iP++) probability[iP] = getpid.GetProbability((AliPID::EParticleType)iP);
-
- if(getpid.GetProbability(sel,fPriors)>fCut) 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;
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; iP<AliPID::kSPECIES; iP++) {
//
// method for the pid-cut selction
//
+ Bool_t sel = kFALSE;
if (!track) return kFALSE ;
TString className(track->ClassName());
- 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<AliESDtrack*>(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<AliAODTrack *>(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
//
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<AliPID::kSPECIES; j++){
for(Int_t i=0; i< kNdets; i++){
if(!fDets[i]) continue;
AliDebug(1,Form("-----> 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; nn<AliPID::kSPECIES; nn++) {
- combpid[nn]=pid[kNdets][nn];
- sumesdpid+=fPriors[nn]*combpid[nn];
- if(fIsQAOn) {
- if(!fhCombResp[nn]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() 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; k<AliPID::kSPECIES; k++){
- if(fIsQAOn) {
- if(!fhCombResp[k]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() 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; nn<AliPID::kSPECIES; nn++) sumesdpid+=pid[kNdets][nn];
+ if(sumesdpid<=0) {
+ AliDebug(1,"priors or ESDpid are inconsistent, please check them");
+ return;
+ } else {
+ for(Int_t k=0; k<AliPID::kSPECIES; k++){
+ combpid[k] = pid[kNdets][k]/sumesdpid;
+ if(fIsQAOn) {
+ if(!fhCombResp[k]) AliDebug(1,Form("\n no fhCombResp[%i], check if pidcut->Init() 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; n<AliPID::kSPECIES; n++) {
- combpid[n]=fPriors[n]*prod[n]/sum;
- if(fIsQAOn) {
- if(!fhCombProb[n]) AliDebug(1,Form("no fhCombProb[%i] defined, check if pidcut->Init() 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]));
+
}
//__________________________________________
//
//
//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()));
+ }
+ }
+ }
+
}
//___________________________________________________
if(!fIsQAOn) return;
DefineHistograms();
- if(fgIsComb){
+ if(fgIsComb || fgIsAOD){
for(Int_t iPart =0; iPart<AliPID::kSPECIES; iPart++){
qalist->Add(fhCombResp[iPart]);
qalist->Add(fhCombProb[iPart]);
}
}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff