]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
MC implementations for D0-HFE correlation (Matthias)
authorprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 11 Sep 2012 13:23:15 +0000 (13:23 +0000)
committerprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 11 Sep 2012 13:23:15 +0000 (13:23 +0000)
17 files changed:
PWGHF/CMakelibPWGHFcorrelationHF.pkg
PWGHF/correlationHF/AliAnalysisTaskDxHFECorrelation.cxx
PWGHF/correlationHF/AliAnalysisTaskDxHFEParticleSelection.cxx
PWGHF/correlationHF/AliAnalysisTaskDxHFEParticleSelection.h
PWGHF/correlationHF/AliDxHFECorrelation.cxx
PWGHF/correlationHF/AliDxHFEParticleSelection.cxx
PWGHF/correlationHF/AliDxHFEParticleSelection.h
PWGHF/correlationHF/AliDxHFEParticleSelectionD0.cxx
PWGHF/correlationHF/AliDxHFEParticleSelectionD0.h
PWGHF/correlationHF/AliDxHFEParticleSelectionEl.cxx
PWGHF/correlationHF/AliDxHFEParticleSelectionEl.h
PWGHF/correlationHF/AliDxHFEParticleSelectionMCD0.cxx [new file with mode: 0644]
PWGHF/correlationHF/AliDxHFEParticleSelectionMCD0.h [new file with mode: 0644]
PWGHF/correlationHF/AliDxHFEParticleSelectionMCEl.cxx [new file with mode: 0644]
PWGHF/correlationHF/AliDxHFEParticleSelectionMCEl.h [new file with mode: 0644]
PWGHF/correlationHF/AliDxHFEToolsMC.cxx [new file with mode: 0644]
PWGHF/correlationHF/AliDxHFEToolsMC.h [new file with mode: 0644]

index 1751c7249d6ff5691e7901d8699c5c8be63445e6..ba0893b322284a8947ab095ed00c12c0b493ebfe 100644 (file)
@@ -30,6 +30,9 @@ set ( CLASS_HDRS
     AliDxHFEParticleSelection.h
     AliDxHFEParticleSelectionD0.h
     AliDxHFEParticleSelectionEl.h
+    AliDxHFEParticleSelectionMCD0.h
+    AliDxHFEParticleSelectionMCEl.h
+    AliDxHFEToolsMC.h
     AliDxHFECorrelation.h
     AliAnalysisTaskDxHFEParticleSelection.h
     AliAnalysisTaskDxHFECorrelation.h
index 260154f70461fa7b529bdd27a780c7e9b22ecffc..7203f1ce256b73cfc130741282cbbd9932646736 100644 (file)
@@ -26,7 +26,9 @@
 #include "AliAnalysisTaskDxHFECorrelation.h"
 #include "AliDxHFECorrelation.h"
 #include "AliDxHFEParticleSelectionD0.h"
+#include "AliDxHFEParticleSelectionMCD0.h"
 #include "AliDxHFEParticleSelectionEl.h"
+#include "AliDxHFEParticleSelectionMCEl.h"
 #include "AliAnalysisManager.h"
 #include "AliLog.h"
 #include "AliESDInputHandler.h"
@@ -143,15 +145,17 @@ void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
   default: selectionD0Options+="FillD0D0bar ";
   }
 
-  fD0s=new AliDxHFEParticleSelectionD0(selectionD0Options);
+  if(fUseMC) fD0s=new AliDxHFEParticleSelectionMCD0(selectionD0Options);
+  else fD0s=new AliDxHFEParticleSelectionD0(selectionD0Options);
   fD0s->SetCuts(fCutsD0);
-  fD0s->InitControlObjects();
+  fD0s->Init();
 
   //Electrons
-  fElectrons=new AliDxHFEParticleSelectionEl;
+  if(fUseMC) fElectrons=new AliDxHFEParticleSelectionMCEl;
+  else fElectrons=new AliDxHFEParticleSelectionEl;
   fElectrons->SetCuts(fPID, AliDxHFEParticleSelectionEl::kCutPID);
   fElectrons->SetCuts(fCutsHFE, AliDxHFEParticleSelectionEl::kCutHFE);
-  fElectrons->InitControlObjects();
+  fElectrons->Init();
 
   //Correlation
   fCorrelation=new AliDxHFECorrelation;
@@ -160,7 +164,7 @@ void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
 
   // Fix for merging:
   // Retrieving the individual objects created
-  // and storing them instead of fD0s,fElectrons etc.. 
+  // and storing them instead of fD0s, fElectrons etc.. 
   TList *list =(TList*)fD0s->GetControlObjects();
   TObject *obj=NULL;
 
index 71e4662c39e3bf38322013f1159ee19dc557db27..7d39fd65a6b5808c8364a4319e148814df6b0a8d 100644 (file)
@@ -26,6 +26,7 @@
 #include "AliAnalysisTaskDxHFEParticleSelection.h"
 #include "AliDxHFEParticleSelection.h"
 #include "AliDxHFEParticleSelectionD0.h"
+#include "AliDxHFEParticleSelectionMCD0.h"
 #include "AliAnalysisManager.h"
 #include "AliAnalysisCuts.h"
 #include "AliLog.h"
