Option to use Bayesian PID added for D0 (Jeremy)
authorfprino <fprino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 7 Apr 2013 20:58:52 +0000 (20:58 +0000)
committerfprino <fprino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 7 Apr 2013 20:58:52 +0000 (20:58 +0000)
PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.cxx
PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.h
PWGHF/vertexingHF/AliRDHFCutsD0toKpi.cxx
PWGHF/vertexingHF/AliRDHFCutsD0toKpi.h
PWGHF/vertexingHF/macros/AddTaskD0Mass.C

index 41d0ad7..e64ac6a 100644 (file)
@@ -23,6 +23,7 @@
 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
 // Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
 // Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
+// Jeremy Wilkinson, jwilkinson@physi.uni-heidelberg.de (weighted Bayesian)
 /////////////////////////////////////////////////////////////
 
 #include <Riostream.h>
@@ -85,7 +86,10 @@ AliAnalysisTaskSE(),
   fUseSelectionBit(kTRUE),
   fWriteVariableTree(kFALSE),
   fVariablesTree(0),
-  fCandidateVariables()
+  fCandidateVariables(),
+  fPIDCheck(kFALSE),
+  fDrawDetSignal(kFALSE),
+  fDetSignal(0)
 {
   // Default constructor
   for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
@@ -118,7 +122,10 @@ AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0t
   fUseSelectionBit(kTRUE),
   fWriteVariableTree(kFALSE),
   fVariablesTree(0),
-  fCandidateVariables()
+  fCandidateVariables(),
+  fPIDCheck(kFALSE),
+  fDrawDetSignal(kFALSE),
+  fDetSignal(0)
 {
   // Default constructor
 
@@ -141,6 +148,8 @@ AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0t
   DefineOutput(6,TList::Class());  //My private output
   // Output slot #7 keeps a tree of the candidate variables after track selection
   DefineOutput(7,TTree::Class());  //My private outpu
+  // Output slot #8 writes into a TList container (Detector signals)
+  DefineOutput(8, TList::Class()); //My private output
 }
 
 //________________________________________________________________________
@@ -177,6 +186,10 @@ AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
     delete fVariablesTree;
     fVariablesTree = 0;
   }
+  if (fDetSignal) {
+    delete fDetSignal;
+    fDetSignal = 0;
+  }
  
 }  
 
@@ -219,6 +232,10 @@ void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
   fDistr->SetOwner();
   fDistr->SetName("distributionslist");
 
+  fDetSignal = new TList();
+  fDetSignal->SetOwner();
+  fDetSignal->SetName("detectorsignals");
+
   TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
   TString nameMassPt="", nameSgnPt="", nameBkgPt="", nameRflPt="";
   Int_t nbins2dPt=60; Float_t binInPt=0., binFinPt=30.;
@@ -708,7 +725,7 @@ void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
   fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
   fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
   if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
