519210ce303cb0e9479effa9ebafa1bff28231ec
[u/mrichter/AliRoot.git] / 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 }