/************************************************************************* * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /////////////////////////////////////////////////////////////////////////// // Dielectron PID // // // // // /* Detailed description */ // // /////////////////////////////////////////////////////////////////////////// #include #include #include #include #include #include #include #include #include #include //!!!!! Remove once Eta correction is treated in the tender #include "AliDielectronVarManager.h" #include "AliDielectronPID.h" ClassImp(AliDielectronPID) TGraph *AliDielectronPID::fgFitCorr=0x0; Double_t AliDielectronPID::fgCorr=0.0; Double_t AliDielectronPID::fgCorrdEdx=1.0; TF1 *AliDielectronPID::fgFunEtaCorr=0x0; TGraph *AliDielectronPID::fgdEdxRunCorr=0x0; AliDielectronPID::AliDielectronPID() : AliAnalysisCuts(), fNcuts(0), fPIDResponse(0x0) { // // Default Constructor // for (Int_t icut=0; icut=kNmaxPID){ AliError(Form("only %d pid cut ranges allowed",kNmaxPID)); return; } if (TMath::Abs(nSigmaUp+99999.)<1e-20){ nSigmaUp=TMath::Abs(nSigmaLow); nSigmaLow=-1.*nSigmaUp; } fDetType[fNcuts]=det; fPartType[fNcuts]=type; fNsigmaLow[fNcuts]=nSigmaLow; fNsigmaUp[fNcuts]=nSigmaUp; fmin[fNcuts]=min; fmax[fNcuts]=max; fExclude[fNcuts]=exclude; fRequirePIDbit[fNcuts]=pidBitType; fActiveCuts[fNcuts]=(var==-1 ? AliDielectronVarManager::kP : var); AliInfo(Form("Add PID cut %d: sigma [% .1f,% .1f] \t cut [% .1f,% .f] \t var %d->%s \n", fNcuts,nSigmaLow,nSigmaUp,min,max,fActiveCuts[fNcuts],AliDielectronVarManager::GetValueName(fActiveCuts[fNcuts]))); ++fNcuts; } //______________________________________________ void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, TF1 * const funUp, Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/, UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/) { // // cut using a TF1 as upper cut // if (funUp==0x0){ AliError("A valid function is required for the upper cut. Not adding the cut!"); return; } fFunUpperCut[fNcuts]=funUp; AddCut(det,type,nSigmaLow,0.,min,max,exclude,pidBitType,var); } //______________________________________________ void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, Double_t nSigmaUp, Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/, UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/) { // // cut using a TF1 as lower cut // if (funLow==0x0){ AliError("A valid function is required for the lower cut. Not adding the cut!"); return; } fFunLowerCut[fNcuts]=funLow; AddCut(det,type,0.,nSigmaUp,min,max,exclude,pidBitType,var); } //______________________________________________ void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, TF1 * const funUp, Double_t min/*=0*/, Double_t max/*=0*/, Bool_t exclude/*=kFALSE*/, UInt_t pidBitType/*=AliDielectronPID::kRequire*/, Int_t var/*=-1*/) { // // cut using a TF1 as lower and upper cut // if ( (funUp==0x0) || (funLow==0x0) ){ AliError("A valid function is required for upper and lower cut. Not adding the cut!"); return; } fFunUpperCut[fNcuts]=funUp; fFunLowerCut[fNcuts]=funLow; AddCut(det,type,0.,0.,min,max,exclude,pidBitType,var); } //______________________________________________ void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp, Double_t min, Double_t max, Bool_t exclude, UInt_t pidBitType, TF1 * const funSigma) { // // cut using a TF1 as lower cut // if (funSigma==0x0){ AliError("A valid function is required for the sigma cut. Not adding the cut!"); return; } fFunSigma[fNcuts]=funSigma; fSigmaFunLow[fNcuts]=min; fSigmaFunUp[fNcuts]=max; AddCut(det,type,nSigmaLow,nSigmaUp,0,0,exclude,pidBitType,-1); } //______________________________________________ Bool_t AliDielectronPID::IsSelected(TObject* track) { // // perform PID cuts // //loop over all cuts AliVTrack *part=static_cast(track); AliESDtrack *esdTrack=0x0; Double_t origdEdx=-1; // apply ETa correction, remove once this is in the tender if( (part->IsA() == AliESDtrack::Class()) && fgFunEtaCorr ){ esdTrack=static_cast(part); origdEdx=esdTrack->GetTPCsignal(); esdTrack->SetTPCsignal(origdEdx/GetEtaCorr(esdTrack)/fgCorrdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN()); } //Fill values Double_t values[AliDielectronVarManager::kNMaxValues]; AliDielectronVarManager::Fill(track,values); Bool_t selected=kFALSE; fPIDResponse=AliDielectronVarManager::GetPIDResponse(); for (UChar_t icut=0; icut1e-20) && (val=max) ) continue; // check if fFunSigma is set, then check if 'part' is in sigma range of the function if(fFunSigma[icut]){ val= fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]); if (fPartType[icut]==AliPID::kElectron){ val-=fgCorr; } min= fFunSigma[icut]->Eval(part->GetTPCmomentum())+fSigmaFunLow[icut]; max= fFunSigma[icut]->Eval(part->GetTPCmomentum())+fSigmaFunUp[icut]; if(val=max) continue; } switch (fDetType[icut]){ case kITS: selected = IsSelectedITS(part,icut); break; case kTPC: selected = IsSelectedTPC(part,icut); break; case kTRD: selected = IsSelectedTRD(part,icut); break; case kTRDeleEff: selected = IsSelectedTRDeleEff(part,icut); break; case kTRDeleEff2D: selected = IsSelectedTRDeleEff(part,icut,AliTRDPIDResponse::kLQ2D); break; case kTOF: selected = IsSelectedTOF(part,icut); break; case kEMCAL: selected = IsSelectedEMCAL(part,icut); break; } if (!selected) { if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN()); return kFALSE; } } if (esdTrack) esdTrack->SetTPCsignal(origdEdx,esdTrack->GetTPCsignalSigma(),esdTrack->GetTPCsignalN()); return selected; } //______________________________________________ Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut) { // // ITS part of the PID check // Don't accept the track if there was no pid bit set // if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kITSpid)) return kFALSE; if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kITSpid)) return kTRUE; Double_t mom=part->P(); Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasITS(part, fPartType[icut]); // test if we are supposed to use a function for the cut if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom); if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom); Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut]; return selected; } //______________________________________________ Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut) { // // TPC part of the PID check // Don't accept the track if there was no pid bit set // if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTPCpid)) return kFALSE; if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTPCpid)) return kTRUE; Double_t mom=part->GetTPCmomentum(); Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTPC(part, fPartType[icut]); if (fPartType[icut]==AliPID::kElectron){ numberOfSigmas-=fgCorr; } // test if we are supposed to use a function for the cut if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom); if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom); Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut]; return selected; } //______________________________________________ Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const part, Int_t icut) { // // TRD part of the pid check // the TRD checks on the probabilities. // if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kFALSE; if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kTRUE; if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable && (part->GetTRDntrackletsPID()<4)) return kTRUE; Double_t p[AliPID::kSPECIES]={0.}; fPIDResponse->ComputeTRDProbability(part,AliPID::kSPECIES,p); Float_t particleProb=p[fPartType[icut]]; Bool_t selected=((particleProb>=fNsigmaLow[icut])&&(particleProb<=fNsigmaUp[icut]))^fExclude[icut]; return selected; } //______________________________________________ Bool_t AliDielectronPID::IsSelectedTRDeleEff(AliVTrack * const part, Int_t icut, AliTRDPIDResponse::ETRDPIDMethod PIDmethod) { // // TRD part of the pid check using electron efficiency requirement // in this case the upper limit as well as the particle specie is ignored // and the lower limit regarded as the requested electron efficiency // if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kFALSE; if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTRDpid)) return kTRUE; Double_t centrality = -1.; if(part->IsA() == AliESDtrack::Class()) centrality=(const_cast( (static_cast(part))->GetESDEvent()) )->GetCentrality()->GetCentralityPercentile("V0M"); if(part->IsA() == AliAODTrack::Class()) centrality=(const_cast( (static_cast(part))->GetAODEvent()) )->GetCentrality()->GetCentralityPercentile("V0M"); Bool_t selected=fPIDResponse->IdentifiedAsElectronTRD(part,fNsigmaLow[icut], centrality, PIDmethod)^fExclude[icut]; return selected; } //______________________________________________ Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut) { // // TOF part of the PID check // Don't accept the track if there was no pid bit set // if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliVTrack::kTOFpid)) return kFALSE; if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliVTrack::kTOFpid)) return kTRUE; Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasTOF(part, fPartType[icut]); Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut]; return selected; } //______________________________________________ Bool_t AliDielectronPID::IsSelectedEMCAL(AliVTrack * const part, Int_t icut) { // // emcal pid selecttion // //TODO: correct way to check for emcal pid? Float_t numberOfSigmas=fPIDResponse->NumberOfSigmasEMCAL(part, fPartType[icut]); Bool_t hasPID=numberOfSigmas>-998.; if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!hasPID) return kFALSE; if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!hasPID) return kTRUE; Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut]; return selected; } //______________________________________________ void AliDielectronPID::SetDefaults(Int_t def){ // // initialise default pid strategies // if (def==0){ // 2sigma bands TPC: // - include e // - exclude mu,K,pi,p // -complete p range AddCut(kTPC,AliPID::kElectron,2); AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,0.,kTRUE); AddCut(kTPC,AliPID::kPion,-2.,2.,0.,0.,kTRUE); AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,0.,kTRUE); AddCut(kTPC,AliPID::kProton,-2.,2.,0.,0.,kTRUE); } else if (def==1) { // 2sigma bands TPC: // - include e 0SetParameters(-2.7,-0.4357); AddCut(kTPC,AliPID::kElectron,lowerCut,3.); AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5); } else if (def==7) { // wide TPC cut // TOF ele band 3sigma 0SetParameters(-2.65,-0.6757); AddCut(kTPC,AliPID::kElectron,lowerCut,4.); AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable); } else if (def==10) { AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable); AddCut(kTPC,AliPID::kElectron,3.); AddCut(kTPC,AliPID::kPion,-3.,3.,0.,0.,kTRUE); AddCut(kTPC,AliPID::kProton,-3.,3.,0.,0.,kTRUE); } else if (def==11) { // lower cut TPC: parametrisation by HFE // only use from period d on !! // upper cut TPC: 3 sigma // TOF ele band 3sigma 0SetParameters(-3.718,-0.4,0.27); AddCut(kTPC,AliPID::kElectron,lowerCut,3.); AddCut(kTOF,AliPID::kElectron,-3,3,0,2.); } else if (def==12) { // lower cut TPC: parametrisation by HFE // only use from period d on !! // upper cut TPC: 3 sigma // TOF 5 sigma inclusion if TOFpid available // this should reduce K,p,Pi to a large extent TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)+[2]", 0, 100); lowerCut->SetParameters(-3.718,-0.4,0.27); AddCut(kTPC,AliPID::kElectron,lowerCut,4.); AddCut(kTOF,AliPID::kElectron,-5,5,0,200,kFALSE,AliDielectronPID::kIfAvailable); } } //______________________________________________ void AliDielectronPID::SetCorrVal(Double_t run) { // // set correction value for run // fgCorr=0.; fgCorrdEdx=1.; if (fgFitCorr){ fgCorr=fgFitCorr->Eval(run); if (runGetX()[0]) fgCorr=fgFitCorr->GetY()[0]; if (run>fgFitCorr->GetX()[fgFitCorr->GetN()-1]) fgCorr=fgFitCorr->GetY()[fgFitCorr->GetN()-1]; } if (fgdEdxRunCorr){ fgCorrdEdx=fgdEdxRunCorr->Eval(run); if (runGetX()[0]) fgCorrdEdx=fgdEdxRunCorr->GetY()[0]; if (run>fgdEdxRunCorr->GetX()[fgFitCorr->GetN()-1]) fgCorrdEdx=fgdEdxRunCorr->GetY()[fgdEdxRunCorr->GetN()-1]; } } //______________________________________________ Double_t AliDielectronPID::GetEtaCorr(AliVTrack *track) { // // return eta correction // if (!fgFunEtaCorr) return 1; return fgFunEtaCorr->Eval(track->Eta()); }