]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Update in D-electon correlation code (Hege, Matthias)
authorprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Jul 2012 08:51:45 +0000 (08:51 +0000)
committerprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Jul 2012 08:51:45 +0000 (08:51 +0000)
14 files changed:
PWGHF/CMakelibPWGHFcorrelationHF.pkg
PWGHF/correlationHF/AliAnalysisTaskDxHFECorrelation.cxx
PWGHF/correlationHF/AliAnalysisTaskDxHFECorrelation.h
PWGHF/correlationHF/AliAnalysisTaskDxHFEParticleSelection.cxx
PWGHF/correlationHF/AliAnalysisTaskDxHFEParticleSelection.h
PWGHF/correlationHF/AliDxHFECorrelation.cxx
PWGHF/correlationHF/AliDxHFECorrelation.h
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/compile-DxHFE.C

index 47d4c1c1d49ac6f3d275bea84fc1a2db961bf881..1751c7249d6ff5691e7901d8699c5c8be63445e6 100644 (file)
@@ -49,8 +49,8 @@ set ( MODULE_DHDR )
 set ( EINCLUDE  PWGHF/vertexingHF PWGHF/hfe include)
 
 set ( ELIBS  
-#   PWGHFhfe            # likely to be linked later
-#   PWGHFvertexingHF    # likely to be linked later
+    PWGHFhfe
+    PWGHFvertexingHF
     ANALYSIS
     ANALYSISalice
     )
index fdef14e3b8b5f4693879847c2805c9283ec8205d..260154f70461fa7b529bdd27a780c7e9b22ecffc 100644 (file)
 #include "AliAnalysisManager.h"
 #include "AliLog.h"
 #include "AliESDInputHandler.h"
+#include "AliAnalysisDataSlot.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliAnalysisManager.h"
+#include "AliVertexerTracks.h"
+#include "AliAODHandler.h"
+#include "AliInputEventHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODRecoCascadeHF.h"
+#include "AliRDHFCutsD0toKpi.h"
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+#include "AliHFEcontainer.h"
+#include "AliHFEpid.h"
+#include "AliHFEpidBase.h"
+#include "AliHFEcuts.h"
+#include "AliHFEtools.h"
+#include "TObject.h"
 #include "TChain.h"
 #include "TSystem.h"
 #include "TFile.h"
 #include <memory>
 
+using namespace std;
+
 /// ROOT macro for the implementation of ROOT specific class methods
 ClassImp(AliAnalysisTaskDxHFECorrelation)
 
@@ -45,12 +69,18 @@ AliAnalysisTaskDxHFECorrelation::AliAnalysisTaskDxHFECorrelation(const char* opt
   , fCorrelation(NULL)
   , fD0s(NULL)
   , fElectrons(NULL)
+  , fCutsD0(NULL)
+  , fCutsHFE(NULL)
+  , fPID(NULL)
+  , fFillOnlyD0D0bar(0)
+  , fUseMC(kFALSE)
 {
   // constructor
   //
   //
-
   DefineSlots();
+  fPID = new AliHFEpid("hfePid");
+
 }
 
 int AliAnalysisTaskDxHFECorrelation::DefineSlots()
@@ -58,6 +88,7 @@ int AliAnalysisTaskDxHFECorrelation::DefineSlots()
   // define the data slots
   DefineInput(0, TChain::Class());
   DefineOutput(1, TList::Class());
+  DefineOutput(2,AliRDHFCutsD0toKpi::Class());
   return 0;
 }
 
@@ -80,56 +111,209 @@ AliAnalysisTaskDxHFECorrelation::~AliAnalysisTaskDxHFECorrelation()
   fElectrons=NULL;
   if (fCorrelation) delete fCorrelation;
   fCorrelation=NULL;
+  if (fCutsD0) delete fCutsD0;
+  fCutsD0=NULL;
+  if (fCutsHFE) delete fCutsHFE;
+  fCutsHFE=NULL;
+  if(fPID) delete fPID;
+  fPID=NULL;
+
+
 }
 
 void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
 {
   // create result objects and add to output list
 
-  std::auto_ptr<TList> Output(new TList);
-  std::auto_ptr<AliDxHFEParticleSelection> D0s(new AliDxHFEParticleSelectionD0);
-  std::auto_ptr<AliDxHFEParticleSelection> Electrons(new AliDxHFEParticleSelectionEl);
-  std::auto_ptr<AliDxHFECorrelation> Correlation(new AliDxHFECorrelation);
+  //Initialize PID for electron selection
+  if(!fPID->GetNumberOfPIDdetectors()) { 
+    fPID->AddDetector("TOF",0);
+    fPID->AddDetector("TPC",1);
+  }
+  fPID->InitializePID();
 
-  if (!Output     .get() ||
-      !D0s        .get() ||
-      !Electrons  .get() ||
-      !Correlation.get()) {
-    AliFatal("allocation of worker classes failed");
-    return;
+  fOutput = new TList;
+  fOutput->SetOwner();
+
+  // setting up for D0s
+  TString selectionD0Options;
+  switch (fFillOnlyD0D0bar) {
+  case 1: selectionD0Options+="FillOnlyD0 "; break;
+  case 2: selectionD0Options+="FillOnlyD0bar "; break;
+  default: selectionD0Options+="FillD0D0bar ";
   }
 
-  fOutput      = Output     .release();
-  fD0s         = D0s        .release();
-  fElectrons   = Electrons  .release();
-  fCorrelation = Correlation.release();
+  fD0s=new AliDxHFEParticleSelectionD0(selectionD0Options);
+  fD0s->SetCuts(fCutsD0);
+  fD0s->InitControlObjects();
 
-  fOutput->SetOwner();
+  //Electrons
+  fElectrons=new AliDxHFEParticleSelectionEl;
+  fElectrons->SetCuts(fPID, AliDxHFEParticleSelectionEl::kCutPID);
+  fElectrons->SetCuts(fCutsHFE, AliDxHFEParticleSelectionEl::kCutHFE);
+  fElectrons->InitControlObjects();
+
+  //Correlation
+  fCorrelation=new AliDxHFECorrelation;
+  fCorrelation->SetCuts(dynamic_cast<AliRDHFCutsD0toKpi*>(fCutsD0));
+  fCorrelation->Init();
+
+  // Fix for merging:
+  // Retrieving the individual objects created
+  // and storing them instead of fD0s,fElectrons etc.. 
+  TList *list =(TList*)fD0s->GetControlObjects();
+  TObject *obj=NULL;
+
+  TIter next(list);
+  while((obj = next())){
+    fOutput->Add(obj);
+  }
+
+  list=(TList*)fCorrelation->GetControlObjects();
+  next=TIter(list);
+  while((obj= next())){
+    fOutput->Add(obj);
+  }
+
+  list=(TList*)fElectrons->GetControlObjects();
+  next=TIter(list);
+  while((obj = next()))
+    fOutput->Add(obj);
+
+  if (fCutsD0) {
+    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(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);
+
 }
 
 void AliAnalysisTaskDxHFECorrelation::UserExec(Option_t* /*option*/)
 {
   // process the event
-
-  // TODO: implement correct input, this is likely not to be the
-  // ESD
   TObject* pInput=InputEvent();
   if (!pInput) {
     AliError("failed to get input");
     return;
   }
   AliVEvent *pEvent = dynamic_cast<AliVEvent*>(pInput);
-  if(!pEvent){
-    AliError(Form("input of wrong class type %s, expecting AliVEvent", pInput->ClassName()));
+  TClonesArray *inputArray=0;
+
+  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.    
+    pEvent = dynamic_cast<AliAODEvent*> (AODEvent());
+    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
+    // have to taken from the AOD event hold by the AliAODExtension
+    AliAODHandler* aodHandler = (AliAODHandler*) 
+      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+    
+    if(aodHandler->GetExtensions()) {
+      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
+      AliAODEvent* aodFromExt = ext->GetAOD();
+      inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
+    }
+  } else if(pEvent) {
+    inputArray=(TClonesArray*)pEvent->GetList()->FindObject("D0toKpi");
+  }
+  if(!inputArray || !pEvent) {
+    AliError("Input branch not found!\n");
+    return;
+  }
+  // fix for temporary bug in ESDfilter
+  // the AODs with null vertex pointer didn't pass the PhysSel
+  if(!pEvent->GetPrimaryVertex() || TMath::Abs(pEvent->GetMagneticField())<0.001){
+    AliDebug(2,"Rejected at GetPrimaryvertex");
+    return;
+  }
+
+  fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsAll);
+  AliRDHFCuts* cutsd0=dynamic_cast<AliRDHFCuts*>(fCutsD0);
+  if (!cutsd0) return; // Fatal thrown already in initialization
+
+  if(!cutsd0->IsEventSelected(pEvent)) {
+    AliDebug(2,"rejected at IsEventSelected");
+    return;
+  }
+
+  if(!fPID->IsInitialized()){ 
+    // Initialize PID with the given run number
+    AliWarning("PID not initialised, get from Run no");
+    fPID->InitializePID(pEvent->GetRunNumber());
+  }
+
+  AliPIDResponse *pidResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler())->GetPIDResponse();
+  if(!pidResponse){
+    AliDebug(1, "Using default PID Response");
+    pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class()); 
+  }
+  fPID->SetPIDResponse(pidResponse);
+
+
+  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()) {
+    return;
+  }
+  Int_t nD0Selected = pSelectedD0s->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");
     return;
   }
+  fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsD0);
+
 
-  std::auto_ptr<TObjArray> pSelectedD0s(fD0s->Select(pEvent));
   std::auto_ptr<TObjArray> pSelectedElectrons(fElectrons->Select(pEvent));
+
+  // TODO: use the array of selected track for something, right now
+  // only the control histograms of the selection class are filled
+  // 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");
+    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());
+
   if (iResult<0) {
     AliError(Form("%s processing failed with error %d", fCorrelation->GetName(), iResult));
   }
index d43cbcd856e74e7fbf95f4cf122db29e23bc1f04..f4ebdd3bf9dba6a77fc15b526f5ee4e4612a2233 100644 (file)
 
 #include "AliAnalysisTaskSE.h"
 #include "TString.h"
+
+class AliPID;
+class AliPIDResponse;
 class TList;
 class AliDxHFEParticleSelection;
+class AliDxHFEParticleSelectionD0;
+class AliDxHFEParticleSelectionEl;
 class AliDxHFECorrelation;
+class AliAnalysisCuts;
+class AliHFEpid;
+class AliHFEcuts;
 
 /**
  * @class AliAnalysisTaskDxHFECorrelation
@@ -47,8 +55,15 @@ class AliAnalysisTaskDxHFECorrelation : public AliAnalysisTaskSE {
 
   /// set options
   void SetOption(const char* opt) { fOption = opt; }
+  void SetFillOnlyD0D0bar(Int_t flagfill){fFillOnlyD0D0bar=flagfill;}
+  virtual void SetUseMC(Bool_t useMC){fUseMC=useMC;}
+  virtual void SetCutsD0(AliAnalysisCuts* cuts){fCutsD0=cuts;}
+  virtual void SetCutsHFE(AliHFEcuts* cuts){fCutsHFE=cuts;}
+
   /// overloaded from TObject: get option
   virtual Option_t* GetOption() const { return fOption;}
+  Int_t  GetFillOnlyD0D0bar() const {return fFillOnlyD0D0bar;}
+  Bool_t GetUseMC() const {return fUseMC;}
 
  protected:
 
@@ -61,12 +76,18 @@ class AliAnalysisTaskDxHFECorrelation : public AliAnalysisTaskSE {
   int DefineSlots();
 
   TList* fOutput;                        //! list send on output slot 1
-  TString fOption;                       // option string
-  AliDxHFECorrelation* fCorrelation;     // correlation worker class
-  AliDxHFEParticleSelection* fD0s;       // selection of D0s
-  AliDxHFEParticleSelection* fElectrons; // selection of electrons
+  TString fOption;                       //  option string
+  AliDxHFECorrelation* fCorrelation;     //  correlation worker class
+  AliDxHFEParticleSelection* fD0s;       //  selection of D0s
+  AliDxHFEParticleSelection* fElectrons; //  selection of electrons
+  AliAnalysisCuts *fCutsD0;              //  Cuts D0 
+  AliHFEcuts *fCutsHFE;                  //  Cuts HFE
+  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
+
 
-  ClassDef(AliAnalysisTaskDxHFECorrelation, 1);
+  ClassDef(AliAnalysisTaskDxHFECorrelation, 2);
 };
 
 #endif
index 689983354708828872c4ca65534cd5c5b3d2c037..71e4662c39e3bf38322013f1159ee19dc557db27 100644 (file)
 #include "AliDxHFEParticleSelection.h"
 #include "AliDxHFEParticleSelectionD0.h"
 #include "AliAnalysisManager.h"
+#include "AliAnalysisCuts.h"
 #include "AliLog.h"
+#include "TH1F.h"
 #include "AliESDInputHandler.h"
+#include "AliAODHandler.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliRDHFCutsD0toKpi.h"
 #include "TChain.h"
 #include "TSystem.h"
 #include "TFile.h"
+#include "TObjArray.h"
 #include <memory>
 
+using namespace std;
+
 /// ROOT macro for the implementation of ROOT specific class methods
 ClassImp(AliAnalysisTaskDxHFEParticleSelection)
 
@@ -41,7 +49,9 @@ AliAnalysisTaskDxHFEParticleSelection::AliAnalysisTaskDxHFEParticleSelection(con
   : AliAnalysisTaskSE("AliAnalysisTaskDxHFEParticleSelection")
   , fOutput(0)
   , fOption(opt)
+  , fCuts(NULL)
   , fSelector(NULL)
+  , fUseMC(kFALSE)
 {
   // constructor
   //
@@ -77,6 +87,11 @@ AliAnalysisTaskDxHFEParticleSelection::~AliAnalysisTaskDxHFEParticleSelection()
     delete fSelector;
     fSelector=NULL;
   }
+
+  if (fCuts) {
+    delete fCuts;
+    fCuts=NULL;
+  }
 }
 
 void AliAnalysisTaskDxHFEParticleSelection::UserCreateOutputObjects()
@@ -86,10 +101,34 @@ void AliAnalysisTaskDxHFEParticleSelection::UserCreateOutputObjects()
   fOutput = new TList;
   fOutput->SetOwner();
 
-  fSelector=new AliDxHFEParticleSelectionD0;
-  fSelector->InitControlObjects();
-  fOutput->Add(fSelector);
-  
+  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");
+  }
+  selector->InitControlObjects();
+
+  fSelector=selector.release();
+
+  // Retrieving the list containing histos and THnSparse
+  // and storing them instead of fSelector
+  // Fix to be able to merge
+  TList *list =(TList*)fSelector->GetControlObjects();
+  TObject *obj=NULL;
+
+  TIter next(list);
+  while((obj = next())){
+    fOutput->Add(obj);
+  }
+
   // all tasks must post data once for all outputs
   PostData(1, fOutput);
 }
@@ -105,13 +144,53 @@ void AliAnalysisTaskDxHFEParticleSelection::UserExec(Option_t* /*option*/)
     AliError("failed to get input");
     return;
   }
+
+  // check if input is an ESD
   AliVEvent *pEvent = dynamic_cast<AliVEvent*>(pInput);
-  if(!pEvent){
-    AliError(Form("input of wrong class type %s, expecting AliVEvent", pInput->ClassName()));
+  TClonesArray *inputArray=0;
+
+  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.    
+    pEvent = AODEvent();
+    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
+    // have to taken from the AOD event hold by the AliAODExtension
+    AliAODHandler* aodHandler = (AliAODHandler*) 
+      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+    
+    if(aodHandler->GetExtensions()) {
+      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
+      AliAODEvent* aodFromExt = ext->GetAOD();
+      inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
+    }
+  } else if(pEvent) {
+    inputArray=(TClonesArray*)pEvent->GetList()->FindObject("D0toKpi");
+  }
+  if(!inputArray || !pEvent) {
+    AliError("Input branch not found!\n");
+    return;
+  }
+  // fix for temporary bug in ESDfilter
+  // the AODs with null vertex pointer didn't pass the PhysSel
+  if(!pEvent->GetPrimaryVertex() || TMath::Abs(pEvent->GetMagneticField())<0.001){
+    AliDebug(2,"Rejected at GetPrimaryvertex");
+    return;
+  }
+
+  fSelector->HistogramEventProperties(AliDxHFEParticleSelection::kEventsAll);
+  AliRDHFCuts* cuts=dynamic_cast<AliRDHFCuts*>(fCuts);
+  if (!cuts) return; // Fatal thrown already in initialization
+
+  if(!cuts->IsEventSelected(pEvent)) {
+    AliDebug(2,"rejected at IsEventSelected");
     return;
   }
 
-  std::auto_ptr<TObjArray> selectedTracks(fSelector->Select(pEvent));
+  Int_t nInD0toKpi = inputArray->GetEntriesFast();
+
+  fSelector->HistogramEventProperties(AliDxHFEParticleSelection::kEventsSel);
+             
+  std::auto_ptr<TObjArray> selectedTracks(fSelector->Select(inputArray,pEvent));
   // TODO: use the array of selected track for something, right now
   // only the control histograms of the selection class are filled
   // note: the pointer is deleted automatically once the scope is left
@@ -119,7 +198,22 @@ void AliAnalysisTaskDxHFEParticleSelection::UserExec(Option_t* /*option*/)
   // first, however some other cleanup will be necessary in that case
   // probably a clone with a reduced AliVParticle implementation is
   // appropriate.
-  if (selectedTracks.get()) {
+
+  if(! selectedTracks.get()) {
+    cout << "No selected D0s in this event" << endl; 
+    return;
+  }
+
+  //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);
+
+  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 d069a85c2fe952cc708f647fab180f1a90057f06..1bb1bdd8a43563fbb64c7904d4a352be28655243 100644 (file)
 
 #include "AliAnalysisTaskSE.h"
 #include "TString.h"
-class TList;
+
 class AliDxHFEParticleSelection;
+class AliAnalysisCuts;
+class TList;
 
 /**
  * @class AliAnalysisTaskDxHFEParticleSelection
@@ -47,6 +49,9 @@ class AliAnalysisTaskDxHFEParticleSelection : public AliAnalysisTaskSE {
   virtual void Terminate(Option_t*);
 
   void SetOption(const char* opt) { fOption = opt; }
+  virtual void SetUseMC(Bool_t useMC){fUseMC=useMC;}
+  virtual void SetCuts(AliAnalysisCuts* cuts){fCuts=cuts;}
+  Bool_t GetUseMC() const {return fUseMC;}
 
  protected:
 
@@ -58,11 +63,13 @@ class AliAnalysisTaskDxHFEParticleSelection : public AliAnalysisTaskSE {
 
   int DefineSlots();
 
-  TList* fOutput;                  //! list send on output slot 1
-  TString fOption;                 // option string
+  TList* fOutput;                       // list send on output slot 1
+  TString fOption;                      //  option string
+  AliAnalysisCuts* fCuts;               //  Cuts 
   AliDxHFEParticleSelection* fSelector; // selector instance
+  bool fUseMC;                          // use MC info
 
-  ClassDef(AliAnalysisTaskDxHFEParticleSelection, 1);
+  ClassDef(AliAnalysisTaskDxHFEParticleSelection, 2);
 };
 
 #endif
index 601b81b10c559a1a06766fcfc510913681124ccf..c0fd9ab53587cbde46acbb1e2304c239c625ff33 100644 (file)
 #include "AliDxHFECorrelation.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 "AliRDHFCutsD0toKpi.h"
 #include "TObjArray.h"
 #include "TH1D.h"
 #include "TH2D.h"
+#include "THnSparse.h"
 #include "TMath.h"
 #include "TFile.h"
 #include "TCanvas.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
 #include <iostream>
 #include <cerrno>
 #include <memory>
@@ -42,7 +51,12 @@ ClassImp(AliDxHFECorrelation)
 
 AliDxHFECorrelation::AliDxHFECorrelation(const char* name)
   : TNamed(name?name:"AliDxHFECorrelation", "")
-  , fHistograms(NULL)
+  , fHistograms(NULL)  
+  , fControlObjects(NULL)
+  , fCorrProperties(NULL)
+  , fhEventControlCorr(NULL)
+  , fCuts(NULL)
+  , fUseMC(kFALSE)
 {
   // default constructor
   // 
@@ -50,93 +64,200 @@ AliDxHFECorrelation::AliDxHFECorrelation(const char* name)
 
 }
 
+const char* AliDxHFECorrelation::fgkEventControlBinNames[]={
+  "nEventsAll",
+  "nEventsSelected",
+  "nEventsD0"
+  "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
   //
   //
-  if (fHistograms) {
-    delete fHistograms;
-    fHistograms=NULL;
-  }
+  if (fHistograms) delete fHistograms;
+  fHistograms=NULL;
+
+  // NOTE: fControlObjects owns the object, and they are deleted in the
+  // destructor of TList
+  if (fControlObjects) delete fControlObjects;
+  fControlObjects=NULL;
+  fCorrProperties=NULL;
+  fhEventControlCorr=NULL;
+
+  // NOTE: the external object is deleted elsewhere
+  fCuts=NULL;
 }
 
 int AliDxHFECorrelation::Init()
 {
-  /// init class and create histograms
-  if (fHistograms) delete fHistograms;
-  fHistograms=new TObjArray;
-  if (!fHistograms) return -ENOMEM;
-  fHistograms->SetOwner(kTRUE);
-  AliInfo(Form("initializing %s ", GetName()));
+  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
 
-  // TODO: This is just a mockup to illustrate the functionality
-  // the specific histograms need to be defined
+  /*
+    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;
+
+  for (iLabel=0; iLabel<sizeEventdphi; iLabel++)
+    CorrProperties->GetAxis(iLabel)->SetTitle(fgkCorrControlBinNames[iLabel]);
+
+  //----------------------------------------------
+  // Histogram for storing event information
+
+  std::auto_ptr<TH1D> hEventControl(new TH1D("hEventControlCorr", "hEventControlCorr", 10, 0, 10));
+  if (!hEventControl.get()) {
+    return -ENOMEM;
+  }
+
+  for (iLabel=0; iLabel<kNEventControlLabels; iLabel++)
+    hEventControl->GetXaxis()->SetBinLabel(iLabel, fgkEventControlBinNames[iLabel]);
+
+  fCorrProperties=CorrProperties.release();
+  AddControlObject(fCorrProperties);
+  fhEventControlCorr=hEventControl.release();
+  AddControlObject(fhEventControlCorr);
+
+  return 0;
+}
+
+int AliDxHFECorrelation::AddControlObject(TObject* pObj)
+{
+  AliInfo("Adding object");
+  /// add control object to list, the base class becomes owner of the object
+  if (!pObj) return -EINVAL;
+  if (!fControlObjects) {
+    fControlObjects=new TList;
+    if (!fControlObjects) return -ENOMEM;
+    fControlObjects->SetOwner();
+  }
+  if (fControlObjects->FindObject(pObj->GetName())) {
+    AliError(Form("ignoring duplicate object '%s' of type %s", pObj->GetName(), pObj->ClassName()));
+    return -EEXIST;
+  }
+  fControlObjects->Add(pObj);
+  return 0;
+}
+
+int AliDxHFECorrelation::HistogramEventProperties(int bin)
+{
+  /// histogram event properties
+  if (!fhEventControlCorr) return 0;
+  fhEventControlCorr->Fill(bin);
 
-  // avoid the objects to be added to the global directory 
-  bool statusAddDirectory=TH1::AddDirectoryStatus();
-  TH1::AddDirectory(false);
-  const double Pii=TMath::Pi();
-  TObjArray* a=fHistograms;
-  a->AddAt(new TH1D("hD0pT"       , "D0 pT"              ,100,0.0,10.0)   , khD0pT        );
-  a->AddAt(new TH1D("hD0Phi"      , "D0 phi"             ,180,0.0,2*Pii)  , khD0Phi       );
-  a->AddAt(new TH1D("hD0Eta"      , "D0 eta"             ,100,-10.0,10.0) , khD0Eta       );
-  a->AddAt(new TH1D("hElectronpT" , "Electron pT"        ,100,0.0,10.0)   , khElectronpT  );
-  a->AddAt(new TH1D("hElectronPhi", "Electron phi"       ,180,0.0,2*Pii)  , khElectronPhi );
-  a->AddAt(new TH1D("hElectronEta", "Electron eta"       ,100,-10.0,10.0) , khElectronEta );
-  a->AddAt(new TH1D("hDeltaPhi"   , "#Delta#Phi D0 - HFE",180,0.0,2*Pii)  , khDeltaPhi    );
-
-  TH1::AddDirectory(statusAddDirectory);
   return 0;
 }
 
-int AliDxHFECorrelation::Fill(const TObjArray* candidatesD0, const TObjArray* candidatesElectron)
+int AliDxHFECorrelation::Fill(const TObjArray* triggerCandidates, const TObjArray* associatedTracks)
 {
-  /// fill histograms from array of AliVParticle objects
-  if (!candidatesD0 || !candidatesElectron) return -EINVAL;
-  if (!fHistograms) {
+  /// fill ThnSparse from array of AliVParticle objects
+  if (!triggerCandidates || !associatedTracks) return -EINVAL;
+  if (!fControlObjects) {
     Init();
   }
-  if (!fHistograms) {
-    AliError("Initialisation failed, can not fill histograms");
+  if (!fControlObjects) {
+    AliError("Initialisation failed, can not fill THnSparse");
   }
 
   const double Pii=TMath::Pi();
 
-  TIter itrigger(candidatesD0);
+  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. 
   while ((otrigger=itrigger())!=NULL) {
     // loop over trigger D0 particle
     ctrigger++;
     AliVParticle* ptrigger=dynamic_cast<AliVParticle*>(otrigger);
+    AliAODRecoDecayHF2Prong *d0 = dynamic_cast<AliAODRecoDecayHF2Prong*>(ptrigger);
+
     if (!ptrigger) continue;
-    ((TH1D*)fHistograms->At(khD0pT)) ->Fill(ptrigger->Pt());
-    ((TH1D*)fHistograms->At(khD0Phi))->Fill(ptrigger->Phi());
-    ((TH1D*)fHistograms->At(khD0Eta))->Fill(ptrigger->Eta());
-    // TODO: add further correlation specific cuts here, e.g acceptance
-    // which are no primarily part of the particle selection
 
-    TIter iElectron(candidatesElectron);
+    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;
-      ((TH1D*)fHistograms->At(khElectronpT)) ->Fill(pElectron->Pt());
-      ((TH1D*)fHistograms->At(khElectronPhi))->Fill(pElectron->Phi());
-      ((TH1D*)fHistograms->At(khElectronEta))->Fill(pElectron->Eta());
 
-      // phi difference to trigger particle
-      Double_t DeltaPhi = ptrigger->Phi() - pElectron->Phi();
-      if (DeltaPhi<-0.5*Pii) DeltaPhi += 2*Pii;
-      if (DeltaPhi>1.5*Pii)  DeltaPhi -= 2*Pii;
-      ((TH1D*)fHistograms->At(khDeltaPhi))->Fill(DeltaPhi);
+      //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 electrons
-  } // loop over D0 trigger particle
+    } // loop over associated tracks
+  } // loop over trigger particle
 
   return 0;
 }
@@ -162,33 +283,13 @@ void AliDxHFECorrelation::Print(Option_t */*option*/) const
 void AliDxHFECorrelation::Draw(Option_t */*option*/)
 {
   /// overloaded from TObject: draw histograms
-  if (fHistograms) {
-    TString name;
-    int canvasno=1;
-    int padno=1;
-    const char* drawoption="";
-    name.Form("%s_%d", GetName(), canvasno++);
-    TCanvas* c=new TCanvas(name);
-    c->SetWindowSize(1600,800);
-    c->SetTitle(Form("%s: particle properties", GetName()));
-    c->Divide(3,2);
-    padno=1;
-    c->cd(padno++); fHistograms->At(khD0pT)       ->Draw(drawoption);
-    c->cd(padno++); fHistograms->At(khD0Phi)      ->Draw(drawoption);
-    c->cd(padno++); fHistograms->At(khD0Eta)      ->Draw(drawoption);
-    c->cd(padno++); fHistograms->At(khDeltaPhi)   ->Draw(drawoption);
-    c->cd(padno++); fHistograms->At(khElectronpT) ->Draw(drawoption);
-    c->cd(padno++); fHistograms->At(khElectronPhi)->Draw(drawoption);
-    c->cd(padno++); fHistograms->At(khElectronEta)->Draw(drawoption);
-    c->Print(".png");
-  }
 }
 
 TObject* AliDxHFECorrelation::FindObject(const char *name) const
 {
   /// overloaded from TObject: find object by name
-  if (fHistograms) {
-    return fHistograms->FindObject(name);
+  if (fControlObjects) {
+    return fControlObjects->FindObject(name);
   }
   return NULL;
 }
@@ -196,13 +297,13 @@ TObject* AliDxHFECorrelation::FindObject(const char *name) const
 TObject* AliDxHFECorrelation::FindObject(const TObject *obj) const
 {
   /// overloaded from TObject: find object by pointer
-  if (fHistograms) {
-    return fHistograms->FindObject(obj);
+  if (fControlObjects) {
+    return fControlObjects->FindObject(obj);
   }
   return NULL;
 }
 
-void     AliDxHFECorrelation::SaveAs(const char *filename,Option_t */*option*/) const
+void AliDxHFECorrelation::SaveAs(const char *filename, Option_t */*option*/) const
 {
   /// overloaded from TObject: save to file
   std::auto_ptr<TFile> output(TFile::Open(filename, "RECREATE"));
@@ -211,13 +312,14 @@ void     AliDxHFECorrelation::SaveAs(const char *filename,Option_t */*option*/)
     return;
   }
   output->cd();
-  if (fHistograms) fHistograms->Write();
+  if (fControlObjects) fControlObjects->Write();
   output->Close();
 }
 
 AliDxHFECorrelation& AliDxHFECorrelation::operator+=(const AliDxHFECorrelation& other)
 {
   /// add histograms from another instance
+  // TODO - need to change this to ThnSparse?
   if (!fHistograms || !other.fHistograms) return *this;
   
   for (int i=0; i<kNofHistograms; i++) {
index 74efd33d76471ad3c9a2c3e28ff6e3d5f752ce79..8ddd5131c85c37ac427d159fe30c9b52ec0c33ff 100644 (file)
 
 #include "TNamed.h"
 
-class TH1F;
-class TH2F;
+class AliRDHFCutsD0toKpi;
+class TH1;
+class THnSparse;
+class TObject;
+class TList;
 
 class AliDxHFECorrelation : public TNamed {
  public:
@@ -26,12 +29,24 @@ class AliDxHFECorrelation : public TNamed {
   /// destructor
   virtual ~AliDxHFECorrelation();
 
+  // event control histogram
+  enum {
+    kEventsAll = 0, // all events
+    kEventsSel,     // selected events
+    kEventsD0 ,     // events with D0s
+    kEventsD0e,     // events with correlated D0s
+    kNEventControlLabels
+  };
+
   // init
   int Init();
 
   /// fill histograms from particles
   int Fill(const TObjArray* candidatesD0, const TObjArray* candidatesElectron);
 
+  /// histogram event properties
+  virtual int HistogramEventProperties(int bin);
+
   /// overloaded from TObject: cleanup
   virtual void Clear(Option_t * option ="");
   /// overloaded from TObject: print info
@@ -45,8 +60,18 @@ 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 SetUseMC(Bool_t useMC){fUseMC=useMC;}
+
+  Bool_t GetUseMC() const {return fUseMC;}
+  const TList* GetControlObjects() const {return fControlObjects;}
+
+
   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 {
     khD0pT,         // TH1F
     khD0Phi,        // TH1F
@@ -54,19 +79,29 @@ class AliDxHFECorrelation : public TNamed {
     khElectronpT,   // TH1F
     khElectronPhi,  // TH1F
     khElectronEta,  // TH1F
-    khDeltaPhi,        // TH1F
     kNofHistograms
   };
 
  protected:
+  /// add control object to list, the base class becomes owner of the object
+  int AddControlObject(TObject* pObj);
+
  private:
   /// copy constructor
   AliDxHFECorrelation(const AliDxHFECorrelation& other);
   /// assignment operator
   AliDxHFECorrelation& operator=(const AliDxHFECorrelation& other);
 
-  TObjArray* fHistograms; // the histograms
+  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
+  AliRDHFCutsD0toKpi *fCuts;  //  Cuts 
+  Bool_t fUseMC;              // use MC info
+
+  static const char* fgkEventControlBinNames[];
+  static const char* fgkCorrControlBinNames[];
 
-  ClassDef(AliDxHFECorrelation, 1)
+  ClassDef(AliDxHFECorrelation, 2)
 };
 #endif
index 899da3c1cd773fc78e87f09d42cbba94b0349f0f..2a5201ca1a24947c6264362c3a4ecef1c6e61748 100644 (file)
@@ -49,6 +49,8 @@ AliDxHFEParticleSelection::AliDxHFEParticleSelection(const char* name, const cha
   , fControlObjects(NULL)
   , fhEventControl(NULL)
   , fhTrackControl(NULL)
+  , fUseMC(false)
+  , fVerbosity(0)
 {
   // constructor
   // 
@@ -57,6 +59,12 @@ AliDxHFEParticleSelection::AliDxHFEParticleSelection(const char* name, const cha
   // 
 }
 
+const char* AliDxHFEParticleSelection::fgkEventControlBinNames[]={
+  "nEventsAll",
+  "nEventsSelected",
+  "nEventsD0"
+};
+
 AliDxHFEParticleSelection::~AliDxHFEParticleSelection()
 {
   // destructor
@@ -64,16 +72,25 @@ AliDxHFEParticleSelection::~AliDxHFEParticleSelection()
   fSelectedTracks=NULL;
   if (fControlObjects) delete fControlObjects;
   fControlObjects=NULL;
+  fhEventControl=NULL;
+  fhTrackControl=NULL;
 }
 
 int AliDxHFEParticleSelection::InitControlObjects()
 {
+  // init control objects
+  if (fVerbosity>0) {
+    AliInfo("Setting up control objects");
+  }
+
   /// init the control objects, can be overloaded by childs which should
   /// call AliDxHFEParticleSelection::InitControlObjects() explicitly
   std::auto_ptr<TH1D> hEventControl(new TH1D("hEventControl", "hEventControl", 10, 0, 10));
   std::auto_ptr<TH1D> hTrackControl(new TH1D("hTrackControl", "hTrackControl", 10, 0, 10));
 
   fhEventControl=hEventControl.release();
+  for (int iLabel=0; iLabel<kNEventPropertyLabels; iLabel++)
+    fhEventControl->GetXaxis()->SetBinLabel(iLabel, fgkEventControlBinNames[iLabel]);
   AddControlObject(fhEventControl);
   fhTrackControl=hTrackControl.release();
   AddControlObject(fhTrackControl);
@@ -94,11 +111,25 @@ int AliDxHFEParticleSelection::AddControlObject(TObject* pObj)
     AliError(Form("ignoring duplicate object '%s' of type %s", pObj->GetName(), pObj->ClassName()));
     return -EEXIST;
   }
+  if (GetVerbosity()>0) {
+    AliInfo(Form("Adding object '%s' of type %s",pObj->GetName(),pObj->ClassName()));
+  }
   fControlObjects->Add(pObj);
   return 0;
 }
 
-int AliDxHFEParticleSelection::HistogramParticleProperties(AliVParticle* p, bool selected)
+int AliDxHFEParticleSelection::HistogramEventProperties(int bin)
+{
+  /// histogram event properties
+  if (!fControlObjects) return 0;
+
+  // TODO: use enums for the bins of the control histogram
+  // for now: 0=all, 1=events with D0s, 2=events with correlated D0s
+  fhEventControl->Fill(bin);
+  return 0;
+}
+
+int AliDxHFEParticleSelection::HistogramParticleProperties(AliVParticle* p, int selected)
 {
   /// histogram particle properties
   if (!p) return -EINVAL;
@@ -112,7 +143,8 @@ int AliDxHFEParticleSelection::HistogramParticleProperties(AliVParticle* p, bool
 
 TObjArray* AliDxHFEParticleSelection::Select(const AliVEvent* pEvent)
 {
-  /// create selection, array contains only pointers but does not own the objects
+  /// create selection from 'Tracks' member of the event,
+  /// array contains only pointers but does not own the objects
   /// object array needs to be deleted by caller
   if (!pEvent) return NULL;
   TObjArray* selectedTracks=new TObjArray;
@@ -120,29 +152,30 @@ TObjArray* AliDxHFEParticleSelection::Select(const AliVEvent* pEvent)
   int nofTracks=pEvent->GetNumberOfTracks();
   for (int itrack=0; itrack<nofTracks; itrack++) {
     AliVParticle* track=pEvent->GetTrack(itrack);
-    bool selected=IsSelected(track);
-    HistogramParticleProperties(track, selected);
-    if (!selected) continue;
+    int selectionCode=IsSelected(track);
+    HistogramParticleProperties(track, selectionCode);
+    if (selectionCode==0) continue;
     selectedTracks->Add(track);
   }
   return selectedTracks;
 }
 
-TObjArray* AliDxHFEParticleSelection::Select(TObjArray* pTracks)
+TObjArray* AliDxHFEParticleSelection::Select(TObjArray* pParticles, const AliVEvent* pEvent)
 {
-  /// create selection, array contains only pointers but does not own the objects
+  /// create selection from the array of particles,
+  /// array contains only pointers but does not own the objects
   /// object array needs to be deleted by caller
-  if (!pTracks) return NULL;
+  if (!pParticles) return NULL;
   TObjArray* selectedTracks=new TObjArray;
   if (!selectedTracks) return NULL;
-  TIter itrack(pTracks);
+  TIter next(pParticles);
   TObject* pObj=NULL;
-  while ((pObj=itrack())!=NULL) {
+  while ((pObj=next())) {
     AliVParticle* track=dynamic_cast<AliVParticle*>(pObj);
     if (!track) continue;
-    bool selected=IsSelected(track);
-    HistogramParticleProperties(track, selected);
-    if (!selected) continue;
+    int selectionCode=IsSelected(track, pEvent);
+    HistogramParticleProperties(track, selectionCode);
+    if (selectionCode ==0) continue;
     selectedTracks->Add(track);
   }
   return selectedTracks;
@@ -155,11 +188,11 @@ int AliDxHFEParticleSelection::CheckAndAdd(AliVParticle* /*p*/)
   return -ENOSYS;
 }
 
-bool AliDxHFEParticleSelection::IsSelected(AliVParticle* /*p*/)
+int AliDxHFEParticleSelection::IsSelected(AliVParticle* /*p*/, const AliVEvent* /*e*/)
 {
   /// check particle if it passes the selection criteria
   /// childs can overload, by default all tracks are selected
-  return true;
+  return 1;
 }
 
 void AliDxHFEParticleSelection::AliDxHFEParticleSelection::Clear(Option_t * /*option*/)
@@ -175,10 +208,15 @@ void AliDxHFEParticleSelection::Print(Option_t */*option*/) const
   if (fControlObjects) fControlObjects->Print();
 }
  
-void AliDxHFEParticleSelection::SaveAs(const char* filename,Option_t */*option*/) const
+void AliDxHFEParticleSelection::SaveAs(const char* filename, Option_t */*option*/) const
 {
   /// inherited from TObject: save selection criteria
-  std::auto_ptr<TFile> output(TFile::Open(filename, "RECREATE"));
+  TString fileoption;
+  // TODO: options recreate
+  fileoption="RECREATE";
+  //else fileoption="UPDATE";
+
+  std::auto_ptr<TFile> output(TFile::Open(filename,fileoption));
   if (!output.get() || output->IsZombie()) {
     AliError(Form("can not open file %s from writing", filename));
     return;
index c2fd51598003dc9ab739bdbdce72b3d42d761920..c6712e7e11defc85d01973909e261dc6842c0e3f 100644 (file)
@@ -42,6 +42,13 @@ class AliDxHFEParticleSelection : public TNamed {
   /// destructor
   virtual ~AliDxHFEParticleSelection();
 
+  enum {
+    kEventsAll = 0,
+    kEventsSel,
+    kEventsD0,
+    kNEventPropertyLabels
+  };
+
   /// set options
   void SetOption(const char* opt) { fOption = opt; }
   /// overloaded from TObject: get option
@@ -50,28 +57,42 @@ class AliDxHFEParticleSelection : public TNamed {
   /// init the control objects
   virtual int InitControlObjects();
 
-  /// create selection, array contains only pointers but does not own the objects
+  /// create selection from 'Tracks' member of the event,
+  /// array contains only pointers but does not own the objects
   /// object array needs to be deleted by caller
   virtual TObjArray* Select(const AliVEvent* pEvent);
-
-  /// create selection, array contains only pointers but does not own the objects
+  /// create selection from the array of particles,
+  /// array contains only pointers but does not own the objects
   /// object array needs to be deleted by caller
-  virtual TObjArray* Select(TObjArray* pTracks);
+  virtual TObjArray* Select(TObjArray* particles, const AliVEvent* pEvent);
+
+  // Get the list fControlObjects. 
+  const TList* GetControlObjects() const {return fControlObjects;}
+
+  /// histogram event properties
+  virtual int HistogramEventProperties(int bin);
 
   /// check and add track to internal array
   int CheckAndAdd(AliVParticle* p);
+
+  /// set cuts object: general TObject pointer is used as argument to support
+  // different types; a type cast check is implemented in the method
+  virtual void SetCuts(TObject* /*cuts*/, int /*level*/=0) {}
+
+  Bool_t GetUseMC() const {return fUseMC;}
+
   /// get selected tracks
   const TObjArray* GetSelected() const {return fSelectedTracks;}
 
   /// check particle if it passes the selection criteria
-  virtual bool IsSelected(AliVParticle* p);
+  virtual int IsSelected(AliVParticle* p, const AliVEvent *pEvent=NULL);
 
   /// inherited from TObject: cleanup
   virtual void Clear(Option_t * option ="");
   /// inherited from TObject: print info
   virtual void Print(Option_t *option="") const;
   /// inherited from TObject: safe selection criteria
-  virtual void SaveAs(const char *filename="",Option_t *option="") const;
+  virtual void SaveAs(const char *filename="", Option_t *option="") const;
   /// inherited from TObject: draw content
   virtual void Draw(Option_t *option="");
   /// inherited from TObject: find object by name
@@ -79,12 +100,23 @@ class AliDxHFEParticleSelection : public TNamed {
   /// inherited from TObject: find object by pointer
   virtual TObject* FindObject(const TObject *obj) const;
 
+  /// set verbosity
+  void SetVerbosity(int verbosity) {fVerbosity=verbosity;}
+  /// get verbosity
+  int GetVerbosity() const {return fVerbosity;}
+
  protected:
   /// add control object to list, the base class becomes owner of the object
   int AddControlObject(TObject* pObj);
 
   /// histogram particle properties
-  virtual int HistogramParticleProperties(AliVParticle* p, bool selected=true);
+  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
@@ -99,8 +131,13 @@ class AliDxHFEParticleSelection : public TNamed {
   TList* fControlObjects; // list of control objects
   TH1* fhEventControl; //! event control histogram
   TH1* fhTrackControl; //! track control histogram
+  bool fUseMC;         // specific implementation for MC selection
+  int fVerbosity;      //! verbosity
+
+  static const char* fgkEventControlBinNames[]; //! bin labels for event control histogram
+
+  ClassDef(AliDxHFEParticleSelection, 2);
 
-  ClassDef(AliDxHFEParticleSelection, 1);
 };
 
 #endif
index d09c62faf0eccd68cc84f603bcbd4ac2ae80a66c..c6c36bd142f0aee3eed6779fcb70eacc0681f862 100644 (file)
 
 #include "AliDxHFEParticleSelectionD0.h"
 #include "AliVParticle.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 "AliRDHFCutsD0toKpi.h"
 #include "TObjArray.h"
 #include "THnSparse.h"
 #include "TAxis.h"
+#include "TString.h"
 #include <iostream>
 #include <cerrno>
 #include <memory>
 
+using namespace std;
+
 /// ROOT macro for the implementation of ROOT specific class methods
 ClassImp(AliDxHFEParticleSelectionD0)
 
 AliDxHFEParticleSelectionD0::AliDxHFEParticleSelectionD0(const char* opt)
   : AliDxHFEParticleSelection("D0", opt)
   , fD0Properties(NULL)
+  , fD0Daughter0(NULL)
+  , fD0Daughter1(NULL)
+  , fCuts(NULL)
+  , fFillOnlyD0D0bar(0)
 {
   // constructor
   // 
   // 
   // 
   // 
+  TString strOption(opt);
+  // TODO: one might need a proper argument parsing including
+  // chopping whole string into individual arguments
+  if (strOption.Contains("FillD0D0bar")) fFillOnlyD0D0bar=0;
+  else if (strOption.Contains("FillOnlyD0")) fFillOnlyD0D0bar=1;
+  else if (strOption.Contains("FillOnlyD0bar")) fFillOnlyD0D0bar=2;
 }
 
 AliDxHFEParticleSelectionD0::~AliDxHFEParticleSelectionD0()
 {
   // destructor
+  if (fD0Properties) {
+    delete fD0Properties;
+    fD0Properties=NULL;
+  }
+  if (fD0Daughter0) {
+    delete fD0Daughter0;
+    fD0Daughter0=NULL;
+  }
+  if (fD0Daughter1) {
+    delete fD0Daughter1;
+    fD0Daughter1=NULL;
+  }
+
+  // Note: external object deleted elsewhere  
+  fCuts=NULL;
 }
 
 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: very specific D0 for the moment, sort out later
   // 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};
+  //                            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));
@@ -75,7 +108,7 @@ int AliDxHFEParticleSelectionD0::InitControlObjects()
     return -ENOMEM;
   }
   int axis=0;
-  D0Properties->GetAxis(axis++)->SetTitle("mass");
+  D0Properties->GetAxis(axis++)->SetTitle("D0 inv mass");
   D0Properties->GetAxis(axis++)->SetTitle("Pt");
   D0Properties->GetAxis(axis++)->SetTitle("Phi"); 
   D0Properties->GetAxis(axis++)->SetTitle("Ptbin"); 
@@ -83,35 +116,169 @@ int AliDxHFEParticleSelectionD0::InitControlObjects()
   fD0Properties=D0Properties.release();
   AddControlObject(fD0Properties);
 
+  //Adding control objects for the daughters
+  InitControlObjectsDaughters("pi information",0);
+  InitControlObjectsDaughters("K information",1);
+
   return AliDxHFEParticleSelection::InitControlObjects();
 }
 
-int AliDxHFEParticleSelectionD0::HistogramParticleProperties(AliVParticle* p, bool selected)
+int AliDxHFEParticleSelectionD0::InitControlObjectsDaughters(TString name, int daughter)
+{
+  //Setting up Control objects for the daughters.
+  AliInfo("Setting up daughter THnSparse");
+
+  const int thnSize2 = 5;
+  const double Pi=TMath::Pi();
+  //                          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.};
+
+  std::auto_ptr<THnSparseF> DaughterProperties(new THnSparseF(name, name, thnSize2, thnBins, thnMin, thnMax));
+
+  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"); 
+
+  if(daughter==0){ 
+    fD0Daughter0=DaughterProperties.release();
+    AddControlObject(fD0Daughter0);
+  }
+  
+  if(daughter==1){
+    fD0Daughter1=DaughterProperties.release();
+    AddControlObject(fD0Daughter1);
+  }
+  return 0;
+}
+
+int AliDxHFEParticleSelectionD0::HistogramParticleProperties(AliVParticle* p, int selectionCode)
 {
+
   /// histogram particle properties
   if (!p) return -EINVAL;
 
   // fill the common histograms
-  AliDxHFEParticleSelection::HistogramParticleProperties(p, selected);
+  AliDxHFEParticleSelection::HistogramParticleProperties(p, selectionCode);
 
-  // TODO: histograms for all and selected particles
-  if (!selected) return 0;
+  // no daughters to fill if 0 (= no candidate)
+  if (selectionCode==0) return 0;
 
-  // TODO: find out which type is necessary
   AliAODRecoDecayHF2Prong* part=dynamic_cast<AliAODRecoDecayHF2Prong*>(p);
-  if (part) {
+
+  if(!part) return 0;
+  // Convention: 1. daughter is postive track, 2. = negative
+  AliAODTrack *prongpos=(AliAODTrack*)part->GetDaughter(0);
+  AliAODTrack *prongneg=(AliAODTrack*)part->GetDaughter(1);
+
+  if(!prongpos || !prongneg) {
+    return 0;
+  }
+  // Only D0s are filled 
+  // TODO: Also include D0bar
+  if ((selectionCode==1 || selectionCode==3) && fFillOnlyD0D0bar<2) {
     Double_t invmassD0 = part->InvMassD0();
-    // TODO: use cut object to define pt bin
-    Int_t ptbin=0;//cuts->PtBin(part->Pt());
-    Double_t D0Stuff[] = {invmassD0,part->Pt(),part->Phi(),ptbin};
-    if (fD0Properties) fD0Properties->Fill(D0Stuff);
+    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);
+    if(fD0Daughter0) fD0Daughter0->Fill(piProperties);
+    if(fD0Daughter1) fD0Daughter1->Fill(KProperties);
   }
 
   return 0;
 }
 
-bool AliDxHFEParticleSelectionD0::IsSelected(AliVParticle* /*p*/)
+TObjArray* AliDxHFEParticleSelectionD0::Select(TObjArray* pTracks, const AliVEvent *pEvent)
+{
+  /// create selection, array contains only pointers but does not own the objects
+  /// object array needs to be deleted by caller
+  if (!pTracks) return NULL;
+  TObjArray* selectedTracks=new TObjArray;
+  if (!selectedTracks) return NULL;
+  TIter itrack(pTracks);
+  TObject* pObj=NULL;
+  while ((pObj=itrack())!=NULL) {
+    AliVParticle* track=dynamic_cast<AliVParticle*>(pObj);
+    if (!track) continue;
+    int selectionCode=IsSelected(track,pEvent);
+    HistogramParticleProperties(track, selectionCode);
+    //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;
+    selectedTracks->Add(track);
+  }
+  return selectedTracks;
+}
+
+int AliDxHFEParticleSelectionD0::IsSelected(AliVParticle* p, const AliVEvent* pEvent)
 {
   /// TODO: implement specific selection of D0 candidates
-  return true;
+  /// Could also return values based on where where selection "failed"
+  int selectionCode=0;
+
+  AliAODRecoDecayHF2Prong *d0 = dynamic_cast<AliAODRecoDecayHF2Prong*>(p);
+  if(d0->GetSelectionMap()) if(!d0->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
+      AliDebug(1,"Skip D0 from Dstar");
+      return 0; //skip the D0 from Dstar
+    }
+
+  // TODO: the cuts instance should be const but the function definition of
+  // AliRDHFCuts::IsSelected does not allow this
+  AliRDHFCuts* cuts=const_cast<AliRDHFCuts*>(fCuts);
+  if (!cuts) {
+    selectionCode=1;
+  } 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");
+      return 0;
+    } //out of bounds
+
+    // TODO: the aod pointer should also be const but the function definition of
+    // AliRDHFCuts::IsSelected does not allow this
+    AliAODEvent* aod=NULL;
+    if (pEvent) aod=dynamic_cast<AliAODEvent*>(const_cast<AliVEvent*>(pEvent));
+  
+    // 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));
+    TObjArray daughters;
+    daughters.AddAt((AliAODTrack*)d0->GetDaughter(0),0);
+    daughters.AddAt((AliAODTrack*)d0->GetDaughter(1),1);
+
+    //check daughters
+    if(!daughters.UncheckedAt(0) || !daughters.UncheckedAt(1)) {
+      AliDebug(1,"at least one daughter not found!");
+      daughters.Clear();
+      return 0;
+    }
+  }
+
+  return selectionCode;
+}
+
+void AliDxHFEParticleSelectionD0::SetCuts(TObject* cuts, int /*level*/)
+{
+  /// set cuts objects
+  fCuts=dynamic_cast<AliRDHFCuts*>(cuts);
+  if (!fCuts && cuts) {
+    AliError(Form("cuts object is not of required type AliRDHFCuts but %s", cuts->ClassName()));
+  }
 }
index fc80e82b338b918d1a941be75e7f76775135fa01..93c43f151a37187e44c906bf722253c386ed1919 100644 (file)
@@ -16,6 +16,9 @@
 
 #include "AliDxHFEParticleSelection.h"
 
+class AliVEvent;
+class AliRDHFCuts;
+
 /**
  * @class AliDxHFEParticleSelectionD0
  * D0 selection for D-HFE correlations, implements the specific
@@ -31,12 +34,25 @@ class AliDxHFEParticleSelectionD0 : public AliDxHFEParticleSelection {
   /// overloaded from AliDxHFEParticleSelection: init the control objects
   virtual int InitControlObjects();
 
+  //Function for daughter control objects
+  //TODO: move to AliDxHFEParticleSelection to be used for several particles?
+  virtual int InitControlObjectsDaughters(TString name, int daughter);
+
+  //Overloaded from AliDxHFEParticleSelection
+  virtual TObjArray* Select(TObjArray* particles, const AliVEvent* pEvent);
+  using AliDxHFEParticleSelection::Select;
+
   /// overloaded from AliDxHFEParticleSelection: check particle
-  virtual bool IsSelected(AliVParticle* p);
+  virtual int IsSelected(AliVParticle* p, const AliVEvent *pEvent=NULL);
+  /// overloaded from AliDxHFEParticleSelection: set cuts
+  virtual void SetCuts(TObject* /*cuts*/, int level=0);
+
+  //AliRDHFCutsD0toKpi GetCuts()const {return fCuts;}
+  Int_t  GetFillOnlyD0D0bar() const {return fFillOnlyD0D0bar;}
 
  protected:
   /// overloaded from AliDxHFEParticleSelection: histogram particle properties
-  virtual int HistogramParticleProperties(AliVParticle* p, bool selected=true);
+  virtual int HistogramParticleProperties(AliVParticle* p, int selected=1);
 
  private:
   /// copy contructor prohibited
@@ -44,12 +60,17 @@ class AliDxHFEParticleSelectionD0 : public AliDxHFEParticleSelection {
   /// assignment operator prohibited
   AliDxHFEParticleSelectionD0& operator=(const AliDxHFEParticleSelectionD0&);
 
-  THnSparse* fD0Properties; //! the particle properties of selected particles
+  THnSparse* fD0Properties;     //  the particle properties of selected particles
+  THnSparse* fD0Daughter0;      //  the particle properties of selected particles
+  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)
+
   // TODO: at the moment the dimensions of the different THnSparse objects are different
   // needs to be consolidated
   // TODO: one might need particle properties of all and/or at different cut stages
 
-  ClassDef(AliDxHFEParticleSelectionD0, 1);
+  ClassDef(AliDxHFEParticleSelectionD0, 2);
 };
 
 #endif
index ed136cf63590ae6ccce8c43c57cfb6dc84eee291..0835d9740f6ec68bab8123e9b5055c7d17fadf02 100644 (file)
 /// @date   2012-03-19
 /// @brief  D0 selection for D0-HFE correlation
 ///
-
 #include "AliDxHFEParticleSelectionEl.h"
 #include "AliVParticle.h"
+#include "AliVEvent.h"
+#include "AliPID.h"
+#include "AliPIDResponse.h"
+#include "AliHFEcontainer.h"
+#include "AliHFEpid.h"
+#include "AliHFEpidBase.h"
+#include "AliHFEtools.h"
+#include "AliHFEcuts.h"
+#include "AliAODTrack.h"
+#include "AliAnalysisDataSlot.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliAnalysisManager.h"
+#include "AliCFManager.h"
+#include "THnSparse.h"
+#include "TH1F.h"
+#include "TAxis.h"
 #include "TObjArray.h"
+#include <iostream>
+#include <cerrno>
+#include <memory>
+
+using namespace std;
 
 /// ROOT macro for the implementation of ROOT specific class methods
 ClassImp(AliDxHFEParticleSelectionEl)
 
 AliDxHFEParticleSelectionEl::AliDxHFEParticleSelectionEl(const char* opt)
-  : AliDxHFEParticleSelection(opt)
+  : AliDxHFEParticleSelection("electron", opt)
+  , fPID(NULL)
+  , fElectronProperties(NULL)
+  , fWhichCut(NULL)  
+  , fCuts(NULL)
+  , fCFM(NULL)
 {
   // constructor
   // 
@@ -43,10 +68,179 @@ AliDxHFEParticleSelectionEl::AliDxHFEParticleSelectionEl(const char* opt)
 AliDxHFEParticleSelectionEl::~AliDxHFEParticleSelectionEl()
 {
   // destructor
+  if (fElectronProperties) {
+    delete fElectronProperties;
+    fElectronProperties=NULL;
+  }
+  if(fCFM){
+    delete fCFM;
+    fCFM=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));
+
+  if (ElectronProperties.get()==NULL) {
+    return -ENOMEM;
+  }
+  int axis=0;
+  ElectronProperties->GetAxis(axis++)->SetTitle("Pt");
+  ElectronProperties->GetAxis(axis++)->SetTitle("Phi");
+  ElectronProperties->GetAxis(axis++)->SetTitle("Eta"); 
+
+  fElectronProperties=ElectronProperties.release();
+
+  AddControlObject(fElectronProperties);
+
+  fWhichCut= new TH1F("fWhichCut","effective cut for a rejected particle",6,-0.5,5.5);
+  AddControlObject(fWhichCut);
+
+  //--------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
+  AliInfo("Setting up CFM");
+  fCFM = new AliCFManager;
+  // the setup of cut objects is done in AliHFEcuts::Initialize
+  // the ids used by this class must be the same, the code might be broken if
+  // the sequence in AliHFEcuts::Initialize is changed
+  const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack;
+  // reset pointers in the CF manager
+  fCFM->SetNStepParticle(kNcutSteps);
+  for(Int_t istep = 0; istep < kNcutSteps; istep++) {
+    fCFM->SetParticleCutsList(istep, NULL);
+  }
+  if(!fCuts) {
+    AliWarning("Cuts not available. Default cuts will be used");
+    fCuts = new AliHFEcuts;
+    fCuts->CreateStandardCuts();
+  }
+  // TODO: error handling?
+  fCuts->Initialize(fCFM);
+
+  return 0;
+}
+
+int AliDxHFEParticleSelectionEl::HistogramParticleProperties(AliVParticle* p, int selectionCode)
+{
+  /// histogram particle properties
+  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()};
+  if(fElectronProperties) fElectronProperties->Fill(eProperties);
+  
+  return 0;
+}
+
+int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent*)
+{
+  /// select El candidates
+  AliAODTrack *track=(AliAODTrack*)pEl;
+
+  //--------track cut selection-----------------------
+  //Using AliHFECuts:
+  // RecKine: ITSTPC cuts  
+  if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)){
+    fWhichCut->Fill(0);
+    return 0;
+  }
+  
+  // RecPrim
+  if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) {
+    fWhichCut->Fill(1);
+    return 0;
+  }
+  
+  // HFEcuts: ITS layers cuts
+  if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) {
+    fWhichCut->Fill(2);
+    return 0;
+  }
+  
+  // HFE cuts: TOF PID and mismatch flag
+  if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) {
+    fWhichCut->Fill(3);
+    return 0;
+  }
+  
+  // HFE cuts: TPC PID cleanup
+  if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)){
+    fWhichCut->Fill(4);
+    return 0;
+  } 
+  // HFEcuts: Nb of tracklets TRD0
+  //if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
+
+
+  //--------PID selection-----------------------
+  AliHFEpidObject hfetrack;
+  hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
+  hfetrack.SetRecTrack(track);
+
+  // TODO: configurable colliding system
+  //if(IsPbPb()) hfetrack.SetPbPb();
+  hfetrack.SetPP();
+
+  if(fPID && fPID->IsSelected(&hfetrack)) {
+    AliDebug(3,"Inside FilldPhi, electron is selected");
+
+    return 1;
+  }
+  else{
+    fWhichCut->Fill(5);
+    return 0;
+  }
+}
+
+void AliDxHFEParticleSelectionEl::SetCuts(TObject* cuts, int level)
+{
+  /// set cut objects
+  if (level==kCutHFE) {
+    fCuts=dynamic_cast<AliHFEcuts*>(cuts);
+    if (!fCuts && cuts) {
+      AliError(Form("Cut object is not of required type AliHFEcuts but %s", cuts->ClassName()));
+    }
+    return;
+  }
+
+  if (level==kCutPID) {
+    fPID=dynamic_cast<AliHFEpid*>(cuts);
+    if (!fPID && cuts) {
+      AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
+    }
+    return;
+  }
 }
 
-bool AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* /*p*/)
+//________________________________________________________________________
+Bool_t AliDxHFEParticleSelectionEl::ProcessCutStep(Int_t cutStep, AliVParticle *track)
 {
-  /// TODO: implement specific selection of D0 candidates
-  return true;
+  // Check single track cuts for a given cut step
+  const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
+  if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
+  //if(!fCuts->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
+  return kTRUE;
 }
index 5c3b627464b13bafeb1672889766a41c869c1c37..e3a91b665cb0609d6ad3c9e9b85bec547e5165e8 100644 (file)
 #define ALIDXHFEPARTICLESELECTIONEL_H
 
 #include "AliDxHFEParticleSelection.h"
+class AliPID;
+class AliPIDResponse;
+class AliHFEcuts;
+class AliHFEvarManager;
+class AliHFEpid;
+class AliHFEpidBase;
+class AliHFEtools;
+class AliVEvent;
+class AliCFManager;
+class TH1;
 
 /**
  * @class AliDxHFEParticleSelectionEl
@@ -28,8 +38,25 @@ class AliDxHFEParticleSelectionEl : public AliDxHFEParticleSelection {
   /// destructor
   virtual ~AliDxHFEParticleSelectionEl();
 
+  enum {
+    kCutHFE = 0,
+    kCutPID = 1,
+    kNCuts
+  };
+
+  /// overloaded from AliDxHFEParticleSelection: init the control objects
+  virtual int InitControlObjects();
+
   /// overloaded from AliDxHFEParticleSelection: check particle
-  virtual bool IsSelected(AliVParticle* p);
+  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);
+
+  /// set cuts object: a type cast check is implemented in the method
+  virtual void SetCuts(TObject* /*cuts*/, int /*level*/=0);
 
  protected:
 
@@ -39,7 +66,17 @@ class AliDxHFEParticleSelectionEl : public AliDxHFEParticleSelection {
   /// assignment operator prohibited
   AliDxHFEParticleSelectionEl& operator=(const AliDxHFEParticleSelectionEl&);
 
-  ClassDef(AliDxHFEParticleSelectionEl, 1);
+  /// check cut of specified step, e.g.
+  bool ProcessCutStep(Int_t cutStep, AliVParticle *track);
+
+  AliHFEpid*    fPID;                //! the PID object
+  THnSparse*    fElectronProperties; // the particle properties of selected particles
+  TH1*          fWhichCut;           // effective cut for a rejected particle
+  AliHFEcuts*   fCuts;               //! Cuts for HF electrons
+  AliCFManager* fCFM;                //! Correction Framework Manager
+
+
+  ClassDef(AliDxHFEParticleSelectionEl, 2);
 };
 
 #endif
index 29cddf18e0b1d0ee2a9711c5b96df20ce219719d..e0e1e030cda4277c19e807dba3a876cec3c11cab 100644 (file)
 #if defined(__CINT__) && !defined(__MAKECINT__)
 {
   //----------- Loading the required libraries ---------//
+  gInterpreter->ExecuteMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
 
   gSystem->Load("libSTEERBase");
   gSystem->Load("libESD");
   gSystem->Load("libAOD");
   gSystem->Load("libANALYSIS");
   gSystem->Load("libANALYSISalice");
-  gSystem->AddIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include");
+  gSystem->Load("libPWGHFhfe.so");
+  gSystem->Load("libCORRFW");
+  gSystem->AddIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT/PWGHF/vertexingHF -I$ALICE_ROOT/PWGHF/base -I$ALICE_ROOT/PWGHF/hfe ");
   gROOT->LoadMacro("$ALICE_ROOT/PWGHF/correlationHF/AliDxHFEParticleSelection.cxx+");
   gROOT->LoadMacro("$ALICE_ROOT/PWGHF/correlationHF/AliDxHFEParticleSelectionD0.cxx+");
   gROOT->LoadMacro("$ALICE_ROOT/PWGHF/correlationHF/AliDxHFEParticleSelectionEl.cxx+");