-  if(fFillVarHists){
+  if(fFillVarHists || fPIDCheck){
     fNentries->GetXaxis()->SetBinLabel(8,"PID=0");
     fNentries->GetXaxis()->SetBinLabel(9,"PID=1");
     fNentries->GetXaxis()->SetBinLabel(10,"PID=2");
@@ -736,16 +753,42 @@ void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
   Int_t nVar = 15;
   fCandidateVariables = new Double_t [nVar];
   TString * fCandidateVariableNames = new TString[nVar];
-  fCandidateVariableNames[0]="massD0"; fCandidateVariableNames[1]="massD0bar"; fCandidateVariableNames[2]="pt";
-  fCandidateVariableNames[3]="dca"; fCandidateVariableNames[4]="costhsD0"; fCandidateVariableNames[5]="costhsD0bar";
-  fCandidateVariableNames[6]="ptk"; fCandidateVariableNames[7]="ptpi";
-  fCandidateVariableNames[8]="d0k"; fCandidateVariableNames[9]="d0pi"; fCandidateVariableNames[10]="d0xd0";
-  fCandidateVariableNames[11]="costhp"; fCandidateVariableNames[12]="costhpxy"; 
-  fCandidateVariableNames[13]="lxy";  fCandidateVariableNames[14]="specialcuts";
+  fCandidateVariableNames[0] = "massD0";
+  fCandidateVariableNames[1] = "massD0bar";
+  fCandidateVariableNames[2] = "pt";
+  fCandidateVariableNames[3] = "dca";
+  fCandidateVariableNames[4] = "costhsD0";
+  fCandidateVariableNames[5] = "costhsD0bar";
+  fCandidateVariableNames[6] = "ptk";
+  fCandidateVariableNames[7] = "ptpi";
+  fCandidateVariableNames[8] = "d0k";
+  fCandidateVariableNames[9] = "d0pi";
+  fCandidateVariableNames[10] = "d0xd0";
+  fCandidateVariableNames[11] = "costhp";
+  fCandidateVariableNames[12] = "costhpxy";
+  fCandidateVariableNames[13] = "lxy";
+  fCandidateVariableNames[14] = "specialcuts";
   for(Int_t ivar=0; ivar<nVar; ivar++){
     fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/d",fCandidateVariableNames[ivar].Data()));
   }
 
+  //
+  // Output slot 8 : List for detector response histograms
+  //
+  fDetSignal = new TList();
+  if (fDrawDetSignal) {
+    TH2F *TOFSigBefPID = new TH2F("TOFSigBefPID", "TOF signal of daughters before PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
+    TH2F *TOFSigAftPID = new TH2F("TOFSigAftPID", "TOF signal after PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
+
+    TH2F *TPCSigBefPID = new TH2F("TPCSigBefPID", "TPC dE/dx before PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
+    TH2F *TPCSigAftPID = new TH2F("TPCSigAftPID", "TPC dE/dx after PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
+
+    fDetSignal->Add(TOFSigBefPID);
+    fDetSignal->Add(TOFSigAftPID);
+    fDetSignal->Add(TPCSigBefPID);
+    fDetSignal->Add(TPCSigAftPID);
+  }
+
   // Post the data
   PostData(1,fOutputMass);
   PostData(2,fDistr);
@@ -753,6 +796,7 @@ void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
   PostData(5,fCounter);  
   PostData(6,fOutputMassPt);
   PostData(7,fVariablesTree);
+  PostData(8, fDetSignal);
 
   return;
 }
@@ -934,7 +978,21 @@ void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
        FillVarHists(aod,d,mcArray,fCuts,fDistr);
       }
 
+      if (fDrawDetSignal) {
+       DrawDetSignal(d, fDetSignal);
+      }
+
       FillMassHists(d,mcArray,mcHeader,fCuts,fOutputMass);
+      if (fPIDCheck) {
+       Int_t isSelectedPIDfill = 3;
+       if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(d); //0 rejected,1 D0,2 Dobar, 3 both
+
+       if (isSelectedPIDfill == 0)fNentries->Fill(7);
+       if (isSelectedPIDfill == 1)fNentries->Fill(8);
+       if (isSelectedPIDfill == 2)fNentries->Fill(9);
+       if (isSelectedPIDfill == 3)fNentries->Fill(10);
+      }
+
     }
   
     fDaughterTracks.Clear();
@@ -949,6 +1007,45 @@ void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
   PostData(5,fCounter);
   PostData(6,fOutputMassPt);
   PostData(7,fVariablesTree);
+  PostData(8, fDetSignal);
+
+  return;
+}
+
+//____________________________________________________________________________
+void AliAnalysisTaskSED0Mass::DrawDetSignal(AliAODRecoDecayHF2Prong *part, TList *ListDetSignal)
+{
+  //
+  // Function called in UserExec for drawing detector signal histograms:
+  //
+  fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(0), 0);
+  fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(1), 1);
+
+  AliESDtrack *esdtrack1 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
+  AliESDtrack *esdtrack2 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
+
+
+  //For filling detector histograms
+  Int_t isSelectedPIDfill = 3;
+  if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
+
+
+  //fill "before PID" histos with every daughter
+  ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
+  ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
+  ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
+  ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
+
+  if (isSelectedPIDfill != 0)  { //fill "After PID" with everything that's not rejected
+    ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
+    ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
+    ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
+    ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
+
+  }
+
+  delete esdtrack1;
+  delete esdtrack2;
 
   return;
 }
@@ -966,11 +1063,13 @@ void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Pr
   //Double_t pt = d->Pt(); //mother pt
   Int_t isSelectedPID=3;
   if(!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
-  if (isSelectedPID==0)fNentries->Fill(7);
-  if (isSelectedPID==1)fNentries->Fill(8);
-  if (isSelectedPID==2)fNentries->Fill(9);
-  if (isSelectedPID==3)fNentries->Fill(10);
-  //fNentries->Fill(8+isSelectedPID);
+  if (!fPIDCheck) {  //if fPIDCheck, this is already filled elsewhere
+    if (isSelectedPID==0)fNentries->Fill(7);
+    if (isSelectedPID==1)fNentries->Fill(8);
+    if (isSelectedPID==2)fNentries->Fill(9);
+    if (isSelectedPID==3)fNentries->Fill(10);
+    //fNentries->Fill(8+isSelectedPID);
+  }
 
   if(fCutOnDistr && !fIsSelectedCandidate) return; 
   //printf("\nif no cuts or cuts passed\n");
@@ -1575,6 +1674,15 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon
   Bool_t isPrimary=kTRUE;
   if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
 
+  //Define weights for Bayesian (if appropriate)
+
+  Double_t weigD0=1.;
+  Double_t weigD0bar=1.;
+  if (fCuts->GetCombPID() && (fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeight || fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeightNoFilter)) {
+    weigD0=fCuts->GetWeightsNegative()[AliPID::kKaon] * fCuts->GetWeightsPositive()[AliPID::kPion];
+    weigD0bar=fCuts->GetWeightsPositive()[AliPID::kKaon] * fCuts->GetWeightsNegative()[AliPID::kPion];
+  }
+  
   //count candidates selected by cuts
   fNentries->Fill(1);
   //count true D0 selected by cuts
@@ -1601,17 +1709,17 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon
          //      cout<<"Fill S with D0"<<endl;
          fillthis="histSgn_";
          fillthis+=ptbin;
-         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
 
          if(fFillPtHist){ 
            fillthismasspt="histSgnPt";
-           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt);
+           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
          }
          if(fFillImpParHist){ 
-           if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
+           if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0);
            else {
-             fHistMassPtImpParTC[2]->Fill(arrayForSparse);
-             fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
+             fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0);
+             fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0);
            }
          }
 
@@ -1619,32 +1727,32 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon
            if(TMath::Abs(invmassD0 - mPDG) < 0.027 && fFillVarHists){
              fillthis="histSgn27_";
              fillthis+=ptbin;
-             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
            }
          }
        } else{ //it was a D0bar
          fillthis="histRfl_";
          fillthis+=ptbin;
-         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
 
          if(fFillPtHist){ 
            fillthismasspt="histRflPt";
            //      cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;            
-           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt);
+           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
          }
 
        }
       } else {//background
        fillthis="histBkg_";
        fillthis+=ptbin;
-       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
 
        if(fFillPtHist){ 
          fillthismasspt="histBkgPt";
          //      cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
-         ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt);
+         ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
        }
-       if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse);
+       if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0);
 
       }
 
@@ -1654,16 +1762,17 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon
       //      cout<<"Filling "<<fillthis<<endl;
 
       //      printf("Fill mass with D0");
-      ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
+       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
+      
 
       if(fFillPtHist){ 
        fillthismasspt="histMassPt";
        //      cout<<"Filling "<<fillthismasspt<<endl;
-       ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt);
+       ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
       }
       if(fFillImpParHist) {
        //      cout << "Filling fHistMassPtImpParTC[0]"<<endl;
-       fHistMassPtImpParTC[0]->Fill(arrayForSparse);
+       fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0);
       }
 
     }
@@ -1688,7 +1797,7 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon
        if (pdgD0==-421){ //D0bar
          fillthis="histSgn_";
          fillthis+=ptbin;
-         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
+         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
          // if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
          //   fillthis="histSgn27_";
          //   fillthis+=ptbin;
@@ -1698,52 +1807,54 @@ void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClon
          if(fFillPtHist){ 
            fillthismasspt="histSgnPt";
            //      cout<<" Filling "<< fillthismasspt << endl;
-           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt);
+           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
          }
          if(fFillImpParHist){ 
            //      cout << " Filling impact parameter thnsparse"<<endl;
-           if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
+           if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0bar);
            else {
-             fHistMassPtImpParTC[2]->Fill(arrayForSparse);
-             fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
+             fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0bar);
+             fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0bar);
            }
          }
          
        } else{
          fillthis="histRfl_";
          fillthis+=ptbin;
-         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
+         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
          if(fFillPtHist){ 
            fillthismasspt="histRflPt";
            //      cout << " Filling "<<fillthismasspt<<endl;
-           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt);
+           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
          }
        }
       } else {//background or LS
        fillthis="histBkg_";
        fillthis+=ptbin;
-       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
+       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
 
        if(fFillPtHist){ 
          fillthismasspt="histBkgPt";
          //      cout<<" Filling "<< fillthismasspt << endl;
-         ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt);
+         ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
        }
-       if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse);
+       if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0bar);
 
       }
     }else{
       fillthis="histMass_";
       fillthis+=ptbin;
       //      printf("Fill mass with D0bar");
-      ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
+
+       ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar,weigD0bar);
+      
 
       if(fFillPtHist){ 
        fillthismasspt="histMassPt";
        //      cout<<" Filling "<< fillthismasspt << endl;
-       ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt);
+       ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
       }
-      if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse);
+      if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0bar);
 
     }
   }
@@ -1844,6 +1955,13 @@ void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
     printf("ERROR: fCounter not available\n");
     return;
   }
+  if (fDrawDetSignal) {
+    fDetSignal = dynamic_cast<TList*>(GetOutputData(8));
+    if (!fDetSignal) {
+      printf("ERROR: fDetSignal not available\n");
+      return;
+    }
+  }
   Int_t nptbins=fCuts->GetNPtBins();
   for(Int_t ipt=0;ipt<nptbins;ipt++){ 
 
index 1e2caf4..d0b127a 100644 (file)
@@ -58,6 +58,9 @@ class AliAnalysisTaskSED0Mass : public AliAnalysisTaskSE
   void SetRejectSDDClusters(Bool_t flag) { fIsRejectSDDClusters=flag; }
   void SetUseSelectionBit(Bool_t flag) { fUseSelectionBit=flag; }
   void SetWriteVariableTree(Bool_t flag) { fWriteVariableTree=flag; }
+  void SetDrawDetSignal(Bool_t flag) { fDrawDetSignal=flag; }
+  void SetPIDCheck(Bool_t flag) { fPIDCheck=flag; }
+
 
   Bool_t GetCutOnDistr() const {return fCutOnDistr;}
   Bool_t GetUsePid4Distr() const {return fUsePid4Distr;}
@@ -69,11 +72,15 @@ class AliAnalysisTaskSED0Mass : public AliAnalysisTaskSE
   Bool_t GetRejectSDDClusters() const { return fIsRejectSDDClusters; }
   Bool_t GetUseSelectionBit() const { return fUseSelectionBit; }
   Bool_t GetWriteVariableTree() const {return fWriteVariableTree;}
+  Bool_t GetDrawDetSignal() const {return fDrawDetSignal;}
+  Bool_t GetPIDCheck() const {return fPIDCheck;}
 
  private:
 
   AliAnalysisTaskSED0Mass(const AliAnalysisTaskSED0Mass &source);
   AliAnalysisTaskSED0Mass& operator=(const AliAnalysisTaskSED0Mass& source); 
+  void    DrawDetSignal(AliAODRecoDecayHF2Prong *part, TList *ListDetSignal);
+
   void     FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi *cuts, TList *listout);
   void     FillVarHists(AliAODEvent *aodev,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout);
   AliAODVertex* GetPrimaryVtxSkipped(AliAODEvent *aodev);
@@ -107,9 +114,11 @@ class AliAnalysisTaskSED0Mass : public AliAnalysisTaskSE
   Bool_t    fWriteVariableTree;       // flag to decide whether to write the candidate variables on a tree variables
   TTree    *fVariablesTree;           //! tree of the candidate variables after track selection on output slot 7
   Double_t *fCandidateVariables;      //!  variables to be written to the tree
+  Bool_t       fPIDCheck;                      // flag to decide whether to fill "PID = x" bins in fNentrie
+  Bool_t    fDrawDetSignal;            // flag to decide whether to draw the TPC dE/dx and TOF signal before/after PID
+  TList           *fDetSignal;         //!Detector signal histograms (on output slot 8)
 
-
-  ClassDef(AliAnalysisTaskSED0Mass,17); // AliAnalysisTaskSE for D0->Kpi
+  ClassDef(AliAnalysisTaskSED0Mass,18); // AliAnalysisTaskSE for D0->Kpi
 };
 
 #endif
index b04d450..340bdeb 100644 (file)
@@ -49,7 +49,11 @@ fDefaultPID(kFALSE),
 fUseKF(kFALSE),
 fPtLowPID(2.),
 fPtMaxSpecialCuts(9999.),
-fmaxPtrackForPID(9999)
+fmaxPtrackForPID(9999),
+fCombPID(kFALSE),
+fWeightsPositive(new Double_t[AliPID::kSPECIES]),
+fWeightsNegative(new Double_t[AliPID::kSPECIES]),
+fBayesianStrategy(0)
 {
   //
   // Default Constructor
@@ -104,7 +108,11 @@ AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
   fUseKF(source.fUseKF),
   fPtLowPID(source.fPtLowPID),
   fPtMaxSpecialCuts(source.fPtMaxSpecialCuts),
-  fmaxPtrackForPID(source.fmaxPtrackForPID)
+  fmaxPtrackForPID(source.fmaxPtrackForPID),
+  fCombPID(source.fCombPID),
+  fWeightsPositive(source.fWeightsPositive),
+  fWeightsNegative(source.fWeightsNegative),
+  fBayesianStrategy(source.fBayesianStrategy)
 {
   //
   // Copy constructor
@@ -127,10 +135,31 @@ AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &sour
   fPtLowPID=source.fPtLowPID;
   fPtMaxSpecialCuts=source.fPtMaxSpecialCuts;
   fmaxPtrackForPID=source.fmaxPtrackForPID;
+  fCombPID = source.fCombPID;
+  fWeightsPositive = source.fWeightsPositive;
+  fWeightsNegative = source.fWeightsNegative;
+  fBayesianStrategy = source.fBayesianStrategy;
   return *this;
 }
 
 
+
+//---------------------------------------------------------------------------
+AliRDHFCutsD0toKpi::~AliRDHFCutsD0toKpi()
+{
+//
+// Destructor
+//
+  if (fWeightsNegative) {
+    delete[] fWeightsNegative;
+    fWeightsNegative = 0;
+  }
+  if (fWeightsPositive) {
+    delete[] fWeightsNegative;
+    fWeightsNegative = 0;
+  }
+}
+
 //---------------------------------------------------------------------------
 void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {
   // 
@@ -371,11 +400,38 @@ Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEve
   if(selectionLevel==AliRDHFCuts::kAll || 
      selectionLevel==AliRDHFCuts::kCandidate ||
      selectionLevel==AliRDHFCuts::kPID) {
-    returnvaluePID = IsSelectedPID(d);
-    fIsSelectedPID=returnvaluePID;
-    if(!returnvaluePID) return 0;
+    if (!fCombPID) {
+      returnvaluePID = IsSelectedPID(d);
+      fIsSelectedPID=returnvaluePID;
+      if(!returnvaluePID) return 0;
+   } else {
+      switch (fBayesianStrategy) {
+      case kBayesSimple: {
+       returnvaluePID = IsSelectedSimpleBayesianPID(d);
+       break;
+      }
+      case kBayesWeightNoFilter: {
+       CalculateBayesianWeights(d);
+       returnvaluePID = 3;
+       break;
+      }
+      case(kBayesWeight): {
+       returnvaluePID = IsSelectedCombPID(d);
+       break;
+      }
+      case(kBayesMomentum): {
+       returnvaluePID = IsSelectedCombPID(d);
+       break;
+      }
+      default: {
+       returnvaluePID = IsSelectedCombPID(d);
+       break;
+      }
+      }
+    }
   }
 
+
   Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
 
   if(!returnvalueComb) return 0;
@@ -682,6 +738,7 @@ Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d)
 
   if(!fUsePID) return 3;
   if(fDefaultPID) return IsSelectedPIDdefault(d);
+  if (fCombPID) return IsSelectedCombPID(d); //to use Bayesian method if applicable
   fWhyRejection=0;
   Int_t isD0D0barPID[2]={1,2};
   Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value; 
@@ -1147,6 +1204,7 @@ void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)
 }
 
 
+//---------------------------------------------------------------------------
 void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
   //
   //STANDARD CUTS USED FOR 2010 pp analysis 
@@ -1272,6 +1330,7 @@ void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
 }
 
 
+//---------------------------------------------------------------------------
 void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() {
   //
   // STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV
@@ -1397,6 +1456,7 @@ void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() {
 }
 
 
+//---------------------------------------------------------------------------
 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
   //
   // Final CUTS USED FOR 2010 PbPb analysis of central events
@@ -1537,6 +1597,7 @@ void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
 
 }
 
+//---------------------------------------------------------------------------
 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {
   // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%
 
@@ -1673,6 +1734,7 @@ void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {
 }
 
 
+//---------------------------------------------------------------------------
 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {
   
   // Default 2010 PbPb cut object
@@ -1690,3 +1752,319 @@ void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {
   EnableSemiCentralTrigger();
 
 }
+
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedCombPID(AliAODRecoDecayHF* d)
+{
+  // ############################################################
+  //
+  // Apply Bayesian PID selection
+  // Implementation of Bayesian PID: Jeremy Wilkinson
+  //
+  // ############################################################
+
+  if (!fUsePID || !d) return 3;
+
+  switch (fBayesianStrategy) {
+
+  case kBayesWeightNoFilter: {
+    CalculateBayesianWeights(d);   //If weighted, no filtering: Calculate weights, return as unidentified
+    return 3;
+    break;
+  }
+  case kBayesSimple: {
+    return IsSelectedSimpleBayesianPID(d);   // If simple, go to simple method
+    break;
+  }
+  case(kBayesWeight): {
+    break;
+  }
+  case(kBayesMomentum): {
+    break;  //If weighted or momentum method, stay in this function
+  }
+
+  }
+
+  //  Int_t isD0D0barPID[2]={1,2};
+
+
+
+
+  Int_t isD0 = 0;
+  Int_t isD0bar = 0;
+  Int_t returnvalue = 0;
+
+  Int_t isPosKaon = 0, isNegKaon = 0;
+  //   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  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: unknown (but compatible)
+
+
+  //   combinedPID[0][0] = 0; // First daughter, kaon
+  //   combinedPID[0][1] = 0; // First daughter, pion
+  //   combinedPID[1][0] = 0; // Second daughter, kaon
+  //   combinedPID[1][1] = 0; // Second daughter, pion
+
+  Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
+  AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
+
+  if ((aodtrack[0]->Charge()*aodtrack[1]->Charge()) != -1) {
+    return 0;  //reject if charges not opposite
+  }
+  for (Int_t daught = 0; daught < 2; daught++) {
+    //Loop over prongs
+
+    if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
+      checkPIDInfo[daught] = kFALSE;
+      continue;
+    }
+    CalculateBayesianWeights(d);
+    //      Double_t prob0[AliPID::kSPECIES];
+    //      fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
+    ///Three possible Bayesian methods: Pick one (only first implemented so far. TODO: implement other two methods)
+    ///A: Standard max. probability method   
+    if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
+      isPosKaon = 1;  //flag [daught] as a kaon
+    }
+
+    if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
+      isNegKaon = 1;  //flag [daught] as a kaon
+    }
+    ///B: Method based on probability greater than prior
+    /* Double_t kaonpriorprob;
+       TH1D *kaonprior = fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon);
+       kaonpriorprob =   kaonprior->GetBinContent(kaonprior->FindBin(aodtrack[daught]->P()));
+
+       if (prob0[AliPID::kKaon] > kaonpriorprob){combinedPID[daught][0] = 1;}
+
+    */
+    ///C: Probability threshold (TODO: Add necessary threshold variables; currently unused)
+    /*
+      if(prob0[AliPID::kKaon] >= fMinProbKaon) {combinedPID[daught][0] = 1;}
+
+    */
+
+  }
+   
+   
+  //Loop over daughters ends here
+
+  if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
+    return 0;      //Reject if both daughters lack both TPC and TOF info
+  }
+  //Momentum-based selection
+
+
+  if (isNegKaon && isPosKaon) { // If both are kaons, reject
+    isD0 = 0;
+    isD0bar = 0;
+  } else if (isNegKaon) {       //If negative kaon present, D0
+    isD0 = 1;
+  } else if (isPosKaon) {       //If positive kaon present, D0bar
+    isD0bar = 1;
+  } else {                      //If neither ID'd as kaon, subject to extra tests
+    isD0 = 1;
+    isD0bar = 1;
+    if (aodtrack[0]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[0], "TOF"))) { //If it's low momentum and not ID'd as a kaon, we assume it must be a pion
+      if (aodtrack[0]->Charge() == -1) {
+       isD0 = 0;  //Check charge and reject based on pion hypothesis
+      }
+      if (aodtrack[0]->Charge() == 1) {
+       isD0bar = 0;
+      }
+    }
+    if (aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
+      if (aodtrack[1]->Charge() == -1) {
+       isD0 = 0;
+      }
+      if (aodtrack[1]->Charge() == 1) {
+       isD0bar = 0;
+      }
+    }
+  }
+
+   
+
+  if (isD0 && isD0bar) {
+    returnvalue = 3;
+  }
+  if (isD0 && !isD0bar) {
+    returnvalue = 1;
+  }
+  if (!isD0 && isD0bar) {
+    returnvalue = 2;
+  }
+  if (!isD0 && !isD0bar) {
+    returnvalue = 0;
+  }
+
+  return returnvalue;
+
+
+
+  /*
+
+  //OLD CODE
+  if (combinedPID[0][0] == 1 && combinedPID[1][1] == 1) {          //First track is kaon, second is pion
+  if (aodtrack[0]->Charge() == -1) isD0 = 1;            //Kaon negative: particle is D0
+  else isD0bar = 1;                        //Kaon positive: particle is anti
+  } else if (combinedPID[1][0] == 1 && combinedPID[0][1] == 1) {         //First track is pion, second is kaon
+  if (aodtrack[1]->Charge() == -1) isD0 = 1;              //Kaon negative: particle is D0
+  else isD0bar = 1;                        //Kaon positive: particle is anti
+  }
+
+  else if (combinedPID[0][0] == 1 && combinedPID[1][0] == 1) {
+  isD0 = 1;   //If both are kaons, leave to top. cuts
+  isD0bar = 1;
+  }
+
+
+  else {
+  isD0 = 1;
+  isD0bar = 1;
+  if (combinedPID[0][0] == 0 && aodtrack[0]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[0], "TOF"))) { //If it's low momentum and not ID'd as a kaon, we assume it must be a pion
+  if (aodtrack[0]->Charge() == -1) {
+  isD0 = 0;  //Check charge and reject depending on pion hypothesis
+  }
+  if (aodtrack[0]->Charge() == 1) {
+  isD0bar = 0;
+  }
+  }
+  if (combinedPID[1][0] == 0 && aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
+  if (aodtrack[1]->Charge() == -1) {
+  isD0 = 0;
+  }
+  if (aodtrack[1]->Charge() == 1) {
+  isD0bar = 0;
+  }
+  }
+  }*/
+  /*    if(combinedPID[0][1] == 1 && combinedPID[1][1] == 1) {         //If both are ID'd as pions:
+       if ((aodtrack[0]->P()) <1.0 && (aodtrack[1]->P() < 1.0)) {     //Check both momenta. If both low, we require one kaon.
+       isD0 == 0; isD0bar == 0;}              //Both pions, both low momenta: reject
+       else {isD0 == 1; isD0bar == 1;}              //Both pions, both high momenta: Allow as both (filter with topological cuts)
+       }
+  */
+
+
+}
+
+
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedSimpleBayesianPID(AliAODRecoDecayHF* d)
+
+{
+  //Apply Bayesian PID selection; "simple" method.
+  //Looks for a kaon via max. probability. If neither daughter is a kaon, candidate is rejected.
+
+  Int_t isPosKaon = 0, isNegKaon = 0;
+  Int_t returnvalue;          //0 = rejected, 1 = D0, 2 = D0bar. This version will not output 3 (both).
+  Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
+
+  AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
+  if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
+    return 0;  //Reject if charges do not oppose
+  }
+  for (Int_t daught = 0; daught < 2; daught++) {
+    //Loop over prongs to check PID information
+
+    if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
+      checkPIDInfo[daught] = kFALSE;
+    }
+
+  }
+  //Loop over daughters ends here.
+
+  if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
+    return 0;
+  }
+
+  CalculateBayesianWeights(d);
+
+   
+  // identify kaon
+  ///Three possible Bayesian methods: Pick one (only first implemented so far. TODO: implement other two methods)
+  ///A: Standard max. probability method
+
+   
+  if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
+    isPosKaon = 1;  //flag [daught] as a kaon
+  }
+
+  if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
+    isNegKaon = 1;  //flag [daught] as a kaon
+  }
+
+
+  ///B: Method based on probability greater than prior
+  /* Double_t kaonpriorprob;
+     TH1D *kaonprior = fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon);
+     kaonpriorprob =   kaonprior->GetBinContent(kaonprior->FindBin(aodtrack[daught]->P()));
+
+     if (prob0[AliPID::kKaon] > kaonpriorprob){combinedPID[daught][0] = 1;}
+
+  */
+  ///C: Probability threshold (TODO: Add necessary threshold variables; currently unused)
+  /*
+    if(prob0[AliPID::kKaon] >= fMinProbKaon) {combinedPID[daught][0] = 1;}
+
+  */
+
+
+
+
+  if (isPosKaon && isNegKaon)   {  //If both are ID'd as kaons, reject
+    returnvalue = 0;
+  } else if (isNegKaon)   {     //If negative kaon, D0
+    returnvalue = 1;
+  } else if (isPosKaon)   {     //If positive kaon, D0-bar
+    returnvalue = 2;
+  } else {
+    returnvalue = 0;  //If neither kaon, reject
+  }
+
+  return returnvalue;
+}
+
+
+
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::CalculateBayesianWeights(AliAODRecoDecayHF* d)
+
+{
+  //Function to compute weights for Bayesian method
+  Double_t *prob0;
+  prob0 = new Double_t[AliPID::kSPECIES];
+
+
+  AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
+  if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
+    return;  //Reject if charges do not oppose
+  }
+  for (Int_t daught = 0; daught < 2; daught++) {
+    //Loop over prongs
+
+    if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
+      continue;
+    }
+
+
+    // identify kaon
+    fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
+
+    if (aodtrack[daught]->Charge() == +1) {
+      SetWeightsPositive(prob0);
+    }
+    if (aodtrack[daught]->Charge() == -1) {
+      SetWeightsNegative(prob0);
+    }
+
+  }
+}
+
+
+
index 57d4961..149e214 100644 (file)
@@ -23,7 +23,7 @@ class AliRDHFCutsD0toKpi : public AliRDHFCuts
 
   AliRDHFCutsD0toKpi(const char* name="CutsD0toKpi");
   
