]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/correlationHF/AliDxHFEParticleSelectionMCEl.cxx
Correct backward incompatible change (R.Romita)
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliDxHFEParticleSelectionMCEl.cxx
index 0874d841c31e0318ceda579ddc18a0bada331dd4..57405192f483696aeadf09e7a887cfa5582ee126 100644 (file)
@@ -30,6 +30,7 @@
 #include "TH1F.h"
 #include "TAxis.h"
 #include "AliAODTrack.h"
+#include "AliReducedParticle.h"
 #include <iostream>
 #include <cerrno>
 #include <memory>
@@ -44,10 +45,15 @@ ClassImp(AliDxHFEParticleSelectionMCEl)
 AliDxHFEParticleSelectionMCEl::AliDxHFEParticleSelectionMCEl(const char* opt)
   : AliDxHFEParticleSelectionEl(opt)
   , fMCTools()
-  , fElectronProperties(NULL)
-  , fWhichCut(NULL)
+  , fPDGnotMCElectron(NULL)
+  , fPDGNotHFMother(NULL)
   , fOriginMother(0)
   , fResultMC(0)
+  , fUseKine(kFALSE)
+  , fMotherPDGs(0)
+  , fUseMCReco(kFALSE)
+  , fSelectionStep(AliDxHFEParticleSelectionEl::kNoCuts)
+  , fStoreCutStepInfo(kFALSE)
 {
   // constructor
   // 
@@ -57,9 +63,15 @@ AliDxHFEParticleSelectionMCEl::AliDxHFEParticleSelectionMCEl(const char* opt)
 
   fMCTools.~AliDxHFEToolsMC();
 
+  ParseArguments(opt);
+
+  // This is also checked in AddTasks, but just for safety! 
+  if(fUseMCReco && fUseKine) AliFatal("CAN'T SET BOTH usekine AND elmcreco AT THE SAME TIME");
+
   // TODO: argument scan, build tool options accordingly
   // e.g. set mc mode first/last, skip control histograms
   TString toolopt("pdg=11 mc-last");
+  if(fUseKine) toolopt+=" usekine";
   new (&fMCTools) AliDxHFEToolsMC(toolopt);
 }
 
@@ -78,9 +90,55 @@ const char* AliDxHFEParticleSelectionMCEl::fgkPDGMotherBinLabels[]={
   "others"
 };
 
+                             
+const char*  AliDxHFEParticleSelectionMCEl::fgkPDGBinLabels[]={
+  "positron",
+  "electron",
+  "#mu+",
+  "#mu-",
+  "#pi+",
+  "#pi-",
+  "K+",
+  "K-",
+  "proton",
+  "antiproton",
+  "others"
+};
+
 AliDxHFEParticleSelectionMCEl::~AliDxHFEParticleSelectionMCEl()
 {
   // destructor
+
+  if(fPDGnotMCElectron){
+    delete fPDGnotMCElectron;
+    fPDGnotMCElectron=NULL;
+  }
+
+}
+
+int AliDxHFEParticleSelectionMCEl::Init()
+{
+  //
+  // init function
+  // 
+
+  // Particles considered HFE background. can be expanded
+  fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGpi0); 
+  fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGeta);
+  fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGgamma);
+  fMotherPDGs.push_back(AliDxHFEToolsMC::kPDGJpsi);
+
+  int iResult=0;
+  iResult=AliDxHFEParticleSelectionEl::Init();
+  if (iResult<0) return iResult;
+
+  // Histo containing PDG of track which was not MC truth electron
+  // TODO: Add to list in baseclass
+  fPDGnotMCElectron=CreateControlHistogram("fPDGnotMCElectron","PDG of track not MC truth electron",AliDxHFEToolsMC::kNofPDGLabels,fgkPDGBinLabels);  
+  fPDGNotHFMother=CreateControlHistogram("fPDGNotHFMother","PDG of mother not HF",5000);  
+  AddControlObject(fPDGnotMCElectron);
+  AddControlObject(fPDGNotHFMother);
+  return 0;
 }
 
 THnSparse* AliDxHFEParticleSelectionMCEl::DefineTHnSparse()
@@ -88,47 +146,77 @@ THnSparse* AliDxHFEParticleSelectionMCEl::DefineTHnSparse()
   //
   // Defines the THnSparse. 
 
-  const int thnSize = 6;
-  InitTHnSparseArray(thnSize);
   const double Pi=TMath::Pi();
   TString name;
