moving filling of V0 histogram later
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / AliAnalysisTaskPhiCorrelations.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id:$ */
17
18 #include <TROOT.h>
19 #include <TChain.h>
20 #include <TFile.h>
21 #include <TList.h>
22 #include <TMath.h>
23 #include <TTree.h>
24 #include <TH2F.h>
25 #include <TH3F.h>
26 #include <TRandom.h>
27 #include <TParameter.h>
28
29 #include "AliAnalysisTaskPhiCorrelations.h"
30 #include "AliAnalyseLeadingTrackUE.h"
31 #include "AliUEHistograms.h"
32 #include "AliUEHist.h"
33
34 #include "AliAnalysisManager.h"
35 #include "AliAODHandler.h"
36 #include "AliAODInputHandler.h"
37 #include "AliAODMCParticle.h"
38 #include "AliInputEventHandler.h"
39 #include "AliLog.h"
40 #include "AliMCEventHandler.h"
41 #include "AliVParticle.h"
42 #include "AliCFContainer.h"
43
44 #include "AliESDEvent.h"
45 #include "AliESDInputHandler.h"
46 #include "AliMultiplicity.h"
47 #include "AliCentrality.h"
48 #include "AliStack.h"
49 #include "AliAODMCHeader.h"
50 #include "AliGenCocktailEventHeader.h"
51 #include "AliGenEventHeader.h"
52
53 #include "AliEventPoolManager.h"
54
55 #include "AliESDZDC.h"
56 #include "AliESDtrackCuts.h"
57
58 #include "AliHelperPID.h"
59 #include "AliAnalysisUtils.h"
60 #include "TMap.h"
61
62 ////////////////////////////////////////////////////////////////////////
63 //
64 // Analysis class for azimuthal correlation studies
65 // Based on the UE task from Sara Vallero and Jan Fiete
66 //
67 // This class needs input AODs.
68 // The output is a list of analysis-specific containers.
69 //
70 // The AOD can be either connected to the InputEventHandler  
71 // for a chain of AOD files 
72 // or 
73 // to the OutputEventHandler
74 // for a chain of ESD files,
75 // in this case the class should be in the train after the jet-finder
76 //
77 //    Authors:
78 //    Jan Fiete Grosse-Oetringhaus
79 // 
80 ////////////////////////////////////////////////////////////////////////
81
82
83 ClassImp( AliAnalysisTaskPhiCorrelations )
84 ClassImp( AliDPhiBasicParticle )
85
86 //____________________________________________________________________
87 AliAnalysisTaskPhiCorrelations:: AliAnalysisTaskPhiCorrelations(const char* name):
88 AliAnalysisTask(name,""),
89 // general configuration
90 fDebug(0),
91 fMode(0),
92 fReduceMemoryFootprint(kFALSE),
93 fFillMixed(kTRUE),
94 fMixingTracks(50000),
95 fTwoTrackEfficiencyStudy(kFALSE),
96 fTwoTrackEfficiencyCut(0),
97 fUseVtxAxis(kFALSE),
98 fCourseCentralityBinning(kFALSE),
99 fSkipTrigger(kFALSE),
100 fInjectedSignals(kFALSE),
101 fHelperPID(0x0),
102 fAnalysisUtils(0x0),
103 fMap(0x0),
104 // pointers to UE classes
105 fAnalyseUE(0x0),
106 fHistos(0x0),
107 fHistosMixed(0),
108 fEfficiencyCorrectionTriggers(0),
109 fEfficiencyCorrectionAssociated(0),
110 // handlers and events
111 fAOD(0x0),
112 fESD(0x0),
113 fArrayMC(0x0),
114 fInputHandler(0x0),
115 fMcEvent(0x0),
116 fMcHandler(0x0),
117 fPoolMgr(0x0),
118 // histogram settings
119 fListOfHistos(0x0), 
120 // event QA
121 fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default) 
122 fZVertex(7.),
123 fCentralityMethod("V0M"),
124 // track cuts
125 fTrackEtaCut(0.8),
126 fTrackEtaCutMin(-1.),
127 fOnlyOneEtaSide(0),
128 fPtMin(0.5),
129 fDCAXYCut(0),
130 fFilterBit(0xFF),
131 fTrackStatus(0),
132 fSelectBit(AliVEvent::kMB|AliVEvent::kUserDefined),
133 fUseChargeHadrons(kFALSE),
134 fParticleSpeciesTrigger(-1),
135 fParticleSpeciesAssociated(-1),
136 fCheckMotherPDG(kTRUE),
137 fSelectCharge(0),
138 fTriggerSelectCharge(0),
139 fAssociatedSelectCharge(0),
140 fTriggerRestrictEta(-1),
141 fEtaOrdering(kFALSE),
142 fCutConversions(kFALSE),
143 fCutResonances(kFALSE),
144 fRejectResonanceDaughters(-1),
145 fFillOnlyStep0(kFALSE),
146 fSkipStep6(kFALSE),
147 fRejectCentralityOutliers(kFALSE),
148 fRejectZeroTrackEvents(kFALSE),
149 fRemoveWeakDecays(kFALSE),
150 fRemoveDuplicates(kFALSE),
151 fSkipFastCluster(kFALSE),
152 fWeightPerEvent(kFALSE),
153 fCustomBinning(),
154 fPtOrder(kTRUE),
155 fFillpT(kFALSE)
156 {
157   // Default constructor
158   // Define input and output slots here
159   // Input slot #0 works with a TChain
160   DefineInput(0, TChain::Class());
161   // Output slot #0 writes into a TList container
162   DefineOutput(0, TList::Class());
163 }
164
165 AliAnalysisTaskPhiCorrelations::~AliAnalysisTaskPhiCorrelations() 
166
167   // destructor
168   
169   if (fListOfHistos  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) 
170     delete fListOfHistos;
171 }
172
173 //____________________________________________________________________
174 void AliAnalysisTaskPhiCorrelations::ConnectInputData(Option_t* /*option*/)
175 {
176   
177   // Connect the input data  
178   if (fDebug > 1) AliInfo("ConnectInputData() ");
179   
180   // Since AODs can either be connected to the InputEventHandler
181   // or to the OutputEventHandler ( the AOD is created by a previus task in the train )
182   // we need to get the pointer to the AODEvent correctly.
183   
184   // Delta AODs are also accepted.
185   
186   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
187   
188   if( handler && handler->InheritsFrom("AliAODInputHandler") ) 
189   { // input AOD
190     fAOD = ((AliAODInputHandler*)handler)->GetEvent();
191     if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
192   } 
193   else 
194   {  //output AOD
195     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
196     if (handler && handler->InheritsFrom("AliAODHandler") ) 
197     {
198       fAOD = ((AliAODHandler*)handler)->GetAOD();
199       if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
200     } 
201     else 
202     {  // no AOD
203       AliWarning("I can't get any AOD Event Handler");
204     }
205   }
206   
207   if (handler && handler->InheritsFrom("AliESDInputHandler") ) 
208   { // input ESD
209     // pointer received per event in ::Exec
210     if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliESDInputHandler");
211   } 
212   
213   // Initialize common pointers
214   Initialize();
215 }
216
217 //____________________________________________________________________
218 void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
219 {
220   // Create the output container
221   
222   if (fDebug > 1) AliInfo("CreateOutputObjects()");
223    
224   // Initialize class with main algorithms, event and track selection. 
225   fAnalyseUE = new AliAnalyseLeadingTrackUE();
226   fAnalyseUE->SetParticleSelectionCriteria(fFilterBit, fUseChargeHadrons, fTrackEtaCut, fTrackEtaCutMin, fPtMin);
227   fAnalyseUE->SetDCAXYCut(fDCAXYCut);
228   fAnalyseUE->SetTrackStatus(fTrackStatus);
229   fAnalyseUE->SetCheckMotherPDG(fCheckMotherPDG);
230   fAnalyseUE->SetDebug(fDebug); 
231   fAnalyseUE->DefineESDCuts(fFilterBit);
232   fAnalyseUE->SetEventSelection(fSelectBit);
233   fAnalyseUE->SetHelperPID(fHelperPID);
234   if ((fParticleSpeciesTrigger != -1 || fParticleSpeciesAssociated != -1) && !fHelperPID)
235     AliFatal("HelperPID object should be set in the steering macro");
236
237   // Initialize output list of containers
238   if (fListOfHistos != NULL){
239         delete fListOfHistos;
240         fListOfHistos = NULL;
241         }
242   if (!fListOfHistos){
243         fListOfHistos = new TList();
244         fListOfHistos->SetOwner(kTRUE); 
245         }
246
247   // Initialize class to handle histograms 
248   TString histType = "4R";
249   if (fUseVtxAxis == 1)
250     histType = "5R";
251   else if (fUseVtxAxis == 2)
252     histType = "6R";
253   if (fCourseCentralityBinning)
254     histType += "C";
255   fHistos = new AliUEHistograms("AliUEHistogramsSame", histType, fCustomBinning);
256   fHistosMixed = new AliUEHistograms("AliUEHistogramsMixed", histType, fCustomBinning);
257   
258   fHistos->SetSelectCharge(fSelectCharge);
259   fHistosMixed->SetSelectCharge(fSelectCharge);
260   
261   fHistos->SetSelectTriggerCharge(fTriggerSelectCharge);
262   fHistosMixed->SetSelectTriggerCharge(fTriggerSelectCharge);
263   
264   fHistos->SetSelectAssociatedCharge(fAssociatedSelectCharge);
265   fHistosMixed->SetSelectAssociatedCharge(fAssociatedSelectCharge);
266
267   fHistos->SetTriggerRestrictEta(fTriggerRestrictEta);
268   fHistosMixed->SetTriggerRestrictEta(fTriggerRestrictEta);
269   
270   fHistos->SetOnlyOneEtaSide(fOnlyOneEtaSide);
271   fHistosMixed->SetOnlyOneEtaSide(fOnlyOneEtaSide);
272   
273   fHistos->SetEtaOrdering(fEtaOrdering);
274   fHistosMixed->SetEtaOrdering(fEtaOrdering);
275
276   fHistos->SetPairCuts(fCutConversions, fCutResonances);
277   fHistosMixed->SetPairCuts(fCutConversions, fCutResonances);
278   
279   fHistos->SetRejectResonanceDaughters(fRejectResonanceDaughters);
280   fHistosMixed->SetRejectResonanceDaughters(fRejectResonanceDaughters);
281   
282   fHistos->SetTrackEtaCut(fTrackEtaCut);
283   fHistosMixed->SetTrackEtaCut(fTrackEtaCut);
284   
285   fHistos->SetWeightPerEvent(fWeightPerEvent);
286   fHistosMixed->SetWeightPerEvent(fWeightPerEvent);
287   
288   fHistos->SetPtOrder(fPtOrder);
289   fHistosMixed->SetPtOrder(fPtOrder);
290   
291   if (fEfficiencyCorrectionTriggers)
292    {
293     fHistos->SetEfficiencyCorrectionTriggers(fEfficiencyCorrectionTriggers);
294     fHistosMixed->SetEfficiencyCorrectionTriggers((THnF*) fEfficiencyCorrectionTriggers->Clone());
295    }
296   if (fEfficiencyCorrectionAssociated)
297   {
298     fHistos->SetEfficiencyCorrectionAssociated(fEfficiencyCorrectionAssociated);
299     fHistosMixed->SetEfficiencyCorrectionAssociated((THnF*) fEfficiencyCorrectionAssociated->Clone());
300   }
301   
302   // add histograms to list
303   fListOfHistos->Add(fHistos);
304   fListOfHistos->Add(fHistosMixed);
305   // add HelperPID to list
306   if (fHelperPID)
307     fListOfHistos->Add(fHelperPID);
308   // add TMap to list
309   if (fMap)
310     fListOfHistos->Add(fMap);
311   
312   fListOfHistos->Add(new TH2F("trackletsVsV0Cent", ";L1 clusters;v0 centrality", 100, -0.5, 9999.5, 101, 0, 101));
313   fListOfHistos->Add(new TH2F("processIDs", ";#Delta#phi;process id", 100, -0.5 * TMath::Pi(), 1.5 * TMath::Pi(), kPNoProcess + 1, -0.5, kPNoProcess + 0.5));
314   fListOfHistos->Add(new TH1F("eventStat", ";;events", 4, -0.5, 3.5));
315   fListOfHistos->Add(new TH2F("mixedDist", ";centrality;tracks;events", 101, 0, 101, 200, 0, fMixingTracks * 1.5));
316   fListOfHistos->Add(new TH1F("pids", ";pdg;tracks", 2001, -1000.5, 1000.5));
317   fListOfHistos->Add(new TH2F("referenceMultiplicity", ";centrality;tracks;events", 101, 0, 101, 200, 0, 200));
318   fListOfHistos->Add(new TH2F("V0AMult", "V0A multiplicity;V0A multiplicity;V0A multiplicity (scaled)", 1000, -.5, 999.5, 1000, -.5, 999.5));
319   fListOfHistos->Add(new TH2F("V0AMultCorrelation", "V0A multiplicity;V0A multiplicity;SPD tracklets", 1000, -.5, 999.5, 1000, -.5, 999.5));
320   
321   PostData(0,fListOfHistos);
322   
323   // Add task configuration to output list 
324   AddSettingsTree();
325
326   // event mixing
327   Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemention of AliEventPoolManager
328    
329   Int_t nCentralityBins  = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(1);
330   Double_t* centralityBins = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray();
331   
332   const Int_t kNZvtxBins  = 10+(1+10)*4;
333   // bins for further buffers are shifted by 100 cm
334   Double_t vertexBins[kNZvtxBins+1] = { -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10, 
335                                        90,  92,  94,  96,  98, 100, 102, 104, 106, 108, 110, 
336                                       190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210, 
337                                       290, 292, 294, 296, 298, 300, 302, 304, 306, 308, 310, 
338                                       390, 392, 394, 396, 398, 400, 402, 404, 406, 408, 410 };
339
340   Int_t nZvtxBins  = kNZvtxBins;
341   Double_t* zvtxbin = vertexBins;
342
343   if (fMode == 0 && fHistos->GetUEHist(2)->GetEventHist()->GetNVar() > 2)
344   {
345     nZvtxBins = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(2);
346     zvtxbin = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetXbins()->GetArray();
347   }
348
349   fPoolMgr = new AliEventPoolManager(poolsize, fMixingTracks, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
350   fPoolMgr->SetTargetValues(fMixingTracks, 0.1, 5);
351 }
352
353 //____________________________________________________________________
354 void  AliAnalysisTaskPhiCorrelations::Exec(Option_t */*option*/)
355 {
356   // exec (per event)
357   fAnalyseUE->NextEvent();
358   
359   // receive ESD pointer if we are not running AOD analysis
360   if (!fAOD)
361   {
362     AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
363     if (handler && handler->InheritsFrom("AliESDInputHandler"))
364       fESD = ((AliESDInputHandler*)handler)->GetEvent();
365   }
366
367   if (fMode)
368   {
369     // correction mode
370     
371     if (fMcHandler)
372       fMcEvent = fMcHandler->MCEvent();
373
374     if (fAOD)
375     {
376       // array of MC particles
377       fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
378       if (!fArrayMC)
379         AliFatal("No array of MC particles found !!!");
380     }
381     
382     AnalyseCorrectionMode();
383   }
384   else AnalyseDataMode();
385 }
386
387 /******************** ANALYSIS METHODS *****************************/
388
389 //____________________________________________________________________
390 void  AliAnalysisTaskPhiCorrelations::AddSettingsTree()
391 {
392   //Write settings to output list
393   TTree *settingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
394   settingsTree->Branch("fnTracksVertex", &fnTracksVertex,"nTracksVertex/I");
395   settingsTree->Branch("fZVertex", &fZVertex,"ZVertex/D");
396   //settingsTree->Branch("fCentralityMethod", fCentralityMethod.Data(),"CentralityMethod/C");
397   settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
398   settingsTree->Branch("fTrackEtaCutMin", &fTrackEtaCutMin, "TrackEtaCutMin/D");
399   settingsTree->Branch("fOnlyOneEtaSide", &fOnlyOneEtaSide,"OnlyOneEtaSide/I");
400   settingsTree->Branch("fPtMin", &fPtMin, "PtMin/D");
401   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
402   settingsTree->Branch("fTrackStatus", &fTrackStatus,"TrackStatus/I");
403   settingsTree->Branch("fSelectBit", &fSelectBit,"EventSelectionBit/I");
404   settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
405   settingsTree->Branch("fParticleSpeciesTrigger", &fParticleSpeciesTrigger,"ParticleSpeciesTrigger/I");
406   settingsTree->Branch("fParticleSpeciesAssociated", &fParticleSpeciesAssociated,"ParticleSpeciesAssociated/I");
407   settingsTree->Branch("fCheckMotherPDG", &fCheckMotherPDG,"CheckMotherPDG/I");
408   settingsTree->Branch("fSelectCharge", &fSelectCharge,"SelectCharge/I");
409   settingsTree->Branch("fTriggerSelectCharge", &fTriggerSelectCharge,"TriggerSelectCharge/I");
410   settingsTree->Branch("fAssociatedSelectCharge", &fAssociatedSelectCharge,"fAssociatedSelectCharge/I");
411   settingsTree->Branch("fTriggerRestrictEta", &fTriggerRestrictEta,"TriggerRestrictEta/D");
412   settingsTree->Branch("fEtaOrdering", &fEtaOrdering,"EtaOrdering/O");
413   settingsTree->Branch("fCutConversions", &fCutConversions,"CutConversions/O");
414   settingsTree->Branch("fCutResonances", &fCutResonances,"CutResonances/O");
415   settingsTree->Branch("fRejectResonanceDaughters", &fRejectResonanceDaughters,"RejectResonanceDaughters/I");
416   settingsTree->Branch("fFillpT", &fFillpT,"FillpT/O");
417   settingsTree->Branch("fMixingTracks", &fMixingTracks,"MixingTracks/I");
418   settingsTree->Branch("fSkipTrigger", &fSkipTrigger,"SkipTrigger/O");
419   settingsTree->Branch("fInjectedSignals", &fInjectedSignals,"SkipTrigger/O");
420   settingsTree->Branch("fRejectCentralityOutliers", &fRejectCentralityOutliers,"RejectCentralityOutliers/O");
421   settingsTree->Branch("fRejectZeroTrackEvents", &fRejectZeroTrackEvents,"RejectZeroTrackEvents/O");
422   settingsTree->Branch("fRemoveWeakDecays", &fRemoveWeakDecays,"RemoveWeakDecays/O");
423   settingsTree->Branch("fRemoveDuplicates", &fRemoveDuplicates,"RemoveDuplicates/O");
424   settingsTree->Branch("fSkipFastCluster", &fSkipFastCluster,"SkipFastCluster/O");
425   settingsTree->Branch("fWeightPerEvent", &fWeightPerEvent,"WeightPerEvent/O");
426   settingsTree->Branch("fPtOrder", &fPtOrder,"PtOrder/O");
427   //fCustomBinning
428   
429   settingsTree->Fill();
430   fListOfHistos->Add(settingsTree);
431 }  
432
433 //____________________________________________________________________
434 void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
435 {
436   // Run the analysis on MC to get the correction maps
437   //
438  
439   if ( fDebug > 3 ) AliInfo( " Processing event in Corrections mode ..." );
440   
441   Double_t centrality = 0;
442   
443   if (fCentralityMethod.Length() > 0)
444   {
445     AliCentrality *centralityObj = 0;
446     if (fAOD)
447       centralityObj = fAOD->GetHeader()->GetCentralityP();
448     else if (fESD)
449       centralityObj = fESD->GetCentrality();
450     
451     if (centralityObj)
452     {
453       centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
454       AliInfo(Form("Centrality is %f", centrality));
455     }
456     else
457     {
458       Printf("WARNING: Centrality object is 0");
459       centrality = -1;
460      }
461   }
462   
463   // Support for ESD and AOD based analysis
464   AliVEvent* inputEvent = fAOD;
465   if (!inputEvent)
466     inputEvent = fESD;
467
468   Float_t bSign = 0;
469   
470   if (inputEvent)
471   {
472     fHistos->SetRunNumber(inputEvent->GetRunNumber());
473     bSign = (inputEvent->GetMagneticField() > 0) ? 1 : -1;
474   }
475     
476   TObject* mc = fArrayMC;
477   if (!mc)
478     mc = fMcEvent;
479   
480   // count all events
481   fHistos->FillEvent(centrality, -1);
482   
483   if (centrality < 0)
484     return;
485
486   // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
487   TObject* vertexSupplier = fMcEvent;
488   if (fAOD) // AOD
489     vertexSupplier = fAOD->FindListObject(AliAODMCHeader::StdBranchName());
490     
491   if (!fAnalyseUE->VertexSelection(vertexSupplier, 0, fZVertex)) 
492     return;
493   
494   Float_t zVtx = 0;
495   if (fAOD)
496     zVtx = ((AliAODMCHeader*) vertexSupplier)->GetVtxZ();
497   else
498     zVtx = fMcEvent->GetPrimaryVertex()->GetZ();
499   Float_t weight = 1;
500   if (fFillpT)
501     weight = -1;
502   
503   // For productions with injected signals, figure out above which label to skip particles/tracks
504   Int_t skipParticlesAbove = 0;
505   if (fInjectedSignals)
506   {
507     AliGenEventHeader* eventHeader = 0;
508     Int_t headers = 0;
509
510     if (fMcEvent)
511     {
512       // ESD
513       AliHeader* header = (AliHeader*) fMcEvent->Header();
514       if (!header)
515         AliFatal("fInjectedSignals set but no MC header found");
516         
517       AliGenCocktailEventHeader* cocktailHeader = dynamic_cast<AliGenCocktailEventHeader*> (header->GenEventHeader());
518       if (!cocktailHeader)
519       {
520         header->Dump();
521         AliFatal("fInjectedSignals set but no MC cocktail header found");
522       }
523
524       headers = cocktailHeader->GetHeaders()->GetEntries();
525       eventHeader = dynamic_cast<AliGenEventHeader*> (cocktailHeader->GetHeaders()->First());
526     }
527     else
528     {
529       // AOD
530       AliAODMCHeader* header = (AliAODMCHeader*) fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
531       if (!header)
532         AliFatal("fInjectedSignals set but no MC header found");
533       
534       headers = header->GetNCocktailHeaders();
535       eventHeader = header->GetCocktailHeader(0);
536     }
537     
538     if (!eventHeader)
539     {
540       // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing 
541       // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
542       AliError("First event header not found. Skipping this event.");
543       fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
544       return;
545     }
546     
547     skipParticlesAbove = eventHeader->NProduced();
548     AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
549   }
550   
551   // Get MC primaries
552   // triggers
553   TObjArray* tmpList = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
554   CleanUp(tmpList, mc, skipParticlesAbove);
555   TObjArray* tracksMC = CloneAndReduceTrackList(tmpList);
556   delete tmpList;
557   
558   // associated
559   TObjArray* tracksCorrelateMC = tracksMC;
560   if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
561   {
562     tmpList = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
563     CleanUp(tmpList, mc, skipParticlesAbove);
564     tracksCorrelateMC = CloneAndReduceTrackList(tmpList);
565     delete tmpList;
566   }
567   
568   /*
569   if (fAOD)
570   {
571     for (Int_t i=0; i<fArrayMC->GetEntriesFast(); i++)
572       ((TH1F*) fListOfHistos->FindObject("pids"))->Fill(((AliAODMCParticle*) fArrayMC->At(i))->PdgCode());
573   }
574   else
575   {
576     for (Int_t i=0; i<fMcEvent->GetNumberOfTracks(); i++)
577       ((TH1F*) fListOfHistos->FindObject("pids"))->Fill(fMcEvent->GetTrack(i)->PdgCode());
578   }
579   */
580   
581   if (fFillOnlyStep0)
582     zVtx = 0;
583   
584   // (MC-true all particles)
585   // STEP 0
586   fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, tracksCorrelateMC, weight);
587   
588   // mixed event
589   if (fFillMixed)
590   {
591     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
592     if (pool->IsReady())
593       for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
594         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
595     pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateMC));
596   }
597   
598 //   Printf("trigger: %d", ((AliInputEventHandler*)fInputHandler)->IsEventSelected());
599   
600   // Trigger selection ************************************************
601   if (!fFillOnlyStep0 && (fSkipTrigger || fAnalyseUE->TriggerSelection(fInputHandler)))
602   {  
603     // (MC-true all particles)
604     // STEP 1
605     if (!fReduceMemoryFootprint)
606       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTriggered, tracksMC, tracksCorrelateMC, weight);
607     else
608       fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
609       
610     if (!inputEvent) {
611       AliFatal("UNEXPECTED: inputEvent is 0. Trigger selection should have failed");
612       return;
613     }
614     
615     // Vertex selection *************************************************
616     if (fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex))
617     {
618       // fill here for tracking efficiency
619       // loop over particle species
620       
621       for (Int_t particleSpecies = 0; particleSpecies < 4; particleSpecies++)
622       {
623         TObjArray* primMCParticles = fAnalyseUE->GetAcceptedParticles(mc, 0x0, kTRUE, particleSpecies, kTRUE);
624         TObjArray* primRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE, kFALSE);
625         TObjArray* allRecoTracksMatched  = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE, kFALSE);
626         TObjArray* primRecoTracksMatchedPID = 0;
627         TObjArray* allRecoTracksMatchedPID  = 0;
628         
629         if (fHelperPID)
630         {
631           primRecoTracksMatchedPID = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE, kTRUE);
632           allRecoTracksMatchedPID  = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE, kTRUE);
633         }
634         
635         CleanUp(primMCParticles, mc, skipParticlesAbove);
636         CleanUp(primRecoTracksMatched, mc, skipParticlesAbove);
637         CleanUp(allRecoTracksMatched, mc, skipParticlesAbove);
638         CleanUp(primRecoTracksMatchedPID, mc, skipParticlesAbove);
639         CleanUp(allRecoTracksMatchedPID, mc, skipParticlesAbove);
640         
641         // select charges
642         if (fTriggerSelectCharge != 0)
643         {
644           SelectCharge(primMCParticles);
645           SelectCharge(primRecoTracksMatched);
646           SelectCharge(allRecoTracksMatched);
647           SelectCharge(primRecoTracksMatchedPID);
648           SelectCharge(allRecoTracksMatchedPID);
649         }
650       
651         fHistos->FillTrackingEfficiency(primMCParticles, primRecoTracksMatched, allRecoTracksMatched, primRecoTracksMatchedPID, allRecoTracksMatchedPID, 0, particleSpecies, centrality, zVtx);
652         
653 //      Printf("%d --> %d %d %d", particleSpecies, primMCParticles->GetEntries(), primRecoTracksMatched->GetEntries(), allRecoTracksMatched->GetEntries());
654
655         delete primMCParticles;
656         delete primRecoTracksMatched;
657         delete allRecoTracksMatched;
658       }
659       
660       TObjArray* fakeParticles = fAnalyseUE->GetFakeParticles(inputEvent, mc, kFALSE, -1, kTRUE);
661       CleanUp((TObjArray*) fakeParticles->At(0), mc, skipParticlesAbove);
662       CleanUp((TObjArray*) fakeParticles->At(1), mc, skipParticlesAbove);
663
664       fHistos->FillTrackingEfficiency(0, 0, 0, 0, 0, (TObjArray*) fakeParticles->At(2), 0, centrality, zVtx);
665       fHistos->FillFakePt(fakeParticles, centrality);
666 //       Printf(">>>>> %d %d %d fakes", ((TObjArray*) fakeParticles->At(0))->GetEntriesFast(), ((TObjArray*) fakeParticles->At(1))->GetEntriesFast(), ((TObjArray*) fakeParticles->At(2))->GetEntriesFast());
667       delete fakeParticles;
668     
669       // (MC-true all particles)
670       // STEP 2
671       if (!fReduceMemoryFootprint)
672         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepVertex, tracksMC, tracksCorrelateMC, weight);
673       else
674         fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
675       
676       // Get MC primaries that match reconstructed track
677       // triggers
678       tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, fParticleSpeciesTrigger, kTRUE);
679       CleanUp(tmpList, mc, skipParticlesAbove);
680       TObjArray* tracksRecoMatchedPrim = CloneAndReduceTrackList(tmpList);
681       delete tmpList;
682       
683       // associated
684       TObjArray* tracksCorrelateRecoMatchedPrim = tracksRecoMatchedPrim;
685       if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
686       {
687         tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, fParticleSpeciesAssociated, kTRUE);
688         CleanUp(tmpList, mc, skipParticlesAbove);
689         tracksCorrelateRecoMatchedPrim = CloneAndReduceTrackList(tmpList);
690         delete tmpList;
691       }
692
693       // (RECO-matched (quantities from MC particle) primary particles)
694       // STEP 4
695       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, tracksCorrelateRecoMatchedPrim, weight);
696
697       // mixed event
698       if (fFillMixed)
699       {
700         AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx + 200);
701         if (pool->IsReady())
702           for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
703             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
704         pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateRecoMatchedPrim));
705       }
706       
707       // Get MC primaries + secondaries that match reconstructed track
708       // triggers
709       tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, fParticleSpeciesTrigger, kTRUE);
710       CleanUp(tmpList, mc, skipParticlesAbove);
711       TObjArray* tracksRecoMatchedAll = CloneAndReduceTrackList(tmpList);
712       delete tmpList;
713       
714       // associated
715       TObjArray* tracksCorrelateRecoMatchedAll = tracksRecoMatchedAll;
716       if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
717       {
718         tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, fParticleSpeciesAssociated, kTRUE);
719         CleanUp(tmpList, mc, skipParticlesAbove);
720         tracksCorrelateRecoMatchedAll = CloneAndReduceTrackList(tmpList);
721         delete tmpList;
722       }
723       
724       // (RECO-matched (quantities from MC particle) all particles)
725       // STEP 5
726       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, tracksCorrelateRecoMatchedAll, weight);
727       
728       // mixed event
729       if (fFillMixed)
730       {
731         AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx + 300);
732         if (pool->IsReady())
733           for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
734             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
735         pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateRecoMatchedAll));
736       }
737       
738       // Get RECO tracks
739       // triggers
740       tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
741       CleanUp(tmpList, mc, skipParticlesAbove);
742       TObjArray* tracks = CloneAndReduceTrackList(tmpList);
743       delete tmpList;
744       
745       // associated
746       TObjArray* tracksCorrelate = tracks;
747       if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
748       {
749         tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
750         CleanUp(tmpList, mc, skipParticlesAbove);
751         tracksCorrelate = CloneAndReduceTrackList(tmpList);
752         delete tmpList;
753       }
754      
755       // (RECO all tracks)
756       // STEP 6
757       if (!fSkipStep6)
758         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, tracksCorrelate, weight);
759       
760       // two track cut, STEP 8
761       if (fTwoTrackEfficiencyCut > 0)
762         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, tracksCorrelate, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut);
763
764       // apply correction efficiency, STEP 10
765       if (fEfficiencyCorrectionTriggers || fEfficiencyCorrectionAssociated)
766       {
767           // with or without two track efficiency depending on if fTwoTrackEfficiencyCut is set
768         Bool_t twoTrackCut = (fTwoTrackEfficiencyCut > 0);
769         
770         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepCorrected, tracks, tracksCorrelate, weight, kTRUE, twoTrackCut, bSign, fTwoTrackEfficiencyCut, kTRUE);
771       }
772       
773       // mixed event
774       if (fFillMixed)
775       {
776         AliEventPool* pool2 = fPoolMgr->GetEventPool(centrality, zVtx + 100);
777         if (pool2->IsReady())
778         {
779           for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
780           {
781             // STEP 6
782             if (!fSkipStep6)
783               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0));
784             
785             // two track cut, STEP 8
786             if (fTwoTrackEfficiencyCut > 0)
787               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut);
788             
789             // apply correction efficiency, STEP 10
790             if (fEfficiencyCorrectionTriggers || fEfficiencyCorrectionAssociated)
791             {
792               // with or without two track efficiency depending on if fTwoTrackEfficiencyCut is set
793               Bool_t twoTrackCut = (fTwoTrackEfficiencyCut > 0);
794               
795               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepCorrected, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0), twoTrackCut, bSign, fTwoTrackEfficiencyCut, kTRUE);
796             }
797           }
798         }
799         pool2->UpdatePool(CloneAndReduceTrackList(tracksCorrelate));
800       }
801       
802       if (0 && !fReduceMemoryFootprint)
803       {
804         // make list of secondaries (matched with MC)
805         TObjArray* tracksRecoMatchedSecondaries = new TObjArray;
806         for (Int_t i=0; i<tracksRecoMatchedAll->GetEntriesFast(); i++)
807           if (((AliAODMCParticle*)tracksRecoMatchedAll->At(i))->IsPhysicalPrimary() == kFALSE)
808             tracksRecoMatchedSecondaries->Add(tracksRecoMatchedAll->At(i));
809       
810         // Study: Use only secondaries as trigger particles and plot the correlation vs. all particles; store in step 9
811         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy2, tracksRecoMatchedSecondaries, tracksRecoMatchedAll, weight);
812         
813         // Study: Use only primaries as trigger particles and plot the correlation vs. secondaries; store in step 8
814         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksRecoMatchedPrim, tracksRecoMatchedSecondaries, weight);
815       
816         // plot delta phi vs process id of secondaries
817         // trigger particles: primaries in 4 < pT < 10
818         // associated particles: secondaries in 1 < pT < 10
819         
820         for (Int_t i=0; i<tracksRecoMatchedPrim->GetEntriesFast(); i++)
821         {
822           AliVParticle* triggerParticle = (AliVParticle*) tracksRecoMatchedPrim->At(i);
823           
824           if (triggerParticle->Pt() < 4 || triggerParticle->Pt() > 10)
825             continue;
826           
827           for (Int_t j=0; j<tracksRecoMatchedSecondaries->GetEntriesFast(); j++)
828           {
829             AliAODMCParticle* particle = (AliAODMCParticle*) tracksRecoMatchedSecondaries->At(j);
830             
831             if (particle->Pt() < 1 || particle->Pt() > 10)
832               continue;
833             
834             if (particle->Pt() > triggerParticle->Pt())
835               continue;
836               
837             Double_t deltaPhi = triggerParticle->Phi() - particle->Phi();
838             if (deltaPhi > 1.5 * TMath::Pi()) 
839               deltaPhi -= TMath::TwoPi();
840             if (deltaPhi < -0.5 * TMath::Pi())
841               deltaPhi += TMath::TwoPi();
842               
843             Int_t processID = fMcEvent->Stack()->Particle(particle->GetLabel())->GetUniqueID();
844               
845             ((TH2F*) fListOfHistos->FindObject("processIDs"))->Fill(deltaPhi, processID);
846           }
847         }
848         
849         delete tracksRecoMatchedSecondaries;
850       }
851   
852       if (tracksCorrelateRecoMatchedPrim != tracksRecoMatchedPrim)
853         delete tracksCorrelateRecoMatchedPrim;
854       delete tracksRecoMatchedPrim;
855
856       if (tracksCorrelateRecoMatchedAll != tracksRecoMatchedAll)
857         delete tracksCorrelateRecoMatchedAll;
858       delete tracksRecoMatchedAll;
859       
860       if (tracksCorrelate != tracks)
861         delete tracksCorrelate;
862       delete tracks;
863     }
864   }
865   
866   if (tracksMC != tracksCorrelateMC)
867     delete tracksCorrelateMC;
868   delete tracksMC;
869 }
870
871 //____________________________________________________________________
872 void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
873 {
874
875   // Run the analysis on DATA or MC to get raw distributions
876  
877   if ( fDebug > 3 ) AliInfo( " Processing event in Data mode ..." );
878
879   ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(0);
880
881   if (!fInputHandler)
882     return;
883     
884   // skip not selected events here (the AOD is not updated for those)
885   if (!fSkipTrigger && !(fInputHandler->IsEventSelected() & fSelectBit))
886     return;
887   
888   // skip fast cluster events here if requested
889   if (fSkipFastCluster && (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly))
890     return;
891  
892   // Support for ESD and AOD based analysis
893   AliVEvent* inputEvent = fAOD;
894   if (!inputEvent)
895     inputEvent = fESD;
896
897   Double_t centrality = 0;
898   
899   AliCentrality *centralityObj = 0;
900   if (fCentralityMethod.Length() > 0)
901   {
902     if (fCentralityMethod == "ZNA_MANUAL")
903     {
904       Bool_t zna = kFALSE;
905       for(Int_t j = 0; j < 4; ++j) {
906         if (fESD->GetZDCData()->GetZDCTDCData(12,j) != 0) {
907           zna = kTRUE;
908         }
909       }
910
911 //       Printf("%d %f", zna, fZNAtower[0]);
912       if (zna)
913       {
914         // code from Chiara O (23.10.12)
915         const Double_t *fZNAtower = fESD->GetZDCData()->GetZN2TowerEnergy();
916         Float_t znacut[4] = {681., 563., 413., 191.};
917         
918         if(fZNAtower[0]>znacut[0]) centrality = 1;
919         else if(fZNAtower[0]>znacut[1]) centrality = 21;
920         else if(fZNAtower[0]>znacut[2]) centrality = 41;
921         else if(fZNAtower[0]>znacut[3]) centrality = 61;
922         else centrality = 81;
923       }
924       else
925         centrality = -1;
926     }
927     else if (fCentralityMethod == "TRACKS_MANUAL")
928     {
929       // for pp
930       TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
931       centrality = tracks->GetEntriesFast();
932       if (centrality > 40)
933         centrality = 41;
934 //       Printf("%d %f", tracks->GetEntriesFast(), centrality);
935       delete tracks;
936     }
937     else if (fCentralityMethod == "V0A_MANUAL")
938     {
939       // for pp
940       
941       //Total multiplicity in the VZERO A detector
942       Float_t MV0A=inputEvent->GetVZEROData()->GetMTotV0A();
943       Float_t MV0AScaled=0.;
944       if (fMap){
945         TParameter<float>* sf=(TParameter<float>*)fMap->GetValue(Form("%d",inputEvent->GetRunNumber()));
946         if(sf)MV0AScaled=MV0A*sf->GetVal();
947       }
948       
949       if (MV0AScaled > 0)
950         centrality = MV0AScaled;
951       else
952         centrality = -1;
953     }
954     else
955     {
956       if (fAOD)
957         centralityObj = fAOD->GetHeader()->GetCentralityP();
958       else if (fESD)
959         centralityObj = fESD->GetCentrality();
960       
961       if (centralityObj)
962         centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
963         //centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
964       else
965         centrality = -1;
966
967       if (fAOD)
968       {
969         // remove outliers
970         if (centrality == 0)
971         {
972           if (fAOD->GetVZEROData())
973           {
974             Float_t multV0 = 0;
975             for (Int_t i=0; i<64; i++)
976               multV0 += fAOD->GetVZEROData()->GetMultiplicity(i);
977             if (multV0 < 19500)
978             {
979               centrality = -1;
980               AliInfo("Rejecting event due to too small V0 multiplicity");
981             }
982           }
983         }
984       }
985     }
986     
987     AliInfo(Form("Centrality is %f", centrality));
988   }
989   
990   Float_t bSign = (inputEvent->GetMagneticField() > 0) ? 1 : -1;
991
992   fHistos->SetRunNumber(inputEvent->GetRunNumber());
993   
994   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
995   fHistos->FillEvent(centrality, AliUEHist::kCFStepAll);
996   
997   // Trigger selection ************************************************
998   if (!fSkipTrigger && !fAnalyseUE->TriggerSelection(fInputHandler)) return;
999
1000   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
1001   fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
1002   
1003   // Pileup selection ************************************************
1004   if (fAnalysisUtils && fAnalysisUtils->IsPileUpEvent(inputEvent)) 
1005   {
1006     // count the removed events
1007     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1008
1009     return;
1010   }
1011   
1012   // Vertex selection *************************************************
1013   if(!fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex)) return;
1014   
1015   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
1016   fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
1017  
1018   // fill V0 control histograms
1019   if (fCentralityMethod == "V0A_MANUAL")
1020   {
1021     ((TH2F*) fListOfHistos->FindObject("V0AMult"))->Fill(inputEvent->GetVZEROData()->GetMTotV0A(), centrality);
1022     if (fAOD)
1023       ((TH2F*) fListOfHistos->FindObject("V0AMultCorrelation"))->Fill(inputEvent->GetVZEROData()->GetMTotV0A(), fAOD->GetTracklets()->GetNumberOfTracklets());
1024   }
1025     
1026   // optimization
1027   if (centrality < 0)
1028     return;
1029
1030   TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
1031   //Printf("Accepted %d tracks", tracks->GetEntries());
1032   
1033   // check for outlier in centrality vs number of tracks (rough constants extracted from correlation histgram)
1034   Bool_t reject = kFALSE;
1035   if (fRejectCentralityOutliers)
1036   {
1037     if (centrality > 40 && centrality <= 50 && tracks->GetEntriesFast() > 1160)
1038       reject = kTRUE;
1039     if (centrality > 50 && centrality <= 60 && tracks->GetEntriesFast() > 650)
1040       reject = kTRUE;
1041     if (centrality > 60 && centrality <= 70 && tracks->GetEntriesFast() > 370)
1042       reject = kTRUE;
1043     if (centrality > 70 && centrality <= 80 && tracks->GetEntriesFast() > 220)
1044       reject = kTRUE;
1045     if (centrality > 80 && centrality <= 90 && tracks->GetEntriesFast() > 130)
1046       reject = kTRUE;
1047     if (centrality > 90 && tracks->GetEntriesFast() > 75)
1048       reject = kTRUE;
1049   }
1050   
1051   if (reject)
1052   {
1053     AliInfo(Form("Rejecting event due to centrality vs tracks correlation: %f %d", centrality, tracks->GetEntriesFast()));
1054     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1055     delete tracks;
1056     return;
1057   }
1058   
1059   if (fRejectZeroTrackEvents && tracks->GetEntriesFast() == 0)
1060   {
1061     AliInfo(Form("Rejecting event because it has no tracks: %f %d", centrality, tracks->GetEntriesFast()));
1062     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1063     delete tracks;
1064     return;
1065   }
1066
1067   // correlate particles with...
1068   TObjArray* tracksCorrelate = 0;
1069   if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
1070     tracksCorrelate = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
1071   
1072   // reference multiplicity
1073   Int_t referenceMultiplicity = -1;
1074   if (fESD)
1075     referenceMultiplicity = AliESDtrackCuts::GetReferenceMultiplicity(fESD);
1076   else if (fAOD)
1077     referenceMultiplicity = tracks->GetEntriesFast(); // TODO to be replaced by the estimator once available in the AOD
1078 //    referenceMultiplicity = fAOD->GetHeader()->GetRefMultiplicityComb05();
1079
1080   ((TH2F*) fListOfHistos->FindObject("referenceMultiplicity"))->Fill(centrality, referenceMultiplicity);
1081   
1082   const AliVVertex* vertex = inputEvent->GetPrimaryVertex();
1083   Double_t zVtx = vertex->GetZ();
1084   
1085   Float_t weight = 1;
1086   if (fFillpT)
1087     weight = -1;
1088
1089   // Fill containers at STEP 6 (reconstructed)
1090   if (centrality >= 0)
1091   {
1092     if (!fSkipStep6)
1093       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, tracksCorrelate, weight, kTRUE, kFALSE, 0, 0.02, kTRUE);
1094     
1095     ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(1);
1096     
1097     if (fTwoTrackEfficiencyCut > 0)
1098       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, tracksCorrelate, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut, kTRUE);
1099   }
1100
1101   // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
1102   TObjArray* tracksClone = CloneAndReduceTrackList(tracks);
1103   delete tracks;
1104   
1105   if (fFillMixed)
1106   {
1107     // event mixing
1108     
1109     // 1. First get an event pool corresponding in mult (cent) and
1110     //    zvertex to the current event. Once initialized, the pool
1111     //    should contain nMix (reduced) events. This routine does not
1112     //    pre-scan the chain. The first several events of every chain
1113     //    will be skipped until the needed pools are filled to the
1114     //    specified depth. If the pool categories are not too rare, this
1115     //    should not be a problem. If they are rare, you could lose
1116     //    statistics.
1117
1118     // 2. Collect the whole pool's content of tracks into one TObjArray
1119     //    (bgTracks), which is effectively a single background super-event.
1120
1121     // 3. The reduced and bgTracks arrays must both be passed into
1122     //    FillCorrelations(). Also nMix should be passed in, so a weight
1123     //    of 1./nMix can be applied.
1124
1125     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
1126     
1127     if (!pool)
1128       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, zVtx));
1129     
1130 //     pool->SetDebug(1);
1131      
1132     if (pool->IsReady()) 
1133     {
1134       
1135       Int_t nMix = pool->GetCurrentNEvents();
1136 //       cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
1137       
1138       ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
1139       ((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
1140       if (pool->IsReady())
1141         ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
1142     
1143       // Fill mixed-event histos here  
1144       for (Int_t jMix=0; jMix<nMix; jMix++) 
1145       {
1146         TObjArray* bgTracks = pool->GetEvent(jMix);
1147         
1148         if (!fSkipStep6)
1149           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kFALSE, 0, 0.02, kTRUE);
1150
1151         if (fTwoTrackEfficiencyCut > 0)
1152           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut, kTRUE);
1153       }
1154     }
1155     
1156     // ownership is with the pool now
1157     if (tracksCorrelate)
1158     {
1159       pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelate));
1160       delete tracksClone;
1161     }
1162     else
1163       pool->UpdatePool(tracksClone);
1164     //pool->PrintInfo();
1165   }
1166   else
1167   {
1168     delete tracksClone;
1169     if (tracksCorrelate)
1170       delete tracksCorrelate;
1171   }
1172 }
1173
1174 TObjArray* AliAnalysisTaskPhiCorrelations::CloneAndReduceTrackList(TObjArray* tracks)
1175 {
1176   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
1177   
1178   TObjArray* tracksClone = new TObjArray;
1179   tracksClone->SetOwner(kTRUE);
1180   
1181   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
1182   {
1183     AliVParticle* particle = (AliVParticle*) tracks->UncheckedAt(i);
1184     AliDPhiBasicParticle* copy = new AliDPhiBasicParticle(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge());
1185     copy->SetUniqueID(particle->GetUniqueID());
1186     tracksClone->Add(copy);
1187   }
1188   
1189   return tracksClone;
1190 }
1191
1192 //____________________________________________________________________
1193 void  AliAnalysisTaskPhiCorrelations::Initialize()
1194 {
1195   // input handler
1196   fInputHandler = (AliInputEventHandler*)
1197          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1198   // MC handler
1199   fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1200 }
1201
1202 //____________________________________________________________________
1203 void AliAnalysisTaskPhiCorrelations::RemoveDuplicates(TObjArray* tracks)
1204 {
1205   // remove particles with the same label
1206   
1207   Int_t before = tracks->GetEntriesFast();
1208
1209   for (Int_t i=0; i<before; ++i) 
1210   {
1211     AliVParticle* part = (AliVParticle*) tracks->At(i);
1212     
1213     for (Int_t j=i+1; j<before; ++j) 
1214     {
1215       AliVParticle* part2 = (AliVParticle*) tracks->At(j);
1216       
1217       if (part->GetLabel() == part2->GetLabel())
1218       {
1219         Printf("Removing %d with label %d (duplicated in %d)", i, part->GetLabel(), j); part->Dump(); part2->Dump();
1220         TObject* object = tracks->RemoveAt(i);
1221         if (tracks->IsOwner())
1222           delete object;
1223         break;
1224       }
1225     }
1226   }
1227  
1228   tracks->Compress();
1229   
1230   if (before > tracks->GetEntriesFast())
1231     AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast())); 
1232 }
1233
1234 void AliAnalysisTaskPhiCorrelations::CleanUp(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
1235 {
1236   // calls RemoveInjectedSignals, RemoveWeakDecays and RemoveDuplicates
1237   
1238   if (!tracks)
1239     return;
1240   
1241   if (fInjectedSignals)
1242     fAnalyseUE->RemoveInjectedSignals(tracks, mcObj, maxLabel);
1243   if (fRemoveWeakDecays)
1244     fAnalyseUE->RemoveWeakDecays(tracks, mcObj);
1245   if (fRemoveDuplicates)
1246     RemoveDuplicates(tracks);
1247 }
1248
1249 //____________________________________________________________________
1250 void AliAnalysisTaskPhiCorrelations::SelectCharge(TObjArray* tracks)
1251 {
1252   // remove particles with charge not selected (depending on fTriggerSelectCharge)
1253   
1254   if (!tracks)
1255     return;
1256   
1257   Int_t before = tracks->GetEntriesFast();
1258
1259   for (Int_t i=0; i<before; ++i) 
1260   {
1261     AliVParticle* part = (AliVParticle*) tracks->At(i);
1262     
1263     if (part->Charge() * fTriggerSelectCharge < -1)
1264     {
1265 //       Printf("Removing %d with charge %d", i, part->Charge());
1266       TObject* object = tracks->RemoveAt(i);
1267       if (tracks->IsOwner())
1268         delete object;
1269     }
1270   }
1271  
1272   tracks->Compress();
1273   
1274   if (before > tracks->GetEntriesFast())
1275     AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast())); 
1276 }