-  virtual ~AliRDHFCutsD0toKpi(){}
+  virtual ~AliRDHFCutsD0toKpi();
 
   AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi& source);
   AliRDHFCutsD0toKpi& operator=(const AliRDHFCutsD0toKpi& source); 
@@ -39,6 +39,10 @@ class AliRDHFCutsD0toKpi : public AliRDHFCuts
                          {return IsSelected(obj,selectionLevel,0);}
   virtual Int_t IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod);
 
+  virtual Int_t IsSelectedCombPID(AliAODRecoDecayHF* d); 
+  virtual Int_t IsSelectedSimpleBayesianPID(AliAODRecoDecayHF* d);
+  virtual void CalculateBayesianWeights(AliAODRecoDecayHF* d);
+
   Float_t GetMassCut(Int_t iPtBin=0) const { return (GetCuts() ? fCutsRD[GetGlobalIndex(0,iPtBin)] : 1.e6);}
   Float_t GetDCACut(Int_t iPtBin=0) const { return (GetCuts() ? fCutsRD[GetGlobalIndex(1,iPtBin)] : 1.e6);}
   Int_t CombineSelectionLevels(Int_t selectionvalTrack,Int_t selectionvalCand,Int_t selectionvalPID)const;
@@ -63,6 +67,31 @@ class AliRDHFCutsD0toKpi : public AliRDHFCuts
   Double_t GetPtForPIDtight()const {return fPtLowPID;}
   void SetUseKF(Bool_t useKF);
   Bool_t GetIsUsedKF() const {return fUseKF;}