+  THnSparse* thn=NULL;
   name.Form("%s info", GetName());
 
-  //                          0    1      2     3       4         5
-  //                          Pt   Phi   Eta  'stat e'  mother "pdg e"
-  int    thnBins[thnSize] = { 1000,  200, 500,    2,       14,      10 };
-  double thnMin [thnSize] = {    0,    0, -1., -0.5,     -1.5,    -0.5 };
-  double thnMax [thnSize] = {  100, 2*Pi,  1.,  1.5,     12.5,     9.5 };
-  const char* thnNames[thnSize]={
-    "Pt",
-    "Phi",
-    "Eta", 
-    "MC statistics electron",
-    "Mother", 
-    "PDG of selected electron"
-  };
-  return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
+  if(fStoreCutStepInfo){
+    const int thnSizeExt =5;
+
+    InitTHnSparseArray(thnSizeExt);
+
+    // TODO: Redo binning of distributions more?
+    //                        0    1      2     3     
+    //                                Pt   Phi   Eta   mother 
+    int    thnBinsExt[thnSizeExt] = { 100,  100, 100,    15, kNCutLabels };
+    double thnMinExt [thnSizeExt] = {   0,    0, -1.,  -1.5, kRecKineITSTPC };
+    double thnMaxExt [thnSizeExt] = {  10, 2*Pi,  1.,  13.5, kSelected };
+    const char* thnNamesExt[thnSizeExt]={
+      "Pt",
+      "Phi",
+      "Eta", 
+      "Mother", //bin==-1: Not MC truth electron
+      "Last survived cut step"
+    };
+    thn=(THnSparse*)CreateControlTHnSparse(name,thnSizeExt,thnBinsExt,thnMinExt,thnMaxExt,thnNamesExt);
+  }
+  else{
+
+    const int thnSize =4;
+    InitTHnSparseArray(thnSize);
+
+    // TODO: Redo binning of distributions more?
+    //                        0    1      2     3     
+    //                                Pt   Phi   Eta   mother 
+    int    thnBins[thnSize] = { 100,  100, 100,    15  };
+    double thnMin [thnSize] = {   0,    0, -1.,  -1.5  };
+    double thnMax [thnSize] = {  10, 2*Pi,  1.,  13.5  };
+    const char* thnNames[thnSize]={
+      "Pt",
+      "Phi",
+      "Eta", 
+      "Mother", //bin==-1: Not MC truth electron
+    };
+    thn=(THnSparse*)CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
+
+  }
+  return thn;
 }
 
 int AliDxHFEParticleSelectionMCEl::FillParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
 {
   // fill the data array from the particle data
   if (!data) return -EINVAL;
-  AliAODTrack *track=(AliAODTrack*)p;
-  if (!track) return -ENODATA;
+  if (!p) return -ENODATA;
+  // handle different types of tracks, can be extended
+  AliReducedParticle *trRP=dynamic_cast<AliReducedParticle*>(p);
   int i=0;
   if (dimension!=GetDimTHnSparse()) {
     // TODO: think about filling only the available data and throwing a warning
     return -ENOSPC;
   }
-  data[i++]=track->Pt();
-  data[i++]=track->Phi();
-  data[i++]=track->Eta();
-  data[i++]=fResultMC;     // stat electron (MC electron or not)
-  data[i++]=fOriginMother; // at the moment not included background. Should expand
-  data[i++]=1;             // PDG e - not sure if needed, maybe only needed as separate histo? 
-  
+  memset(data, 0, dimension*sizeof(data[0]));
+  data[i++]=p->Pt();
+  data[i++]=p->Phi();
+  data[i++]=p->Eta();
+  if (trRP) data[i]=trRP->GetOriginMother();
+  else data[i]=fOriginMother;
+  i++; // take out of conditionals to be save
+  if (i<dimension) {
+    if(fStoreCutStepInfo) data[i]=GetLastSurvivedCutsStep();
+    i++; // take out of conditionals to be save
+  }
   return i;
 }
 
@@ -145,6 +233,9 @@ int AliDxHFEParticleSelectionMCEl::IsSelected(AliVParticle* p, const AliVEvent*
   if (!p || !pEvent){
     return -EINVAL;
   }
+  fOriginMother=-1;
+
+  if(!fUseKine && !fUseMCReco){
   // step 1:
   // optional MC selection before the particle selection
   if (fMCTools.MCFirst() && (iResult=CheckMC(p, pEvent))==0) {
@@ -157,15 +248,32 @@ int AliDxHFEParticleSelectionMCEl::IsSelected(AliVParticle* p, const AliVEvent*
   iResult=AliDxHFEParticleSelectionEl::IsSelected(p, pEvent);
   if (fMCTools.MCFirst() || iResult==0) return iResult;
 
+
+  }
+
   // step 2, only executed if MC check is last
   // optional MC selection after the particle selection
   // result stored to be filled into THnSparse
   // TODO: strictly speaken the particles should be rejected
   // if not mc selected, however skip this for the moment, because of
   // the logic outside
+  // This line will run always when running directly over the stack or over MC reconstructed tracks
+  // For MC reconstructed tracks, should consider adding more constrictions on the tracks,
+  // maybe do the track selection first (which means could call AliDxHFEParticleSelectionEl::IsSelected()
+  // and don't do PID
+  
+  if(fUseMCReco) {// && ( fSelectionStep > AliDxHFEParticleSelectionEl::kNoCuts )){
+    // always check the base class method in mode fUseMCReco
+    iResult=AliDxHFEParticleSelectionEl::IsSelected(p, pEvent);
+    if(iResult == 0) return iResult;
+  }
   fResultMC=CheckMC(p, pEvent);
+  
+  if(fUseKine || fUseMCReco)
+    return fResultMC;
 
   return iResult;
+   
 }
 
 int AliDxHFEParticleSelectionMCEl::CheckMC(AliVParticle* p, const AliVEvent* pEvent)
@@ -174,39 +282,44 @@ int AliDxHFEParticleSelectionMCEl::CheckMC(AliVParticle* p, const AliVEvent* pEv
   if (!p || !pEvent){
     return -EINVAL;
   }
-  fOriginMother=-1;
   int iResult=0;
 
   if (!fMCTools.IsInitialized() && (iResult=fMCTools.InitMCParticles(pEvent))<0) {
     // TODO: message? but has to be filtered in order to avoid message flood
     return 0; // no meaningful filtering on mc possible
   }
-  if (fMCTools.RejectByPDG(p,false)) {
+
+  if(fUseKine){
+    // If run on kinematical level, check if particle is electron (only ones who are relevant) 
+    // and if pass IsPhysicalPrimary()
+    Int_t test = fMCTools.CheckMCParticle(p);
+    if(test==1) return 0;
+    // Add additional contraints on the electrons? 
+    // remove delta e? also remove E<300MeV?
+    // Should also mark dalitz decay and gamma conversion..
+    AliAODMCParticle *mcp=dynamic_cast<AliAODMCParticle*>(p);
+    if(!mcp || !mcp->IsPhysicalPrimary()) return 0; 
+
+  }
+
+  int pdgParticle=-1;
+  if (!fUseKine && fMCTools.RejectByPDG(p,false, &pdgParticle)) {
     // rejected by pdg
-    // would also like to use this info? or not meaningful?
-    // NEED TO CHANGE THIS!
+    // TODO: Move this to fMCTools???? Can this be part of the statistics in the MC class?
+    fPDGnotMCElectron->Fill(fMCTools.MapPDGLabel(pdgParticle));
     return 0;
   }
 
   int pdgMother=0;
+  // Find PDG of first mother
   pdgMother=fMCTools.FindMotherPDG(p,AliDxHFEToolsMC::kGetFirstMother);
 
-  // Particles considered HFE background. can be expanded
-  // Should be created only once 
-  // TODO: that needs to be configured once to avoid performance penalty
-  vector<int> motherPDGs;
-  motherPDGs.push_back(AliDxHFEToolsMC::kPDGpi0); 
-  motherPDGs.push_back(AliDxHFEToolsMC::kPDGeta);
-  motherPDGs.push_back(AliDxHFEToolsMC::kPDGgamma);
-  motherPDGs.push_back(AliDxHFEToolsMC::kPDGJpsi);
-
-  if(fMCTools.RejectByPDG(pdgMother,motherPDGs)){
-    pdgMother=fMCTools.FindMotherPDG(p,AliDxHFEToolsMC::kGetOriginMother);
-    fOriginMother=fMCTools.GetOriginMother();
-  }
-  else{
-    //Could this be done in a more elegant way?
+  // Check if first mother is counted as background
+  Bool_t isNotBackground=fMCTools.RejectByPDG(pdgMother,fMotherPDGs);
+
+  if(!isNotBackground){
+    // Set fOriginMother if mother is counted as background
+    // TODO: Could this be done in a more elegant way?
     switch(pdgMother){
     case(AliDxHFEToolsMC::kPDGpi0): fOriginMother=AliDxHFEToolsMC::kNrOrginMother; break;
     case(AliDxHFEToolsMC::kPDGeta): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+1; break;
@@ -214,14 +327,24 @@ int AliDxHFEParticleSelectionMCEl::CheckMC(AliVParticle* p, const AliVEvent* pEv
     case(AliDxHFEToolsMC::kPDGJpsi): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+3;break;
     }
   }
+  else{
+    // If loop over Stack, also checks if First mother is a HF meson
+    Bool_t isHFmeson =fMCTools.TestMotherHFMeson(TMath::Abs(pdgMother));
 
-  /*if (fMCTools.RejectByMotherPDG(p)) {
-    // rejected by pdg of original mother
-    // H: want pdg of origin process to be stored in THnSparse
-    // Not sure if this is needed... Need to check, using AliDxHFEToolsMC, who
-    // first mother are, and also what origin is. Use this info here. 
-    return 0;
-    }*/
+    if(isHFmeson){
+      // If first mother is a HF meson, loops back to find 
+      // original quark + if there was a gluon. Result is 
+      // stored in fOriginMother
+      pdgMother=fMCTools.FindMotherPDG(p,AliDxHFEToolsMC::kGetOriginMother);
+      fOriginMother=fMCTools.GetOriginMother();
+    }
+    else{
+      //NotHFmother
+      fPDGNotHFMother->Fill(pdgMother);
+      fOriginMother=AliDxHFEToolsMC::kNrOrginMother+4;
+    }
+
+  }
 
   return 1;
 }
@@ -231,3 +354,62 @@ void AliDxHFEParticleSelectionMCEl::Clear(const char* option)
   /// clear internal memory
   fMCTools.Clear(option);
 }
+
+AliVParticle *AliDxHFEParticleSelectionMCEl::CreateParticle(AliVParticle* track)
+{
+
+  AliReducedParticle *part = new AliReducedParticle(track->Eta(), track->Phi(), track->Pt(),track->Charge(),fOriginMother);
+
+  return part;
+
+}
+
+int AliDxHFEParticleSelectionMCEl::ParseArguments(const char* arguments)
+{
+  // parse arguments and set internal flags
+  TString strArguments(arguments);
+  auto_ptr<TObjArray> tokens(strArguments.Tokenize(" "));
+  if (!tokens.get()) return 0;
+
+  AliInfo(strArguments);
+  TIter next(tokens.get());
+  TObject* token;
+  while ((token=next())) {
+    TString argument=token->GetName();
+    if (argument.BeginsWith("usekine") ){
+      fUseKine=kTRUE;
+      continue;
+    }
+    if (argument.BeginsWith("elmcreco")){
+      fUseMCReco=kTRUE;
+      if(argument.BeginsWith("elmcreco=")){
+       argument.ReplaceAll("elmcreco=", "");
+       if(argument.CompareTo("alltracks")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kNoCuts;
+       if(argument.CompareTo("afterreckineitstpc")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kRecKineITSTPC;
+       if(argument.CompareTo("afterrecprim")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kRecPrim;
+       if(argument.CompareTo("afterhfeits")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kHFEcutsITS;
+       if(argument.CompareTo("afterhfetof")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kHFEcutsTOF;
+       if(argument.CompareTo("afterhfetpc")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kHFEcutsTPC;
+       if(argument.CompareTo("aftertrackcuts")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kHFEcutsTPC;
+       if(argument.CompareTo("aftertofpid")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kPIDTOF;
+       if(argument.CompareTo("aftertpcpid")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kPIDTPC;
+       if(argument.CompareTo("afterfullpid")==0) fSelectionStep=AliDxHFEParticleSelectionEl::kPIDTOFTPC;
+
+       AliDxHFEParticleSelectionEl::SetFinalCutStep(fSelectionStep);
+      }
+       continue;
+    }
+    if(argument.BeginsWith("storelastcutstep")){
+      AliInfo("Stores the last cut step");
+      fUseMCReco=kTRUE;
+      fStoreCutStepInfo=kTRUE;
+      AliDxHFEParticleSelectionEl::SetStoreLastCutStep(kTRUE);
+      continue;
+    }
+    // forwarding of single argument works, unless key-option pairs separated
+    // by blanks are introduced
+    AliDxHFEParticleSelection::ParseArguments(argument);
+  }
+  
+  return 0;
+}