#include "AliAODRecoDecayHF2Prong.h"
#include "AliAODTrack.h"
#include "AliESDtrack.h"
+#include "AliAODPid.h"
+#include "AliTPCPIDResponse.h"
+#include "AliAODVertex.h"
ClassImp(AliRDHFCutsD0toKpi)
//--------------------------------------------------------------------------
AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) :
-AliRDHFCuts(name)
+AliRDHFCuts(name),
+fRemoveDaughtersFromPrimary(kFALSE),
+fUseSpecialCuts(kFALSE)
{
//
// Default Constructor
SetVarsForOpt(4,forOpt);
Float_t limits[2]={0,999999999.};
SetPtBins(2,limits);
+
}
//--------------------------------------------------------------------------
AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
- AliRDHFCuts(source)
+ AliRDHFCuts(source),
+ fRemoveDaughtersFromPrimary(source.fRemoveDaughtersFromPrimary),
+ fUseSpecialCuts(source.fUseSpecialCuts)
{
//
// Copy constructor
if(&source == this) return *this;
AliRDHFCuts::operator=(source);
+ fRemoveDaughtersFromPrimary=source.fRemoveDaughtersFromPrimary;
+ fUseSpecialCuts=source.fUseSpecialCuts;
return *this;
}
return;
}
//---------------------------------------------------------------------------
-Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel) {
+Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod) {
//
// Apply selection
//
return 0;
}
-
// selection on daughter tracks
if(selectionLevel==AliRDHFCuts::kAll ||
selectionLevel==AliRDHFCuts::kTracks) {
}
+
// returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
Int_t returnvaluePID=3;
Int_t returnvalueCuts=3;
- Int_t returnvalue=3;
// selection on PID
if(selectionLevel==AliRDHFCuts::kAll ||
+ selectionLevel==AliRDHFCuts::kCandidate ||
selectionLevel==AliRDHFCuts::kPID) {
returnvaluePID = IsSelectedPID(d);
}
// selection on candidate
if(selectionLevel==AliRDHFCuts::kAll ||
selectionLevel==AliRDHFCuts::kCandidate) {
+
+ //recalculate vertex w/o daughters
+ AliAODVertex *origownvtx=0x0;
+ AliAODVertex *recvtx=0x0;
+ if(aod && fRemoveDaughtersFromPrimary) {
+ if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
+ recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
+ if(!recvtx){
+ AliDebug(2,"Removal of daughter tracks failed");
+ //recvtx=d->GetPrimaryVtx();
+ if(origownvtx){
+ delete origownvtx;
+ origownvtx=NULL;
+ }
+ return 0;
+ }
+ //set recalculed primary vertex
+ d->SetOwnPrimaryVtx(recvtx);
+ delete recvtx; recvtx=NULL;
+ }
+
Double_t pt=d->Pt();
Int_t okD0=0,okD0bar=0;
Int_t ptbin=PtBin(pt);
-
+ if (ptbin==-1) {
+ if(origownvtx){
+ d->SetOwnPrimaryVtx(origownvtx);
+ delete origownvtx;
+ origownvtx=NULL;
+ }
+ else d->UnsetOwnPrimaryVtx();
+ return 0;
+ }
Double_t mD0,mD0bar,ctsD0,ctsD0bar;
okD0=1; okD0bar=1;
if(d->PtProng(1) < fCutsRD[GetGlobalIndex(3,ptbin)] || d->PtProng(0) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
if(d->PtProng(0) < fCutsRD[GetGlobalIndex(3,ptbin)] || d->PtProng(1) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
-
- if(!okD0 && !okD0bar) return 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
- if(!okD0 && !okD0bar) return 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
- if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0;
+ if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) returnvalueCuts=0;
d->InvMassD0(mD0,mD0bar);
if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
- if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
- if(!okD0 && !okD0bar) return 0;
+ if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
d->CosThetaStarD0(ctsD0,ctsD0bar);
- if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
+ if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
- if(!okD0 && !okD0bar) return 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
- if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) return 0;
+ if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) returnvalueCuts=0;
- if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) return 0;
+ if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) returnvalueCuts=0;
- if (okD0) returnvalueCuts=1; //cuts passed as D0
- if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
- if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
- }
+ if (returnvalueCuts!=0) {
+ if (okD0) returnvalueCuts=1; //cuts passed as D0
+ if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
+ if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
+ }
+ // call special cuts
+ Int_t special=1;
+ if(fUseSpecialCuts) special=IsSelectedSpecialCuts(d);
+ if(!special) returnvalueCuts=0;
+
+ // unset recalculated primary vertex when not needed any more
+ if(origownvtx) {
+ d->SetOwnPrimaryVtx(origownvtx);
+ delete origownvtx;
+ origownvtx=NULL;
+ } else if(fRemoveDaughtersFromPrimary) {
+ d->UnsetOwnPrimaryVtx();
+ AliDebug(3,"delete new vertex\n");
+ }
- // combine PID and Cuts results
- switch(returnvaluePID) {
- case 0:
- returnvalue=0;
- break;
- case 1:
- returnvalue=((returnvalueCuts==1 || returnvalueCuts==3) ? 1 : 0);
- break;
- case 2:
- returnvalue=((returnvalueCuts==2 || returnvalueCuts==3) ? 2 : 0);
- break;
- case 3:
- returnvalue=returnvalueCuts;
- break;
- default:
- returnvalue=0;
- break;
}
- return returnvalue;
+
+
+ // cout<<"Pid = "<<returnvaluePID<<endl;
+ return CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
}
//---------------------------------------------------------------------------
Bool_t AliRDHFCutsD0toKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
if(pt > 5.) {
// applying cut for pt > 5 GeV
- AliDebug(4,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
+ AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
if (TMath::Abs(y) > 0.8){
return kFALSE;
}
// appliying smooth cut for pt < 5 GeV
Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
- AliDebug(4,Form("pt of D0 = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
+ AliDebug(2,Form("pt of D0 = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
if (y < minFiducialY || y > maxFiducialY){
return kFALSE;
}
return kTRUE;
}
//---------------------------------------------------------------------------
-Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* rd) const {
+Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d)
+{
+ // ############################################################
//
// Apply PID selection
- // (return: 0 not sel, 1 only D0, 2 only D0bar, 3 both)
//
+ //
+ // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)
+ //
+ // d must be a AliAODRecoDecayHF2Prong object
+ // returns 0 if both D0 and D0bar are rejectecd
+ // 1 if D0 is accepted while D0bar is rejected
+ // 2 if D0bar is accepted while D0 is rejected
+ // 3 if both are accepted
+ // fWhyRejection variable (not returned for the moment, print it if needed)
+ // keeps some information on why a candidate has been
+ // rejected according to the following (unfriendly?) scheme
+ // if more rejection cases are considered interesting, just add numbers
+ //
+ // TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant)
+ // from 20 to 30: "detector" selection (PID acceptance)
+ // 26: TPC refit
+ // 27: ITS refit
+ // 28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)
+ //
+ // from 30 to 40: PID selection
+ // 31: no Kaon compatible tracks found between daughters
+ // 32: no Kaon identified tracks found (strong sel. at low momenta)
+ // 33: both mass hypotheses are rejected
+ //
+ // ############################################################
+
+ if(!fUsePID) return 3;
+ fWhyRejection=0;
+ Int_t isD0D0barPID[2]={1,2};
+ Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid
+ Double_t tofSig,times[5];// used fot TOF pid
+ Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters
+ Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];
+ Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
+ // same for prong 2
+ // values convention -1 = discarded
+ // 0 = not identified (but compatible) || No PID (->hasPID flag)
+ // 1 = identified
+ // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
+ // Initial hypothesis: unknwon (but compatible)
+ isKaonPionTOF[0][0]=0;
+ isKaonPionTOF[0][1]=0;
+ isKaonPionTOF[1][0]=0;
+ isKaonPionTOF[1][1]=0;
+
+ isKaonPionTPC[0][0]=0;
+ isKaonPionTPC[0][1]=0;
+ isKaonPionTPC[1][0]=0;
+ isKaonPionTPC[1][1]=0;
+
+ combinedPID[0][0]=0;
+ combinedPID[0][1]=0;
+ combinedPID[1][0]=0;
+ combinedPID[1][1]=0;
+
+
+
+ for(Int_t daught=0;daught<2;daught++){
+ //Loop con prongs
+
+ // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################
+
+ AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
+
+ if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
+ fWhyRejection=26;
+ return 0;
+ }
+ if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
+ fWhyRejection=27;
+ return 0;
+ }
+
+ AliAODPid *pid=aodtrack->GetDetPid();
+ if(!pid) {
+ //delete esdtrack;
+ hasPID[daught]--;
+ continue;
+ }
+
+ // ########### Step 1- Check of TPC and TOF response ####################
+
+ Double_t ptrack=aodtrack->P();
+ //#################### TPC PID #######################
+ if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){
+ // NO TPC PID INFO FOR THIS TRACK
+ hasPID[daught]--;
+ }
+ else {
+ static AliTPCPIDResponse theTPCpid;
+ nsigmaTPCpi = theTPCpid.GetNumberOfSigmas((Float_t)aodtrack->P(),(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);
+ nsigmaTPCK = theTPCpid.GetNumberOfSigmas((Float_t)aodtrack->P(),(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);
+ if(ptrack<0.6){
+ if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
+ else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;
+ else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ else if(ptrack<.8){
+ if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;
+ else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;
+ else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ else {
+ // if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;
+ if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;
+ if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ }
+
+
+ // ##### TOF PID: do not ask nothing for pion/protons ############
+ if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){
+ // NO TOF PID INFO FOR THIS TRACK
+ hasPID[daught]--;
+ }
+ else{
+ tofSig=pid->GetTOFsignal();
+ pid->GetIntegratedTimes(times);
+ if(TMath::Abs(tofSig-times[3])>3.*160.){
+ isKaonPionTOF[daught][0]=-1;
+ }
+ else {
+ if(ptrack<1.5){
+ isKaonPionTOF[daught][0]=1;
+ }
+ }
+ }
+
+ //######### Step 2: COMBINE TOF and TPC PID ###############
+ // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown
+ combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];
+ combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];
+
+
+ //######### Step 3: USE PID INFO
+
+ if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded
+ isD0D0barPID[0]=0;
+ isD0D0barPID[1]=0;
+ }
+ else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-), if k for both TPC and TOF -> is K
+ if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
+ else isD0D0barPID[0]=0;// if K+ D0 excluded
+ }
+ else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-) and k- only for TPC or TOF -> reject
+ isD0D0barPID[0]=0;
+ isD0D0barPID[1]=0;
+ }
+ else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
+ if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
+ else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
+ }
+ else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
+ if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
+ else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
+ }
+
+ // ########## ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification ###############################
+ // ########## more tolerant criteria for single particle ID-> more selective criteria for D0 ##############################
+ // ############### NOT OPTIMIZED YET ###################################
+ if(d->Pt()<2.){
+ isKaonPionTPC[daught][0]=0;
+ isKaonPionTPC[daught][1]=0;
+ if(ptrack<0.6){
+ if(TMath::Abs(nsigmaTPCK)<3.)isKaonPionTPC[daught][0]=1;
+ else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
+ else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ else if(ptrack<.8){
+ if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
+ else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
+ else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ else {
+ if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ }
+
+ }// END OF LOOP ON DAUGHTERS
+
+ // FURTHER PID REQUEST (both daughter info is needed)
+ if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
+ fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
+ return 0;
+ }
+ else if(hasPID[0]==0&&hasPID[1]==0){
+ fWhyRejection=28;// reject cases in which no PID info is available
+ return 0;
+ }
+ if(d->Pt()<2.){
+ // request of K identification at low D0 pt
+ combinedPID[0][0]=0;
+ combinedPID[0][1]=0;
+ combinedPID[1][0]=0;
+ combinedPID[1][1]=0;
+
+ combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];
+ combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];
+ combinedPID[0][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];
+ combinedPID[0][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];
+
+ if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
+ fWhyRejection=32;// reject cases where the Kaon is not identified
+ return 0;
+ }
+ }
+
+ // cout<<"Why? "<<fWhyRejection<<endl;
+ return isD0D0barPID[0]+isD0D0barPID[1];
+}
+
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::CombineSelectionLevels(Int_t selectionvalTrack,
+ Int_t selectionvalCand,
+ Int_t selectionvalPID) const
+{
+ //
+ // This method combines the tracks, PID and cuts selection results
+ //
+ if(selectionvalTrack==0) return 0;
+
+ Int_t returnvalue;
- if(!fUsePID || !rd) return 3;
+ switch(selectionvalPID) {
+ case 0:
+ returnvalue=0;
+ break;
+ case 1:
+ returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);
+ break;
+ case 2:
+ returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);
+ break;
+ case 3:
+ returnvalue=selectionvalCand;
+ break;
+ default:
+ returnvalue=0;
+ break;
+ }
+
+ return returnvalue;
+}
+//----------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
+{
+ //
+ // Note: this method is temporary
+ // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters
+ //
- Int_t returnvalue=0;
+ //apply cuts
- // HERE LOOP ON DAUGHTERS, USE AliAODpidHF CLASS, ETC...
+ Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;
+ // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066
+ // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution),
+ Int_t returnvalue=3; //cut passed
+ for(Int_t i=0;i<2/*prongs*/;i++){
+ if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
+ }
+ if(d->DecayLength()<decLengthCut) return 0; //decLengthCut not passed
+ if(d->NormalizedDecayLength()<normDecLengthCut) return 0; //decLengthCut not passed
+
return returnvalue;
}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff