]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliRDHFCutsLctopKpi.cxx
Fix in PID application for LambdaC (Rossella)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCutsLctopKpi.cxx
index b27a60098a6fcd632c3343d54419e687c257e30d..150bc6534991e801a9c392c992e343cf43362031 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 /////////////////////////////////////////////////////////////
 //
 // Class for cuts on AOD reconstructed Lc->pKpi
@@ -43,12 +45,12 @@ fRecoKF(kFALSE)
   //
   // Default Constructor
   //
-  Int_t nvars=12;
+  Int_t nvars=13;
   SetNVars(nvars);
-  TString varNames[12]={"inv. mass [GeV]",
+  TString varNames[13]={"inv. mass [GeV]",
+                       "pTK [GeV/c]",
                        "pTP [GeV/c]",
-                       "pTPi [GeV/c]",
-                       "d0P [cm]   lower limit!",
+                       "d0K [cm]   lower limit!",
                        "d0Pi [cm]  lower limit!",
                        "dist12 (cm)",
                        "sigmavert (cm)",
@@ -56,8 +58,9 @@ fRecoKF(kFALSE)
                        "pM=Max{pT1,pT2,pT3} (GeV/c)",
                        "cosThetaPoint",
                        "Sum d0^2 (cm^2)",
-                       "dca cut (cm)"};
-  Bool_t isUpperCut[12]={kTRUE,
+                       "dca cut (cm)",
+                       "cut on pTpion [GeV/c]"};
+  Bool_t isUpperCut[13]={kTRUE,
                         kFALSE,
                         kFALSE,
                         kFALSE,
@@ -68,21 +71,24 @@ fRecoKF(kFALSE)
                         kFALSE,
                         kFALSE,
                         kFALSE,
-                        kTRUE};
+                        kTRUE,
+                        kFALSE
+                        };
   SetVarNames(nvars,varNames,isUpperCut);
-  Bool_t forOpt[12]={kFALSE,
-                    kFALSE,
+  Bool_t forOpt[13]={kFALSE,
+                    kTRUE,
+                    kTRUE,
                     kFALSE,
                     kFALSE,
                     kFALSE,
                     kFALSE,
                     kTRUE,
-                    kTRUE,
-                    kTRUE,
-                    kTRUE,
-                    kTRUE,
-                    kFALSE};
-  SetVarsForOpt(5,forOpt);
+                    kFALSE,
+                    kFALSE,
+                    kFALSE,
+                    kFALSE,
+                    kTRUE};
+  SetVarsForOpt(4,forOpt);
   Float_t limits[2]={0,999999999.};
   SetPtBins(2,limits);
 }
@@ -150,7 +156,7 @@ void AliRDHFCutsLctopKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,In
   if(fVarsForOpt[1]){
     iter++;
     for(Int_t iprong=0;iprong<3;iprong++){
-      if(TMath::Abs(pdgdaughters[iprong])==2212) {
+      if(TMath::Abs(pdgdaughters[iprong])==321) {
        vars[iter]=dd->PtProng(iprong);
       }
     }
@@ -158,7 +164,7 @@ void AliRDHFCutsLctopKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,In
   if(fVarsForOpt[2]){
     iter++;
     for(Int_t iprong=0;iprong<3;iprong++){
-      if(TMath::Abs(pdgdaughters[iprong])==211) {
+      if(TMath::Abs(pdgdaughters[iprong])==2212) {
        vars[iter]=dd->PtProng(iprong);
       }
     }
@@ -211,6 +217,14 @@ void AliRDHFCutsLctopKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,In
     iter++;
     vars[iter]=dd->GetDCA();
   }
+  if(fVarsForOpt[12]){
+    iter++;
+    for(Int_t iprong=0;iprong<3;iprong++){
+      if(TMath::Abs(pdgdaughters[iprong])==211) {
+       vars[iter]=dd->PtProng(iprong);
+      }
+    }
+  }
 
   return;
 }
@@ -233,18 +247,13 @@ Int_t AliRDHFCutsLctopKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEv
   }
 
 
-  // selection on daughter tracks 
-  if(selectionLevel==AliRDHFCuts::kAll || 
-     selectionLevel==AliRDHFCuts::kTracks) {
-    if(!AreDaughtersSelected(d)) return 0;
-  }
-
+  if(fKeepSignalMC) if(IsSignalMC(d,aod,4122)) return 3;
 
   Int_t returnvalue=3;
   Int_t returnvaluePID=3;
 
-     if(fUsePID || selectionLevel==AliRDHFCuts::kPID) returnvaluePID = IsSelectedPID(d);
-     if(returnvaluePID==0) return 0;
+  if(d->Pt()<fMinPtCand) return 0;
+  if(d->Pt()>fMaxPtCand) return 0;
 
   // selection on candidate
   if(selectionLevel==AliRDHFCuts::kAll || 
@@ -267,10 +276,8 @@ Int_t AliRDHFCutsLctopKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEv
     if(!okLcpKpi && !okLcpiKp) return 0;
 
     if(TMath::Abs(d->PtProng(1)) < fCutsRD[GetGlobalIndex(1,ptbin)] || TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]) return 0;//Kaon
-    //if(TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)] || TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)]) return 0;//Proton
-    //if(TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)] || TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)]) return 0;//Pion
-    if((TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)]) || (TMath::Abs(d->PtProng(2)) < 0.4)) okLcpKpi=0;
-    if((TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)]) || (TMath::Abs(d->PtProng(0)) < 0.4))okLcpiKp=0;
+    if((TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)]) || (TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(12,ptbin)])) okLcpKpi=0;
+    if((TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)]) || (TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(12,ptbin)]))okLcpiKp=0;
     if(!okLcpKpi && !okLcpiKp) return 0;
 
     
@@ -317,8 +324,25 @@ Int_t AliRDHFCutsLctopKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEv
     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
-    
+   
   }
+
+  if(selectionLevel==AliRDHFCuts::kAll ||
+     selectionLevel==AliRDHFCuts::kCandidate|| 
+     selectionLevel==AliRDHFCuts::kPID) {
+     returnvaluePID = IsSelectedPID(d);
+     fIsSelectedPID=returnvaluePID;
+     }
+  //  if(fUsePID || selectionLevel==AliRDHFCuts::kPID) returnvaluePID = IsSelectedCombinedPID(d);   // to test!!
+  if(returnvaluePID==0) return 0;
+
+  // selection on daughter tracks 
+  if(selectionLevel==AliRDHFCuts::kAll || 
+     selectionLevel==AliRDHFCuts::kTracks) {
+    if(!AreDaughtersSelected(d)) return 0;
+  }
+
+
   Int_t returnvalueTot=CombinePIDCuts(returnvalue,returnvaluePID);
   return returnvalueTot;
 }
@@ -326,19 +350,22 @@ Int_t AliRDHFCutsLctopKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEv
 Int_t AliRDHFCutsLctopKpi::IsSelectedPID(AliAODRecoDecayHF* obj) {
 
 
-    if(!obj) {return 3;}
+    if(!fUsePID || !obj) return 3;
     Int_t okLcpKpi=0,okLcpiKp=0;
     Int_t returnvalue=0;
     Bool_t isPeriodd=fPidHF->GetOnePad();
-    Bool_t isPbPb=fPidHF->GetPbPb();
+    Bool_t isMC=fPidHF->GetMC();
     Bool_t ispion0=kTRUE,ispion2=kTRUE;
     Bool_t isproton0=kFALSE,isproton2=kFALSE;
     Bool_t iskaon1=kFALSE;
-    Double_t sigmaTOF=120.;
-    if(isPeriodd) sigmaTOF=160.;
-    if(isPbPb) sigmaTOF=160.;
-    fPidObjprot->SetTofSigma(sigmaTOF);
-    fPidHF->SetTofSigma(sigmaTOF);
+    if(isPeriodd) {
+     fPidObjprot->SetOnePad(kTRUE);
+     fPidObjpion->SetOnePad(kTRUE);
+    }
+    if(isMC) {
+     fPidObjprot->SetMC(kTRUE);
+     fPidObjpion->SetMC(kTRUE);
+    }
 
     for(Int_t i=0;i<3;i++){
      AliAODTrack *track=(AliAODTrack*)obj->GetDaughter(i);
@@ -350,6 +377,7 @@ Int_t AliRDHFCutsLctopKpi::IsSelectedPID(AliAODRecoDecayHF* obj) {
        iskaon1=kTRUE;
       if(fPidHF->MakeRawPid(track,2)>=1) iskaon1=kFALSE;
       }
+      if(!iskaon1) return 0;
      
      }else{
      //pion or proton
@@ -370,6 +398,7 @@ Int_t AliRDHFCutsLctopKpi::IsSelectedPID(AliAODRecoDecayHF* obj) {
       if(isProton>=1) isproton0=kTRUE;
 
      }
+      if(!ispion0 && !isproton0) return 0;
      if(i==2) {
       if(isPion<0) ispion2=kFALSE;
       if(isProton>=1) isproton2=kTRUE;
@@ -386,6 +415,67 @@ Int_t AliRDHFCutsLctopKpi::IsSelectedPID(AliAODRecoDecayHF* obj) {
 
  return returnvalue;
 }
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsLctopKpi::IsSelectedCombinedPID(AliAODRecoDecayHF* obj) {
+
+  //  Printf(" -------- IsSelectedCombinedPID --------------");
+
+    if(!obj) {return 3;}
+    Int_t okLcpKpi=0,okLcpiKp=0;
+    Int_t returnvalue=0;
+    Bool_t isPeriodd=fPidHF->GetOnePad();
+    Bool_t isMC=fPidHF->GetMC();
+    Bool_t isKaon = kFALSE;
+    Bool_t isPion = kFALSE;
+    Bool_t isProton = kFALSE;
+
+    if(isPeriodd) {
+           fPidObjprot->SetOnePad(kTRUE);
+           fPidObjpion->SetOnePad(kTRUE);
+    }
+    if(isMC) {
+           fPidObjprot->SetMC(kTRUE);
+           fPidObjpion->SetMC(kTRUE);
+    }
+
+    Double_t probComb[AliPID::kSPECIES]={0.};  // array to store the info for the combined ITS|TPC|TOF probabilities
+    fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetITS|AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);  // the mask could become a member of the cut object
+
+    for(Int_t i=0;i<3;i++){
+           AliAODTrack *track=(AliAODTrack*)obj->GetDaughter(i);
+           if(!track) return 0;
+           // identify kaon
+           /*UInt_t detUsed =*/ fPidHF->GetPidCombined()->ComputeProbabilities(track, fPidHF->GetPidResponse(), probComb); // for the moment we don't check detUsed...
+           //      Printf("[%d] %x %f %f %f %f %f",i,detUsed,probComb[0],probComb[1],probComb[2],probComb[3],probComb[4]);
+           if(i==1) {
+                   if(TMath::MaxElement(AliPID::kSPECIES,probComb) == probComb[3]) {  // the probability to be a Kaon is the highest
+                     isKaon=kTRUE;
+                     //                      if ( probComb[3] > 0.8 ) isKaon=kTRUE;
+                     //                      else return 0;
+                   }
+                   else return 0; // prong at position 1 is not a Kaon, returning 0                
+           }
+           else {
+                   //pion or proton
+                   
+                   if(TMath::MaxElement(AliPID::kSPECIES,probComb) == probComb[4]) {  // the probability to be a proton is the highest
+                           isProton=kTRUE;
+                           if (isPion) okLcpiKp = 1;  // the pion was already identified, so here we must be at i == 2 --> Lc --> pi K p (K already found)
+                   }
+                   if(TMath::MaxElement(AliPID::kSPECIES,probComb) == probComb[2]) {  // the probability to be a pion is the highest
+                           isPion=kTRUE;
+                           if (isProton) okLcpKpi = 1;  // the proton was already identified, so here we must be at i == 2 --> Lc --> p K pi (K already found)
+                   }
+                   
+           }
+    }
+    
+    if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
+    if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
+    if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
+    
+    return returnvalue;
+}
 //-----------------------
 Int_t AliRDHFCutsLctopKpi::CombinePIDCuts(Int_t returnvalue, Int_t returnvaluePID) const {
 
@@ -427,7 +517,7 @@ void AliRDHFCutsLctopKpi::SetStandardCutsPP2010() {
  AddTrackCuts(esdTrackCuts);
 
  const Int_t nptbins=4;
- const Int_t nvars=12;
+ const Int_t nvars=13;
  Float_t* ptbins;
  ptbins=new Float_t[nptbins+1];
  
@@ -459,6 +549,7 @@ void AliRDHFCutsLctopKpi::SetStandardCutsPP2010() {
   prodcutsval[9][ipt]=0.;
   prodcutsval[10][ipt]=0.;
   prodcutsval[11][ipt]=0.05;
+  prodcutsval[12][ipt]=0.4;
  }
  SetCuts(nvars,nptbins,prodcutsval);
 
@@ -472,6 +563,7 @@ void AliRDHFCutsLctopKpi::SetStandardCutsPP2010() {
  pidObjK->SetITS(kTRUE);
  Double_t plimK[2]={0.5,0.8};
  pidObjK->SetPLimit(plimK,2);
+ pidObjK->SetTOFdecide(kTRUE);
 
  SetPidHF(pidObjK);
 
@@ -479,6 +571,7 @@ void AliRDHFCutsLctopKpi::SetStandardCutsPP2010() {
  pidObjpi->SetTPC(kTRUE);
  Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
  pidObjpi->SetSigma(sigmaspi);
+ pidObjpi->SetTOFdecide(kTRUE);
  SetPidpion(pidObjpi);
 
  AliAODPidHF* pidObjp=new AliAODPidHF();
@@ -491,6 +584,7 @@ void AliRDHFCutsLctopKpi::SetStandardCutsPP2010() {
  pidObjp->SetITS(kTRUE);
  Double_t plimp[2]={1.,2.};
  pidObjp->SetPLimit(plimp,2);
+ pidObjp->SetTOFdecide(kTRUE);
 
  SetPidprot(pidObjp);
 
@@ -511,7 +605,113 @@ void AliRDHFCutsLctopKpi::SetStandardCutsPP2010() {
 
  return;
 }
+//------------------
+void AliRDHFCutsLctopKpi::SetStandardCutsPbPb2010() {
+
+ SetName("LctopKpiProdCuts");
+ SetTitle("Production cuts for Lc analysis");
 
+ AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
+
+ esdTrackCuts->SetRequireTPCRefit(kTRUE);
+ esdTrackCuts->SetMinNClustersTPC(70);
+ esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
+                                          AliESDtrackCuts::kAny);
+ esdTrackCuts->SetRequireITSRefit(kTRUE);
+ esdTrackCuts->SetMinNClustersITS(4);
+ esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0100*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
+ esdTrackCuts->SetEtaRange(-0.8,0.8);
+ esdTrackCuts->SetMaxDCAToVertexXY(1.);
+ esdTrackCuts->SetMaxDCAToVertexZ(1.);
+ esdTrackCuts->SetPtRange(0.8,1.e10);
+ AddTrackCuts(esdTrackCuts);
+
+ const Int_t nptbins=4;
+ const Int_t nvars=13;
+ Float_t* ptbins;
+ ptbins=new Float_t[nptbins+1];
+ ptbins[0]=0.;
+ ptbins[1]=2.;
+ ptbins[2]=3.;
+ ptbins[3]=4.;
+ ptbins[4]=9999.;
+
+ SetGlobalIndex(nvars,nptbins);
+ SetPtBins(nptbins+1,ptbins);
+
+ Float_t** prodcutsval;
+ prodcutsval=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++){
+  prodcutsval[iv]=new Float_t[nptbins];
+ }
+
+ for(Int_t ipt=0;ipt<nptbins;ipt++){
+  prodcutsval[0][ipt]=0.13;
+  prodcutsval[1][ipt]=0.5;
+  prodcutsval[2][ipt]=0.6;
+  prodcutsval[3][ipt]=0.;
+  prodcutsval[4][ipt]=0.;
+  prodcutsval[5][ipt]=0.01;
+  prodcutsval[6][ipt]=0.04;
+  prodcutsval[7][ipt]=0.006;
+  prodcutsval[8][ipt]=0.8;
+  prodcutsval[9][ipt]=0.3;
+  prodcutsval[10][ipt]=0.;
+  prodcutsval[11][ipt]=0.05;
+  prodcutsval[12][ipt]=0.4;
+ }
+ SetCuts(nvars,nptbins,prodcutsval);
+
+ AliAODPidHF* pidObjK=new AliAODPidHF();
+ Double_t sigmasK[5]={3.,1.,1.,3.,2.};
+ pidObjK->SetSigma(sigmasK);
+ pidObjK->SetAsym(kTRUE);
+ pidObjK->SetMatch(1);
+ pidObjK->SetTPC(kTRUE);
+ pidObjK->SetTOF(kTRUE);
+ pidObjK->SetITS(kTRUE);
+ Double_t plimK[2]={0.5,0.8};
+ pidObjK->SetPLimit(plimK,2);
+
+ SetPidHF(pidObjK);
+
+ AliAODPidHF* pidObjpi=new AliAODPidHF();
+ pidObjpi->SetTPC(kTRUE);
+ Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
+ pidObjpi->SetSigma(sigmaspi);
+ SetPidpion(pidObjpi);
+
+ AliAODPidHF* pidObjp=new AliAODPidHF();
+ Double_t sigmasp[5]={3.,1.,1.,3.,2.};
+ pidObjp->SetSigma(sigmasp);
+ pidObjp->SetAsym(kTRUE);
+ pidObjp->SetMatch(1);
+ pidObjp->SetTPC(kTRUE);
+ pidObjp->SetTOF(kTRUE);
+ pidObjp->SetITS(kTRUE);
+ Double_t plimp[2]={1.,2.};
+ pidObjp->SetPLimit(plimp,2);
+
+ SetPidprot(pidObjp);
+
+ SetUsePID(kTRUE);
+
+ PrintAll();
+
+ for(Int_t iiv=0;iiv<nvars;iiv++){
+  delete [] prodcutsval[iiv];
+ }
+ delete [] prodcutsval;
+ prodcutsval=NULL;
+ delete [] ptbins;
+ ptbins=NULL;
+
+ delete pidObjp;
+ pidObjp=NULL;
+
+ return;
+}
 //------------------
 Bool_t AliRDHFCutsLctopKpi::ReconstructKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field) const{
 
@@ -536,3 +736,4 @@ Bool_t AliRDHFCutsLctopKpi::ReconstructKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs
  d->SetSecondaryVtx(vertexAOD);
  return kTRUE;
 }
+