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   , fCuts(NULL)
  81   , fUseMC(kFALSE)
  82   , fUseEventMixing(kFALSE)
  83   , fSystem(0)
  84   , fSelectedD0s(NULL)
  85   , fSelectedElectrons(NULL)
  86   , fListHFE(NULL)
  87   , fTriggerParticle(AliDxHFECorrelation::kD)
  88   , fUseKine(kFALSE)
  89   , fMCArray(NULL)
  90   , fCorrelationArguments("")
  91 {
  92   // constructor
  93   //
  94   //
  95   DefineSlots();
  96 }
  97
  98 int AliAnalysisTaskDxHFECorrelation::DefineSlots()
  99 {
 100   // define the data slots
 101   DefineInput(0, TChain::Class());
 102   DefineOutput(1, TList::Class());
 103   DefineOutput(2,AliRDHFCutsD0toKpi::Class());
 104   DefineOutput(3,TList::Class());
 105   DefineOutput(4,AliHFAssociatedTrackCuts::Class());
 106   return 0;
 107 }
 108
 109 AliAnalysisTaskDxHFECorrelation::~AliAnalysisTaskDxHFECorrelation()
 110 {
 111   // destructor
 112   //
 113   //
 114
 115   // histograms are in the output list and deleted when the output
 116   // list is deleted by the TSelector dtor
 117
 118   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 119     delete fOutput;
 120     fOutput = 0;
 121   }
 122   if (fD0s) delete fD0s;
 123   fD0s=NULL;
 124   if (fElectrons) delete fElectrons;
 125   fElectrons=NULL;
 126   if (fCorrelation) delete fCorrelation;
 127   fCorrelation=NULL;
 128   // external object, do not delete
 129   fCutsD0=NULL;
 130   // external object, do not delete
 131   // TODO: also delete it here? this class is set as owner...
 132   fListHFE=NULL;
 133   if(fSelectedElectrons) delete fSelectedElectrons;
 134   fSelectedElectrons=NULL;
 135   if(fSelectedD0s) delete fSelectedD0s;
 136   fSelectedD0s=NULL;
 137   if(fListHFE) delete fListHFE;
 138   fListHFE=NULL;
 139   if(fMCArray) delete fMCArray;
 140   fMCArray=NULL;
 141
 142
 143 }
 144
 145 void AliAnalysisTaskDxHFECorrelation::UserCreateOutputObjects()
 146 {
 147   // create result objects and add to output list
 148   int iResult=0;
 149   // ParseArguments will also define the strings that are used as input
 150   // for the particle selection classes and the correlation class
 151
 152   ParseArguments(fOption.Data());
 153
 154   fOutput = new TList;
 155   fOutput->SetOwner();
 156
 157   // D0s ===============================================
 158   if(fUseMC) fD0s=new AliDxHFEParticleSelectionMCD0(fOption);
 159   else fD0s=new AliDxHFEParticleSelectionD0(fOption);
 160   fD0s->SetCuts(fCutsD0,AliDxHFEParticleSelectionD0::kCutD0);
 161   iResult=fD0s->Init();
 162   if (iResult<0) {
 163     AliFatal(Form("initialization of worker class instance fD0s failed with error %d", iResult));
 164   }
 165
 166   //Electrons ============================================
 167   fListHFE->SetOwner(); // Not sure if needed
 168   if(fUseMC) fElectrons=new AliDxHFEParticleSelectionMCEl(fOption);
 169   else fElectrons=new AliDxHFEParticleSelectionEl(fOption);
 170   fElectrons->SetCuts(fListHFE, AliDxHFEParticleSelectionEl::kCutList);
 171   iResult=fElectrons->Init();
 172   if (iResult<0) {
 173     AliFatal(Form("initialization of worker class instance fElectrons failed with error %d", iResult));
 174   }
 175
 176   //Correlation ===========================================
 177   if(fUseMC) fCorrelation=new AliDxHFECorrelationMC;
 178   else fCorrelation=new AliDxHFECorrelation;
 179   fCorrelation->SetCuts(fCuts);
 180   iResult=fCorrelation->Init(fOption);
 181   if (iResult<0) {
 182     AliFatal(Form("initialization of worker class instance fCorrelation failed with error %d", iResult));
 183   }
 184
 185   // Fix for merging:
 186   // Retrieving the individual objects created
 187   // and storing them instead of fD0s, fElectrons etc..
 188   TList *list =(TList*)fCorrelation->GetControlObjects();
 189   TObject *obj=NULL;
 190
 191   TIter next(list);
 192   while((obj = next())){
 193     fOutput->Add(obj);
 194   }
 195
 196   list=(TList*)fD0s->GetControlObjects();
 197   next=TIter(list);
 198   while((obj= next())){
 199     fOutput->Add(obj);
 200   }
 201
 202   list=(TList*)fElectrons->GetControlObjects();
 203   next=TIter(list);
 204   while((obj = next()))
 205     fOutput->Add(obj);
 206
 207   if (!fCutsD0) {
 208     AliFatal(Form("cut object for D0 missing"));
 209     return;
 210   }
 211
 212   if (!dynamic_cast<AliRDHFCutsD0toKpi*>(fCutsD0)) {
 213     AliFatal(Form("cut object %s is of incorrect type %s, expecting AliRDHFCutsD0toKpi", fCutsD0->GetName(), fCutsD0->ClassName()));
 214     return;
 215   }
 216   // that's the copy for the output stream
 217   AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(dynamic_cast<AliRDHFCutsD0toKpi&>(*fCutsD0));
 218   const char* nameoutput=GetOutputSlot(2)->GetContainer()->GetName();
 219   copyfCuts->SetName(nameoutput);
 220
 221   // all tasks must post data once for all outputs
 222   PostData(1, fOutput);
 223   PostData(2,copyfCuts);
 224   PostData(3,fListHFE);
 225   PostData(4,fCuts);
 226
 227 }
 228
 229 void AliAnalysisTaskDxHFECorrelation::UserExec(Option_t* /*option*/)
 230 {
 231   // process the event
 232   TObject* pInput=InputEvent();
 233   if (!pInput) {
 234     AliError("failed to get input");
 235     return;
 236   }
 237   AliVEvent *pEvent = dynamic_cast<AliVEvent*>(pInput);
 238   TClonesArray *inputArray=0;
 239
 240   fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsAll);
 241
 242   if(!pEvent && AODEvent() && IsStandardAOD()) { //Not sure if this is needed.. Keep it for now.
 243     // In case there is an AOD handler writing a standard AOD, use the AOD
 244     // event in memory rather than the input (ESD) event.
 245     pEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 246     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 247     // have to taken from the AOD event hold by the AliAODExtension
 248     AliAODHandler* aodHandler = (AliAODHandler*)
 249       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 250
 251     if(aodHandler->GetExtensions()) {
 252       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 253       AliAODEvent* aodFromExt = ext->GetAOD();
 254       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 255     }
 256   } else if(pEvent) {
 257     inputArray=(TClonesArray*)pEvent->GetList()->FindObject("D0toKpi");
 258   }
 259   if(!inputArray || !pEvent) {
 260     AliError("Input branch not found!\n");
 261     return;
 262   }
 263   // fix for temporary bug in ESDfilter
 264   // the AODs with null vertex pointer didn't pass the PhysSel
 265   if(!pEvent->GetPrimaryVertex() || TMath::Abs(pEvent->GetMagneticField())<0.001){
 266     AliDebug(2,"Rejected at GetPrimaryvertex");
 267     return;
 268   }
 269
 270   AliRDHFCuts* cutsd0=dynamic_cast<AliRDHFCuts*>(fCutsD0);
 271   if (!cutsd0) return; // Fatal thrown already in initialization
 272
 273   if(!cutsd0->IsEventSelected(pEvent)) {
 274     // TODO: Fill histograms based on why the event is rejected
 275     AliDebug(2,"rejected at IsEventSelected");
 276     return;
 277   }
 278
 279   if(fSystem==1){
 280     // Not really used anywhere, just as a test of centralityselection
 281     // (Could also be used in histo, so have not removed it)
 282     AliCentrality *centralityObj = 0;
 283     Double_t MultipOrCent = -1;
 284     AliAODEvent* aodEvent=dynamic_cast<AliAODEvent*>(pEvent);
 285     if (aodEvent) {
 286       centralityObj = aodEvent->GetHeader()->GetCentralityP();
 287       if (centralityObj) {
 288         MultipOrCent = centralityObj->GetCentralityPercentileUnchecked("V0M");
 289       }
 290     }
 291     AliInfo(Form("Centrality is %f", MultipOrCent));
 292   }
 293
 294   // Gets the PID response from the analysis manager
 295   AliPIDResponse *pidResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler())->GetPIDResponse();
 296   if(!pidResponse){
 297     // TODO: consider issuing fatal instead of debug in case pidresponse not available
 298     AliDebug(1, "Using default PID Response");
 299     pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class());
 300   }
 301
 302   // Fetching the PID objects from the list, checks if the objects are AliHFEpids
 303   // If so, checks if they are initialized and also sets the pidresponse
 304   TObject *obj=NULL;
 305   TIter next(fListHFE);
 306   while((obj = next())){
 307     AliHFEpid* pidObj=dynamic_cast<AliHFEpid*>(obj);
 308     if(pidObj){
 309       if(!pidObj->IsInitialized()){
 310         AliWarning("PID not initialised, get from Run no");
 311         pidObj->InitializePID(pEvent->GetRunNumber());
 312       }
 313       pidObj->SetPIDResponse(pidResponse);
 314     }
 315   }
 316
 317   // Also sends the pidresponse to the particle selection class for electron
 318   fElectrons->SetPIDResponse(pidResponse);
 319
 320   // Retrieving process from the AODMCHeader.
 321   // TODO: Move it somewhere else? (keep it here for the moment since only need to read once pr event)
 322   if(fUseMC){
 323     AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(pEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 324
 325     if (!mcHeader) {
 326       AliError("Could not find MC Header in AOD");
 327       return;
 328     }
 329     Int_t eventType = mcHeader->GetEventType();
 330     fCorrelation->SetEventType(eventType);
 331     if(fUseKine){
 332       fMCArray = dynamic_cast<TObjArray*>(pEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 333       if(fUseMC && !fMCArray){
 334         AliError("Array of MC particles not found");
 335         return;
 336       }
 337     }
 338   }
 339
 340   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 341
 342   fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsSel);
 343
 344   Int_t nElSelected=0;
 345   // Run electron selection first if trigger particle is an electron
 346   if(fTriggerParticle==AliDxHFECorrelation::kElectron){
 347
 348     if (fSelectedElectrons) delete fSelectedElectrons;
 349     // If run on kinematical level, send in MCarray instead of reconstructed tracks
 350     if(fUseKine) fSelectedElectrons=(fElectrons->Select(fMCArray, pEvent));
 351     else fSelectedElectrons=(fElectrons->Select(pEvent));
 352
 353     if(! fSelectedElectrons) {
 354       return;
 355     }
 356
 357     nElSelected =  fSelectedElectrons->GetEntriesFast();
 358
 359     // No need to go further if no electrons are found, except for event mixing. Will here anyway correlate D0s with electrons from previous events
 360     if(!fUseEventMixing && nElSelected==0){
 361       AliDebug(4,"No electrons found in this event");
 362       return;
 363     }
 364     // Fill bin with triggered events if electrons are the trigger (only for events with nr electrons >0)
 365     if(nElSelected>0) fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsTriggered);
 366
 367   }
 368
 369
 370   // D0 selection
 371   if(fSelectedD0s) delete fSelectedD0s;
 372   if(fUseKine)  fSelectedD0s=(fD0s->Select(fMCArray,pEvent));
 373   else fSelectedD0s=(fD0s->Select(inputArray,pEvent));
 374
 375   if(! fSelectedD0s) {
 376     return;
 377   }
 378   Int_t nD0Selected = fSelectedD0s->GetEntriesFast();
 379
 380   // When not using EventMixing, no need to go further if no D0s are found.
 381   // For Event Mixing, need to store all found electrons in the pool
 382   if(!fUseEventMixing && nD0Selected==0){
 383     AliDebug(4,"No D0s found in this event");
 384     return;
 385   }
 386   //  Run electron selection second if trigger particle is D0
 387   if(fTriggerParticle!=AliDxHFECorrelation::kElectron){
 388
 389     // Fill bin with triggered events here if D0 are the trigger (only for events with nr D0 >0)
 390     if(nD0Selected>0) fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsTriggered);
 391
 392     if (fSelectedElectrons) delete fSelectedElectrons;
 393
 394     // If run on kinematical level, send in MCarray instead of reconstructed tracks
 395     if(fUseKine) fSelectedElectrons=(fElectrons->Select(fMCArray, pEvent));
 396     else fSelectedElectrons=(fElectrons->Select(pEvent));
 397     if(! fSelectedElectrons) {
 398       return;
 399     }
 400
 401     nElSelected =  fSelectedElectrons->GetEntriesFast();
 402
 403     // No need to go further if no electrons are found, except for event mixing. Will here anyway correlate D0s with electrons from previous events
 404     if(!fUseEventMixing && nElSelected==0){
 405       AliDebug(4,"No electrons found in this event");
 406       return;
 407     }
 408   }
 409
 410   // Should not be necessary:
 411   if(!fUseEventMixing && (nD0Selected==0 && nElSelected==0)){
 412     AliDebug(4,"Neither D0 nor electrons in this event");
 413     return;
 414   }
 415
 416   AliDebug(4,Form("Number of D0->Kpi Start: %d , End: %d    Electrons Selected: %d\n", nInD0toKpi, nD0Selected, nElSelected));
 417
 418   if(nD0Selected >0 && nElSelected>0) fCorrelation->HistogramEventProperties(AliDxHFECorrelation::kEventsCorrelated);
 419
 420   int iResult=0;
 421   if(fTriggerParticle==AliDxHFECorrelation::kD)
 422     fCorrelation->Fill(fSelectedD0s, fSelectedElectrons, pEvent);
 423   else
 424     fCorrelation->Fill(fSelectedElectrons, fSelectedD0s, pEvent);
 425
 426   if (iResult<0) {
 427     AliFatal(Form("%s processing failed with error %d", fCorrelation->GetName(), iResult));
 428   }
 429
 430   PostData(1, fOutput);
 431   return;
 432
 433 }
 434
 435 int AliAnalysisTaskDxHFECorrelation::ParseArguments(const char* arguments)
 436 {
 437   // parse arguments and set internal flags
 438   TString strArguments(arguments);
 439   auto_ptr<TObjArray> tokens(strArguments.Tokenize(" "));
 440   if (!tokens.get()) return 0;
 441   AliInfo(strArguments);
 442   TIter next(tokens.get());
 443   TObject* token;
 444   while ((token=next())) {
 445     TString argument=token->GetName();
 446
 447     if (argument.BeginsWith("event-mixing")) {
 448       fUseEventMixing=true;
 449       AliInfo("Running with Event mixing");
 450       continue;
 451     }
 452
 453     if (argument.BeginsWith("mc")) {
 454       fUseMC=true;
 455       AliInfo("Running on MC data");
 456       continue;
 457     }
 458     if (argument.BeginsWith("usekine") || argument.BeginsWith("kine")) {
 459       fUseKine=true;
 460       AliInfo("Running on MC stack");
 461       continue;
 462     }
 463     if (argument.BeginsWith("system=")) {
 464       argument.ReplaceAll("system=", "");
 465       if (argument.CompareTo("pp")==0) {fSystem=0;}
 466       else if (argument.CompareTo("Pb-Pb")==0){ fSystem=1;}
 467       else if (argument.CompareTo("p-Pb")==0){ fSystem=2;}
 468       else {
 469         AliWarning(Form("can not set collision system, unknown parameter '%s'", argument.Data()));
 470         // TODO: check what makes sense
 471         fSystem=0;
 472       }
 473       continue;
 474     }
 475     if (argument.BeginsWith("trigger=")){
 476       argument.ReplaceAll("trigger=", "");
 477       if (argument.CompareTo("D0")==0) {fTriggerParticle=AliDxHFECorrelation::kD; AliInfo("CorrTask: trigger on D0");}
 478       else if (argument.CompareTo("D")==0){ fTriggerParticle=AliDxHFECorrelation::kD; AliInfo("CorrTask: trigger on D0");}
 479       else if (argument.CompareTo("electron")==0) { fTriggerParticle=AliDxHFECorrelation::kElectron; AliInfo("CorrTask: trigger on electron");}
 480       continue;
 481     }
 482     AliWarning(Form("unknown argument '%s'", argument.Data()));
 483
 484   }
 485
 486   return 0;
 487 }
 488
 489 void AliAnalysisTaskDxHFECorrelation::FinishTaskOutput()
 490 {
 491   // end of the processing
 492 }
 493
 494 void AliAnalysisTaskDxHFECorrelation::Terminate(Option_t *)
 495 {
 496   // last action on the client
 497   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 498   if (!fOutput) {
 499     // looks like that is a valid condition if the task is run
 500     // in mode "terminate"
 501     return;
 502   }
 503 }