@@ -52,6 +53,7 @@ AliAnalysisTaskDxHFEParticleSelection::AliAnalysisTaskDxHFEParticleSelection(con
   , fCuts(NULL)
   , fSelector(NULL)
   , fUseMC(kFALSE)
+  , fFillOnlyD0D0bar(0)
 {
   // constructor
   //
@@ -101,22 +103,18 @@ void AliAnalysisTaskDxHFEParticleSelection::UserCreateOutputObjects()
   fOutput = new TList;
   fOutput->SetOwner();
 
-  std::auto_ptr<AliDxHFEParticleSelection> selector(new AliDxHFEParticleSelectionD0);
-  if (!selector.get()) return;
-  if (fCuts) {
-    AliRDHFCutsD0toKpi* cuts=dynamic_cast<AliRDHFCutsD0toKpi*>(fCuts);
-    if (!cuts) {
-      AliFatal(Form("cut object %s is of incorrect type %s, expecting AliRDHFCutsD0toKpi", fCuts->GetName(), fCuts->ClassName()));
-      return;
-    }
-    selector->SetCuts(fCuts);
-    AliInfo(Form("Initializing particle selection %s, using %d pt bins", selector->GetName(), cuts->GetNPtBins()));
-  } else {
-    AliWarning("no cut object available for particle selection");
+  // setting up for D0s
+  TString selectionD0Options;
+  switch (fFillOnlyD0D0bar) {
+  case 1: selectionD0Options+="FillOnlyD0 "; break;
+  case 2: selectionD0Options+="FillOnlyD0bar "; break;
+  default: selectionD0Options+="FillD0D0bar ";
   }
-  selector->InitControlObjects();
 
-  fSelector=selector.release();
+  if(fUseMC) fSelector=new AliDxHFEParticleSelectionMCD0(selectionD0Options);
+  else fSelector=new AliDxHFEParticleSelectionD0(selectionD0Options);
+  fSelector->SetCuts(fCuts);
+  fSelector->Init();
 
   // Retrieving the list containing histos and THnSparse
   // and storing them instead of fSelector
@@ -186,7 +184,7 @@ void AliAnalysisTaskDxHFEParticleSelection::UserExec(Option_t* /*option*/)
     return;
   }
 
-  Int_t nInD0toKpi = inputArray->GetEntriesFast();
+  //  Int_t nInD0toKpi = inputArray->GetEntriesFast();
 
   fSelector->HistogramEventProperties(AliDxHFEParticleSelection::kEventsSel);
              
@@ -206,14 +204,12 @@ void AliAnalysisTaskDxHFEParticleSelection::UserExec(Option_t* /*option*/)
 
   //Test to see if I have read in D0s and retrieved them after selection
   Int_t nD0Selected = selectedTracks->GetEntriesFast();
-  printf("Number of D0->Kpi Start: %d , End: %d\n",nInD0toKpi,nD0Selected);
 
-  fSelector->HistogramEventProperties(AliDxHFEParticleSelection::kEventsD0);
+  fSelector->HistogramEventProperties(AliDxHFEParticleSelection::kEventsWithParticle);
 
   for(Int_t iD0toKpi = 0; iD0toKpi < nD0Selected; iD0toKpi++) {
     AliAODRecoDecayHF2Prong *particle = (AliAODRecoDecayHF2Prong*)selectedTracks->UncheckedAt(iD0toKpi);
     if (!particle) continue;
-    cout << "D0s inv mass: " << particle->InvMassD0() << endl;
   }
 
   PostData(1, fOutput);
index 1bb1bdd8a43563fbb64c7904d4a352be28655243..11ee0a311a6f8d67b8cb80a05f36857aaf93d082 100644 (file)
@@ -49,6 +49,7 @@ class AliAnalysisTaskDxHFEParticleSelection : public AliAnalysisTaskSE {
   virtual void Terminate(Option_t*);
 
   void SetOption(const char* opt) { fOption = opt; }
+  void SetFillOnlyD0D0bar(Int_t flagfill){fFillOnlyD0D0bar=flagfill;}
   virtual void SetUseMC(Bool_t useMC){fUseMC=useMC;}
   virtual void SetCuts(AliAnalysisCuts* cuts){fCuts=cuts;}
   Bool_t GetUseMC() const {return fUseMC;}
@@ -68,6 +69,8 @@ class AliAnalysisTaskDxHFEParticleSelection : public AliAnalysisTaskSE {
   AliAnalysisCuts* fCuts;               //  Cuts 
   AliDxHFEParticleSelection* fSelector; // selector instance
   bool fUseMC;                          // use MC info
+  Int_t     fFillOnlyD0D0bar;            // flag to set what to fill (0 = both, 1 = D0 only, 2 = D0bar only)
+
 
   ClassDef(AliAnalysisTaskDxHFEParticleSelection, 2);
 };
index c0fd9ab53587cbde46acbb1e2304c239c625ff33..446c21648ac78282786754e7558ae3d967c62ed4 100644 (file)
@@ -67,7 +67,7 @@ AliDxHFECorrelation::AliDxHFECorrelation(const char* name)
 const char* AliDxHFECorrelation::fgkEventControlBinNames[]={
   "nEventsAll",
   "nEventsSelected",
-  "nEventsD0"
+  "nEventsD0",
   "nEventsD0e"
 };
 
@@ -162,7 +162,7 @@ int AliDxHFECorrelation::Init()
   }
 
   for (iLabel=0; iLabel<kNEventControlLabels; iLabel++)
-    hEventControl->GetXaxis()->SetBinLabel(iLabel, fgkEventControlBinNames[iLabel]);
+    hEventControl->GetXaxis()->SetBinLabel(iLabel+1, fgkEventControlBinNames[iLabel]);
 
   fCorrProperties=CorrProperties.release();
   AddControlObject(fCorrProperties);
index 2a5201ca1a24947c6264362c3a4ecef1c6e61748..4d72677c14f45d71c7f74564dcded27326e14f11 100644 (file)
@@ -49,7 +49,7 @@ AliDxHFEParticleSelection::AliDxHFEParticleSelection(const char* name, const cha
   , fControlObjects(NULL)
   , fhEventControl(NULL)
   , fhTrackControl(NULL)
-  , fUseMC(false)
+  , fUseMC(0)
   , fVerbosity(0)
 {
   // constructor
@@ -62,7 +62,18 @@ AliDxHFEParticleSelection::AliDxHFEParticleSelection(const char* name, const cha
 const char* AliDxHFEParticleSelection::fgkEventControlBinNames[]={
   "nEventsAll",
   "nEventsSelected",
-  "nEventsD0"
+  "nEventsWParticle"
+};
+
+const char* AliDxHFEParticleSelection::fgkTrackControlBinNames[]={
+  "nTrackAll",
+  "nTrackSelected",
+};
+
+const char* AliDxHFEParticleSelection::fgkTrackControlDimNames[]={
+  "Pt",
+  "Phi",
+  "Eta"
 };
 
 AliDxHFEParticleSelection::~AliDxHFEParticleSelection()
@@ -76,9 +87,23 @@ AliDxHFEParticleSelection::~AliDxHFEParticleSelection()
   fhTrackControl=NULL;
 }
 
+int AliDxHFEParticleSelection::Init()
+{
+  //
+  // Init part sel. Calls InitControlObjects()
+  //
+
+  InitControlObjects();
+  return 0;
+}
+
+
 int AliDxHFEParticleSelection::InitControlObjects()
 {
+  //
   // init control objects
+  // TODO: Change to private now that have Init()?
+  //
   if (fVerbosity>0) {
     AliInfo("Setting up control objects");
   }
@@ -90,14 +115,63 @@ int AliDxHFEParticleSelection::InitControlObjects()
 
   fhEventControl=hEventControl.release();
   for (int iLabel=0; iLabel<kNEventPropertyLabels; iLabel++)
-    fhEventControl->GetXaxis()->SetBinLabel(iLabel, fgkEventControlBinNames[iLabel]);
+    fhEventControl->GetXaxis()->SetBinLabel(iLabel+1, fgkEventControlBinNames[iLabel]);
   AddControlObject(fhEventControl);
   fhTrackControl=hTrackControl.release();
+  for (int iLabel=0; iLabel<kNTrackPropertyLabels; iLabel++)
+    fhTrackControl->GetXaxis()->SetBinLabel(iLabel+1, fgkTrackControlBinNames[iLabel]);
   AddControlObject(fhTrackControl);
 
   return 0;
 }
 
+THnSparse* AliDxHFEParticleSelection::CreateControlTHnSparse(const char* name,
+                                                            int thnSize,
+                                                            int* thnBins,
+                                                            double* thnMin,
+                                                            double* thnMax,
+                                                            const char** binLabels) const
+{
+  //
+  // Creates THnSparse. 
+  //
+
+  AliInfo("Setting up THnSparse");
+
+  std::auto_ptr<THnSparseF> th(new THnSparseF(name, name, thnSize, thnBins, thnMin, thnMax));
+  if (th.get()==NULL) {
+    return NULL;
+  }
+  for (int iLabel=0; iLabel<thnSize; iLabel++) {
+    th->GetAxis(iLabel)->SetTitle(binLabels[iLabel]);    
+   
+  }
+ return th.release();
+
+}
+
+THnSparse* AliDxHFEParticleSelection::DefineTHnSparse() const
+{
+  //
+  // Defines the THnSparse. For now, only calls CreatControlTHnSparse
+
+  //TODO: Should make it more general. Or maybe one can use this here, and skip in PartSelEl?
+
+  const int thnSize = 3;
+  const double Pi=TMath::Pi();
+  TString name;
+  //                          0    1       2
+  //                          Pt   Phi    Eta
+  int    thnBins[thnSize] = { 1000,  200, 500};
+  double thnMin [thnSize] = {    0,    0, -1.};
+  double thnMax [thnSize] = {  100, 2*Pi,  1.};
+
+  name.Form("%s info", GetName());
+
+  return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,fgkTrackControlDimNames);
+
+}
+
 int AliDxHFEParticleSelection::AddControlObject(TObject* pObj)
 {
   /// add control object to list, the base class becomes owner of the object
@@ -132,15 +206,34 @@ int AliDxHFEParticleSelection::HistogramEventProperties(int bin)
 int AliDxHFEParticleSelection::HistogramParticleProperties(AliVParticle* p, int selected)
 {
   /// histogram particle properties
+
   if (!p) return -EINVAL;
   if (!fControlObjects) return 0;
 
   // TODO: use enums for the bins of the control histogram
-  fhTrackControl->Fill(0);
-  if (selected) fhTrackControl->Fill(1);
+  fhTrackControl->Fill(kTrackAll);
+  if (selected) fhTrackControl->Fill(kTrackSel);
   return 0;
 }
 
+int AliDxHFEParticleSelection::DefineParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
+{
+  // fill the data array from the particle data
+  if (!data) return -EINVAL;
+  int i=0;
+  // TODO: this corresponds to the THnSparse dimensions which is available in the same class
+  // use this consistently
+  const int requiredDimension=3;
+  if (dimension!=requiredDimension) {
+    // TODO: think about filling only the available data and throwing a warning
+    return -ENOSPC;
+  }
+  data[i++]=p->Pt();
+  data[i++]=p->Phi();
+  data[i++]=p->Eta();
+  return i;
+}
+
 TObjArray* AliDxHFEParticleSelection::Select(const AliVEvent* pEvent)
 {
   /// create selection from 'Tracks' member of the event,
@@ -152,7 +245,7 @@ TObjArray* AliDxHFEParticleSelection::Select(const AliVEvent* pEvent)
   int nofTracks=pEvent->GetNumberOfTracks();
   for (int itrack=0; itrack<nofTracks; itrack++) {
     AliVParticle* track=pEvent->GetTrack(itrack);
-    int selectionCode=IsSelected(track);
+    int selectionCode=IsSelected(track,pEvent);
     HistogramParticleProperties(track, selectionCode);
     if (selectionCode==0) continue;
     selectedTracks->Add(track);
index c6712e7e11defc85d01973909e261dc6842c0e3f..62f29accb2d37f08156a82523d03132ac01c4d78 100644 (file)
@@ -45,9 +45,15 @@ class AliDxHFEParticleSelection : public TNamed {
   enum {
     kEventsAll = 0,
     kEventsSel,
-    kEventsD0,
+    kEventsWithParticle,
     kNEventPropertyLabels
   };
+  
+  enum {
+    kTrackAll = 0,
+    kTrackSel,
+    kNTrackPropertyLabels
+  };
 
   /// set options
   void SetOption(const char* opt) { fOption = opt; }
@@ -55,6 +61,7 @@ class AliDxHFEParticleSelection : public TNamed {
   virtual Option_t* GetOption() const { return fOption;}
 
   /// init the control objects
+  virtual int Init();
   virtual int InitControlObjects();
 
   /// create selection from 'Tracks' member of the event,
@@ -72,6 +79,9 @@ class AliDxHFEParticleSelection : public TNamed {
   /// histogram event properties
   virtual int HistogramEventProperties(int bin);
 
+  // TODO: function can be renamed to better describe what it's doing
+  virtual int DefineParticleProperties(AliVParticle* p, Double_t* date, int dimension) const;
+
   /// check and add track to internal array
   int CheckAndAdd(AliVParticle* p);
 
@@ -79,6 +89,8 @@ class AliDxHFEParticleSelection : public TNamed {
   // different types; a type cast check is implemented in the method
   virtual void SetCuts(TObject* /*cuts*/, int /*level*/=0) {}
 
+  // TODO: check whether that is needed, should be covered by the specific
+  // child implementation
   Bool_t GetUseMC() const {return fUseMC;}
 
   /// get selected tracks
@@ -104,6 +116,16 @@ class AliDxHFEParticleSelection : public TNamed {
   void SetVerbosity(int verbosity) {fVerbosity=verbosity;}
   /// get verbosity
   int GetVerbosity() const {return fVerbosity;}
+  
+  /// create control THnSparse
+  THnSparse* CreateControlTHnSparse(const char* name,
+                                   int thnSize,
+                                   int* thnBins,
+                                   double* thnMin,
+                                   double* thnMax,
+                                   const char** binLabels) const;
+
+  virtual THnSparse* DefineTHnSparse() const;
 
  protected:
   /// add control object to list, the base class becomes owner of the object
@@ -112,12 +134,6 @@ class AliDxHFEParticleSelection : public TNamed {
   /// histogram particle properties
   virtual int HistogramParticleProperties(AliVParticle* p, int selected=1);
 
-  /// set from the specific child class implementation
-  /// TODO: think about removing this function
-  /// in principle that should not be necessary, the analysis should
-  /// be blind after the particle selection
-  void SetUseMC(bool useMC=true) {fUseMC=useMC;}
-
  private:
   /// copy contructor prohibited
   AliDxHFEParticleSelection(const AliDxHFEParticleSelection&);
@@ -135,6 +151,8 @@ class AliDxHFEParticleSelection : public TNamed {
   int fVerbosity;      //! verbosity
 
   static const char* fgkEventControlBinNames[]; //! bin labels for event control histogram
+  static const char* fgkTrackControlBinNames[]; //! bin labels for track control histogram
+  static const char* fgkTrackControlDimNames[]; //! bin labels fro track thnsparse
 
   ClassDef(AliDxHFEParticleSelection, 2);
 
index c6c36bd142f0aee3eed6779fcb70eacc0681f862..648d7da3477f8efd26eebe07304544caa0ec1ec2 100644 (file)
@@ -51,6 +51,8 @@ AliDxHFEParticleSelectionD0::AliDxHFEParticleSelectionD0(const char* opt)
   , fD0Daughter1(NULL)
   , fCuts(NULL)
   , fFillOnlyD0D0bar(0)
+  , fD0InvMass(0.0)
+  , fPtBin(-1)
 {
   // constructor
   // 
@@ -85,35 +87,28 @@ AliDxHFEParticleSelectionD0::~AliDxHFEParticleSelectionD0()
   fCuts=NULL;
 }
 
+const char* AliDxHFEParticleSelectionD0::fgkTrackControlBinNames[]={
+  "Pt",
+  "Phi",
+  "Ptbin", 
+  "D0InvMass", 
+  "Eta"
+};
+
+const char* AliDxHFEParticleSelectionD0::fgkDgTrackControlBinNames[]={
+  "Pt",
+  "Phi",
+  "Ptbin", 
+  "D0InvMass", 
+  "Eta"
+};
+
 int AliDxHFEParticleSelectionD0::InitControlObjects()
 {
   /// init the control objects, can be overloaded by childs which should
   /// call AliDxHFEParticleSelection::InitControlObjects() explicitly
-  AliInfo("Setting up THnSparse");
-  TString name;
-  const int thnSize = 4;
-  const double pi=TMath::Pi();
-
-  // TODO: theta?
-  //                            0       1     2      3
-  //                         mass      Pt   Phi    Ptbin
-  int    thnBins[thnSize] = {   200,   1000,  100,   100  };
-  double thnMin [thnSize] = {  1.5648,   0,    0,     0   };
-  double thnMax [thnSize] = {  2.1648, 100,  (2*pi), 100  };
-
-  name.Form("%s info", GetName());
-  std::auto_ptr<THnSparseF> D0Properties(new THnSparseF(name, name, thnSize, thnBins, thnMin, thnMax));
-
-  if (D0Properties.get()==NULL) {
-    return -ENOMEM;
-  }
-  int axis=0;
-  D0Properties->GetAxis(axis++)->SetTitle("D0 inv mass");
-  D0Properties->GetAxis(axis++)->SetTitle("Pt");
-  D0Properties->GetAxis(axis++)->SetTitle("Phi"); 
-  D0Properties->GetAxis(axis++)->SetTitle("Ptbin"); 
 
-  fD0Properties=D0Properties.release();
+  fD0Properties=DefineTHnSparse();
   AddControlObject(fD0Properties);
 
   //Adding control objects for the daughters
@@ -123,9 +118,55 @@ int AliDxHFEParticleSelectionD0::InitControlObjects()
   return AliDxHFEParticleSelection::InitControlObjects();
 }
 
+THnSparse* AliDxHFEParticleSelectionD0::DefineTHnSparse() const
+{
+  //
+  // Defines the THnSparse. For now, only calls CreateControlTHnSparse
+  // TODO: remove pt?? (Have ptbin)
+
+  const int thnSize2 = 5;
+  const double Pi=TMath::Pi();
+  TString name;
+  name.Form("%s info", GetName());
+
+  //                          0    1      2      3          4    
+  //                          Pt   Phi   Ptbin  D0InvMass  Eta  
+  int    thnBins[thnSize2] = { 1000,  200, 21,     200,     500 };
+  double thnMin [thnSize2] = {    0,    0,  0,    1.5648,   -1. };
+  double thnMax [thnSize2] = {  100, 2*Pi, 20,    2.1648,    1. };
+
+  return CreateControlTHnSparse(name,thnSize2,thnBins,thnMin,thnMax,fgkTrackControlBinNames);
+
+}
+
+int AliDxHFEParticleSelectionD0::DefineParticleProperties(AliVParticle* p, Double_t* data, int dimension) const
+{
+  // fill the data array from the particle data
+  if (!data) return -EINVAL;
+  //  AliAODTrack *track=(AliAODTrack*)p;
+  AliAODRecoDecayHF2Prong* track=dynamic_cast<AliAODRecoDecayHF2Prong*>(p);
+  if (!track) return -ENODATA;
+  int i=0;
+  // TODO: this corresponds to the THnSparse dimensions which is available in the same class
+  // use this consistently
+  const int requiredDimension=5;
+  if (dimension!=requiredDimension) {
+    // TODO: think about filling only the available data and throwing a warning
+    return -ENOSPC;
+  }
+  data[i++]=track->Pt();
+  data[i++]=track->Phi();
+  data[i++]=fPtBin;
+  data[i++]=fD0InvMass;
+  data[i++]=track->Eta();
+
+  return i;
+}
+
 int AliDxHFEParticleSelectionD0::InitControlObjectsDaughters(TString name, int daughter)
 {
-  //Setting up Control objects for the daughters.
+  // Setting up Control objects for the daughters.
+  // Move to ParticleSelection?? 
   AliInfo("Setting up daughter THnSparse");
 
   const int thnSize2 = 5;
@@ -141,12 +182,9 @@ int AliDxHFEParticleSelectionD0::InitControlObjectsDaughters(TString name, int d
   if (DaughterProperties.get()==NULL) {
     return -ENOMEM;
   }
-  int axis=0;
-  DaughterProperties->GetAxis(axis++)->SetTitle("Pt");
-  DaughterProperties->GetAxis(axis++)->SetTitle("Phi");
-  DaughterProperties->GetAxis(axis++)->SetTitle("Ptbin"); 
-  DaughterProperties->GetAxis(axis++)->SetTitle("D0InvMass"); 
-  DaughterProperties->GetAxis(axis++)->SetTitle("Eta"); 
+
+  for(int iLabel=0; iLabel< 5;iLabel++)
+    DaughterProperties->GetAxis(iLabel)->SetTitle(fgkDgTrackControlBinNames[iLabel]);  
 
   if(daughter==0){ 
     fD0Daughter0=DaughterProperties.release();
@@ -186,14 +224,17 @@ int AliDxHFEParticleSelectionD0::HistogramParticleProperties(AliVParticle* p, in
   // Only D0s are filled 
   // TODO: Also include D0bar
   if ((selectionCode==1 || selectionCode==3) && fFillOnlyD0D0bar<2) {
-    Double_t invmassD0 = part->InvMassD0();
-    AliDebug(3,Form("pt %f,  nr pt bins %d",part->Pt(),fCuts->GetNPtBins()));
-    Int_t ptbin=fCuts->PtBin(part->Pt());
-    Double_t D0Properties[] = {invmassD0,part->Pt(),part->Phi(),ptbin};
-    Double_t KProperties[]={prongneg->Pt(),prongneg->Phi(),ptbin, invmassD0,prongneg->Eta()};
-    Double_t piProperties[]={prongpos->Pt(),prongpos->Phi(),ptbin,invmassD0,prongpos->Eta()};
-
-    if(fD0Properties) fD0Properties->Fill(D0Properties);
+    fD0InvMass= part->InvMassD0();
+    fPtBin=fCuts->PtBin(part->Pt());
+
+    // TODO: avoid repeated allocation of the arrays
+    Double_t KProperties[]={prongneg->Pt(),prongneg->Phi(),fPtBin, fD0InvMass,prongneg->Eta()};
+    Double_t piProperties[]={prongpos->Pt(),prongpos->Phi(),fPtBin,fD0InvMass,prongpos->Eta()};
+
+    // TODO: make array a member, consistent dimensions for THnSparse and array
+    Double_t d0Properties[5]={0.0, 0.0, 0.0, 0.0, 0.0};
+    DefineParticleProperties(p, d0Properties, 5);
+    if(fD0Properties) fD0Properties->Fill(d0Properties);
     if(fD0Daughter0) fD0Daughter0->Fill(piProperties);
     if(fD0Daughter1) fD0Daughter1->Fill(KProperties);
   }
@@ -215,8 +256,8 @@ TObjArray* AliDxHFEParticleSelectionD0::Select(TObjArray* pTracks, const AliVEve
     if (!track) continue;
     int selectionCode=IsSelected(track,pEvent);
     HistogramParticleProperties(track, selectionCode);
-    //TODO: Also add selection for D0bar
 
+    //TODO: Also add selection for D0bar
     // Add track if it is either defined as D0(selectionCode==1) or both 
     // D0bar and a D0 (selectionCode==3)
     if (! ((selectionCode==1 || selectionCode==3) && fFillOnlyD0D0bar<2)) continue;
@@ -228,7 +269,9 @@ TObjArray* AliDxHFEParticleSelectionD0::Select(TObjArray* pTracks, const AliVEve
 int AliDxHFEParticleSelectionD0::IsSelected(AliVParticle* p, const AliVEvent* pEvent)
 {
   /// TODO: implement specific selection of D0 candidates
-  /// Could also return values based on where where selection "failed"
+  /// Could also return values based on where where selection "failed
+  /// Selected. Return 0 (none), 1(D0), 2(D0bar) or 3 (both)
+
   int selectionCode=0;
 
   AliAODRecoDecayHF2Prong *d0 = dynamic_cast<AliAODRecoDecayHF2Prong*>(p);
@@ -241,9 +284,9 @@ int AliDxHFEParticleSelectionD0::IsSelected(AliVParticle* p, const AliVEvent* pE
   // AliRDHFCuts::IsSelected does not allow this
   AliRDHFCuts* cuts=const_cast<AliRDHFCuts*>(fCuts);
   if (!cuts) {
-    selectionCode=1;
+    selectionCode=0;
   } else if(cuts->IsInFiducialAcceptance(d0->Pt(),d0->Y(421)) ) {
-    //    if(cuts->IsSelected(d0,AliRDHFCuts::kTracks,pEvent))fNentries->Fill(6);       
+
     Int_t ptbin=cuts->PtBin(d0->Pt());
     if(ptbin==-1) {
       AliDebug(1,"Pt out of bounds");
@@ -255,7 +298,7 @@ int AliDxHFEParticleSelectionD0::IsSelected(AliVParticle* p, const AliVEvent* pE
     AliAODEvent* aod=NULL;
     if (pEvent) aod=dynamic_cast<AliAODEvent*>(const_cast<AliVEvent*>(pEvent));
   
-    // Selected. Return 0 (none), 1(D0), 2(D0bar) or 3 (both)
+    // Selected. Return 0 (none), 1 (D0), 2 (D0bar) or 3 (both)
     selectionCode=cuts->IsSelected(d0,AliRDHFCuts::kAll,aod); 
 
     AliDebug(1,Form("Candidate is %d \n", selectionCode));
index 93c43f151a37187e44c906bf722253c386ed1919..f0bacb203d62b28f8716cd49a7b464dc923fe85e 100644 (file)
@@ -33,6 +33,9 @@ class AliDxHFEParticleSelectionD0 : public AliDxHFEParticleSelection {
 
   /// overloaded from AliDxHFEParticleSelection: init the control objects
   virtual int InitControlObjects();
+  virtual THnSparse* DefineTHnSparse() const;
+  // TODO: function can be renamed to better describe what it's doing
+  virtual int DefineParticleProperties(AliVParticle* p, Double_t* date, int dimension) const;
 
   //Function for daughter control objects
   //TODO: move to AliDxHFEParticleSelection to be used for several particles?
@@ -49,6 +52,8 @@ class AliDxHFEParticleSelectionD0 : public AliDxHFEParticleSelection {
 
   //AliRDHFCutsD0toKpi GetCuts()const {return fCuts;}
   Int_t  GetFillOnlyD0D0bar() const {return fFillOnlyD0D0bar;}
+  Int_t GetPtBin() const {return fPtBin;}
+  Double_t GetInvMass() const {return fD0InvMass;}
 
  protected:
   /// overloaded from AliDxHFEParticleSelection: histogram particle properties
@@ -65,6 +70,11 @@ class AliDxHFEParticleSelectionD0 : public AliDxHFEParticleSelection {
   THnSparse* fD0Daughter1;      //  the particle properties of selected particles
   AliRDHFCuts* fCuts;           //! pointer to external cuts object 
   Int_t     fFillOnlyD0D0bar;   //  flag to set what to fill (0 = both, 1 = D0 only, 2 = D0bar only)
+  Double_t fD0InvMass;          // D0InvMass
+  Int_t fPtBin;                 // Pt Bin
+
+  static const char* fgkTrackControlBinNames[]; //! bin labels for track control histogram
+  static const char* fgkDgTrackControlBinNames[]; //! bin labels for track control histogram
 
   // TODO: at the moment the dimensions of the different THnSparse objects are different
   // needs to be consolidated
index 0835d9740f6ec68bab8123e9b5055c7d17fadf02..27677052b494160144d26e70ff363b17852705c4 100644 (file)
@@ -76,42 +76,43 @@ AliDxHFEParticleSelectionEl::~AliDxHFEParticleSelectionEl()
     delete fCFM;
     fCFM=NULL;
   }
+  if(fWhichCut){
+    delete fWhichCut;
+    fWhichCut=NULL;
+  }
 
   // NOTE: external objects fPID and fCuts are not deleted here
   fPID=NULL;
   fCuts=NULL;
 }
-int AliDxHFEParticleSelectionEl::InitControlObjects()
-{
-  /// init control and monitoring objects
-  AliInfo("Electron THnSparse");
-  const int thnSize = 3;
-  const double Pi=TMath::Pi();
-  TString name;// ="e information";
-  //                          0    1       2
-  //                          Pt   Phi    Eta
-  int    thnBins[thnSize] = { 1000,  200, 500};
-  double thnMin [thnSize] = {    0,    0, -1.};
-  double thnMax [thnSize] = {  100, 2*Pi,  1.};
 
-  name.Form("%s info", GetName());
-  std::auto_ptr<THnSparseF> ElectronProperties(new THnSparseF(name, name, thnSize, thnBins, thnMin, thnMax));
+const char* AliDxHFEParticleSelectionEl::fgkTrackControlBinNames[]={
+  "Pt",
+  "Phi",
+  "Eta"
+};
 
-  if (ElectronProperties.get()==NULL) {
-    return -ENOMEM;
-  }
-  int axis=0;
-  ElectronProperties->GetAxis(axis++)->SetTitle("Pt");
-  ElectronProperties->GetAxis(axis++)->SetTitle("Phi");
-  ElectronProperties->GetAxis(axis++)->SetTitle("Eta"); 
+const char* AliDxHFEParticleSelectionEl::fgkCutBinNames[]={
+  "kRecKineITSTPC",
+  "kRecPrim",
+  "kHFEcuts",
+  "kHFEcutsTOFPID",
+  "kHFEcutsTPCPID",
+  "kPID",
+  "Selected e"
+};
 
-  fElectronProperties=ElectronProperties.release();
+int AliDxHFEParticleSelectionEl::Init()
+{
+  //
+  // init function
+  // 
+  int iResult=0;
 
-  AddControlObject(fElectronProperties);
-
-  fWhichCut= new TH1F("fWhichCut","effective cut for a rejected particle",6,-0.5,5.5);
-  AddControlObject(fWhichCut);
+  // Implicit call to InitControlObjects() before setting up CFM and fCuts
+  // (if not there)
+  iResult=AliDxHFEParticleSelection::Init();
+  if (iResult<0) return iResult;
 
   //--------Initialize correction Framework and Cuts
   // Consider moving this, either to separate function or
@@ -137,6 +138,43 @@ int AliDxHFEParticleSelectionEl::InitControlObjects()
   fCuts->Initialize(fCFM);
 
   return 0;
+
+}
+
+THnSparse* AliDxHFEParticleSelectionEl::DefineTHnSparse() const
+{
+  //
+  // Defines the THnSparse. For now, only calls CreatControlTHnSparse
+  const int thnSize = 3;
+  const double Pi=TMath::Pi();
+  TString name;
+  name.Form("%s info", GetName());
+
+  //                          0    1       2
+  //                          Pt   Phi    Eta
+  int    thnBins[thnSize] = { 1000,  200, 500};
+  double thnMin [thnSize] = {    0,    0, -1.};
+  double thnMax [thnSize] = {  100, 2*Pi,  1.};
+
+  return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,fgkTrackControlBinNames);
+
+}
+
+int AliDxHFEParticleSelectionEl::InitControlObjects()
+{
+  /// init control and monitoring objects
+  AliInfo("Electron THnSparse");
+
+  fElectronProperties=DefineTHnSparse();
+  AddControlObject(fElectronProperties);
+
+  //
+  fWhichCut= new TH1F("fWhichCut","effective cut for a rejected particle",kNCutLabels,-0.5,kNCutLabels-0.5);
+  for (int iLabel=0; iLabel<kNCutLabels; iLabel++)
+    fWhichCut->GetXaxis()->SetBinLabel(iLabel+1, fgkCutBinNames[iLabel]);
+  AddControlObject(fWhichCut);
+
+  return AliDxHFEParticleSelection::InitControlObjects();
 }
 
 int AliDxHFEParticleSelectionEl::HistogramParticleProperties(AliVParticle* p, int selectionCode)
@@ -145,17 +183,40 @@ int AliDxHFEParticleSelectionEl::HistogramParticleProperties(AliVParticle* p, in
   if (!p) return -EINVAL;
   //if (!fControlObjects) return 0;
   if(selectionCode==0) return  0;
-
-  AliAODTrack *track=(AliAODTrack*)p;
-  Double_t eProperties[]={track->Pt(),track->Phi(),track->Eta()};
+  // TODO: make array a member, consistent dimensions for THnSparse and array
+  Double_t eProperties[3]={0.0, 0.0, 0.0};
+  DefineParticleProperties(p, eProperties, 3);
   if(fElectronProperties) fElectronProperties->Fill(eProperties);
-  
+    
   return 0;
 }
 
+int AliDxHFEParticleSelectionEl::DefineParticleProperties(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;
+  int i=0;
+  // TODO: this corresponds to the THnSparse dimensions which is available in the same class
+  // use this consistently
+  const int requiredDimension=3;
+  if (dimension!=requiredDimension) {
+    // 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();
+  return i;
+}
+
+
 int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent*)
 {
   /// select El candidates
+  // TODO: How to handle MC? would be too much duplicate code if copy entire IsSelected. 
+
   AliAODTrack *track=(AliAODTrack*)pEl;
  
 
@@ -163,31 +224,31 @@ int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent*)
   //Using AliHFECuts:
   // RecKine: ITSTPC cuts  
   if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)){
-    fWhichCut->Fill(0);
+    fWhichCut->Fill(kRecKineITSTPC);
     return 0;
   }
   
   // RecPrim
   if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) {
-    fWhichCut->Fill(1);
+    fWhichCut->Fill(kRecPrim);
     return 0;
   }
   
   // HFEcuts: ITS layers cuts
   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) {
-    fWhichCut->Fill(2);
+    fWhichCut->Fill(kHFEcutsITS);
     return 0;
   }
   
   // HFE cuts: TOF PID and mismatch flag
   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) {
-    fWhichCut->Fill(3);
+    fWhichCut->Fill(kHFEcutsTOF);
     return 0;
   }
   
   // HFE cuts: TPC PID cleanup
   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)){
-    fWhichCut->Fill(4);
+    fWhichCut->Fill(kHFEcutsTPC);
     return 0;
   } 
  
@@ -206,11 +267,11 @@ int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent*)
 
   if(fPID && fPID->IsSelected(&hfetrack)) {
     AliDebug(3,"Inside FilldPhi, electron is selected");
-
+    fWhichCut->Fill(kSelected);
     return 1;
   }
   else{
-    fWhichCut->Fill(5);
+    fWhichCut->Fill(kPID);
     return 0;
   }
 }
index e3a91b665cb0609d6ad3c9e9b85bec547e5165e8..914baf85d5e05144291cef3ee928b04b04509970 100644 (file)
@@ -44,16 +44,31 @@ class AliDxHFEParticleSelectionEl : public AliDxHFEParticleSelection {
     kNCuts
   };
 
+  enum {
+    kRecKineITSTPC=0,
+    kRecPrim,
+    kHFEcutsITS,
+    kHFEcutsTOF,
+    kHFEcutsTPC,
+    kPID,
+    kSelected,
+    kNCutLabels
+  };
+
+  ///overloaded from AliDxHFEParticleSelection: Init
+  virtual int Init();
+
   /// overloaded from AliDxHFEParticleSelection: init the control objects
   virtual int InitControlObjects();
+  virtual THnSparse* DefineTHnSparse() const;
 
   /// overloaded from AliDxHFEParticleSelection: check particle
   virtual int IsSelected(AliVParticle* p, const AliVEvent*);
 
   virtual int HistogramParticleProperties(AliVParticle* p, int selected);
 
-  // overloaded from AliDxHFEParticleSelection: specific for electrons
-  //virtual TObjArray* Select(const AliVEvent *pEvent);
+  // TODO: function can be renamed to better describe what it's doing
+  virtual int DefineParticleProperties(AliVParticle* p, Double_t* date, int dimension) const;
 
   /// set cuts object: a type cast check is implemented in the method
   virtual void SetCuts(TObject* /*cuts*/, int /*level*/=0);
@@ -75,6 +90,8 @@ class AliDxHFEParticleSelectionEl : public AliDxHFEParticleSelection {
   AliHFEcuts*   fCuts;               //! Cuts for HF electrons
   AliCFManager* fCFM;                //! Correction Framework Manager
 
+  static const char* fgkTrackControlBinNames[]; //! bin labels for track control histogram
+  static const char* fgkCutBinNames[]; //! bin labels for cuts histogram
 
   ClassDef(AliDxHFEParticleSelectionEl, 2);
 };
diff --git a/PWGHF/correlationHF/AliDxHFEParticleSelectionMCD0.cxx b/PWGHF/correlationHF/AliDxHFEParticleSelectionMCD0.cxx
new file mode 100644 (file)
index 0000000..1f843fe
--- /dev/null
@@ -0,0 +1,197 @@
+// $Id$
+
+//**************************************************************************
+//* This file is property of and copyright by the ALICE Project            * 
+//* ALICE Experiment at CERN, All rights reserved.                         *
+//*                                                                        *
+//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
+//*                  Hege Erdal       <hege.erdal@gmail.com>               *
+//*                                                                        *
+//* 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.                  *
+//**************************************************************************
+
+/// @file   AliDxHFEParticleSelectionMCD0.cxx
+/// @author Hege Erdal, Matthias Richter
+/// @date   2012-07-19
+/// @brief  MC D0 selection for D0-HFE correlation
+///
+
+#include "AliDxHFEParticleSelectionMCD0.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliVParticle.h"
+#include <iostream>
+#include <cerrno>
+#include <memory>
+
+using namespace std;
+
+/// ROOT macro for the implementation of ROOT specific class methods
+ClassImp(AliDxHFEParticleSelectionMCD0)
+
+// TODO: can that be a common definition
+const char* AliDxHFEParticleSelectionMCD0::fgkTrackControlBinNames[]={
+"Pt",
+"Phi",
+"Ptbin", 
+"D0InvMass", 
+"Eta",
+"Statistics D0",
+"Mother of D0"
+};
+
+AliDxHFEParticleSelectionMCD0::AliDxHFEParticleSelectionMCD0(const char* opt)
+  : AliDxHFEParticleSelectionD0(opt)
+  , fMCTools()
+  , fResultMC(0)
+  , fOriginMother(0)
+{
+  // constructor
+  // 
+  // 
+  // 
+  // 
+
+  // TODO: argument scan, pass only relevant arguments to tools
+  fMCTools.~AliDxHFEToolsMC();
+  TString toolopt("pdg=421 mc-last");
+  new (&fMCTools) AliDxHFEToolsMC(toolopt);
+}
+
+AliDxHFEParticleSelectionMCD0::~AliDxHFEParticleSelectionMCD0()
+{
+  // destructor
+}
+THnSparse* AliDxHFEParticleSelectionMCD0::DefineTHnSparse() const
+{
+  //
+  // Defines the THnSparse.
+  // could/should remove Pt and leave Ptbin
+
+  const int thnSize2 = 7;
+  const double Pi=TMath::Pi();
+  TString name;
+  name.Form("%s info", GetName());
+
+  //                          0    1      2      3          4     5         6      
+  //                          Pt   Phi   Ptbin  D0InvMass  Eta  'stat D0'  mother 
+  int    thnBins[thnSize2] = { 1000,  200, 21,     200,     500,     2,       10  };
+  double thnMin [thnSize2] = {    0,    0,  0,    1.5648,   -1.,  -0.5,     -1.5  };
+  double thnMax [thnSize2] = {  100, 2*Pi, 20,    2.1648,    1.,   1.5,      8.5  };
+
+  return CreateControlTHnSparse(name,thnSize2,thnBins,thnMin,thnMax,fgkTrackControlBinNames);
+
+}
+
+int AliDxHFEParticleSelectionMCD0::DefineParticleProperties(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;
+  int i=0;
+  // TODO: this corresponds to the THnSparse dimensions which is available in the same class
+  // use this consistently
+  const int requiredDimension=7;
+  if (dimension!=requiredDimension) {
+    // TODO: think about filling only the available data and throwing a warning
+    return -ENOSPC;
+  }
+  data[i++]=track->Pt();
+  data[i++]=track->Phi();
+  data[i++]=AliDxHFEParticleSelectionMCD0::GetPtBin(); 
+  data[i++]=AliDxHFEParticleSelectionMCD0::GetInvMass();
+  data[i++]=track->Eta();
+  data[i++]=fResultMC;     // stat electron (MC electron or not)
+  data[i++]=fOriginMother; // at the moment not included background. Should expand
+
+  return i;
+}
+
+int AliDxHFEParticleSelectionMCD0::IsSelected(AliVParticle* p, const AliVEvent* pEvent)
+{
+  /// overloaded from AliDxHFEParticleSelection: check particle
+  /// H: Have changed function. Now doing particle selection first, then run MC over 
+  /// selected tracks. Could configure it to be configurable, but not sure if it
+  /// is needed.  
+  /// result from normal track selection is returned, result from MC is stored in
+  /// THnSparse. 
+
+  int iResult=0;
+  fOriginMother=-1;
+
+  // step 1:
+  // MC selection
+  if (fMCTools.MCFirst() && (iResult=CheckMC(p, pEvent))==0) {
+    // histograming?
+    return iResult;
+  }
+
+  // step 2 or 1, depending on sequence:
+  // normal particle selection
+  iResult=AliDxHFEParticleSelectionD0::IsSelected(p, pEvent);
+  if (fMCTools.MCFirst() || iResult==0) return iResult;
+
+  // step 2, only executed if MC check is last
+  // MC selection  - > Should maybe also distinguish between D0 and D0bar
+  iResult=CheckMC(p, pEvent);
+  // TODO: why do we need to store the result in a member?
+  fResultMC=iResult;
+  return iResult;
+}
+
+int AliDxHFEParticleSelectionMCD0::CheckMC(AliVParticle* p, const AliVEvent* pEvent)
+{
+  /// check if MC criteria are fulfilled
+  // Check both D0 and D0bar (for now only D0)
+
+  if (!p || !pEvent){
+    return -EINVAL;
+  }
+  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
+  }
+
+  AliAODRecoDecayHF2Prong *particle = dynamic_cast<AliAODRecoDecayHF2Prong*>(p);
+
+  if(!particle) return 0;
+
+  Int_t pdgDgD0toKpi[2]={AliDxHFEToolsMC::kPDGkaon,AliDxHFEToolsMC::kPDGpion};
+
+  TClonesArray* fMCArray = dynamic_cast<TClonesArray*>(fMCTools.GetMCArray());
+  if(!fMCArray) {cout << "no array" << endl; return -1;}
+
+  // find associated MC particle for D0->Kpi
+  Int_t MClabel=-9999;
+
+  //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx). Checks both D0s and daughters
+  MClabel=particle->MatchToMC(AliDxHFEToolsMC::kPDGD0,fMCArray,2,pdgDgD0toKpi); 
+  
+  if(MClabel<0){
+    fOriginMother=-1;
+    return 0;
+  }
+
+  fMCTools.SetMClabel(MClabel);
+  fMCTools.FindMotherPDG(p,AliDxHFEToolsMC::kGetOriginMother);
+  fOriginMother=fMCTools.GetOriginMother();
+
+  return 1;
+}
+
+void AliDxHFEParticleSelectionMCD0::Clear(const char* option)
+{
+  /// clear internal memory
+  fMCTools.Clear(option);
+}
diff --git a/PWGHF/correlationHF/AliDxHFEParticleSelectionMCD0.h b/PWGHF/correlationHF/AliDxHFEParticleSelectionMCD0.h
new file mode 100644 (file)
index 0000000..6a3a9cd
--- /dev/null
@@ -0,0 +1,62 @@
+//-*- Mode: C++ -*-
+// $Id$
+
+//* This file is property of and copyright by the ALICE Project        * 
+//* ALICE Experiment at CERN, All rights reserved.                     *
+//* See cxx source for full Copyright notice                           *
+
+/// @file   AliDxHFEParticleSelectionMCD0.h
+/// @author Hege Erdal, Matthias Richter
+/// @date   2012-07-19
+/// @brief  D0 MC selection for D-HFE correlations
+///
+
+#ifndef ALIDXHFEPARTICLESELECTIONMCD0_H
+#define ALIDXHFEPARTICLESELECTIONMCD0_H
+
+#include "AliDxHFEParticleSelectionD0.h"
+#include "AliDxHFEToolsMC.h"
+
+/**
+ * @class AliDxHFEParticleSelectionMCD0
+ * Monte Carlo D0 selection for D-HFE correlations, implements the specific
+ * selection criteria.
+ */
+class AliDxHFEParticleSelectionMCD0 : public AliDxHFEParticleSelectionD0 {
+  public:
+  /// constructor
+  AliDxHFEParticleSelectionMCD0(const char* opt="");
+  /// destructor
+  virtual ~AliDxHFEParticleSelectionMCD0();
+
+  /// overloaded from AliDxHFEParticleSelection: check particle
+  virtual int IsSelected(AliVParticle* p, const AliVEvent *pEvent=NULL);
+
+  virtual THnSparse* DefineTHnSparse() const;
+  // TODO: function can be renamed to better describe what it's doing
+  virtual int DefineParticleProperties(AliVParticle* p, Double_t* date, int dimension) const;
+
+  /// check MC criteria
+  int CheckMC(AliVParticle* p, const AliVEvent* pEvent);
+
+  /// clear internal memory
+  virtual void Clear(const char* option="");
+
+ protected:
+
+ private:
+  /// copy contructor prohibited
+  AliDxHFEParticleSelectionMCD0(const AliDxHFEParticleSelectionMCD0&);
+  /// assignment operator prohibited
+  AliDxHFEParticleSelectionMCD0& operator=(const AliDxHFEParticleSelectionMCD0&);
+
+  AliDxHFEToolsMC fMCTools;  // MC selction tools
+  int fResultMC;             // Result on MC check
+  int fOriginMother;         // Holds info on the original mother particle
+  static const char* fgkTrackControlBinNames[]; //! bin labels for track control histogram
+
+
+  ClassDef(AliDxHFEParticleSelectionMCD0, 1);
+};
+
+#endif
diff --git a/PWGHF/correlationHF/AliDxHFEParticleSelectionMCEl.cxx b/PWGHF/correlationHF/AliDxHFEParticleSelectionMCEl.cxx
new file mode 100644 (file)
index 0000000..2c7418f
--- /dev/null
@@ -0,0 +1,235 @@
+// $Id$
+
+//**************************************************************************
+//* This file is property of and copyright by the ALICE Project            * 
+//* ALICE Experiment at CERN, All rights reserved.                         *
+//*                                                                        *
+//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
+//*                  Hege Erdal       <hege.erdal@gmail.com>               *
+//*                                                                        *
+//* 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.                  *
+//**************************************************************************
+
+/// @file   AliDxHFEParticleSelectionMCEl.cxx
+/// @author Hege Erdal, Matthias Richter
+/// @date   2012-07-19
+/// @brief  MC El selection for D0-HFE correlation
+///
+
+#include "AliDxHFEParticleSelectionMCEl.h"
+#include "AliVParticle.h"
+#include "AliLog.h"
+#include "THnSparse.h"
+#include "AliAODMCParticle.h"
+#include "TH1F.h"
+#include "TAxis.h"
+#include "AliAODTrack.h"
+#include <iostream>
+#include <cerrno>
+#include <memory>
+
+using namespace std;
+using std::vector;
+
+
+/// ROOT macro for the implementation of ROOT specific class methods
+ClassImp(AliDxHFEParticleSelectionMCEl)
+
+AliDxHFEParticleSelectionMCEl::AliDxHFEParticleSelectionMCEl(const char* opt)
+  : AliDxHFEParticleSelectionEl(opt)
+  , fMCTools()
+  , fElectronProperties(NULL)
+  , fWhichCut(NULL)
+  , fOriginMother(0)
+  , fResultMC(0)
+{
+  // constructor
+  // 
+  // 
+  // 
+  // 
+
+  fMCTools.~AliDxHFEToolsMC();
+
+  // TODO: argument scan, build tool options accordingly
+  // e.g. set mc mode first/last, skip control histograms
+  TString toolopt("pdg=11 mc-last");
+  new (&fMCTools) AliDxHFEToolsMC(toolopt);
+}
+
+//at the moment not used, keep for now (to be used when plotting)
+const char* AliDxHFEParticleSelectionMCEl::fgkPDGMotherBinLabels[]={
+  "d",
+  "u",
+  "s",
+  "c",
+  "b",
+  "gluon",
+  "gamma",
+  "#pi^{0}",
+  "#eta",
+  "proton",
+  "others"
+};
+
+const char* AliDxHFEParticleSelectionMCEl::fgkTrackControlBinNames[]={
+  "Pt",
+  "Phi",
+  "Eta", 
+  "MC statistics electron"
+  "Mother", 
+  "PDG of selected electron", 
+};
+
+AliDxHFEParticleSelectionMCEl::~AliDxHFEParticleSelectionMCEl()
+{
+  // destructor
+}
+
+THnSparse* AliDxHFEParticleSelectionMCEl::DefineTHnSparse() const
+{
+  //
+  // Defines the THnSparse. 
+
+  const int thnSize = 6;
+  const double Pi=TMath::Pi();
+  TString name;
+  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 };
+  return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,fgkTrackControlBinNames);
+
+}
+
+int AliDxHFEParticleSelectionMCEl::DefineParticleProperties(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;
+  int i=0;
+  // TODO: this corresponds to the THnSparse dimensions which is available in the same class
+  // use this consistently
+  const int requiredDimension=6;
+  if (dimension!=requiredDimension) {
+    // 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? 
+  
+  return i;
+}
+
+int AliDxHFEParticleSelectionMCEl::IsSelected(AliVParticle* p, const AliVEvent* pEvent)
+{
+  /// overloaded from AliDxHFEParticleSelection: check particle
+  /// H: Have changed function. Now doing particle selection first, then run MC over 
+  /// selected tracks. Could configure it to be configurable, but not sure if it
+  /// is needed.  
+  /// Result from normal track selection is returned, result from MC is stored in
+  /// THnSparse. 
+
+  int iResult=0;
+  if (!p || !pEvent){
+    return -EINVAL;
+  }
+  // step 1:
+  // optional MC selection before the particle selection
+  if (fMCTools.MCFirst() && (iResult=CheckMC(p, pEvent))==0) {
+    // histograming?
+    return iResult;
+  }
+
+  // step 2 or 1, depending on sequence:
+  // normal particle selection
+  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
+  iResult=CheckMC(p, pEvent);
+  // TODO: why do we need to store the result in a member?
+  fResultMC=iResult;
+  return iResult;
+}
+
+int AliDxHFEParticleSelectionMCEl::CheckMC(AliVParticle* p, const AliVEvent* pEvent)
+{
+  /// check if MC criteria are fulfilled
+  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)) {
+    // rejected by pdg
+    // would also like to use this info? or not meaningful?
+    // NEED TO CHANGE THIS!
+    return 0;
+  }
+
+  int pdgMother=0;
+  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?
+    switch(pdgMother){
+    case(AliDxHFEToolsMC::kPDGpi0): fOriginMother=AliDxHFEToolsMC::kNrOrginMother; break;
+    case(AliDxHFEToolsMC::kPDGeta): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+1; break;
+    case(AliDxHFEToolsMC::kPDGgamma): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+2;break;
+    case(AliDxHFEToolsMC::kPDGJpsi): fOriginMother=AliDxHFEToolsMC::kNrOrginMother+3;break;
+    }
+  }
+
+  /*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;
+    }*/
+
+  return 1;
+}
+
+void AliDxHFEParticleSelectionMCEl::Clear(const char* option)
+{
+  /// clear internal memory
+  fMCTools.Clear(option);
+}
diff --git a/PWGHF/correlationHF/AliDxHFEParticleSelectionMCEl.h b/PWGHF/correlationHF/AliDxHFEParticleSelectionMCEl.h
new file mode 100644 (file)
index 0000000..c9780c7
--- /dev/null
@@ -0,0 +1,71 @@
+//-*- Mode: C++ -*-
+// $Id$
+
+//* This file is property of and copyright by the ALICE Project        * 
+//* ALICE Experiment at CERN, All rights reserved.                     *
+//* See cxx source for full Copyright notice                           *
+
+/// @file   AliDxHFEParticleSelectionMCEl.h
+/// @author Hege Erdal, Matthias Richter
+/// @date   2012-07-19
+/// @brief  El MC selection for D-HFE correlations
+///
+
+#ifndef ALIDXHFEPARTICLESELECTIONMCEL_H
+#define ALIDXHFEPARTICLESELECTIONMCEL_H
+
+#include "AliDxHFEParticleSelectionEl.h"
+#include "AliDxHFEToolsMC.h"
+
+class THnSparse;
+class TH1;
+
+/**
+ * @class AliDxHFEParticleSelectionMCEl
+ * Monte Carlo electron selection for D-HFE correlations, implements the specific
+ * selection criteria.
+ */
+class AliDxHFEParticleSelectionMCEl : public AliDxHFEParticleSelectionEl {
+  public:
+  /// constructor
+  AliDxHFEParticleSelectionMCEl(const char* opt="");
+  /// destructor
+  virtual ~AliDxHFEParticleSelectionMCEl();
+
+  // Setting up control objects: overloaded from AliDxHFEParticleSelection
+  virtual THnSparse* DefineTHnSparse() const;
+
+  // TODO: function can be renamed to better describe what it's doing
+  virtual int DefineParticleProperties(AliVParticle* p, Double_t* date, int dimension) const;
+
+  /// overloaded from AliDxHFEParticleSelection: check particle
+  virtual int IsSelected(AliVParticle* p, const AliVEvent *pEvent=NULL);
+
+  /// check MC criteria
+  int CheckMC(AliVParticle* p, const AliVEvent* pEvent);
+
+  /// clear internal memory
+  virtual void Clear(const char* option="");
+
+ protected:
+
+ private:
+  /// copy contructor prohibited
+  AliDxHFEParticleSelectionMCEl(const AliDxHFEParticleSelectionMCEl&);
+  /// assignment operator prohibited
+  AliDxHFEParticleSelectionMCEl& operator=(const AliDxHFEParticleSelectionMCEl&);
+
+  /// TODO: check if the label definitions can be used from the ToolsMC
+  static const char* fgkPDGMotherBinLabels[];
+  static const char* fgkTrackControlBinNames[]; //! bin labels for track control histogram
+
+  AliDxHFEToolsMC fMCTools;            // MC selection tools
+  THnSparse*      fElectronProperties; // the particle properties of selected particles
+  TH1*            fWhichCut;           // effective cut for a rejected particle
+  int fOriginMother;                   //  Holds the origin motherquark (process)
+  int fResultMC;                       // Holds information on check on MC
+
+  ClassDef(AliDxHFEParticleSelectionMCEl, 1);
+};
+
+#endif
diff --git a/PWGHF/correlationHF/AliDxHFEToolsMC.cxx b/PWGHF/correlationHF/AliDxHFEToolsMC.cxx
new file mode 100644 (file)
index 0000000..e8ec3a3
--- /dev/null
@@ -0,0 +1,399 @@
+// $Id$
+
+//**************************************************************************
+//* This file is property of and copyright by the ALICE Project            * 
+//* ALICE Experiment at CERN, All rights reserved.                         *
+//*                                                                        *
+//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
+//*                  Hege Erdal       <hege.erdal@gmail.com>               *
+//*                                                                        *
+//* 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.                  *
+//**************************************************************************
+
+/// @file   AliDxHFEToolsMC.cxx
+/// @author Hege Erdal, Matthias Richter
+/// @date   2012-07-19
+/// @brief  Common Tools for MC particle selection
+///
+
+#include "AliDxHFEToolsMC.h"
+#include "AliAODMCParticle.h"
+#include "AliVEvent.h"
+#include "AliVParticle.h"
+#include "AliLog.h"
+#include "TObjArray.h"
+#include "TString.h"
+#include "TH1D.h"
+#include <iostream>
+#include <cerrno>
+#include <memory>
+
+using namespace std;
+
+/// ROOT macro for the implementation of ROOT specific class methods
+ClassImp(AliDxHFEToolsMC)
+
+AliDxHFEToolsMC::AliDxHFEToolsMC(const char* option)
+  : fSequence(kMCLast)
+  , fMCParticles(NULL)
+  , fPDGs()
+  , fMotherPDGs()
+  , fHistPDG(NULL)
+  , fHistPDGMother(NULL)
+  , fOriginMother(kOriginNone)
+  , fMClabel(-1)
+{
+  // constructor
+  // 
+  // 
+  // 
+  // 
+  Init(option);
+}
+
+const char*  AliDxHFEToolsMC::fgkPDGBinLabels[]={
+  "positron",
+  "electron",
+  "#mu+",
+  "#mu-",
+  "#pi+",
+  "#pi-",
+  "K+",
+  "K-",
+  "proton",
+  "antiproton",
+  "others"
+};
+
+const char*  AliDxHFEToolsMC::fgkPDGMotherBinLabels[]={
+  "d",
+  "u",
+  "s",
+  "c",
+  "b",
+  "gluon",
+  "gamma",
+  "#pi^{0}",
+  "#eta",
+  "proton",
+  "others"
+};
+
+const char*  AliDxHFEToolsMC::fgkStatisticsBinLabels[]={
+  "all",   // all MC particles
+  "found", // selected particles with correct MC
+  "fake"   // selected particles without corresponding MC
+};
+
+AliDxHFEToolsMC::~AliDxHFEToolsMC()
+{
+  // destructor
+  if (fHistPDG) delete fHistPDG;
+  fHistPDG=NULL;
+  if (fHistPDGMother) delete fHistPDGMother;
+  fHistPDGMother=NULL;
+}
+
+int AliDxHFEToolsMC::Init(const char* option)
+{
+  // initialize according to options
+  TString strOption(option);
+  bool bControlHist=true;
+  std::auto_ptr<TObjArray> tokens(strOption.Tokenize(" "));
+  if (tokens.get() && tokens->GetEntriesFast()>0) {
+    for (int itoken=0; itoken<tokens->GetEntriesFast(); itoken++) {
+      if (tokens->At(itoken)==NULL) continue;
+      TString arg=tokens->At(itoken)->GetName();
+      const char* key="";
+      key="pdg=";
+      if (arg.BeginsWith(key)) {
+       arg.ReplaceAll(key, "");
+       fPDGs.push_back(arg.Atoi());
+      }
+      key="mother-pdg=";
+      if (arg.BeginsWith(key)) {
+       arg.ReplaceAll(key, "");
+       fMotherPDGs.push_back(arg.Atoi());
+      }
+      key="control-hist=";
+      if (arg.BeginsWith(key)) {
+       arg.ReplaceAll(key, "");
+       bControlHist=arg.CompareTo("off")!=0;
+      }
+      key="mc-first";
+      if (arg.BeginsWith(key)) {
+       fSequence=kMCFirst;     
+      }
+      key="mc-last";
+      if (arg.BeginsWith(key)) {
+       fSequence=kMCLast;      
+      }
+    }
+  }
+  
+  if (bControlHist) {
+    fHistPDG=CreateControlHistogram("histPDG",
+                                   "pdg code of selected particle",
+                                   sizeof(fgkPDGBinLabels)/sizeof(const char*),
+                                   fgkPDGBinLabels);
+    fHistPDGMother=CreateControlHistogram("histPDGMother",
+                                         "pdg code of first mother of selected particle",
+                                         sizeof(fgkPDGMotherBinLabels)/sizeof(fgkPDGMotherBinLabels[0]),
+                                         fgkPDGMotherBinLabels);
+  }
+  return 0;
+}
+
+int AliDxHFEToolsMC::InitMCParticles(const AliVEvent* pEvent)
+{
+
+  // init MC info from event object
+  if (!pEvent) return -EINVAL;
+
+  // TODO: choose branch name depending on VEvent type; configurable?
+  TString branchname(AliAODMCParticle::StdBranchName());
+  TObject* o=pEvent->FindListObject(branchname);
+  if (!o) {
+    AliWarningClass(Form("can not find MC info '%s' in event of type '%s'", branchname.Data(), pEvent->ClassName()));
+    return -ENOENT;
+  }
+  fMCParticles = dynamic_cast<TObjArray*>(o);
+  if (!fMCParticles) {
+    AliWarningClass(Form("ignoring MC info '%s' of wrong type '%s', expecting TObjArray", branchname.Data(), o->ClassName()));
+    return -ENODATA;
+  }
+
+  return 0;
+}
+
+bool AliDxHFEToolsMC::RejectByPDG(int pdg, const vector<int> &list) const
+{
+  // check if pdg should be rejected, particle is not rejected
+  // if it is in the list, returns always false if list is empty
+  if (list.size()==0) return false;
+  for (vector<int>::const_iterator i=list.begin();
+       i!=list.end(); i++) {
+    if (*i==pdg) return false;
+  }
+  return true;
+}
+
+bool AliDxHFEToolsMC::RejectByPDG(AliVParticle* p, bool doStatistics)
+{
+  // check if pdg should be rejected
+  // always false if not pdg list is initialized
+  if (!p) return false;
+  Int_t MClabel = p->GetLabel();
+  if(MClabel<0) return 0;
+
+  int pdgPart=-1;
+  // TODO: there might be more types of particles to be checked
+  AliAODMCParticle* aodmcp=0;
+  aodmcp=dynamic_cast<AliAODMCParticle*>(fMCParticles->At(MClabel));
+  if (aodmcp)
+    pdgPart=TMath::Abs(aodmcp->GetPdgCode());
+  if (pdgPart<0) return 0;
+
+  bool bReject=RejectByPDG(pdgPart, fPDGs);
+  if (doStatistics && fHistPDG) {
+    // TODO: think about histogramming mode, e.g. histogramming of rejected particles?
+    if (!bReject) {
+      fHistPDG->Fill(MapPDGLabel(p->PdgCode()));
+    }
+  }
+  return bReject;
+}
+
+bool AliDxHFEToolsMC::RejectByMotherPDG(AliVParticle* p, bool doStatistics)
+{
+  // check if pdg should be rejected by mother
+  // always false if not mother pdg list is initialized
+  // H: think maybe this is particle specific, and should be moved to PartSelMCEl
+  if (!p) return false;
+  int motherpdg=FindPdgOriginMother(p);
+  // TODO: This should be tweaked. Want to 
+  bool bReject=RejectByPDG(motherpdg, fMotherPDGs);
+  if (doStatistics && fHistPDGMother) {
+    // TODO: think about histogramming mode, e.g. histogramming of rejected particles?
+    if (!bReject) {
+      fHistPDGMother->Fill(MapPDGMotherLabel(p->PdgCode()));
+    }
+  }
+  return bReject;
+}
+
+int AliDxHFEToolsMC::FindMotherPDG(AliVParticle* p, bool bReturnFirstMother)
+{
+  // Will find and return pdg of the mother, either first or loop down to the
+  // initial quark
+
+  // To reset fOriginMother. 
+  fOriginMother=kOriginNone;
+
+  if (!p) return false;
+  int motherpdg=FindPdgOriginMother(p, bReturnFirstMother);
+
+  return motherpdg;
+}
+
+
+TH1* AliDxHFEToolsMC::CreateControlHistogram(const char* name,
+                                            const char* title,
+                                            int nBins,
+                                            const char** binLabels) const
+{
+  /// create control histogram
+  std::auto_ptr<TH1> h(new TH1D(name, title, nBins, -0.5, nBins-0.5));
+  if (!h.get()) return NULL;
+  for (int iLabel=0; iLabel<nBins; iLabel++) {
+    h->GetXaxis()->SetBinLabel(iLabel+1, binLabels[iLabel]);    
+  }
+  
+  return h.release();
+}
+
+int AliDxHFEToolsMC::MapPDGLabel(int pdg) const
+{
+  /// mapping of pdg code to enum
+  switch (pdg) {
+  case kPDGelectron  : return kPDGLabelElectron;
+  case -kPDGelectron : return kPDGLabelPositron;
+  case kPDGmuon             : return kPDGLabelMuPlus;
+  case -kPDGmuon     : return kPDGLabelMuMinus;
+  case kPDGpion             : return kPDGLabelPiPlus;
+  case -kPDGpion     : return kPDGLabelPiMinus;
+  case kPDGkaon             : return kPDGLabelKPlus;
+  case -kPDGkaon     : return kPDGLabelKMinus;
+  case kPDGproton    : return kPDGLabelProton;
+  case -kPDGproton   : return kPDGLabelAntiproton;
+  default:
+    return kPDGLabelOthers;
+  }
+}
+
+int AliDxHFEToolsMC::MapPDGMotherLabel(int pdg) const
+{
+  /// mapping of pdg code to enum
+  switch (pdg) {
+  case kPDGd     : return kPDGMotherLabelD;
+  case kPDGu     : return kPDGMotherLabelU;
+  case kPDGs     : return kPDGMotherLabelS;
+  case kPDGc     : return kPDGMotherLabelC;
+  case kPDGb     : return kPDGMotherLabelB;
+  case kPDGgluon : return kPDGMotherLabelGluon;
+  case kPDGgamma : return kPDGMotherLabelGamma;
+  case kPDGpi0   : return kPDGMotherLabelPi0;
+  case kPDGeta   : return kPDGMotherLabelEta;
+  case kPDGproton: return kPDGMotherLabelProton;
+  default:
+    return kPDGLabelOthers;
+  }
+}
+
+int AliDxHFEToolsMC::FindPdgOriginMother(AliVParticle* p, bool bReturnFirstMother) 
+{
+  // Return the pgd of original mother particle
+  // TODO: need also to have specific for D0, electron etc
+  // for instance to mark when you have gluon, charm or beauty
+  // among the mothers. Or maybe this will be the same anyway?
+  // TODO: implement tests on origin, if charm/beauty quark and if
+  // they came from gluon. use booleans to set this which can be accessed from 
+  // outside? Something like fSequence. 
+
+  if (!p) return kPDGnone;
+
+  Int_t imother=-1;
+  Int_t MClabel=0;
+
+  // Either using MClabel set from outside (needed for Dmesons), or find it
+  if(fMClabel<0){
+    MClabel = p->GetLabel();
+    if(MClabel<0){
+      return kPDGnone;
+    }
+  }
+  else MClabel=fMClabel;
+
+  // try different classes, unfortunately there is no common base class
+  AliAODMCParticle* aodmcp=0;
+  aodmcp=dynamic_cast<AliAODMCParticle*>(fMCParticles->At(MClabel));
+  if (!aodmcp) {
+    return kPDGnone;
+  }
+  imother = aodmcp->GetMother();
+
+  // Absolute or +/- on pgd ???
+  Int_t pdg=TMath::Abs(aodmcp->GetPdgCode());
+
+  if (imother<0){
+    // also check this particle
+    CheckOriginMother(pdg);
+    return pdg;
+  }
+
+  if (!fMCParticles->At(imother)) {
+    AliErrorClass(Form("no mc particle with label %d", imother));
+    return pdg;
+  }
+
+  AliAODMCParticle * mother=dynamic_cast<AliAODMCParticle*>(fMCParticles->At(imother));
+  if (!mother) {
+    AliErrorClass(Form("mc mother particle of wrong class type %s", fMCParticles->At(imother)->ClassName()));
+    return pdg;
+  }
+
+  if(bReturnFirstMother){
+    return mother->GetPdgCode();
+  }
+
+  CheckOriginMother(pdg);
+  
+  if (mother->GetPdgCode()==kPDGproton){
+    // why? - H: This is the proton, can't get further back
+    // To be discussed whether to do this a different way
+    return pdg;
+  }
+
+  //Reset fMClabel to find pdg of mother if looping on D
+  fMClabel=-1;
+  pdg=FindPdgOriginMother(mother);
+
+  return pdg;  
+}
+
+void AliDxHFEToolsMC::CheckOriginMother(int pdg)
+{
+
+  // Checking if the particle is a quark or gluon and setting fOriginMother accordingly.
+  // Right now only check on quark. Need also check on whether it comes from D or B meson?
+
+  switch(pdg){
+  case(kPDGc):
+    fOriginMother = kOriginCharm; break;
+  case(kPDGb): 
+    fOriginMother = kOriginBeauty; break;
+  case(kPDGgluon): 
+    if(fOriginMother==kOriginCharm) fOriginMother=kOriginGluonCharm;
+    else if(fOriginMother==kOriginBeauty) fOriginMother=kOriginGluonBeauty;
+    else fOriginMother=kOriginGluon;
+    break;
+  case(kPDGd): 
+    fOriginMother=kOriginDown; break;
+  case(kPDGu):
+    fOriginMother=kOriginUp; break;
+  case(kPDGs):
+    fOriginMother=kOriginStrange; break;
+  }
+}
+
+void AliDxHFEToolsMC::Clear(const char* /*option*/)
+{
+  // clear internal memory
+  fMCParticles=NULL;
+}
diff --git a/PWGHF/correlationHF/AliDxHFEToolsMC.h b/PWGHF/correlationHF/AliDxHFEToolsMC.h
new file mode 100644 (file)
index 0000000..5337452
--- /dev/null
@@ -0,0 +1,196 @@
+//-*- Mode: C++ -*-
+// $Id$
+
+//* This file is property of and copyright by the ALICE Project        * 
+//* ALICE Experiment at CERN, All rights reserved.                     *
+//* See cxx source for full Copyright notice                           *
+
+/// @file   AliDxHFEToolsMC.h
+/// @author Hege Erdal, Matthias Richter
+/// @date   2012-07-19
+/// @brief  Common Tools for MC particle selection
+///
+
+#ifndef ALIDXHFETOOLSMC_H
+#define ALIDXHFETOOLSMC_H
+
+#include "TObject.h"
+#include <vector>
+
+class AliVEvent;
+class AliVParticle;
+class TH1;
+
+using std::vector;
+
+/**
+ * @class AliDxHFEToolsMC
+ * Common Tools for MC particle selection.
+ */
+class AliDxHFEToolsMC {
+  public:
+  /// constructor
+  AliDxHFEToolsMC(const char* options="");
+  /// destructor
+  virtual ~AliDxHFEToolsMC();
+
+  // different pdgs
+  enum{
+    kPDGnone=0,
+    kPDGd=1,
+    kPDGu=2,
+    kPDGs=3,
+    kPDGc=4,
+    kPDGb=5,
+    kPDGelectron=11,
+    kPDGmuon=13,
+    kPDGgluon=21,
+    kPDGgamma=22,
+    kPDGpi0=111,
+    kPDGpion=211,
+    kPDGeta=221,
+    kPDGkaon=321,
+    kPDGD0=421,
+    kPDGJpsi=443,
+    kPDGproton=2212
+  };
+
+  enum {
+    kPDGLabelPositron,
+    kPDGLabelElectron,
+    kPDGLabelMuPlus,
+    kPDGLabelMuMinus,
+    kPDGLabelPiPlus,
+    kPDGLabelPiMinus,
+    kPDGLabelKPlus,
+    kPDGLabelKMinus,
+    kPDGLabelProton,
+    kPDGLabelAntiproton,
+    kPDGLabelOthers,
+    kNofPDGLabels
+  };
+
+  enum {
+    kPDGMotherLabelD,
+    kPDGMotherLabelU,
+    kPDGMotherLabelS,
+    kPDGMotherLabelC,
+    kPDGMotherLabelB,
+    kPDGMotherLabelGluon,
+    kPDGMotherLabelGamma,
+    kPDGMotherLabelPi0,
+    kPDGMotherLabelEta,
+    kPDGMotherLabelProton,
+    kPDGMotherLabelOthers,
+    kNofPDGMotherLabels
+  };
+
+  enum {
+    kMCFirst = 0,
+    kMCLast
+  };
+
+  enum {
+    kOriginNone=0,
+    kOriginDown,
+    kOriginUp,
+    kOriginStrange,
+    kOriginCharm,
+    kOriginBeauty,
+    kOriginGluon,
+    kOriginGluonCharm,
+    kOriginGluonBeauty,
+    kNrOrginMother
+  };
+
+  enum {
+    kGetOriginMother=0,
+    kGetFirstMother
+  };
+
+  /// initialize according to options
+  int Init(const char* /*option*/);
+
+  /// init MC info from event object
+  int InitMCParticles(const AliVEvent* pEvent);
+
+  /// Returning MC Array
+  TObjArray* GetMCArray() const {return fMCParticles;} 
+  
+  /// check state
+  bool IsInitialized() const {return fMCParticles!=NULL;}
+
+  /// clear internal memory
+  virtual void Clear(const char* option="");
+
+  /// flag indicating to do the selection on MC first
+  bool MCFirst() const {return fSequence==kMCFirst;}
+  /// flag indicating to do the selection on MC last
+  bool MCLast() const {return fSequence==kMCLast;}
+
+  /// Return the result on check on initial quark
+  int GetOriginMother() const {return fOriginMother;}
+
+  /// check if pdg should be rejected
+  /// always false if pdg list is not initialized
+  bool RejectByPDG(AliVParticle* p, bool doStatistics=true);
+
+  /// check if pdg should be rejected by mother
+  /// always false if mother pdg list is not initialized
+  bool RejectByMotherPDG(AliVParticle* p, bool doStatistics=true);
+
+  /// Finds pdg of first or origin mother and returns value
+  int FindMotherPDG(AliVParticle* p, bool bReturnFirstMother=false);
+
+  /// step through tree and find the original mother particle
+  /// TODO: this can possibly be const, however, a member is set inside
+  /// check whether this is really necessary
+  int FindPdgOriginMother(AliVParticle* p,bool bReturnFirstMother=false);
+
+  // Compare pdg to quark and gluon
+  void CheckOriginMother(int pdg);
+
+  // Setting MC label from outside
+  void SetMClabel(int mclab){fMClabel=mclab;}
+
+  /// TODO: want to have this function to be private again, currently
+  /// used with an external vector
+  /// check if pdg should be rejected, particle is not rejected
+  /// if it is in the list, returns always false if list is empty
+  bool RejectByPDG(int pdg, const vector<int> &list) const;
+
+ protected:
+
+ private:
+  /// copy contructor prohibited
+  AliDxHFEToolsMC(const AliDxHFEToolsMC&);
+  /// assignment operator prohibited
+  AliDxHFEToolsMC& operator=(const AliDxHFEToolsMC&);
+
+  /// create control histogram
+  TH1* CreateControlHistogram(const char* name,
+                             const char* title,
+                             int nBins,
+                             const char** binLabels) const;
+
+  /// mapping of pdg code to enum
+  int MapPDGLabel(int pdg) const;
+  /// mapping of pdg code to enum
+  int MapPDGMotherLabel(int pdg) const;
+
+  static const char* fgkPDGBinLabels[];
+  static const char* fgkPDGMotherBinLabels[];
+  static const char* fgkStatisticsBinLabels[];
+
+  int fSequence;           //  sequence of checks
+  TObjArray* fMCParticles; //! pointer to external array of MC particles
+  vector<int> fPDGs;       //  PDGs to be selected
+  vector<int> fMotherPDGs; //  mother PDGs to be selected
+  TH1* fHistPDG;           //  control histogram pdg of selected particle
+  TH1* fHistPDGMother;     //  control histogram pdg of selected particle
+  int fOriginMother;       //  Holds the origin motherquark (process)
+  int fMClabel;            //  MClabel passed from outside (default =-1)
+
+  ClassDef(AliDxHFEToolsMC, 1);
+};
+#endif