PWGHF/correlationHF/AliAnalysisTaskDxHFECorrelation.cxx

   1 // $Id$
   2
   3 //**************************************************************************
   4 //* This file is property of and copyright by the ALICE Project            *
   5 //* ALICE Experiment at CERN, All rights reserved.                         *
   6 //*                                                                        *
   7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
   8 //*                  Sedat Altinpinar <Sedat.Altinpinar@cern.ch>           *
   9 //*                  Hege Erdal       <hege.erdal@gmail.com>               *
  10 //*                                                                        *
  11 //* Permission to use, copy, modify and distribute this software and its   *
  12 //* documentation strictly for non-commercial purposes is hereby granted   *
  13 //* without fee, provided that the above copyright notice appears in all   *
  14 //* copies and that both the copyright notice and this permission notice   *
  15 //* appear in the supporting documentation. The authors make no claims     *
  16 //* about the suitability of this software for any purpose. It is          *
  17 //* provided "as is" without express or implied warranty.                  *
  18 //**************************************************************************
  19
  20 /// @file   AliAnalysisTaskDxHFECorrelation.cxx
  21 /// @author Sedat Altinpinar, Hege Erdal, Matthias Richter
  22 /// @date   2012-03-19
  23 /// @brief  AnalysisTask D0 - HFE correlation
  24 ///
  25
  26 #include "AliAnalysisTaskDxHFECorrelation.h"
  27 #include "AliDxHFECorrelation.h"
  28 #include "AliDxHFECorrelationMC.h"
  29 #include "AliDxHFEParticleSelectionD0.h"
  30 #include "AliDxHFEParticleSelectionMCD0.h"
  31 #include "AliDxHFEParticleSelectionEl.h"
  32 #include "AliDxHFEParticleSelectionMCEl.h"
  33 #include "AliDxHFEParticleSelection.h"
  34 #include "AliHFCorrelator.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliLog.h"
  37 #include "AliESDInputHandler.h"
  38 #include "AliAnalysisDataSlot.h"
  39 #include "AliAnalysisDataContainer.h"
  40 #include "AliAnalysisManager.h"
  41 #include "AliVertexerTracks.h"
  42 #include "AliAODHandler.h"
  43 #include "AliInputEventHandler.h"
  44 #include "AliAODEvent.h"
  45 #include "AliAODVertex.h"
  46 #include "AliAODTrack.h"
  47 #include "AliAODMCHeader.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliAODRecoDecayHF2Prong.h"
  50 #include "AliAODRecoCascadeHF.h"
  51 #include "AliRDHFCutsD0toKpi.h"
  52 #include "AliPID.h"
  53 #include "AliPIDResponse.h"
  54 #include "AliHFEcontainer.h"
  55 #include "AliHFEpid.h"
  56 #include "AliHFEpidBase.h"
  57 #include "AliHFEcuts.h"
  58 #include "AliHFEtools.h"
  59 #include "TObject.h"
  60 #include "TChain.h"
  61 #include "TSystem.h"
  62 #include "AliReducedParticle.h"
  63 #include "AliHFAssociatedTrackCuts.h" // initialization of event pool
  64 #include "TFile.h"
  65 #include <memory>
  66
  67 using namespace std;
  68
  69 /// ROOT macro for the implementation of ROOT specific class methods
  70 ClassImp(AliAnalysisTaskDxHFECorrelation)
  71
  72 AliAnalysisTaskDxHFECorrelation::AliAnalysisTaskDxHFECorrelation(const char* opt)
  73   : AliAnalysisTaskSE("AliAnalysisTaskDxHFECorrelation")
  74   , fOutput(0)
  75   , fOption(opt)
  76   , fCorrelation(NULL)
  77   , fD0s(NULL)
  78   , fElectrons(NULL)
  79   , fCutsD0(NULL)
  80   , fCutsHFE(NULL)
  81   , fCuts(NULL)
  82   , fPID(NULL)
  83   , fPIDTOF(NULL)
  84   , fFillOnlyD0D0bar(0)
  85   , fUseMC(kFALSE)
  86   , fUseEventMixing(kFALSE)
  87   , fSystem(0)
  88   , fSelectedD0s(NULL)
  89   , fSelectedElectrons(NULL)
  90 {
  91   // constructor
  92   //
  93   //
  94   DefineSlots();
  95   fPID = new AliHFEpid("hfePid");
  96   fPIDTOF = new AliHFEpid("hfePidTOF");
  97
  98 }
  99
 100 int AliAnalysisTaskDxHFECorrelation::DefineSlots()
 101 {
 102   // define the data slots
 103   DefineInput(0, TChain::Class());
 104   DefineOutput(1, TList::Class());
 105   DefineOutput(2,AliRDHFCutsD0toKpi::Class());
 106   DefineOutput(3,AliHFEcuts::Class());
 107   DefineOutput(4,AliHFAssociatedTrackCuts::Class());
 108   return 0;
 109 }
 110
 111 AliAnalysisTaskDxHFECorrelation::~AliAnalysisTaskDxHFECorrelation()
 112 {
 113   // destructor
 114   //
 115   //
 116
 117   // histograms are in the output list and deleted when the output
 118   // list is deleted by the TSelector dtor
 119
 120   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 121     delete fOutput;
 122     fOutput = 0;
 123   }
 124   if (fD0s) delete fD0s;
 125   fD0s=NULL;
 126   if (fElectrons) delete fElectrons;
 127   fElectrons=NULL;
 128   if (fCorrelation) delete fCorrelation;
 129   fCorrelation=NULL;
 130   // external object, do not delete
 131   fCutsD0=NULL;
 132   // external object, do not delete
 133   fCutsHFE=NULL;
 134   if(fPID) delete fPID;
 135   fPID=NULL;
 136   if(fPIDTOF) delete fPIDTOF;
 137   fPIDTOF=NULL;
 138   if(fSelectedElectrons) delete fSelectedElectrons;
 139   fSelectedElectrons=NULL;
 140   if(fSelectedD0s) delete fSelectedD0s;
 141   fSelectedD0s=NULL;
 142
 143
 144 }
 145
 146 void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
 147 {
 148   // create result objects and add to output list
 149   int iResult=0;
 150
 151   //Initialize PID for electron selection
 152   // TODO: Put the initialization of these objects in the AddTask..
 153   // PID for Only TOF
 154   if(!fPIDTOF->GetNumberOfPIDdetectors()) {
 155     fPIDTOF->AddDetector("TOF",0);
 156   }
 157   fPIDTOF->ConfigureTOF(3); // number of sigma TOF
 158
 159   // PID object for TPC and TOF combined
 160   if(!fPID->GetNumberOfPIDdetectors()) {
 161     fPID->AddDetector("TOF",0);
 162     fPID->AddDetector("TPC",1);
 163   }
 164
 165   const int paramSize=4;
 166   Double_t params[paramSize];
 167   memset(params, 0, sizeof(Double_t)*paramSize);
 168   params[0]=-1.;
 169   fPID->ConfigureTPCdefaultCut(NULL, params, 3.);
 170   fPID->InitializePID();
 171
 172   fOutput = new TList;
 173   fOutput->SetOwner();
 174
 175   // setting up for D0s
 176   TString selectionD0Options;
 177   switch (fFillOnlyD0D0bar) {
 178   case 1: selectionD0Options+="FillOnlyD0 "; break;
 179   case 2: selectionD0Options+="FillOnlyD0bar "; break;
 180   default: selectionD0Options+="FillD0D0bar ";
 181   }
 182
 183   if(fUseMC) fD0s=new AliDxHFEParticleSelectionMCD0(selectionD0Options);
 184   else fD0s=new AliDxHFEParticleSelectionD0(selectionD0Options);
 185   fD0s->SetCuts(fCutsD0);
 186   iResult=fD0s->Init();
 187   if (iResult<0) {
 188     AliFatal(Form("initialization of worker class instance fD0s failed with error %d", iResult));
 189   }
 190
 191   //Electrons
 192   if(fUseMC) fElectrons=new AliDxHFEParticleSelectionMCEl;
 193   else fElectrons=new AliDxHFEParticleSelectionEl;
 194   //TODO: Create a TList containing all cut-objects needed for the worker classes
 195   fElectrons->SetCuts(fPID, AliDxHFEParticleSelectionEl::kCutPID);
 196   fElectrons->SetCuts(fPIDTOF, AliDxHFEParticleSelectionEl::kCutPIDTOF);
 197   fElectrons->SetCuts(fCutsHFE, AliDxHFEParticleSelectionEl::kCutHFE);
 198   iResult=fElectrons->Init();
 199   if (iResult<0) {
 200     AliFatal(Form("initialization of worker class instance fElectrons failed with error %d", iResult));
 201   }
 202
 203   //Correlation
 204   if(fUseMC) fCorrelation=new AliDxHFECorrelationMC;
 205   else fCorrelation=new AliDxHFECorrelation;
 206   fCorrelation->SetCuts(fCuts);
 207   // TODO: check if we can get rid of the mc flag in the correlation analysis class
 208   // at the moment this is needed to pass on info to AliHFCorrelator
 209   TString arguments;
 210   if (fUseMC)          arguments+=" use-mc";
 211   if (fUseEventMixing) arguments+=" event-mixing";
 212   // TODO: fSystem is a boolean right now, needs to be changed to fit also p-Pb
 213   if (!fSystem)         arguments+=" system=pp";
 214   else                 arguments+=" system=Pb-Pb";
 215   iResult=fCorrelation->Init(arguments);
 216   if (iResult<0) {
 217     AliFatal(Form("initialization of worker class instance fCorrelation failed with error %d", iResult));
 218   }
 219
 220   // Fix for merging:
 221   // Retrieving the individual objects created
 222   // and storing them instead of fD0s, fElectrons etc..
 223   TList *list =(TList*)fCorrelation->GetControlObjects();
 224   TObject *obj=NULL;
 225
 226   TIter next(list);
 227   while((obj = next())){
 228     fOutput->Add(obj);
 229   }
 230
 231   list=(TList*)fD0s->GetControlObjects();
 232   next=TIter(list);
 233   while((obj= next())){
 234     fOutput->Add(obj);
 235   }
 236
 237   list=(TList*)fElectrons->GetControlObjects();
 238   next=TIter(list);
 239   while((obj = next()))
 240     fOutput->Add(obj);
 241
 242   if (!fCutsD0) {
 243     AliFatal(Form("cut object for D0 missing"));
 244     return;
 245   }
 246
 247   if (!dynamic_cast<AliRDHFCutsD0toKpi*>(fCutsD0)) {
 248     AliFatal(Form("cut object %s is of incorrect type %s, expecting AliRDHFCutsD0toKpi", fCutsD0->GetName(), fCutsD0->ClassName()));
 249     return;
 250   }
 251   // that's the copy for the output stream
 252   AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(dynamic_cast<AliRDHFCutsD0toKpi&>(*fCutsD0));
 253   const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
 254   copyfCuts->SetName(nameoutput);
 255
 256   // all tasks must post data once for all outputs
 257   PostData(1, fOutput);
 258   PostData(2,copyfCuts);
 259   PostData(3,fCutsHFE);
 260   PostData(4,fCuts);
 261
 262 }
 263
 264 void AliAnalysisTaskDxHFECorrelation::UserExec(Option_t* /*option*/)
 265 {
 266   // process the event
 267   TObject* pInput=InputEvent();
 268   if (!pInput) {
 269     AliError("failed to get input");
 270     return;
 271   }
 272   AliVEvent *pEvent = dynamic_cast<AliVEvent*>(pInput);
 273   TClonesArray *inputArray=0;
 274
 275   fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsAll);
 276
 277   if(!pEvent && AODEvent() && IsStandardAOD()) { //Not sure if this is needed.. Keep it for now.
 278     // In case there is an AOD handler writing a standard AOD, use the AOD
 279     // event in memory rather than the input (ESD) event.
 280     pEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 281     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 282     // have to taken from the AOD event hold by the AliAODExtension
 283     AliAODHandler* aodHandler = (AliAODHandler*)
 284       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 285
 286     if(aodHandler->GetExtensions()) {
 287       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 288       AliAODEvent* aodFromExt = ext->GetAOD();
 289       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 290     }
 291   } else if(pEvent) {
 292     inputArray=(TClonesArray*)pEvent->GetList()->FindObject("D0toKpi");
 293   }
 294   if(!inputArray || !pEvent) {
 295     AliError("Input branch not found!\n");
 296     return;
 297   }
 298   // fix for temporary bug in ESDfilter
 299   // the AODs with null vertex pointer didn't pass the PhysSel
 300   if(!pEvent->GetPrimaryVertex() || TMath::Abs(pEvent->GetMagneticField())<0.001){
 301     AliDebug(2,"Rejected at GetPrimaryvertex");
 302     return;
 303   }
 304
 305   AliRDHFCuts* cutsd0=dynamic_cast<AliRDHFCuts*>(fCutsD0);
 306   if (!cutsd0) return; // Fatal thrown already in initialization
 307
 308   if(!cutsd0->IsEventSelected(pEvent)) {
 309     AliDebug(2,"rejected at IsEventSelected");
 310     return;
 311   }
 312
 313   if(!fPID->IsInitialized()){
 314     // Initialize PID with the given run number
 315     AliWarning("PID not initialised, get from Run no");
 316     fPID->InitializePID(pEvent->GetRunNumber());
 317   }
 318   if(!fPIDTOF->IsInitialized()){
 319     // Initialize PID with the given run number
 320     AliWarning("PIDTOF not initialised, get from Run no");
 321     fPIDTOF->InitializePID(pEvent->GetRunNumber());
 322   }
 323
 324   AliPIDResponse *pidResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler())->GetPIDResponse();
 325   if(!pidResponse){
 326     // TODO: consider issuing fatal instead of debug in case pidresponse not available
 327     AliDebug(1, "Using default PID Response");
 328     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 329   }
 330
 331   fPID->SetPIDResponse(pidResponse);
 332   fPIDTOF->SetPIDResponse(pidResponse);
 333   fElectrons->SetPIDResponse(pidResponse);
 334
 335   // Retrieving process from the AODMCHeader.
 336   // TODO: Move it somewhere else? (keep it here for the moment since only need to read once pr event)
 337   if(fUseMC){
 338     AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(pEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 339
 340     if (!mcHeader) {
 341       AliError("Could not find MC Header in AOD");
 342       return;
 343     }
 344     Int_t eventType = mcHeader->GetEventType();
 345     fCorrelation->SetEventType(eventType);
 346   }
 347
 348   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 349
 350   fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsSel);
 351
 352   if(fSelectedD0s) delete fSelectedD0s;
 353   fSelectedD0s=(fD0s->Select(inputArray,pEvent));
 354
 355   if(! fSelectedD0s) {
 356     return;
 357   }
 358   Int_t nD0Selected = fSelectedD0s->GetEntriesFast();
 359
 360
 361   /*std::auto_ptr<TObjArray> pSelectedD0s(fD0s->Select(inputArray,pEvent));
 362   if(! pSelectedD0s.get()) {
 363     return;
 364   }
 365   Int_t nD0Selected = pSelectedD0s->GetEntriesFast();*/
 366
 367   // When not using EventMixing, no need to go further if no D0s are found.
 368   // For Event Mixing, need to store all found electrons in the pool
 369   if(!fUseEventMixing && nD0Selected==0){
 370     AliDebug(4,"No D0s found in this event");
 371     return;
 372   }
 373
 374   fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsD0);
 375
 376   /* std::auto_ptr<TObjArray> pSelectedElectrons(fElectrons->Select(pEvent));
 377   // note: the pointer is deleted automatically once the scope is left
 378   // if the array should be published, the auto pointer must be released
 379   // first, however some other cleanup will be necessary in that case
 380   // probably a clone with a reduced AliVParticle implementation is
 381   // appropriate.
 382
 383   if(! pSelectedElectrons.get()) {
 384     return;
 385   }
 386
 387   Int_t nElSelected =  pSelectedElectrons->GetEntriesFast();*/
 388   if (fSelectedElectrons) delete fSelectedElectrons;
 389   fSelectedElectrons=(fElectrons->Select(pEvent));
 390
 391   if(! fSelectedElectrons) {
 392     return;
 393   }
 394
 395   Int_t nElSelected =  fSelectedElectrons->GetEntriesFast();
 396
 397
 398   // No need to go further if no electrons are found, except for event mixing. Will here anyway correlate D0s with electrons from previous events
 399   if(!fUseEventMixing && nElSelected==0){
 400     AliDebug(4,"No electrons found in this event");
 401     return;
 402   }
 403   if(nD0Selected==0 && nElSelected==0){
 404     AliDebug(4,"Neither D0 nor electrons in this event");
 405     return;
 406   }
 407
 408   AliDebug(4,Form("Number of D0->Kpi Start: %d , End: %d    Electrons Selected: %d\n", nInD0toKpi, nD0Selected, nElSelected));
 409
 410   fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsD0e);
 411
 412   //int iResult=fCorrelation->Fill(pSelectedD0s.get(), pSelectedElectrons.get(), pEvent);
 413   int iResult=fCorrelation->Fill(fSelectedD0s, fSelectedElectrons, pEvent);
 414
 415   if (iResult<0) {
 416     AliFatal(Form("%s processing failed with error %d", fCorrelation->GetName(), iResult));
 417   }
 418
 419   PostData(1, fOutput);
 420   return;
 421
 422 }
 423
 424 void AliAnalysisTaskDxHFECorrelation::FinishTaskOutput()
 425 {
 426   // end of the processing
 427 }
 428
 429 void AliAnalysisTaskDxHFECorrelation::Terminate(Option_t *)
 430 {
 431   // last action on the client
 432   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 433   if (!fOutput) {
 434     // looks like that is a valid condition if the task is run
 435     // in mode "terminate"
 436     return;
 437   }
 438 }