+  void SetWeightsPositive(Double_t* weights){
+     for (Int_t i = 0; i<AliPID::kSPECIES; i++) {
+         fWeightsPositive[i] = weights[i];
+     }
+}
+  Double_t *GetWeightsPositive() const {return fWeightsPositive;}
+  void SetWeightsNegative(Double_t* weights){
+     for (Int_t i = 0; i<AliPID::kSPECIES; i++) {
+         fWeightsNegative[i] = weights[i];
+     }
+     }
+  Double_t *GetWeightsNegative() const {return fWeightsNegative;}
+  void SetBayesianStrategy(Int_t strat) {fBayesianStrategy=strat;}
+  Int_t GetBayesianStrategy() const {return fBayesianStrategy;}
+  
+  enum EBayesianStrategy {
+     kBayesMomentum,
+     kBayesWeight,
+     kBayesWeightNoFilter,
+     kBayesSimple
+  };
+  
+  void SetCombPID(Bool_t CombPID){fCombPID=CombPID;}
+  Bool_t GetCombPID() const {return fCombPID;}
+  
 
  
  protected:
@@ -71,13 +100,21 @@ class AliRDHFCutsD0toKpi : public AliRDHFCuts
   Bool_t fUseSpecialCuts;  // flag to switch on/off special cuts
   Bool_t fLowPt;           // flag to switch on/off different pid for low pt D0
   Bool_t fDefaultPID;      // flag to switch on/off the default pid
+  
   Bool_t fUseKF;           // flag to switch on/off D0 selection via KF 
   Double_t fPtLowPID;      // transverse momentum below which the strong PID is applied
   Double_t fPtMaxSpecialCuts; // transverse momentum below which the special cuts are applied
+  
                               //  if set to zero, used for all pt
   Double_t  fmaxPtrackForPID; // max momentum for applying PID
 
-  ClassDef(AliRDHFCutsD0toKpi,8);  // class for cuts on AOD reconstructed D0->Kpi
+  Bool_t fCombPID;             //switch for Bayesian
+  
+  Double_t* fWeightsPositive;  //Bayesian weights for positive track
+  Double_t* fWeightsNegative;  //Bayesian weights for negative track
+  Int_t fBayesianStrategy;
+
+  ClassDef(AliRDHFCutsD0toKpi,9);  // class for cuts on AOD reconstructed D0->Kpi
 };
 
 #endif
index a4d969f..1d9d45d 100644 (file)
@@ -4,7 +4,8 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read
                                       Float_t minC=0, Float_t maxC=0,
                                       TString finDirname="Loose",
                                       TString finname="",TString finObjname="D0toKpiCuts", Bool_t flagAOD049=kFALSE,
-                                      Bool_t FillMassPt=false, Bool_t FillImpPar=false)
+                                      Bool_t FillMassPt=false, Bool_t FillImpPar=false,
+                                          Bool_t DrawDetSignal=false, Bool_t PIDCheck=false)
 {
   //
   // AddTask for the AliAnalysisTaskSE for D0 candidates
@@ -22,7 +23,7 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read
     return NULL;
   }   
 
-  TString filename="",out1name="",out2name="",out3name="",out4name="",out5name="",out6name="",out7name="",inname="";
+  TString filename="",out1name="",out2name="",out3name="",out4name="",out5name="",out6name="",out7name="",out8name="", inname="";
   filename = AliAnalysisManager::GetCommonFileName();
   filename += ":PWG3_D2H_";
   if(flag==0){
@@ -65,6 +66,16 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read
 
     out7name ="coutputVarTree";
 
+   //detector signal hists
+
+    out8name="detectorSignals";
+    if(cutOnDistr) out8name+="C"; 
+    if(flagD0D0bar==1)out8name+="D0";
+    if(flagD0D0bar==2)out8name+="D0bar";
+
+    
+
+
     inname="cinputmassD0_0";
     if(cutOnDistr) inname+="C"; 
     if(flagD0D0bar==1)inname+="D0";
@@ -109,6 +120,7 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read
   out5name += finDirname.Data(); 
   out6name += finDirname.Data();
   out7name += finDirname.Data();
+  out8name += finDirname.Data();
   inname += finDirname.Data();
 
    //setting my cut values
@@ -172,6 +184,7 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read
   out5name+=centr;
   out6name+=centr;
   out7name+=centr;
+  out8name+=centr;
   inname+=centr;
 
   // Aanalysis task    
@@ -190,7 +203,8 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read
 
   massD0Task->SetFillPtHistos(FillMassPt);
   massD0Task->SetFillImpactParameterHistos(FillImpPar);
-  
+  massD0Task->SetDrawDetSignal(DrawDetSignal);
+  massD0Task->SetPIDCheck(PIDCheck);
   //  massD0Task->SetRejectSDDClusters(kTRUE);
 
   //   massD0Task->SetWriteVariableTree(kTRUE);
@@ -209,6 +223,9 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read
   AliAnalysisDataContainer *coutputmassD05 = mgr->CreateContainer(out5name,AliNormalizationCounter::Class(),AliAnalysisManager::kOutputContainer, filename.Data()); //counter
   AliAnalysisDataContainer *coutputmassD06 = mgr->CreateContainer(out6name,TList::Class(),AliAnalysisManager::kOutputContainer, filename.Data()); //mass vs pt vs impt par
   AliAnalysisDataContainer *coutputmassD07 = mgr->CreateContainer(out7name,TTree::Class(),AliAnalysisManager::kOutputContainer, filename.Data()); //mass vs pt vs impt par
+  AliAnalysisDataContainer *coutputmassD08 = mgr->CreateContainer(out8name,TList::Class(),AliAnalysisManager::kOutputContainer, filename.Data()); //dedx
+  
+
   
   mgr->ConnectInput(massD0Task,0,mgr->GetCommonInputContainer());
 
@@ -219,6 +236,7 @@ AliAnalysisTaskSED0Mass *AddTaskD0Mass(Int_t flag=0/*0 = D0,1 = LS*/,Bool_t read
   mgr->ConnectOutput(massD0Task,5,coutputmassD05);
   mgr->ConnectOutput(massD0Task,6,coutputmassD06);
   mgr->ConnectOutput(massD0Task,7,coutputmassD07);
+  mgr->ConnectOutput(massD0Task,8,coutputmassD08);
 
 
   return massD0Task;