//
+
#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
//
for(Int_t j=0; j< AliPID::kSPECIES; j++) {
fPriors[j]=0.2;
+ }
+ for(Int_t j=0; j< AliPID::kSPECIES; j++) {
fPriorsFunc[j]=0x0;
}
-
for(Int_t jDet=0; jDet< kNdets; jDet++) {
fDets[jDet]=kFALSE;
fDetsInAnd[jDet]=kFALSE;
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
//
for(Int_t j=0; j< AliPID::kSPECIES; j++) {
fPriors[j]=0.2;
+ }
+ for(Int_t j=0; j< AliPID::kSPECIES; j++) {
fPriorsFunc[j]=0x0;
}
-
for(Int_t jDet=0; jDet< kNdets; jDet++) {
fDets[jDet]=kFALSE;
fDetsInAnd[jDet]=kFALSE;
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
}
for(Int_t j=0; j< AliPID::kSPECIES; j++){
fPriors[j]=c.fPriors[j];
+ }
+ for(Int_t j=0; j< AliPID::kSPECIES; j++){
fPriorsFunc[j]=c.fPriorsFunc[j];
fhCombResp[j]=c.fhCombResp[j];
fhCombProb[j]=c.fhCombProb[j];
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];
for(Int_t j=0; j< AliPID::kSPECIES; j++){
this->fPriors[j]=c.fPriors[j];
+ }
+ for(Int_t j=0; j< AliPID::kSPECIES; j++){
this->fhCombResp[j]=c.fhCombResp[j];
this->fhCombProb[j]=c.fhCombProb[j];
this-> fPriorsFunc[j]=c.fPriorsFunc[j];
//
//dtor
//
+
for(Int_t i=0; i< kNdets ; i++ ) {
for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
if(fhResp[i][iP])delete fhResp[i][iP];
if(fhCombProb[j])delete fhCombProb[j];
}
-
}
//__________________________________
void AliCFTrackCutPid::SetDetectors(TString dets)
if(dets.Contains("TRD")) {fDets[kTRD]=kTRUE;}
if(dets.Contains("TOF")) {fDets[kTOF]=kTRUE;}
if(dets.Contains("HMPID")) {fDets[kHMPID]=kTRUE;}
+
+ if(dets.Contains("ALL")) for(Int_t i=0; i< kNdets ; i++) fDets[i]=kTRUE;
}
//__________________________________
void AliCFTrackCutPid::SetPriors(Double_t r[AliPID::kSPECIES])
for(Int_t i=0; i< AliPID::kSPECIES; i++) {
if(pTrk->P()>fPriorsFunc[i]->GetXmin() && pTrk->P() < fPriorsFunc[i]->GetXmax()) fPriors[i]=fPriorsFunc[i]->Eval(pTrk->P());
- else {AliInfo("the track momentum is not in the function range. Priors are equal") fPriors[i] = 0.2;}
+ else {AliInfo("the track momentum is not in the function range. Priors are equal"); fPriors[i] = 0.2;}
}
}
//______________________________________
}
}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
{
// is higher than a lower limit.
// Returns: kTRUE= is acceptable
- AliDebug(1,Form("input particle: %i",iPsel));
+ AliDebug(2,Form("input particle: %i",iPsel));
Bool_t ck=kTRUE;
if(iPsel<0) ck=kFALSE;
// 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(1,Form("resp : %f %f %f %f %f",pid[0],pid[1],pid[2],pid[3],pid[4]));
- AliDebug(1,Form("probab : %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;
+
return iPart;
}
//___________________________________________
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++) {
fhProb[idets][iP]->Fill(probability[iP]);
}
- AliPID toresp(pid,kTRUE); Double_t qapriors[5]={0.2,0.2,0.2,0.2,0.2};
+ AliPID toresp(pid,kTRUE);
+ Double_t qapriors[10]={0.2,0.2,0.2,0.2,0.2,0,0,0,0,0};
toresp.SetPriors(qapriors);
for(Int_t iPr=0; iPr<AliPID::kSPECIES; iPr++) fhResp[idets][iPr]->Fill(toresp.GetProbability((AliPID::EParticleType)iPr));
//
// 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);
+ if (!esdTrack) return kFALSE;
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 (!aodtrack) return kFALSE ;
+ 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;
prod[j]*=pid[i][j];
isdet = kTRUE;
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 loop %i -> particle %i response is %f",i,j,pid[i][j]));
+ 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(1,Form("In det loop %i -> particle %i response is %f",i,j,pid[i][j]));
+ 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
-
+ //no detectors found, then go to ESD pid...
if(!isdet) {
- AliDebug(1,"No proper status for the combined pid -> probabilities are set to zero");
- for(Int_t nn=0; nn<AliPID::kSPECIES; nn++) {combpid[nn]=0;}
- }
-
- else{
- for(Int_t k=0; k<AliPID::kSPECIES; k++){
- if(fIsQAOn) {
- if(!fhCombResp[k]) AliDebug(1,Form("no fhCombResp[%i], check if pidcut->Init() was called",k));
- else fhCombResp[k]->Fill(prod[k]);
- }
- AliDebug(1,Form("species %i priors %f and prod %f",k,fPriors[k],prod[k]));
- comb[k]=fPriors[k]*prod[k];
- AliDebug(1,Form("comb %i %f",k,comb[k]));
- sum+=comb[k];
- }
+ 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++) 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]) Printf("no fhCombResp[%i], check if pidcut->Init() was called",n);
- fhCombProb[n]->Fill(combpid[n]);
- }
- }
- }
+ 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]));
+
}
//__________________________________________
//
}
}
//______________________________________________
-void AliCFTrackCutPid::Init()
-{
- //
- // initialises QA histograms
- //
-
- if(fIsQAOn) DefineHistograms();
-}
-//_________________________________________________
void AliCFTrackCutPid::DefineHistograms()
{
//
//QA histo booking
//
- char *detect[5]={"ITS","TPC","TRD","TOF","HMPID"};
- char *partic[5]={"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);
- fhCombProb[iPart] = new TH1F(Form("pCombPart%i",iPart),Form("%s combined probability ",partic[iPart]),fNbins,fXmin,fXmax);
- }
- }
+ 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()));
+ }
+ }
+ }
+
}
//___________________________________________________
-void AliCFTrackCutPid::AddQAHistograms(TList *qalist) const
+void AliCFTrackCutPid::AddQAHistograms(TList *qalist)
{
//
// adds QA histograms in a TList
//
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]);
}
for(Int_t iDet=0; iDet<kNdets; iDet++){
+ if(!fDets[iDet]) continue;
for(Int_t iP =0; iP<AliPID::kSPECIES; iP++){
- if(!fgIsComb)qalist->Add(fhResp[iDet][iP]);
- if(!fgIsComb)qalist->Add(fhProb[iDet][iP]);
+ qalist->Add(fhResp[iDet][iP]);
+ qalist->Add(fhProb[iDet][iP]);
}
}
+
}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff