]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/correlationHF/AliAnalysisTaskDxHFECorrelation.cxx
Filling histograms separately for the SPD and Tracks vertices
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskDxHFECorrelation.cxx
index 005ae4beb974a5a088489025b629e4925e55ddec..519210ce303cb0e9479effa9ebafa1bff28231ec 100644 (file)
 
 #include "AliAnalysisTaskDxHFECorrelation.h"
 #include "AliDxHFECorrelation.h"
+#include "AliDxHFECorrelationMC.h"
 #include "AliDxHFEParticleSelectionD0.h"
 #include "AliDxHFEParticleSelectionMCD0.h"
 #include "AliDxHFEParticleSelectionEl.h"
 #include "AliDxHFEParticleSelectionMCEl.h"
 #include "AliDxHFEParticleSelection.h"
+#include "AliHFCorrelator.h"
 #include "AliAnalysisManager.h"
 #include "AliLog.h"
 #include "AliESDInputHandler.h"
@@ -57,6 +59,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>
 
@@ -73,17 +77,22 @@ AliAnalysisTaskDxHFECorrelation::AliAnalysisTaskDxHFECorrelation(const char* opt
   , fD0s(NULL)
   , fElectrons(NULL)
   , fCutsD0(NULL)
-  , fCutsHFE(NULL)
-  , fPID(NULL)
-  , fFillOnlyD0D0bar(0)
+  , fCuts(NULL)
   , fUseMC(kFALSE)
+  , fUseEventMixing(kFALSE)
+  , fSystem(0)
+  , fSelectedD0s(NULL)
+  , fSelectedElectrons(NULL)
+  , fListHFE(NULL)
+  , fTriggerParticle(AliDxHFECorrelation::kD)
+  , fUseKine(kFALSE)
+  , fMCArray(NULL)
+  , fCorrelationArguments("")
 {
   // constructor
   //
   //
   DefineSlots();
-  fPID = new AliHFEpid("hfePid");
-
 }
 
 int AliAnalysisTaskDxHFECorrelation::DefineSlots()
@@ -92,6 +101,8 @@ int AliAnalysisTaskDxHFECorrelation::DefineSlots()
   DefineInput(0, TChain::Class());
   DefineOutput(1, TList::Class());
   DefineOutput(2,AliRDHFCutsD0toKpi::Class());
+  DefineOutput(3,TList::Class());
+  DefineOutput(4,AliHFAssociatedTrackCuts::Class());
   return 0;
 }
 
@@ -117,9 +128,16 @@ AliAnalysisTaskDxHFECorrelation::~AliAnalysisTaskDxHFECorrelation()
   // external object, do not delete
   fCutsD0=NULL;
   // external object, do not delete
-  fCutsHFE=NULL;
-  if(fPID) delete fPID;
-  fPID=NULL;
+  // TODO: also delete it here? this class is set as owner...
+  fListHFE=NULL;
+  if(fSelectedElectrons) delete fSelectedElectrons;
+  fSelectedElectrons=NULL;
+  if(fSelectedD0s) delete fSelectedD0s;
+  fSelectedD0s=NULL;
+  if(fListHFE) delete fListHFE;
+  fListHFE=NULL;
+  if(fMCArray) delete fMCArray;
+  fMCArray=NULL;
 
 
 }
@@ -127,48 +145,47 @@ AliAnalysisTaskDxHFECorrelation::~AliAnalysisTaskDxHFECorrelation()
 void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
 {
   // create result objects and add to output list
+  int iResult=0;
+  // ParseArguments will also define the strings that are used as input
+  // for the particle selection classes and the correlation class
 
-  //Initialize PID for electron selection
-  if(!fPID->GetNumberOfPIDdetectors()) { 
-    fPID->AddDetector("TOF",0);
-    fPID->AddDetector("TPC",1);
-  }
-  fPID->InitializePID();
+  ParseArguments(fOption.Data());
 
   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 ";
+  // D0s ===============================================
+  if(fUseMC) fD0s=new AliDxHFEParticleSelectionMCD0(fOption);
+  else fD0s=new AliDxHFEParticleSelectionD0(fOption);
+  fD0s->SetCuts(fCutsD0,AliDxHFEParticleSelectionD0::kCutD0);
+  iResult=fD0s->Init();
+  if (iResult<0) {
+    AliFatal(Form("initialization of worker class instance fD0s failed with error %d", iResult));
+  }
+
+  //Electrons ============================================
+  fListHFE->SetOwner(); // Not sure if needed
+  if(fUseMC) fElectrons=new AliDxHFEParticleSelectionMCEl(fOption);
+  else fElectrons=new AliDxHFEParticleSelectionEl(fOption);
+  fElectrons->SetCuts(fListHFE, AliDxHFEParticleSelectionEl::kCutList);
+  iResult=fElectrons->Init();
+  if (iResult<0) {
+    AliFatal(Form("initialization of worker class instance fElectrons failed with error %d", iResult));
   }
 
-  if(fUseMC) fD0s=new AliDxHFEParticleSelectionMCD0(selectionD0Options);
-  else fD0s=new AliDxHFEParticleSelectionD0(selectionD0Options);
-  fD0s->SetCuts(fCutsD0);
-  fD0s->Init();
-
-  //Electrons
-  if(fUseMC) fElectrons=new AliDxHFEParticleSelectionMCEl;
-  else fElectrons=new AliDxHFEParticleSelectionEl;
-  fElectrons->SetCuts(fPID, AliDxHFEParticleSelectionEl::kCutPID);
-  fElectrons->SetCuts(fCutsHFE, AliDxHFEParticleSelectionEl::kCutHFE);
-  fElectrons->Init();
-
-  //Correlation
-  fCorrelation=new AliDxHFECorrelation;
-  fCorrelation->SetCuts(dynamic_cast<AliRDHFCutsD0toKpi*>(fCutsD0));
-  // TODO: check if we can get rid of the mc flag in the correlation analysis class
-  fCorrelation->SetUseMC(fUseMC);
-  fCorrelation->Init();
+  //Correlation ===========================================
+  if(fUseMC) fCorrelation=new AliDxHFECorrelationMC;
+  else fCorrelation=new AliDxHFECorrelation;
+  fCorrelation->SetCuts(fCuts);
+  iResult=fCorrelation->Init(fOption);
+  if (iResult<0) {
+    AliFatal(Form("initialization of worker class instance fCorrelation failed with error %d", iResult));
+  }
 
   // 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 +193,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,15 +204,16 @@ void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
   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;
-    }
+  if (!fCutsD0) {
+    AliFatal(Form("cut object for D0 missing"));
+    return;
   }
-
-  // TODO: why copy? cleanup?
+  if (!dynamic_cast<AliRDHFCutsD0toKpi*>(fCutsD0)) {
+    AliFatal(Form("cut object %s is of incorrect type %s, expecting AliRDHFCutsD0toKpi", fCutsD0->GetName(), fCutsD0->ClassName()));
+    return;
+  }
+  // that's the copy for the output stream
   AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(dynamic_cast<AliRDHFCutsD0toKpi&>(*fCutsD0));
   const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
   copyfCuts->SetName(nameoutput);
@@ -203,6 +221,8 @@ void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
   // all tasks must post data once for all outputs
   PostData(1, fOutput);
   PostData(2,copyfCuts);
+  PostData(3,fListHFE);
+  PostData(4,fCuts);
 
 }
 
@@ -217,6 +237,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.    
@@ -245,87 +267,223 @@ 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
 
   if(!cutsd0->IsEventSelected(pEvent)) {
+    // TODO: Fill histograms based on why the event is rejected
     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());
+  if(fSystem==1){
+    // Not really used anywhere, just as a test of centralityselection
+    // (Could also be used in histo, so have not removed it)
+    AliCentrality *centralityObj = 0;
+    Double_t MultipOrCent = -1;
+    AliAODEvent* aodEvent=dynamic_cast<AliAODEvent*>(pEvent);
+    if (aodEvent) {
+      centralityObj = aodEvent->GetHeader()->GetCentralityP();
+      if (centralityObj) {
+       MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
+      }
+    }
+    AliInfo(Form("Centrality is %f", MultipOrCent));
   }
 
+  // Gets the PID response from the analysis manager
   AliPIDResponse *pidResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler())->GetPIDResponse();
   if(!pidResponse){
+    // TODO: consider issuing fatal instead of debug in case pidresponse not available
     AliDebug(1, "Using default PID Response");
     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class()); 
   }
-  fPID->SetPIDResponse(pidResponse);
+  
+  // Fetching the PID objects from the list, checks if the objects are AliHFEpids
+  // If so, checks if they are initialized and also sets the pidresponse
+  TObject *obj=NULL;
+  TIter next(fListHFE);
+  while((obj = next())){
+    AliHFEpid* pidObj=dynamic_cast<AliHFEpid*>(obj);
+    if(pidObj){
+      if(!pidObj->IsInitialized()){
+       AliWarning("PID not initialised, get from Run no");
+       pidObj->InitializePID(pEvent->GetRunNumber());
+      }
+      pidObj->SetPIDResponse(pidResponse);
+    }
+  }
 
+  // Also sends the pidresponse to the particle selection class for electron
+  fElectrons->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);
+    if(fUseKine){
+      fMCArray = dynamic_cast<TObjArray*>(pEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+      if(fUseMC && !fMCArray){
+       AliError("Array of MC particles not found");
+       return;
+      }
+    }
+  }
 
   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()) {
+  Int_t nElSelected=0;
+  // Run electron selection first if trigger particle is an electron
+  if(fTriggerParticle==AliDxHFECorrelation::kElectron){
+
+    if (fSelectedElectrons) delete fSelectedElectrons;
+    // If run on kinematical level, send in MCarray instead of reconstructed tracks
+    if(fUseKine) fSelectedElectrons=(fElectrons->Select(fMCArray, pEvent));
+    else fSelectedElectrons=(fElectrons->Select(pEvent));
+  
+    if(! fSelectedElectrons) {
+      return;
+    }
+
+    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;
+    }
+    // Fill bin with triggered events if electrons are the trigger (only for events with nr electrons >0)
+    if(nElSelected>0) fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsTriggered);
+
+  }
+
+
+  // D0 selection 
+  if(fSelectedD0s) delete fSelectedD0s;
+  if(fUseKine)  fSelectedD0s=(fD0s->Select(fMCArray,pEvent));
+  else 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");
+  // 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;
   }
-  fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsD0);
+  //  Run electron selection second if trigger particle is D0
+  if(fTriggerParticle!=AliDxHFECorrelation::kElectron){
 
+    // Fill bin with triggered events here if D0 are the trigger (only for events with nr D0 >0)
+    if(nD0Selected>0) fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsTriggered);
 
-  std::auto_ptr<TObjArray> pSelectedElectrons(fElectrons->Select(pEvent));
+    if (fSelectedElectrons) delete fSelectedElectrons;
 
-  // 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 run on kinematical level, send in MCarray instead of reconstructed tracks
+    if(fUseKine) fSelectedElectrons=(fElectrons->Select(fMCArray, pEvent));
+    else fSelectedElectrons=(fElectrons->Select(pEvent));
+    if(! fSelectedElectrons) {
+      return;
+    }
 
+    nElSelected =  fSelectedElectrons->GetEntriesFast();
 
-  if(! pSelectedElectrons.get()) {
-    return;
+    // 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;
+    }
   }
 
-  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");
+  // Should not be necessary:
+  if(!fUseEventMixing && (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);
+  if(nD0Selected >0 && nElSelected>0) fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsCorrelated);
 
-  //This is only called if there are electrons and D0s present.
-  int iResult=fCorrelation->Fill(pSelectedD0s.get(), pSelectedElectrons.get());
+  int iResult=0;
+  if(fTriggerParticle==AliDxHFECorrelation::kD) 
+    fCorrelation->Fill(fSelectedD0s, fSelectedElectrons, pEvent);
+  else 
+    fCorrelation->Fill(fSelectedElectrons, fSelectedD0s, pEvent);
 
   if (iResult<0) {
-    AliError(Form("%s processing failed with error %d", fCorrelation->GetName(), iResult));
+    AliFatal(Form("%s processing failed with error %d", fCorrelation->GetName(), iResult));
   }
 
   PostData(1, fOutput);
+  return;
+
+}
+
+int AliAnalysisTaskDxHFECorrelation::ParseArguments(const char* arguments)
+{
+  // parse arguments and set internal flags
+  TString strArguments(arguments);
+  auto_ptr<TObjArray> tokens(strArguments.Tokenize(" "));
+  if (!tokens.get()) return 0;
+  AliInfo(strArguments);
+  TIter next(tokens.get());
+  TObject* token;
+  while ((token=next())) {
+    TString argument=token->GetName();
+   
+    if (argument.BeginsWith("event-mixing")) {
+      fUseEventMixing=true;
+      AliInfo("Running with Event mixing");
+      continue;
+    }
+      
+    if (argument.BeginsWith("mc")) {
+      fUseMC=true;
+      AliInfo("Running on MC data");
+      continue;
+    }
+    if (argument.BeginsWith("usekine") || argument.BeginsWith("kine")) {
+      fUseKine=true;
+      AliInfo("Running on MC stack");
+      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 if (argument.CompareTo("p-Pb")==0){ fSystem=2;}
+      else {
+       AliWarning(Form("can not set collision system, unknown parameter '%s'", argument.Data()));
+       // TODO: check what makes sense
+       fSystem=0;
+      }
+      continue;
+    }
+    if (argument.BeginsWith("trigger=")){
+      argument.ReplaceAll("trigger=", "");
+      if (argument.CompareTo("D0")==0) {fTriggerParticle=AliDxHFECorrelation::kD; AliInfo("CorrTask: trigger on D0");}
+      else if (argument.CompareTo("D")==0){ fTriggerParticle=AliDxHFECorrelation::kD; AliInfo("CorrTask: trigger on D0");}
+      else if (argument.CompareTo("electron")==0) { fTriggerParticle=AliDxHFECorrelation::kElectron; AliInfo("CorrTask: trigger on electron");}
+      continue;
+    }
+    AliWarning(Form("unknown argument '%s'", argument.Data()));
+      
+  }
+
+  return 0;
 }
 
 void AliAnalysisTaskDxHFECorrelation::FinishTaskOutput()
@@ -338,8 +496,8 @@ void AliAnalysisTaskDxHFECorrelation::Terminate(Option_t *)
   // last action on the client
   fOutput = dynamic_cast<TList*> (GetOutputData(1));
   if (!fOutput) {
-    AliFatal("failed to get output container");
+    // looks like that is a valid condition if the task is run
+    // in mode "terminate"
     return;
   }
-
 }