Implementation of event mixing (Hege)
authorfprino <fprino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Nov 2012 22:18:41 +0000 (22:18 +0000)
committerfprino <fprino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Nov 2012 22:18:41 +0000 (22:18 +0000)
17 files changed:
PWGHF/correlationHF/AliAnalysisTaskDxHFECorrelation.cxx
PWGHF/correlationHF/AliAnalysisTaskDxHFECorrelation.h
PWGHF/correlationHF/AliDxHFECorrelation.cxx
PWGHF/correlationHF/AliDxHFECorrelation.h
PWGHF/correlationHF/AliDxHFECorrelationMC.cxx [new file with mode: 0644]
PWGHF/correlationHF/AliDxHFECorrelationMC.h [new file with mode: 0644]
PWGHF/correlationHF/AliDxHFEParticleSelection.cxx
PWGHF/correlationHF/AliDxHFEParticleSelection.h
PWGHF/correlationHF/AliDxHFEParticleSelectionD0.cxx
PWGHF/correlationHF/AliDxHFEParticleSelectionD0.h
PWGHF/correlationHF/AliDxHFEParticleSelectionEl.cxx
PWGHF/correlationHF/AliDxHFEParticleSelectionMCD0.cxx
PWGHF/correlationHF/AliDxHFEParticleSelectionMCD0.h
PWGHF/correlationHF/AliDxHFEParticleSelectionMCEl.cxx
PWGHF/correlationHF/AliDxHFEParticleSelectionMCEl.h
PWGHF/correlationHF/AliDxHFEToolsMC.cxx
PWGHF/correlationHF/AliDxHFEToolsMC.h

index 7cedf9b..8b1237b 100644 (file)
@@ -25,6 +25,7 @@
 
 #include "AliAnalysisTaskDxHFECorrelation.h"
 #include "AliDxHFECorrelation.h"
+#include "AliDxHFECorrelationMC.h"
 #include "AliDxHFEParticleSelectionD0.h"
 #include "AliDxHFEParticleSelectionMCD0.h"
 #include "AliDxHFEParticleSelectionEl.h"
@@ -57,6 +58,8 @@
 #include "TObject.h"
 #include "TChain.h"
 #include "TSystem.h"
+#include "AliReducedParticle.h"
+#include "AliHFAssociatedTrackCuts.h" // initialization of event pool
 #include "TFile.h"
 #include <memory>
 
@@ -74,9 +77,14 @@ AliAnalysisTaskDxHFECorrelation::AliAnalysisTaskDxHFECorrelation(const char* opt
   , fElectrons(NULL)
   , fCutsD0(NULL)
   , fCutsHFE(NULL)
+  , fCuts(NULL)
   , fPID(NULL)
   , fFillOnlyD0D0bar(0)
   , fUseMC(kFALSE)
+  , fUseEventMixing(kFALSE)
+  , fSystem(0)
+  , fSelectedD0s(NULL)
+  , fSelectedElectrons(NULL)
 {
   // constructor
   //
@@ -92,6 +100,8 @@ int AliAnalysisTaskDxHFECorrelation::DefineSlots()
   DefineInput(0, TChain::Class());
   DefineOutput(1, TList::Class());
   DefineOutput(2,AliRDHFCutsD0toKpi::Class());
+  DefineOutput(3,AliHFEcuts::Class());
+  DefineOutput(4,AliHFAssociatedTrackCuts::Class());
   return 0;
 }
 
@@ -120,6 +130,10 @@ AliAnalysisTaskDxHFECorrelation::~AliAnalysisTaskDxHFECorrelation()
   fCutsHFE=NULL;
   if(fPID) delete fPID;
   fPID=NULL;
+  if(fSelectedElectrons) delete fSelectedElectrons;
+  fSelectedElectrons=NULL;
+  if(fSelectedD0s) delete fSelectedD0s;
+  fSelectedD0s=NULL;
 
 
 }
@@ -159,16 +173,23 @@ void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
   fElectrons->Init();
 
   //Correlation
-  fCorrelation=new AliDxHFECorrelation;
-  fCorrelation->SetCuts(dynamic_cast<AliRDHFCutsD0toKpi*>(fCutsD0));
+  if(fUseMC) fCorrelation=new AliDxHFECorrelationMC;
+  else fCorrelation=new AliDxHFECorrelation;
+  fCorrelation->SetCuts(fCuts);
   // TODO: check if we can get rid of the mc flag in the correlation analysis class
-  fCorrelation->SetUseMC(fUseMC);
-  fCorrelation->Init();
+  // at the moment this is needed to pass on info to AliHFCorrelator
+  TString arguments;
+  if (fUseMC)          arguments+=" use-mc";
+  if (fUseEventMixing) arguments+=" event-mixing";
+  // TODO: fSystem is a boolean right now, needs to be changed to fit also p-Pb
+  if (!fSystem)         arguments+=" system=pp";
+  else                 arguments+=" system=Pb-Pb";
+  fCorrelation->Init(arguments);
 
   // Fix for merging:
   // Retrieving the individual objects created
   // and storing them instead of fD0s, fElectrons etc.. 
-  TList *list =(TList*)fD0s->GetControlObjects();
+  TList *list =(TList*)fCorrelation->GetControlObjects();
   TObject *obj=NULL;
 
   TIter next(list);
@@ -176,7 +197,7 @@ void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
     fOutput->Add(obj);
   }
 
-  list=(TList*)fCorrelation->GetControlObjects();
+  list=(TList*)fD0s->GetControlObjects();
   next=TIter(list);
   while((obj= next())){
     fOutput->Add(obj);
@@ -187,20 +208,26 @@ void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
   while((obj = next()))
     fOutput->Add(obj);
 
-  PostData(1, fOutput);
   if (fCutsD0) {
+    // TODO: eliminate this copy
     AliRDHFCutsD0toKpi* cuts=dynamic_cast<AliRDHFCutsD0toKpi*>(fCutsD0);
     if (!cuts) {
       AliFatal(Form("cut object %s is of incorrect type %s, expecting AliRDHFCutsD0toKpi", fCutsD0->GetName(), fCutsD0->ClassName()));
       return;
     }
-    // TODO: why copy? cleanup?
-    AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*cuts);
-    const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
-    copyfCuts->SetName(nameoutput);
-
-    PostData(2,copyfCuts);
   }
+
+  // TODO: why copy? cleanup?
+  AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(dynamic_cast<AliRDHFCutsD0toKpi&>(*fCutsD0));
+  const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
+  copyfCuts->SetName(nameoutput);
+
+  // all tasks must post data once for all outputs
+  PostData(1, fOutput);
+  PostData(2,copyfCuts);
+  PostData(3,fCutsHFE);
+  PostData(4,fCuts);
+
 }
 
 void AliAnalysisTaskDxHFECorrelation::UserExec(Option_t* /*option*/)
@@ -214,6 +241,8 @@ void AliAnalysisTaskDxHFECorrelation::UserExec(Option_t* /*option*/)
   AliVEvent *pEvent = dynamic_cast<AliVEvent*>(pInput);
   TClonesArray *inputArray=0;
 
+  fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsAll);
+
   if(!pEvent && AODEvent() && IsStandardAOD()) { //Not sure if this is needed.. Keep it for now. 
     // In case there is an AOD handler writing a standard AOD, use the AOD 
     // event in memory rather than the input (ESD) event.    
@@ -242,7 +271,6 @@ void AliAnalysisTaskDxHFECorrelation::UserExec(Option_t* /*option*/)
     return;
   }
 
-  fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsAll);
   AliRDHFCuts* cutsd0=dynamic_cast<AliRDHFCuts*>(fCutsD0);
   if (!cutsd0) return; // Fatal thrown already in initialization
 
@@ -265,64 +293,93 @@ void AliAnalysisTaskDxHFECorrelation::UserExec(Option_t* /*option*/)
  
   fPID->SetPIDResponse(pidResponse);
 
+  // Retrieving process from the AODMCHeader. 
+  // TODO: Move it somewhere else? (keep it here for the moment since only need to read once pr event)
+  if(fUseMC){
+    AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(pEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+    
+    if (!mcHeader) {
+      AliError("Could not find MC Header in AOD");
+      return;
+    }
+    Int_t eventType = mcHeader->GetEventType();
+    fCorrelation->SetEventType(eventType);
+  }
 
   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 
   fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsSel);
 
-  //TODO: pSelectedD0s will now contain both D0 and D0bar. Probably add argument to Select
-  // which determines which ones to return. 
-  std::auto_ptr<TObjArray> pSelectedD0s(fD0s->Select(inputArray,pEvent));
-  if(! pSelectedD0s.get()) {
+  if(fSelectedD0s) delete fSelectedD0s;
+  fSelectedD0s=(fD0s->Select(inputArray,pEvent));
+  
+  if(! fSelectedD0s) {
     return;
   }
-  Int_t nD0Selected = pSelectedD0s->GetEntriesFast();
+  Int_t nD0Selected = fSelectedD0s->GetEntriesFast();
+
 
-  // No need to go further if no D0s are found. Not sure if this is the best way though..
-  // At the moment this means no control histos can be implemented in AliDxHFECorrelation
-  if(nD0Selected==0){
-    //AliInfo("No D0s found in this event");
+  /*std::auto_ptr<TObjArray> pSelectedD0s(fD0s->Select(inputArray,pEvent));
+  if(! pSelectedD0s.get()) {
     return;
   }
-  fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsD0);
+  Int_t nD0Selected = pSelectedD0s->GetEntriesFast();*/
 
+  // When not using EventMixing, no need to go further if no D0s are found.
+  // For Event Mixing, need to store all found electrons in the pool
+  if(!fUseEventMixing && nD0Selected==0){
+    AliDebug(4,"No D0s found in this event");
+    return;
+  }
 
-  std::auto_ptr<TObjArray> pSelectedElectrons(fElectrons->Select(pEvent));
+  fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsD0);
 
-  // TODO: use the array of selected track for something, right now
-  // only the control histograms of the selection class are filled
+  /* std::auto_ptr<TObjArray> pSelectedElectrons(fElectrons->Select(pEvent));
   // note: the pointer is deleted automatically once the scope is left
   // if the array should be published, the auto pointer must be released
   // first, however some other cleanup will be necessary in that case
   // probably a clone with a reduced AliVParticle implementation is
   // appropriate.
 
-
   if(! pSelectedElectrons.get()) {
     return;
   }
 
-  Int_t nElSelected =  pSelectedElectrons->GetEntriesFast();
-
-  // No need to go further if no electrons are found. Not sure if this is the best way though..
-  // At the moment this means no control histos can be implemented in AliDxHFECorrelation
-  if(nElSelected==0){
-    //AliInfo("No electrons found in this event");
+  Int_t nElSelected =  pSelectedElectrons->GetEntriesFast();*/
+  if (fSelectedElectrons) delete fSelectedElectrons;
+  fSelectedElectrons=(fElectrons->Select(pEvent));
+  
+  if(! fSelectedElectrons) {
     return;
   }
 
+  Int_t nElSelected =  fSelectedElectrons->GetEntriesFast();
+
+
+  // No need to go further if no electrons are found, except for event mixing. Will here anyway correlate D0s with electrons from previous events
+  if(!fUseEventMixing && nElSelected==0){
+    AliDebug(4,"No electrons found in this event");
+    return;
+  }
+  if(nD0Selected==0 && nElSelected==0){
+    AliDebug(4,"Neither D0 nor electrons in this event");
+    return;
+  }
+  
   AliDebug(4,Form("Number of D0->Kpi Start: %d , End: %d    Electrons Selected: %d\n", nInD0toKpi, nD0Selected, nElSelected));
 
   fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsD0e);
 
-  //This is only called if there are electrons and D0s present.
-  int iResult=fCorrelation->Fill(pSelectedD0s.get(), pSelectedElectrons.get());
+  //int iResult=fCorrelation->Fill(pSelectedD0s.get(), pSelectedElectrons.get(), pEvent);
+  int iResult=fCorrelation->Fill(fSelectedD0s, fSelectedElectrons, pEvent);
 
   if (iResult<0) {
     AliError(Form("%s processing failed with error %d", fCorrelation->GetName(), iResult));
   }
 
   PostData(1, fOutput);
+  return;
+
 }
 
 void AliAnalysisTaskDxHFECorrelation::FinishTaskOutput()
index f4ebdd3..adbd9f3 100644 (file)
@@ -27,6 +27,7 @@ class AliDxHFECorrelation;
 class AliAnalysisCuts;
 class AliHFEpid;
 class AliHFEcuts;
+class AliHFAssociatedTrackCuts;
 
 /**
  * @class AliAnalysisTaskDxHFECorrelation
@@ -59,6 +60,9 @@ class AliAnalysisTaskDxHFECorrelation : public AliAnalysisTaskSE {
   virtual void SetUseMC(Bool_t useMC){fUseMC=useMC;}
   virtual void SetCutsD0(AliAnalysisCuts* cuts){fCutsD0=cuts;}
   virtual void SetCutsHFE(AliHFEcuts* cuts){fCutsHFE=cuts;}
+  void SetCuts(AliAnalysisCuts* cuts){fCuts=cuts;}
+  void SetUseEventMixing(Bool_t useMixing) {fUseEventMixing=useMixing;}
+  void SetSystem(Bool_t system){fSystem=system;}
 
   /// overloaded from TObject: get option
   virtual Option_t* GetOption() const { return fOption;}
@@ -82,12 +86,17 @@ class AliAnalysisTaskDxHFECorrelation : public AliAnalysisTaskSE {
   AliDxHFEParticleSelection* fElectrons; //  selection of electrons
   AliAnalysisCuts *fCutsD0;              //  Cuts D0 
   AliHFEcuts *fCutsHFE;                  //  Cuts HFE
+  AliAnalysisCuts *fCuts;                // Cuts which holds info for AliHFCorrelator 
   AliHFEpid *fPID;                       //  dummy
   Int_t     fFillOnlyD0D0bar;            // flag to set what to fill (0 = both, 1 = D0 only, 2 = D0bar only)
   Bool_t fUseMC;                 // use MC info
+  Bool_t fUseEventMixing;        // Run Event Mixing analysis
+  Int_t fSystem;                 // Which system pp/PbPb
+  TObjArray *fSelectedD0s; // Array for selected D0s
+  TObjArray *fSelectedElectrons; // Array for selected Electrons
 
 
-  ClassDef(AliAnalysisTaskDxHFECorrelation, 2);
+  ClassDef(AliAnalysisTaskDxHFECorrelation, 3);
 };
 
 #endif
index 63ec94f..f9d47b6 100644 (file)
 //#include "AliFlowTrackSimple.h"      // required dependency libPWGflowBase.so
 //#include "AliFlowCandidateTrack.h"   // required dependency libPWGflowTasks.so
 //#include "AliCFContainer.h"          // required dependency libCORRFW.so
-#include "AliAODRecoDecayHF2Prong.h"   // libPWGHFvertexingHF
-#include "AliRDHFCutsD0toKpi.h"
 #include "TObjArray.h"
+#include "AliHFCorrelator.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
 #include "TH1D.h"
 #include "TH2D.h"
+#include "TH3D.h"
 #include "THnSparse.h"
 #include "TMath.h"
 #include "TFile.h"
 #include "TCanvas.h"
 #include "TDatabasePDG.h"
 #include "TLorentzVector.h"
+#include "AliReducedParticle.h"
+#include "AliDxHFEParticleSelection.h"
 #include <iostream>
 #include <cerrno>
 #include <memory>
@@ -57,6 +61,16 @@ AliDxHFECorrelation::AliDxHFECorrelation(const char* name)
   , fhEventControlCorr(NULL)
   , fCuts(NULL)
   , fUseMC(kFALSE)
+  , fCorrelator(NULL)
+  , fUseEventMixing(kFALSE)
+  , fSystem(0)
+  , fMinPhi(-TMath::Pi()/2)
+  , fMaxPhi(3*TMath::Pi()/2)
+  , fDeltaPhi(0)
+  , fDeltaEta(0)
+  , fDimThn(0)
+  , fCorrArray(NULL)
+  , fEventType(0)
 {
   // default constructor
   // 
@@ -71,16 +85,6 @@ const char* AliDxHFECorrelation::fgkEventControlBinNames[]={
   "nEventsD0e"
 };
 
-// TODO: maybe delete PtD0 in the future, note: there are more places!
-const char* AliDxHFECorrelation::fgkCorrControlBinNames[]={
-  "D0InvMass",
-  "PtD0",
-  "PhiD0",
-  "PtBinD0",
-  "dPhi",
-  "Pt electron"
-};
-
 AliDxHFECorrelation::~AliDxHFECorrelation()
 {
   // destructor
@@ -95,63 +99,27 @@ AliDxHFECorrelation::~AliDxHFECorrelation()
   fControlObjects=NULL;
   fCorrProperties=NULL;
   fhEventControlCorr=NULL;
+  if(fCorrelator) delete fCorrelator;
+  fCorrelator=NULL;
+  if(fCorrArray) delete fCorrArray;
+  fCorrArray=NULL;
 
   // NOTE: the external object is deleted elsewhere
   fCuts=NULL;
 }
 
-int AliDxHFECorrelation::Init()
+int AliDxHFECorrelation::Init(const char* arguments)
 {
-  AliInfo("Setting up THnSparse");
-  // Are using THnSparse instead of histograms to store information
-  // At the moment use old setup for THnSparse, change to new later
-
-  /*
-    TString name;
-    const int thnSize = 7;
-    const double pi=TMath::Pi();
-    //                                  0      1     2      3      4      5    6
-    //                               mass      Pt   Phi    ePt    ePhi  eEta  DeltaPhi
-    int    thnBins[thnSize] = {   200,   1000,  100,   100,   100,  100,  180  };
-    double thnMin [thnSize] = {  1.5648,   0,    0,     0,    0.0, -1.0,  0.0  };
-    double thnMax [thnSize] = {  2.1648, 100,  (2*pi), 100, (2*pi), 1.0, (2*pi)};
-
-    name.Form("%s info", GetName());
-    std::auto_ptr<THnSparseF> CorrProperties(new THnSparseF(name, name, thnSize, thnBins, thnMin, thnMax));
-
-    if (CorrProperties.get()==NULL) {
-    return -ENOMEM;
-    }
-    int axis=0;
-    CorrProperties->GetAxis(axis++)->SetTitle("D0 Inv Mass");
-    CorrProperties->GetAxis(axis++)->SetTitle("D0 Pt");
-    CorrProperties->GetAxis(axis++)->SetTitle("D0 Phi"); 
-    CorrProperties->GetAxis(axis++)->SetTitle("electron Pt"); 
-    CorrProperties->GetAxis(axis++)->SetTitle("electron Phi"); 
-    CorrProperties->GetAxis(axis++)->SetTitle("electron Eta"); 
-    CorrProperties->GetAxis(axis++)->SetTitle("#Delta#Phi D0 -HFE"); */
-
-
-  // TODO: think about removing PtD0, but remember to change this in all places!
-  static const int sizeEventdphi = 6;  
-  Double_t minPhi= -TMath::Pi()/2;
-  Double_t maxPhi= 3*TMath::Pi()/2;
-  //                                       0         1      2      3       4       5
-  //                                     D0invmass PtD0   PhiD0  PtbinD0  dphi   Pte
-  int    binsEventdphi[sizeEventdphi] = {   200,    1000,    100,     21,  100,    1000};
-  double minEventdphi [sizeEventdphi] = { 1.5648,    0,      0,      0,    minPhi,   0 };
-  double maxEventdphi [sizeEventdphi] = {  2.1648,  100, 2*(TMath::Pi()), 20,  maxPhi, 100 };
-
-  TString name;
-  name.Form("%s info", GetName());
-  std::auto_ptr<THnSparseF> CorrProperties(new THnSparseF(name, name, sizeEventdphi, binsEventdphi,minEventdphi , maxEventdphi));
-  if (CorrProperties.get()==NULL) {
-    return -ENOMEM;
-  }
-  int iLabel=0;
+  //
+  // Will initialize thnsparse, histogram and AliHFCorrelator
+  //
+  AliInfo("Initializing correlation objects");
+  ParseArguments(arguments);
 
-  for (iLabel=0; iLabel<sizeEventdphi; iLabel++)
-    CorrProperties->GetAxis(iLabel)->SetTitle(fgkCorrControlBinNames[iLabel]);
+  //----------------------------------------------
+  // Setting up THnSparse 
+  fCorrProperties=DefineTHnSparse();
+  AddControlObject(fCorrProperties);
 
   //----------------------------------------------
   // Histogram for storing event information
@@ -160,18 +128,192 @@ int AliDxHFECorrelation::Init()
   if (!hEventControl.get()) {
     return -ENOMEM;
   }
-
+  int iLabel=0;
   for (iLabel=0; iLabel<kNEventControlLabels; iLabel++)
     hEventControl->GetXaxis()->SetBinLabel(iLabel+1, fgkEventControlBinNames[iLabel]);
 
-  fCorrProperties=CorrProperties.release();
-  AddControlObject(fCorrProperties);
   fhEventControlCorr=hEventControl.release();
   AddControlObject(fhEventControlCorr);
 
+  //----------------------------------------------
+  // AliHFCorrelator for Event Mixing and correlation
+  // 
+  // fCuts is the hadron cut object, fSystem to switch between pp or PbPb
+  AliHFAssociatedTrackCuts* cuts=dynamic_cast<AliHFAssociatedTrackCuts*>(fCuts);
+  if (!cuts) {
+    if (fCuts)
+      AliError(Form("cuts object of wrong type %s, required AliHFAssociatedTrackCuts", fCuts->ClassName()));
+    else
+      AliError("mandatory cuts object missing");
+    return -EINVAL;
+  }
+  fCorrelator = new AliHFCorrelator("Correlator", cuts, fSystem); 
+  fCorrelator->SetDeltaPhiInterval(fMinPhi,fMaxPhi); //Correct Phi Interval
+  fCorrelator->SetEventMixing(fUseEventMixing);      // mixing Off/On 
+  fCorrelator->SetAssociatedParticleType(AliHFCorrelator::kElectron);
+  // 0: don't calculate d0; 1: return d0; 2: return d0/d0err
+  fCorrelator->SetApplyDisplacementCut(kFALSE); 
+  fCorrelator->SetUseMC(fUseMC);
+  fCorrelator->SetUseReco(kTRUE); // Reco/MCTruth
+  Bool_t pooldef = fCorrelator->DefineEventPool();
+       
+  if(!pooldef) AliInfo("Warning:: Event pool not defined properly");
+
+
+    // ============================= EVENT MIXING CHECKS ======================================
+  // TODO: Not sure if all 4 histos are needed. Keep for now.. 
+  Int_t MaxNofEvents = cuts->GetMaxNEventsInPool();
+  Int_t MinNofTracks = cuts->GetMinNTracksInPool();
+  Int_t NofCentBins = cuts->GetNCentPoolBins();
+  Double_t * CentBins = cuts->GetCentPoolBins();
+  Int_t NofZVrtxBins = cuts->GetNZvtxPoolBins();
+  Double_t *ZVrtxBins = cuts->GetZvtxPoolBins();
+       
+  Int_t k =0;
+       
+  if(fSystem) k = 100; // PbPb centrality
+  if(!fSystem) k = NofCentBins; // pp multiplicity
+       
+       
+  Double_t minvalue = CentBins[0];
+  Double_t maxvalue = CentBins[NofCentBins+1];
+  Double_t Zminvalue = ZVrtxBins[0];
+  Double_t Zmaxvalue = ZVrtxBins[NofCentBins+1];
+
+  Double_t Nevents[]={0,2*MaxNofEvents/10,4*MaxNofEvents/10,6*MaxNofEvents/10,8*MaxNofEvents/10,MaxNofEvents};
+  Double_t * events = Nevents;
+       
+  TH3D * EventsPerPoolBin = new TH3D("EventsPerPoolBin","Number of events in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,5,events);
+  EventsPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");
+  EventsPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");
+  EventsPerPoolBin->GetZaxis()->SetTitle("Number of events in pool bin");
+  if(fUseEventMixing) AddControlObject(EventsPerPoolBin);
+       
+  Int_t MaxNofTracks = (MaxNofEvents+1)*MinNofTracks;
+  Int_t Diff = MaxNofTracks-MinNofTracks;
+       
+  Double_t Ntracks[]={MinNofTracks,MinNofTracks+Diff/5,MinNofTracks+2*Diff/5,MinNofTracks+3*Diff/5,MinNofTracks+4*Diff/5,MaxNofTracks};
+  Double_t  * trackN = Ntracks;
+       
+  TH3D * NofTracksPerPoolBin = new TH3D("NofTracksPerPoolBin","Number of tracks in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins,5,trackN);
+  NofTracksPerPoolBin->GetXaxis()->SetTitle("Centrality/multiplicity ");
+  NofTracksPerPoolBin->GetYaxis()->SetTitle("Z vertex [cm]");
+  NofTracksPerPoolBin->GetZaxis()->SetTitle("Number of tracks per bin");
+       
+  if(fUseEventMixing) AddControlObject(NofTracksPerPoolBin);
+       
+  TH2D * NofPoolBinCalls = new TH2D("NofPoolBinCalls","Number of tracks in bin pool",NofCentBins,CentBins,NofZVrtxBins,ZVrtxBins);
+  NofPoolBinCalls->GetXaxis()->SetTitle("Centrality/multiplicity ");
+  NofPoolBinCalls->GetYaxis()->SetTitle("Z vertex [cm]");
+  if(fUseEventMixing) AddControlObject(NofPoolBinCalls);
+
+  TH2D * EventProps = new TH2D("EventProps","Number of tracks in bin pool",k,minvalue,maxvalue,100,Zminvalue,Zmaxvalue);
+  EventProps->GetXaxis()->SetTitle("Centrality/multiplicity ");
+  EventProps->GetYaxis()->SetTitle("Z vertex [cm]");
+  if(fUseEventMixing) AddControlObject(EventProps);
+
+
   return 0;
 }
 
+int AliDxHFECorrelation::ParseArguments(const char* arguments)
+{
+  // parse arguments and set internal flags
+  TString strArguments(arguments);
+  auto_ptr<TObjArray> tokens(strArguments.Tokenize(" "));
+  if (!tokens.get()) return -ENOMEM;
+
+  TIter next(tokens.get());
+  TObject* token;
+  while ((token=next())) {
+    TString argument=token->GetName();
+   
+    if (argument.BeginsWith("event-mixing")) {
+      fUseEventMixing=true;
+      continue;
+    }
+      
+    if (argument.BeginsWith("use-mc")) {
+      fUseMC=true;
+      continue;
+    }
+    if (argument.BeginsWith("system=")) {
+      argument.ReplaceAll("system=", "");
+      if (argument.CompareTo("pp")==0) fSystem=0;
+      else if (argument.CompareTo("Pb-Pb")==0) fSystem=1;
+      else {
+       AliWarning(Form("can not set collision system, unknown parameter '%s'", argument.Data()));
+       // TODO: check what makes sense
+       fSystem=0;
+      }
+      continue;
+    }
+    AliWarning(Form("unknown argument '%s'", argument.Data()));
+      
+  }
+
+  return 0;
+}
+
+THnSparse* AliDxHFECorrelation::DefineTHnSparse()
+{
+  //
+  //Defines the THnSparse. For now, only calls CreateControlTHnSparse
+
+  // here is the only place to change the dimension
+  static const int sizeEventdphi = 7;  
+  InitTHnSparseArray(sizeEventdphi);
+  const double pi=TMath::Pi();
+
+  //TODO: add phi for electron??
+  //                                           0           1       2      3         4     5      6   
+  //                                         D0invmass   PtD0    PhiD0  PtbinD0    Pte   dphi   deta   
+  int         binsEventdphi[sizeEventdphi] = {   200,      1000,   100,    21,     1000,   100,    100};
+  double      minEventdphi [sizeEventdphi] = { 1.5648,      0,       0,     0,       0,  fMinPhi, -2};
+  double      maxEventdphi [sizeEventdphi] = { 2.1648,     100,   2*pi,    20,      100, fMaxPhi, 2};
+  const char* nameEventdphi[sizeEventdphi] = {
+    "D0InvMass",
+    "PtD0",
+    "PhiD0",
+    "PtBinD0",
+    "PtEl",
+    "#Delta#Phi", 
+    "#Delta#eta"
+  };
+
+  TString name;
+  name.Form("%s info", GetName());
+
+
+  return CreateControlTHnSparse(name,sizeEventdphi,binsEventdphi,minEventdphi,maxEventdphi,nameEventdphi);
+
+}
+
+THnSparse* AliDxHFECorrelation::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<THnSparseD> th(new THnSparseD(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();
+
+}
+
 int AliDxHFECorrelation::AddControlObject(TObject* pObj)
 {
   AliInfo("Adding object");
@@ -199,9 +341,11 @@ int AliDxHFECorrelation::HistogramEventProperties(int bin)
   return 0;
 }
 
-int AliDxHFECorrelation::Fill(const TObjArray* triggerCandidates, const TObjArray* associatedTracks)
+int AliDxHFECorrelation::Fill(const TObjArray* triggerCandidates, const TObjArray* associatedTracks, const AliVEvent* pEvent)
 {
-  /// fill ThnSparse from array of AliVParticle objects
+  //
+  // will use AliHFCorrelator to process D0-electron pair and then fill THnSparse.
+  //
   if (!triggerCandidates || !associatedTracks) return -EINVAL;
   if (!fControlObjects) {
     Init();
@@ -209,57 +353,110 @@ int AliDxHFECorrelation::Fill(const TObjArray* triggerCandidates, const TObjArra
   if (!fControlObjects) {
     AliError("Initialisation failed, can not fill THnSparse");
   }
-
-  const double Pii=TMath::Pi();
+  // set the event to be processed
+  // TODO: change the correlator class to take the const pointer
+  fCorrelator->SetAODEvent(dynamic_cast<AliAODEvent*>(const_cast<AliVEvent*>(pEvent))); 
+
+  Bool_t correlatorON = fCorrelator->Initialize(); //define the pool for mixing
+  if(!correlatorON) {
+    AliInfo("AliHFCorrelator didn't initialize the pool correctly or processed a bad event");
+    return 1;
+  }
 
   TIter itrigger(triggerCandidates);
   TObject* otrigger=NULL;
   int ctrigger=-1;
 
   // For the moment this is very specific to D0-electron correlation. Should be 
-  // changed to be more specific. 
+  // changed to be less specific. 
   while ((otrigger=itrigger())!=NULL) {
     // loop over trigger D0 particle
     ctrigger++;
-    AliAODRecoDecayHF2Prong *d0 = dynamic_cast<AliAODRecoDecayHF2Prong*>(otrigger);
-    if (!d0) continue;
-
-    TIter iElectron(associatedTracks);
-    TObject* oElectron=NULL;
-    int cElectron=-1;
-    while ((oElectron=iElectron())!=NULL) {
-      // loop over electrons
-
-      cElectron++;
-      AliVParticle* pElectron=dynamic_cast<AliVParticle*>(oElectron);
-      if (!pElectron) continue;
-
-      //Calculating dPhi using TLorentzVectors
-      Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
-      TLorentzVector D0vector(0.,0.,0.,0.);
-      TLorentzVector evector(0.,0.,0.,0.);
-      D0vector.SetXYZM(d0->Px(),d0->Py(),d0->Pz(),mPDG);  
-      evector.SetXYZM(pElectron->Px(),pElectron->Py(),pElectron->Pz(),0.000511);
-      Double_t DeltaPhi=D0vector.DeltaPhi(evector);
-      if(DeltaPhi<-TMath::PiOver2()) DeltaPhi=DeltaPhi+(2*Pii);
-
-      /*Double_t CorrelationArray[]={d0->InvMassD0(),d0->Pt(),d0->Phi(),
-       pElectron->Pt(),pElectron->Phi(), pElectron->Eta(),
-       DeltaPhi};*/
-
-      // TODO: think about a method to retrieve the pt bin from the
-      // selection object
-      Int_t ptbin=fCuts->PtBin(d0->Pt());
-      Double_t CorrelationArray[]={d0->InvMassD0(),d0->Pt(),d0->Phi(),ptbin,
-                                  DeltaPhi, pElectron->Pt() };
-      fCorrProperties->Fill(CorrelationArray);
-
-    } // loop over associated tracks
+    AliReducedParticle* ptrigger=dynamic_cast<AliReducedParticle*>(otrigger);
+    if (!ptrigger)  continue;
+
+    Double_t phiTrigger = ptrigger->Phi();
+    Double_t ptTrigger = ptrigger->Pt();
+    Double_t etaTrigger = ptrigger->Eta();
+
+    // set the phi of the D meson in the correct range
+    // TODO: Is this correct to do this??
+    phiTrigger = fCorrelator->SetCorrectPhiRange(phiTrigger);
+    // pass to the object the necessary trigger part parameters
+    fCorrelator->SetTriggerParticleProperties(ptTrigger,phiTrigger,etaTrigger); 
+
+    Bool_t execPool = fCorrelator->ProcessEventPool();
+    if(fUseEventMixing && !execPool) {
+      AliInfo("Mixed event analysis: pool is not ready");
+      continue;
+    }
+    Int_t NofEventsinPool = 1;
+    if(fUseEventMixing) NofEventsinPool = fCorrelator->GetNofEventsInPool();
+               
+    // loop on events in the pool; if it is SE analysis, stops at one
+    for (Int_t jMix =0; jMix < NofEventsinPool; jMix++){
+      Bool_t analyzetracks = fCorrelator->ProcessAssociatedTracks(jMix, associatedTracks);
+                       
+      if(!analyzetracks) {
+       AliInfo("AliHFCorrelator::Cannot process the track array");
+       continue;
+      }
+
+      Int_t NofTracks = fCorrelator->GetNofTracks();
+
+      // looping on track candidates
+      for(Int_t iTrack = 0; iTrack<NofTracks; iTrack++){ 
+       Bool_t runcorrelation = fCorrelator->Correlate(iTrack);
+       if(!runcorrelation) continue;
+                       
+       fDeltaPhi = fCorrelator->GetDeltaPhi();
+       fDeltaEta = fCorrelator->GetDeltaEta();
+       
+       AliReducedParticle *assoc = fCorrelator->GetAssociatedParticle();
+
+       FillParticleProperties(ptrigger,assoc,ParticleProperties(),GetDimTHnSparse());
+       fCorrProperties->Fill(ParticleProperties());
+
+      } // loop over electron tracks in event
+    } // loop over events in pool
   } // loop over trigger particle
 
+  Bool_t updated = fCorrelator->PoolUpdate(associatedTracks);
+  if(fUseEventMixing){
+    if(!updated) AliInfo("Pool was not updated");
+    else {
+      EventMixingChecks(pEvent);
+      AliInfo("Pool was updated");
+    }
+  }
   return 0;
 }
 
+
+
+int AliDxHFECorrelation::FillParticleProperties(AliVParticle* tr, AliVParticle *as, Double_t* data, int dimension) const
+{
+  // fill the data array from the particle data
+  if (!data) return -EINVAL;
+  AliReducedParticle *ptrigger=(AliReducedParticle*)tr;
+  AliReducedParticle *assoc=(AliReducedParticle*)as;
+  if (!ptrigger || !assoc) return -ENODATA;
+  int i=0;
+  if (dimension!=GetDimTHnSparse()) {
+    // TODO: think about filling only the available data and throwing a warning
+    return -ENOSPC;
+  }
+  data[i++]=ptrigger->GetInvMass();
+  data[i++]=ptrigger->Pt();
+  data[i++]=ptrigger->Phi();
+  data[i++]=ptrigger->GetPtBin(); 
+  data[i++]=assoc->Pt();
+  data[i++]=GetDeltaPhi();
+  data[i++]=GetDeltaEta();
+
+  return i;
+}
+
 void AliDxHFECorrelation::Clear(Option_t * /*option*/)
 {
   /// overloaded from TObject: cleanup
@@ -338,3 +535,35 @@ AliDxHFECorrelation& AliDxHFECorrelation::operator+=(const AliDxHFECorrelation&
   }
   return *this;
 }
+
+
+//____________________________  Run checks on event mixing ___________________________________________________
+void AliDxHFECorrelation::EventMixingChecks(const AliVEvent* pEvent){
+       
+  AliAODEvent *AOD= (AliAODEvent*)(pEvent);
+  AliCentrality *centralityObj = 0;
+  Int_t multiplicity = -1;
+  Double_t MultipOrCent = -1;
+       
+  // get the pool for event mixing
+  if(!fSystem){ // pp
+    multiplicity = AOD->GetNTracks();
+    MultipOrCent = multiplicity; // convert from Int_t to Double_t
+  }
+  if(fSystem){ // PbPb         
+    centralityObj = AOD->GetHeader()->GetCentralityP();
+    MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
+    AliInfo(Form("Centrality is %f", MultipOrCent));
+  }
+       
+  AliAODVertex *vtx = AOD->GetPrimaryVertex();
+  Double_t zvertex = vtx->GetZ(); // zvertex
+
+  AliEventPool *pool = fCorrelator->GetPool();
+
+  ((TH2D*)fControlObjects->FindObject("NofPoolBinCalls"))->Fill(MultipOrCent,zvertex); // number of calls of pool
+  ((TH2D*)fControlObjects->FindObject("EventProps"))->Fill(MultipOrCent,zvertex); // event properties
+  ((TH3D*)fControlObjects->FindObject("EventsPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->NTracksInPool()); // number of events in the pool
+  ((TH3D*)fControlObjects->FindObject("NofTracksPerPoolBin"))->Fill(MultipOrCent,zvertex,pool->GetCurrentNEvents()); // number of calls of pool
+}
+       
index 80dd3d2..85cb0a3 100644 (file)
 
 #include "TNamed.h"
 
-class AliRDHFCutsD0toKpi;
 class TH1;
 class THnSparse;
 class TObject;
 class TList;
+class AliHFCorrelator;
+class AliVParticle;
+class TObjArray;
+class AliVEvent;
+class AliAnalysisCuts;
 
 class AliDxHFECorrelation : public TNamed {
  public:
@@ -39,13 +43,26 @@ class AliDxHFECorrelation : public TNamed {
   };
 
   // init
-  int Init();
+  int Init(const char* arguments="");
+
+  // parse argument string
+  int ParseArguments(const char* arguments);
 
   /// fill histograms from particles
-  int Fill(const TObjArray* candidatesD0, const TObjArray* candidatesElectron);
+  int Fill(const TObjArray* candidatesD0, const TObjArray* candidatesElectron, const AliVEvent* pEvent);
 
   /// histogram event properties
   virtual int HistogramEventProperties(int bin);
+  virtual THnSparse* DefineTHnSparse();
+  virtual int FillParticleProperties(AliVParticle* tr, AliVParticle* as, Double_t* data, int dimension) const;
+
+  /// create control THnSparse
+  THnSparse* CreateControlTHnSparse(const char* name,
+                                   int thnSize,
+                                   int* thnBins,
+                                   double* thnMin,
+                                   double* thnMax,
+                                   const char** binLabels) const;
 
   /// overloaded from TObject: cleanup
   virtual void Clear(Option_t * option ="");
@@ -60,16 +77,26 @@ class AliDxHFECorrelation : public TNamed {
   /// overloaded from TObject: save to file
   virtual void     SaveAs(const char *filename="",Option_t *option="") const; // *MENU*
 
-  virtual void SetCuts(AliRDHFCutsD0toKpi* cuts) {fCuts=cuts;}
+  virtual void SetCuts(AliAnalysisCuts* cuts) {fCuts=cuts;}
   virtual void SetUseMC(Bool_t useMC){fUseMC=useMC;}
+  //void SetUseEventMixing(Bool_t useMixing) {fUseEventMixing=useMixing;}
+  //void SetSystem(Bool_t system){fSystem=system;}
+  //void SetPhiRange(Double_t min, Double_t max){fMinPhi=min; fMaxPhi=max;}
+  // TODO: SetEventType only needed for MC. How to avoid this?
+  virtual void SetEventType(int type){fEventType=type;}
 
   Bool_t GetUseMC() const {return fUseMC;}
   const TList* GetControlObjects() const {return fControlObjects;}
+  Double_t GetMinPhi() const {return fMinPhi;}
+  Double_t GetMaxPhi() const {return fMaxPhi;}
+  Double_t GetDeltaPhi() const {return fDeltaPhi;}
+  Double_t GetDeltaEta() const {return fDeltaEta;}
+  inline int GetDimTHnSparse() const {return fDimThn;}
 
+  void EventMixingChecks(const AliVEvent* pEvent);
 
   AliDxHFECorrelation& operator+=(const AliDxHFECorrelation& other);
 
-
   // Probably not needed anymore, since code was changed to THnSparse
   // but keep here in case we need it later
   enum {
@@ -86,6 +113,15 @@ class AliDxHFECorrelation : public TNamed {
   /// add control object to list, the base class becomes owner of the object
   int AddControlObject(TObject* pObj);
 
+ /// set the dimension of THn and allocate filling array
+  void InitTHnSparseArray(int dimension) {
+    fDimThn=dimension; 
+    if (fCorrArray) delete[] fCorrArray; fCorrArray=NULL;
+    if (dimension>0) fCorrArray=new Double_t[dimension];
+  }
+
+  inline Double_t* ParticleProperties() const {return fCorrArray;}
+
  private:
   /// copy constructor
   AliDxHFECorrelation(const AliDxHFECorrelation& other);
@@ -101,16 +137,26 @@ class AliDxHFECorrelation : public TNamed {
   // Solved by marking fhEventControlCorr as transient, the cause, though, is not
   // understood
 
-  TObjArray* fHistograms;     //  the histograms - for the moment not in use. 
-  TList* fControlObjects;     //  list of control objects
-  THnSparse* fCorrProperties; //  the Correlation properties of selected particles
-  TH1* fhEventControlCorr;    //! event control histogram (saved via control object list)
-  AliRDHFCutsD0toKpi *fCuts;  //! Cuts 
-  Bool_t fUseMC;              //! use MC info
+  TObjArray* fHistograms;        //  the histograms - for the moment not in use. 
+  TList* fControlObjects;        //  list of control objects
+  THnSparse* fCorrProperties;    //  the Correlation properties of selected particles
+  TH1* fhEventControlCorr;       //! event control histogram (saved via control object list)
+  AliAnalysisCuts *fCuts;        //! Cuts 
+  Bool_t fUseMC;                 // use MC info
+  AliHFCorrelator *fCorrelator;  //! object for correlations
+  Bool_t fUseEventMixing;        // Run Event Mixing analysis
+  Short_t fSystem;               // Which system pp/PbPb
+  Double_t fMinPhi;              // Holds min phi
+  Double_t fMaxPhi;              // Holds maxa phi
+  Double_t fDeltaPhi;            // Delta Phi  
+  Double_t fDeltaEta;            // Delta Eta
+  int fDimThn;                   // Holds dim of THnSparse
+  Double_t* fCorrArray;          //! filling array for THnSparse
+  Int_t fEventType;              // Event type. Only needed for MC (fix)
+
 
   static const char* fgkEventControlBinNames[];
-  static const char* fgkCorrControlBinNames[];
 
-  ClassDef(AliDxHFECorrelation, 3)
+  ClassDef(AliDxHFECorrelation, 4)
 };
 #endif
diff --git a/PWGHF/correlationHF/AliDxHFECorrelationMC.cxx b/PWGHF/correlationHF/AliDxHFECorrelationMC.cxx
new file mode 100644 (file)
index 0000000..ab04265
--- /dev/null
@@ -0,0 +1,141 @@
+
+//**************************************************************************
+//* 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>        *
+//*                  Sedat Altinpinar <Sedat.Altinpinar@cern.ch>           *
+//*                  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   AliDxHFECorrelationMC.cxx
+/// @author Sedat Altinpinar, Hege Erdal, Matthias Richter
+/// @date   2012-11-02
+/// @brief  Worker class for DxHFE correlation
+///
+
+#include "AliDxHFECorrelationMC.h"
+#include "AliVParticle.h"
+#include "AliLog.h"
+//#include "AliAnalysisCuts.h"         // required dependency libANALYSISalice.so
+//#include "AliFlowTrackSimple.h"      // required dependency libPWGflowBase.so
+//#include "AliFlowCandidateTrack.h"   // required dependency libPWGflowTasks.so
+//#include "AliCFContainer.h"          // required dependency libCORRFW.so
+//#include "AliAODRecoDecayHF2Prong.h"   // libPWGHFvertexingHF
+#include "TObjArray.h"
+#include "AliDxHFECorrelation.h"
+#include "TMath.h"
+#include "TFile.h"
+#include "TCanvas.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+#include "AliReducedParticle.h"
+#include "AliDxHFEParticleSelection.h"
+#include <iostream>
+#include <cerrno>
+#include <memory>
+
+using namespace std;
+
+ClassImp(AliDxHFECorrelationMC)
+
+AliDxHFECorrelationMC::AliDxHFECorrelationMC(const char* name)
+  : AliDxHFECorrelation(name?name:"AliDxHFECorrelationMC")
+  , fMCEventType(0)
+{
+  // default constructor
+  // 
+  //
+
+}
+
+
+AliDxHFECorrelationMC::~AliDxHFECorrelationMC()
+{
+  // destructor
+  //
+  //
+
+}
+
+THnSparse* AliDxHFECorrelationMC::DefineTHnSparse()
+{
+  //
+  // Defines the THnSparse. 
+
+  // here is the only place to change the dimension
+  static const int sizeEventdphi = 10;  
+  InitTHnSparseArray(sizeEventdphi);
+  const double pi=TMath::Pi();
+  Double_t minPhi=GetMinPhi();
+  Double_t maxPhi=GetMaxPhi();
+
+  // TODO: Everything here needed for eventmixing? 
+  //                                           0        1     2      3      4     5       6      7    8      9
+  //                                         D0invmass  PtD0 PhiD0 PtbinD0 Pte  dphi    dEta OrigD0 origEl process
+  int         binsEventdphi[sizeEventdphi] = {   200,   1000,  100,  21,   1000, 100,     100,   10,     14,  100 };
+  double      minEventdphi [sizeEventdphi] = { 1.5648,     0,    0,   0,     0 , minPhi, -0.9, -1.5,   -1.5, -0.5 };
+  double      maxEventdphi [sizeEventdphi] = { 2.1648,   100, 2*pi,  20,    100, maxPhi,  0.9,  8.5,   12.5, 99.5 };
+  const char* nameEventdphi[sizeEventdphi] = {
+    "D0InvMass",
+    "PtD0",
+    "PhiD0",
+    "PtBinD0",
+    "PtEl",
+    "#Delta#Phi",
+    "#Delta#eta", 
+    "Origin D0", 
+    "Origin Electron",
+    "Original Process"
+};
+
+  TString name;
+  name.Form("%s info", GetName());
+
+
+  return CreateControlTHnSparse(name,sizeEventdphi,binsEventdphi,minEventdphi,maxEventdphi,nameEventdphi);
+
+
+}
+
+
+int AliDxHFECorrelationMC::FillParticleProperties(AliVParticle* tr, AliVParticle *as, Double_t* data, int dimension) const
+{
+  // fill the data array from the particle data
+  if (!data) return -EINVAL;
+  AliReducedParticle *ptrigger=(AliReducedParticle*)tr;
+  AliReducedParticle *assoc=(AliReducedParticle*)as;
+  if (!ptrigger || !assoc) return -ENODATA;
+  int i=0;
+  if (dimension!=GetDimTHnSparse()) {
+    // TODO: think about filling only the available data and throwing a warning
+    return -ENOSPC;
+  }
+  data[i++]=ptrigger->GetInvMass();
+  data[i++]=ptrigger->Pt();
+  data[i++]=ptrigger->Phi();
+  data[i++]=ptrigger->GetPtBin(); 
+  data[i++]=assoc->Pt();
+  data[i++]=AliDxHFECorrelation::GetDeltaPhi();
+  data[i++]=AliDxHFECorrelation::GetDeltaEta();
+  data[i++]=ptrigger->GetOriginMother();
+  data[i++]=assoc->GetOriginMother();
+  data[i++]=fMCEventType;
+  
+  return i;
+}
+
+int AliDxHFECorrelationMC::Fill(const TObjArray* triggerCandidates, TObjArray* associatedTracks, const AliVEvent* pEvent)
+{
+  // TODO: Implement more on MC?? (Needed?)
+  return AliDxHFECorrelation::Fill(triggerCandidates,associatedTracks,pEvent);
+}
+
diff --git a/PWGHF/correlationHF/AliDxHFECorrelationMC.h b/PWGHF/correlationHF/AliDxHFECorrelationMC.h
new file mode 100644 (file)
index 0000000..26d508f
--- /dev/null
@@ -0,0 +1,51 @@
+//-*- Mode: C++ -*-
+
+//* 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   AliDxHFECorrelationMC.h
+/// @author Sedat Altinpinar, Hege Erdal, Matthias Richter
+/// @date   2012-11-02
+/// @brief  Worker class for DxHFE correlation on MC
+///
+
+#ifndef ALIDXHFECORRELATIONMC_H
+#define ALIDXHFECORRELATIONMC_H
+
+#include "AliDxHFECorrelation.h"
+
+class AliHFCorrelator;
+class AliVEvent;
+class AliVParticle;
+class TObjArray;
+
+class AliDxHFECorrelationMC : public AliDxHFECorrelation {
+ public:
+  /// default constructor
+  AliDxHFECorrelationMC(const char* name=NULL);
+  /// destructor
+  virtual ~AliDxHFECorrelationMC();
+
+  /// fill histograms from particles
+  virtual int Fill(const TObjArray* candidatesD0, TObjArray* candidatesElectron, const AliVEvent* pEvent);
+
+  /// histogram event properties
+  virtual THnSparse* DefineTHnSparse();
+  virtual int FillParticleProperties(AliVParticle* tr, AliVParticle *as, Double_t* data, int dimension) const;
+
+  virtual void SetEventType(int type){fMCEventType=type;}
+
+ protected:
+
+ private:
+  /// copy constructor
+  AliDxHFECorrelationMC(const AliDxHFECorrelationMC& other);
+  /// assignment operator
+  AliDxHFECorrelationMC& operator=(const AliDxHFECorrelationMC& other);
+
+  int  fMCEventType;  // Holds MC Event type, retrieved from MCHeader
+
+  ClassDef(AliDxHFECorrelationMC, 1)
+};
+#endif
index 71ae886..d5f6612 100644 (file)
@@ -32,6 +32,7 @@
 #include "TMath.h"
 #include "TH1D.h"
 #include "THnSparse.h"
+#include "AliReducedParticle.h"
 #include "TFile.h"
 #include <iostream>
 #include <cerrno>
@@ -240,13 +241,14 @@ TObjArray* AliDxHFEParticleSelection::Select(const AliVEvent* pEvent)
   if (!pEvent) return NULL;
   TObjArray* selectedTracks=new TObjArray;
   if (!selectedTracks) return NULL;
+  selectedTracks->SetOwner(); // creating new track objects below
   int nofTracks=pEvent->GetNumberOfTracks();
   for (int itrack=0; itrack<nofTracks; itrack++) {
     AliVParticle* track=pEvent->GetTrack(itrack);
     int selectionCode=IsSelected(track,pEvent);
     HistogramParticleProperties(track, selectionCode);
     if (selectionCode==0) continue;
-    selectedTracks->Add(track);
+    selectedTracks->Add(CreateParticle(track));
   }
   return selectedTracks;
 }
@@ -259,6 +261,7 @@ TObjArray* AliDxHFEParticleSelection::Select(TObjArray* pParticles, const AliVEv
   if (!pParticles) return NULL;
   TObjArray* selectedTracks=new TObjArray;
   if (!selectedTracks) return NULL;
+  selectedTracks->SetOwner(); // creating new track objects below
   TIter next(pParticles);
   TObject* pObj=NULL;
   while ((pObj=next())) {
@@ -267,7 +270,7 @@ TObjArray* AliDxHFEParticleSelection::Select(TObjArray* pParticles, const AliVEv
     int selectionCode=IsSelected(track, pEvent);
     HistogramParticleProperties(track, selectionCode);
     if (selectionCode ==0) continue;
-    selectedTracks->Add(track);
+    selectedTracks->Add(CreateParticle(track));
   }
   return selectedTracks;
 }
@@ -349,3 +352,12 @@ TObject* AliDxHFEParticleSelection::FindObject(const TObject* obj) const
   }
   return NULL;
 }
+
+AliVParticle *AliDxHFEParticleSelection::CreateParticle(AliVParticle* track)
+{
+  // Creating object with reduced particle properties
+  AliReducedParticle *part = new AliReducedParticle(track->Eta(), track->Phi(), track->Pt(), track->Charge(), 0);
+
+  return part;
+
+}
index b276745..3dccad0 100644 (file)
@@ -80,6 +80,7 @@ class AliDxHFEParticleSelection : public TNamed {
   virtual int HistogramEventProperties(int bin);
 
   virtual int FillParticleProperties(AliVParticle* p, Double_t* date, int dimension) const;
+  virtual AliVParticle* CreateParticle(AliVParticle* track);
 
   /// check and add track to internal array
   int CheckAndAdd(AliVParticle* p);
index dbbfa9b..3c1d8bc 100644 (file)
@@ -33,6 +33,7 @@
 #include "AliRDHFCutsD0toKpi.h"
 #include "TObjArray.h"
 #include "THnSparse.h"
+#include "AliReducedParticle.h"
 #include "TAxis.h"
 #include "TString.h"
 #include <iostream>
@@ -113,8 +114,7 @@ int AliDxHFEParticleSelectionD0::InitControlObjects()
 THnSparse* AliDxHFEParticleSelectionD0::DefineTHnSparse()
 {
   //
-  // Defines the THnSparse. For now, only calls CreateControlTHnSparse
-  // TODO: remove pt?? (Have ptbin)
+  // Defines the THnSparse. 
 
   // here is the only place to change the dimension
   const int thnSize2 = 5;
@@ -124,12 +124,12 @@ THnSparse* AliDxHFEParticleSelectionD0::DefineTHnSparse()
   TString name;
   name.Form("%s info", GetName());
 
-  //                                0    1      2      3          4    
+  //                                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. };
-  const char* thnNames[thnSize2] = { "Pt","Phi","Ptbin","D0InvMass","Eta"};
+  int         thnBins [thnSize2] = {1000, 200,  15,     200,     500 };
+  double      thnMin  [thnSize2] = {  0,    0,   0,    1.5648,   -1. };
+  double      thnMax  [thnSize2] = { 100, 2*Pi, 14,    2.1648,    1. };
+  const char* thnNames[thnSize2] = {"Pt", "Phi","Ptbin","D0InvMass","Eta"};
 
   return CreateControlTHnSparse(name,thnSize2,thnBins,thnMin,thnMax,thnNames);
 }
@@ -138,7 +138,6 @@ int AliDxHFEParticleSelectionD0::FillParticleProperties(AliVParticle* p, Double_
 {
   // 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;
@@ -213,15 +212,16 @@ int AliDxHFEParticleSelectionD0::HistogramParticleProperties(AliVParticle* p, in
     return 0;
   }
  
-  // Only D0s are filled 
-  // TODO: Also include D0bar
-  if ((selectionCode==1 || selectionCode==3) && fFillOnlyD0D0bar<2) {
-    fD0InvMass= part->InvMassD0();
-    fPtBin=fCuts->PtBin(part->Pt());
+  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: 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()};
+
+  // Fills only for D0 or both.. 
+  if ((selectionCode==1 || selectionCode==3) && fFillOnlyD0D0bar<2) {
 
     if(fD0Properties && ParticleProperties()) {
       memset(ParticleProperties(), 0, GetDimTHnSparse()*sizeof(ParticleProperties()[0]));
@@ -231,7 +231,19 @@ int AliDxHFEParticleSelectionD0::HistogramParticleProperties(AliVParticle* p, in
     if(fD0Daughter0) fD0Daughter0->Fill(piProperties);
     if(fD0Daughter1) fD0Daughter1->Fill(KProperties);
   }
-
+  else{
+    // If not D0 (or both), check for D0bar (actually now also checks for both, not sure if this is needed)
+    if ((selectionCode==2 || selectionCode==3) && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) {
+
+      if(fD0Properties && ParticleProperties()) {
+       memset(ParticleProperties(), 0, GetDimTHnSparse()*sizeof(ParticleProperties()[0]));
+       FillParticleProperties(p, ParticleProperties(), GetDimTHnSparse());
+       fD0Properties->Fill(ParticleProperties());
+      }
+      if(fD0Daughter0) fD0Daughter0->Fill(piProperties);
+      if(fD0Daughter1) fD0Daughter1->Fill(KProperties);
+    }
+  }
   return 0;
 }
 
@@ -242,6 +254,7 @@ TObjArray* AliDxHFEParticleSelectionD0::Select(TObjArray* pTracks, const AliVEve
   if (!pTracks) return NULL;
   TObjArray* selectedTracks=new TObjArray;
   if (!selectedTracks) return NULL;
+  selectedTracks->SetOwner();
   TIter itrack(pTracks);
   TObject* pObj=NULL;
   while ((pObj=itrack())!=NULL) {
@@ -250,11 +263,18 @@ TObjArray* AliDxHFEParticleSelectionD0::Select(TObjArray* pTracks, const AliVEve
     int selectionCode=IsSelected(track,pEvent);
     HistogramParticleProperties(track, selectionCode);
 
-    //TODO: Also add selection for D0bar
+    // This should make sure the array only gets filled with D0, D0bar or both:
+
     // 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;
-    selectedTracks->Add(track);
+    if ((selectionCode==1 || selectionCode==3) && fFillOnlyD0D0bar<2) 
+      selectedTracks->Add(CreateParticle(track));
+    else{
+      // Add track if it is either defined as D0bar(selectionCode==2) or both 
+      // D0bar and a D0 (selectionCode==3)
+      if ((selectionCode==2 || selectionCode==3) && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) 
+       selectedTracks->Add(CreateParticle(track));
+    }
   }
   return selectedTracks;
 }
@@ -279,7 +299,8 @@ int AliDxHFEParticleSelectionD0::IsSelected(AliVParticle* p, const AliVEvent* pE
   AliRDHFCuts* cuts=const_cast<AliRDHFCuts*>(fCuts);
   if (!cuts) {
     selectionCode=0;
-  } else if(cuts->IsInFiducialAcceptance(d0->Pt(),d0->Y(421)) ) {
+  } 
+  else if(cuts->IsInFiducialAcceptance(d0->Pt(),d0->Y(421)) ) {
 
     Int_t ptbin=cuts->PtBin(d0->Pt());
     if(ptbin==-1) {
@@ -319,3 +340,15 @@ void AliDxHFEParticleSelectionD0::SetCuts(TObject* cuts, int /*level*/)
     AliError(Form("cuts object is not of required type AliRDHFCuts but %s", cuts->ClassName()));
   }
 }
+
+AliVParticle *AliDxHFEParticleSelectionD0::CreateParticle(AliVParticle* track)
+{
+  //
+  // Creates object containing only the variables needed for correlation
+  // 
+
+  AliReducedParticle *part = new AliReducedParticle(track->Eta(), track->Phi(), track->Pt(),fD0InvMass,fPtBin);
+
+  return part;
+
+}
index eaa1078..fa043c1 100644 (file)
@@ -35,6 +35,7 @@ class AliDxHFEParticleSelectionD0 : public AliDxHFEParticleSelection {
   virtual int InitControlObjects();
   virtual THnSparse* DefineTHnSparse();
   virtual int FillParticleProperties(AliVParticle* p, Double_t* date, int dimension) const;
+  virtual AliVParticle* CreateParticle(AliVParticle* track);
 
   //Function for daughter control objects
   //TODO: move to AliDxHFEParticleSelection to be used for several particles?
index 926d170..686eaed 100644 (file)
@@ -108,7 +108,7 @@ int AliDxHFEParticleSelectionEl::Init()
   iResult=AliDxHFEParticleSelection::Init();
   if (iResult<0) return iResult;
 
-  //--------Initialize correction Framework and Cuts
+  //--------Initialize correction Framework and Cuts-------------------------
   // Consider moving this, either to separate function or
   // add a set function for AliCFManager
   // Do we need this? Can we just call AliHFEcuts::CheckParticleCuts
index c01c5d8..b4e4898 100644 (file)
@@ -27,6 +27,8 @@
 #include "AliAODTrack.h"
 #include "AliAODMCParticle.h"
 #include "AliVParticle.h"
+#include "AliReducedParticle.h"
+#include "TH1F.h"
 #include <iostream>
 #include <cerrno>
 #include <memory>
@@ -39,8 +41,10 @@ ClassImp(AliDxHFEParticleSelectionMCD0)
 AliDxHFEParticleSelectionMCD0::AliDxHFEParticleSelectionMCD0(const char* opt)
   : AliDxHFEParticleSelectionD0(opt)
   , fMCTools()
+  , fPDGnotMCD0(NULL)
   , fResultMC(0)
   , fOriginMother(0)
+  , fUseKine(kFALSE)
 {
   // constructor
   // 
@@ -63,30 +67,32 @@ THnSparse* AliDxHFEParticleSelectionMCD0::DefineTHnSparse()
 {
   //
   // Defines the THnSparse.
-  // could/should remove Pt and leave Ptbin
 
   // here is the only place to change the dimension
-  const int thnSize2 = 7;
+  const int thnSize2 = 6;
   InitTHnSparseArray(thnSize2);
   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  };
+  //                                0     1      2       3        4     5
+  //                                Pt   Phi   Ptbin  D0InvMass  Eta   mother 
+  int         thnBins [thnSize2] = {1000, 200,   15,     200,     500,    10  };
+  double      thnMin  [thnSize2] = {   0,  0,     0,    1.5648,   -1.,  -1.5  };
+  double      thnMax  [thnSize2] = { 100, 2*Pi,  14,    2.1648,    1.,   8.5  };
   const char* thnNames[thnSize2] = {
     "Pt",
     "Phi",
     "Ptbin", 
     "D0InvMass", 
     "Eta",
-    "Statistics D0",
-    "Mother of D0"
+    "Mother of D0"  // Bin -1 = not MC truth D0, rest OK
   };
 
+  // Add Histo displaying pdg of D0 candidates not passing MatchToMC()
+  fPDGnotMCD0= new TH1F("fPDGnotMCD0","PDG of track not MC truth D0",1002,-2.5,999.5);
+  AddControlObject(fPDGnotMCD0);
+
   return CreateControlTHnSparse(name,thnSize2,thnBins,thnMin,thnMax,thnNames);
 }
 
@@ -106,7 +112,6 @@ int AliDxHFEParticleSelectionMCD0::FillParticleProperties(AliVParticle* p, Doubl
   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;
@@ -143,7 +148,7 @@ int AliDxHFEParticleSelectionMCD0::IsSelected(AliVParticle* p, const AliVEvent*
   // if not mc selected, however skip this for the moment, because of
   // the logic outside
   fResultMC=CheckMC(p, pEvent);
+
   return iResult;
 }
 
@@ -178,6 +183,25 @@ int AliDxHFEParticleSelectionMCD0::CheckMC(AliVParticle* p, const AliVEvent* pEv
   MClabel=particle->MatchToMC(AliDxHFEToolsMC::kPDGD0,fMCArray,2,pdgDgD0toKpi); 
   
   if(MClabel<0){
+    // Checking PDG of particle if not MC truth D0
+    // TODO: done the right way??
+    Int_t MCl = p->GetLabel();
+    if(MCl<0) {
+      fPDGnotMCD0->Fill(-2);
+      return 0;
+    }
+    int pdgPart=-1;
+    AliAODMCParticle* aodmcp=0;
+    aodmcp=dynamic_cast<AliAODMCParticle*>(fMCArray->At(MCl));
+    if (aodmcp)
+      pdgPart=TMath::Abs(aodmcp->GetPdgCode());
+    if (pdgPart<0){
+      fPDGnotMCD0->Fill(-1);
+      return 0;
+    }
+    else{
+      fPDGnotMCD0->Fill(pdgPart);
+    }
     fOriginMother=-1;
     return 0;
   }
@@ -194,3 +218,15 @@ void AliDxHFEParticleSelectionMCD0::Clear(const char* option)
   /// clear internal memory
   fMCTools.Clear(option);
 }
+
+AliVParticle *AliDxHFEParticleSelectionMCD0::CreateParticle(AliVParticle* track)
+{
+  //
+  //Created object which contain variables needed for correlation. 
+  //
+
+  AliReducedParticle *part = new AliReducedParticle(track->Eta(), track->Phi(), track->Pt(),AliDxHFEParticleSelectionMCD0::GetInvMass(),AliDxHFEParticleSelectionMCD0::GetPtBin(), fOriginMother);
+
+  return part;
+
+}
index 4664799..93ae087 100644 (file)
@@ -17,6 +17,8 @@
 #include "AliDxHFEParticleSelectionD0.h"
 #include "AliDxHFEToolsMC.h"
 
+class TH1;
+
 /**
  * @class AliDxHFEParticleSelectionMCD0
  * Monte Carlo D0 selection for D-HFE correlations, implements the specific
@@ -34,10 +36,14 @@ class AliDxHFEParticleSelectionMCD0 : public AliDxHFEParticleSelectionD0 {
 
   virtual THnSparse* DefineTHnSparse();
   virtual int FillParticleProperties(AliVParticle* p, Double_t* date, int dimension) const;
+  virtual AliVParticle* CreateParticle(AliVParticle* track);
 
   /// check MC criteria
   int CheckMC(AliVParticle* p, const AliVEvent* pEvent);
 
+  /// Flag to run over MC "stack". Not used at the moment
+  void SetUseKine(bool kine){fUseKine=kine;}
+
   /// clear internal memory
   virtual void Clear(const char* option="");
 
@@ -50,10 +56,12 @@ class AliDxHFEParticleSelectionMCD0 : public AliDxHFEParticleSelectionD0 {
   AliDxHFEParticleSelectionMCD0& operator=(const AliDxHFEParticleSelectionMCD0&);
 
   AliDxHFEToolsMC fMCTools;  // MC selction tools
+  TH1* fPDGnotMCD0;          // holds PDG of not MC truth D0s
   int fResultMC;             // Result on MC check
   int fOriginMother;         // Holds info on the original mother particle
+  bool fUseKine;             // Whether to run over MC particles (true) or Reco (false)
 
-  ClassDef(AliDxHFEParticleSelectionMCD0, 1);
+  ClassDef(AliDxHFEParticleSelectionMCD0, 2);
 };
 
 #endif
index 0874d84..731ecac 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,8 +45,7 @@ ClassImp(AliDxHFEParticleSelectionMCEl)
 AliDxHFEParticleSelectionMCEl::AliDxHFEParticleSelectionMCEl(const char* opt)
   : AliDxHFEParticleSelectionEl(opt)
   , fMCTools()
-  , fElectronProperties(NULL)
-  , fWhichCut(NULL)
+  , fPDGnotMCElectron(NULL)
   , fOriginMother(0)
   , fResultMC(0)
 {
@@ -78,9 +78,47 @@ 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
+  // 
+  int iResult=0;
+  iResult=AliDxHFEParticleSelectionEl::Init();
+  if (iResult<0) return iResult;
+
+  // Histo containing PDG of track which was not MC truth electron
+  fPDGnotMCElectron= new TH1F("fPDGnotMCElectron","PDG of track not MC truth electron",AliDxHFEToolsMC::kNofPDGLabels,-0.5,AliDxHFEToolsMC::kNofPDGLabels-0.5);
+  for (int iLabel=0; iLabel<AliDxHFEToolsMC::kNofPDGLabels; iLabel++)
+    fPDGnotMCElectron->GetXaxis()->SetBinLabel(iLabel+1, fgkPDGBinLabels[iLabel]);
+  AddControlObject(fPDGnotMCElectron);
+  return 0;
 }
 
 THnSparse* AliDxHFEParticleSelectionMCEl::DefineTHnSparse()
@@ -88,24 +126,22 @@ THnSparse* AliDxHFEParticleSelectionMCEl::DefineTHnSparse()
   //
   // Defines the THnSparse. 
 
-  const int thnSize = 6;
+  const int thnSize = 4;
   InitTHnSparseArray(thnSize);
   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 };
+  //                          0    1      2     3     
+  //                          Pt   Phi   Eta   mother 
+  int    thnBins[thnSize] = { 1000,  200, 500,    14  };
+  double thnMin [thnSize] = {    0,    0, -1.,  -1.5  };
+  double thnMax [thnSize] = {  100, 2*Pi,  1.,  12.5  };
   const char* thnNames[thnSize]={
     "Pt",
     "Phi",
     "Eta", 
-    "MC statistics electron",
-    "Mother", 
-    "PDG of selected electron"
+    "Mother", //bin==-1: Not MC truth electron
   };
  
   return CreateControlTHnSparse(name,thnSize,thnBins,thnMin,thnMax,thnNames);
@@ -125,9 +161,7 @@ int AliDxHFEParticleSelectionMCEl::FillParticleProperties(AliVParticle* p, Doubl
   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? 
+  data[i++]=fOriginMother; 
   
   return i;
 }
@@ -182,10 +216,11 @@ int AliDxHFEParticleSelectionMCEl::CheckMC(AliVParticle* p, const AliVEvent* pEv
     return 0; // no meaningful filtering on mc possible
   }
  
-  if (fMCTools.RejectByPDG(p,false)) {
+  int pdgParticle=-1;
+  if (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;
   }
 
@@ -206,7 +241,7 @@ int AliDxHFEParticleSelectionMCEl::CheckMC(AliVParticle* p, const AliVEvent* pEv
     fOriginMother=fMCTools.GetOriginMother();
   }
   else{
-    //Could this be done in a more elegant way?
+    //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;
@@ -231,3 +266,12 @@ 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;
+
+}
index c02dc30..019a8cd 100644 (file)
@@ -32,10 +32,15 @@ class AliDxHFEParticleSelectionMCEl : public AliDxHFEParticleSelectionEl {
   /// destructor
   virtual ~AliDxHFEParticleSelectionMCEl();
 
+  //Overload Init function from electron selection class
+  virtual int Init();
+
   // Setting up control objects: overloaded from AliDxHFEParticleSelection
   virtual THnSparse* DefineTHnSparse();
 
   virtual int FillParticleProperties(AliVParticle* p, Double_t* date, int dimension) const;
+  virtual AliVParticle* CreateParticle(AliVParticle* track);
+
 
   /// overloaded from AliDxHFEParticleSelection: check particle
   virtual int IsSelected(AliVParticle* p, const AliVEvent *pEvent=NULL);
@@ -56,14 +61,14 @@ class AliDxHFEParticleSelectionMCEl : public AliDxHFEParticleSelectionEl {
 
   /// TODO: check if the label definitions can be used from the ToolsMC
   static const char* fgkPDGMotherBinLabels[];
-
+  static const char* fgkPDGBinLabels[];
+                    
   AliDxHFEToolsMC fMCTools;            // MC selection tools
-  THnSparse*      fElectronProperties; // the particle properties of selected particles
-  TH1*            fWhichCut;           // effective cut for a rejected particle
+  TH1*            fPDGnotMCElectron;   //! PDG of track not MC truth electron
   int fOriginMother;                   //  Holds the origin motherquark (process)
   int fResultMC;                       // Holds information on check on MC
 
-  ClassDef(AliDxHFEParticleSelectionMCEl, 1);
+  ClassDef(AliDxHFEParticleSelectionMCEl, 2);
 };
 
 #endif
index e8ec3a3..3e70d3e 100644 (file)
@@ -48,6 +48,7 @@ AliDxHFEToolsMC::AliDxHFEToolsMC(const char* option)
   , fHistPDGMother(NULL)
   , fOriginMother(kOriginNone)
   , fMClabel(-1)
+  , fNrMCParticles(-1)
 {
   // constructor
   // 
@@ -168,7 +169,7 @@ int AliDxHFEToolsMC::InitMCParticles(const AliVEvent* pEvent)
     AliWarningClass(Form("ignoring MC info '%s' of wrong type '%s', expecting TObjArray", branchname.Data(), o->ClassName()));
     return -ENODATA;
   }
-
+  fNrMCParticles=fMCParticles->GetEntriesFast();
   return 0;
 }
 
@@ -184,7 +185,7 @@ bool AliDxHFEToolsMC::RejectByPDG(int pdg, const vector<int> &list) const
   return true;
 }
 
-bool AliDxHFEToolsMC::RejectByPDG(AliVParticle* p, bool doStatistics)
+bool AliDxHFEToolsMC::RejectByPDG(AliVParticle* p, bool doStatistics, int* pdgParticleResult)
 {
   // check if pdg should be rejected
   // always false if not pdg list is initialized
@@ -200,6 +201,9 @@ bool AliDxHFEToolsMC::RejectByPDG(AliVParticle* p, bool doStatistics)
     pdgPart=TMath::Abs(aodmcp->GetPdgCode());
   if (pdgPart<0) return 0;
 
+  if (pdgParticleResult)
+    *pdgParticleResult=pdgPart;
+
   bool bReject=RejectByPDG(pdgPart, fPDGs);
   if (doStatistics && fHistPDG) {
     // TODO: think about histogramming mode, e.g. histogramming of rejected particles?
@@ -383,17 +387,44 @@ void AliDxHFEToolsMC::CheckOriginMother(int pdg)
     else if(fOriginMother==kOriginBeauty) fOriginMother=kOriginGluonBeauty;
     else fOriginMother=kOriginGluon;
     break;
-  case(kPDGd): 
-    fOriginMother=kOriginDown; break;
+  case(kPDGd):
+    if(!TestIfHFquark(fOriginMother))
+      fOriginMother=kOriginDown; 
+    break;
   case(kPDGu):
-    fOriginMother=kOriginUp; break;
+    if(!TestIfHFquark(fOriginMother))
+      fOriginMother=kOriginUp; 
+    break;
   case(kPDGs):
-    fOriginMother=kOriginStrange; break;
+    if(!TestIfHFquark(fOriginMother))
+      fOriginMother=kOriginStrange; 
+    break;
+  }
+}
+
+Bool_t AliDxHFEToolsMC::TestIfHFquark(int origin)
+{
+
+  // Checking if particle has been marked as charm/beauty quark
+  
+  Bool_t test=kFALSE;
+  switch(origin){
+  case(kOriginCharm):
+    test=kTRUE; break;
+  case(kOriginBeauty): 
+    test=kTRUE; break;
+  case(kOriginGluonCharm): 
+    test=kTRUE; break;
+  case(kOriginGluonBeauty): 
+    test=kTRUE; break;
   }
+  return test; 
 }
 
 void AliDxHFEToolsMC::Clear(const char* /*option*/)
 {
   // clear internal memory
   fMCParticles=NULL;
+  fNrMCParticles=-1;
 }
+
index 5337452..0cb85d1 100644 (file)
@@ -133,7 +133,10 @@ class AliDxHFEToolsMC {
 
   /// check if pdg should be rejected
   /// always false if pdg list is not initialized
-  bool RejectByPDG(AliVParticle* p, bool doStatistics=true);
+  bool RejectByPDG(AliVParticle* p, bool doStatistics=true, int* pdgParticleResult=NULL);
+  bool RejectByPDG(AliVParticle* p, int* pdgParticleResult) {
+    return RejectByPDG(p, true, pdgParticleResult);
+  }
 
   /// check if pdg should be rejected by mother
   /// always false if mother pdg list is not initialized
@@ -149,6 +152,8 @@ class AliDxHFEToolsMC {
 
   // Compare pdg to quark and gluon
   void CheckOriginMother(int pdg);
+  // Tests if particle have been marked as HF quark 
+  Bool_t TestIfHFquark(int origin);
 
   // Setting MC label from outside
   void SetMClabel(int mclab){fMClabel=mclab;}
@@ -159,6 +164,13 @@ class AliDxHFEToolsMC {
   /// if it is in the list, returns always false if list is empty
   bool RejectByPDG(int pdg, const vector<int> &list) const;
 
+  int GetNrMCParticles() const {return fNrMCParticles;}
+
+  /// mapping of pdg code to enum
+  int MapPDGLabel(int pdg) const;
+  /// mapping of pdg code to enum
+  int MapPDGMotherLabel(int pdg) const;
+
  protected:
 
  private:
@@ -173,11 +185,6 @@ class AliDxHFEToolsMC {
                              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[];
@@ -190,7 +197,8 @@ class AliDxHFEToolsMC {
   TH1* fHistPDGMother;     //  control histogram pdg of selected particle
   int fOriginMother;       //  Holds the origin motherquark (process)
   int fMClabel;            //  MClabel passed from outside (default =-1)
+  int fNrMCParticles;      //  number of MC particles 
 
-  ClassDef(AliDxHFEToolsMC, 1);
+  ClassDef(AliDxHFEToolsMC, 2);
 };
 #endif