Added Eventplane Dependence in dPhi Correlations code
[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 <TInterpreter.h>
20 #include <TChain.h>
21 #include <TFile.h>
22 #include <TList.h>
23 #include <TMath.h>
24 #include <TTree.h>
25 #include <TH2F.h>
26 #include <TH3F.h>
27 #include <TRandom.h>
28 #include <TParameter.h>
29
30 #include "AliAnalysisTaskPhiCorrelations.h"
31 #include "AliAnalyseLeadingTrackUE.h"
32 #include "AliUEHistograms.h"
33 #include "AliUEHist.h"
34
35 #include "AliAnalysisManager.h"
36 #include "AliAODHandler.h"
37 #include "AliAODInputHandler.h"
38 #include "AliAODMCParticle.h"
39 #include "AliInputEventHandler.h"
40 #include "AliLog.h"
41 #include "AliMCEventHandler.h"
42 #include "AliVParticle.h"
43 #include "AliCFContainer.h"
44 #include "AliEventplane.h"
45
46 #include "AliESDEvent.h"
47 #include "AliESDInputHandler.h"
48 #include "AliMultiplicity.h"
49 #include "AliCentrality.h"
50 #include "AliStack.h"
51 #include "AliAODMCHeader.h"
52 #include "AliGenCocktailEventHeader.h"
53 #include "AliGenEventHeader.h"
54 #include "AliCollisionGeometry.h"
55
56 #include "AliEventPoolManager.h"
57
58 #include "AliESDZDC.h"
59 #include "AliESDtrackCuts.h"
60
61 #include "AliHelperPID.h"
62 #include "AliAnalysisUtils.h"
63 #include "TMap.h"
64
65 ////////////////////////////////////////////////////////////////////////
66 //
67 // Analysis class for azimuthal correlation studies
68 // Based on the UE task from Sara Vallero and Jan Fiete
69 //
70 // This class needs input AODs.
71 // The output is a list of analysis-specific containers.
72 //
73 // The AOD can be either connected to the InputEventHandler  
74 // for a chain of AOD files 
75 // or 
76 // to the OutputEventHandler
77 // for a chain of ESD files,
78 // in this case the class should be in the train after the jet-finder
79 //
80 //    Authors:
81 //    Jan Fiete Grosse-Oetringhaus
82 // 
83 ////////////////////////////////////////////////////////////////////////
84
85
86 ClassImp( AliAnalysisTaskPhiCorrelations )
87 ClassImp( AliDPhiBasicParticle )
88
89 //____________________________________________________________________
90 AliAnalysisTaskPhiCorrelations:: AliAnalysisTaskPhiCorrelations(const char* name):
91 AliAnalysisTask(name,""),
92 // general configuration
93 fDebug(0),
94 fMode(0),
95 fReduceMemoryFootprint(kFALSE),
96 fFillMixed(kTRUE),
97 fMixingTracks(50000),
98 fTwoTrackEfficiencyStudy(kFALSE),
99 fTwoTrackEfficiencyCut(0),
100 fTwoTrackCutMinRadius(0.8),
101 fUseVtxAxis(kFALSE),
102 fCourseCentralityBinning(kFALSE),
103 fSkipTrigger(kFALSE),
104 fInjectedSignals(kFALSE),
105 fRandomizeReactionPlane(kFALSE),
106 fHelperPID(0x0),
107 fAnalysisUtils(0x0),
108 fMap(0x0),
109 // pointers to UE classes
110 fAnalyseUE(0x0),
111 fHistos(0x0),
112 fHistosMixed(0),
113 fEfficiencyCorrectionTriggers(0),
114 fEfficiencyCorrectionAssociated(0),
115 fCentralityWeights(0),
116 // handlers and events
117 fAOD(0x0),
118 fESD(0x0),
119 fArrayMC(0x0),
120 fInputHandler(0x0),
121 fMcEvent(0x0),
122 fMcHandler(0x0),
123 fPoolMgr(0x0),
124 // histogram settings
125 fListOfHistos(0x0), 
126 // event QA
127 fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default) 
128 fZVertex(7.),
129 fAcceptOnlyMuEvents(kFALSE),
130 fCentralityMethod("V0M"),
131 // track cuts
132 fTrackEtaCut(0.8),
133 fTrackEtaCutMin(-1.),
134 fTrackPhiCutEvPlMin(0.),
135 fTrackPhiCutEvPlMax(0.),
136 fOnlyOneEtaSide(0),
137 fPtMin(0.5),
138 fDCAXYCut(0),
139 fSharedClusterCut(-1),
140 fCrossedRowsCut(-1),
141 fFoundFractionCut(-1),
142 fFilterBit(0xFF),
143 fTrackStatus(0),
144 fSelectBit(AliVEvent::kMB|AliVEvent::kUserDefined),
145 fUseChargeHadrons(kFALSE),
146 fParticleSpeciesTrigger(-1),
147 fParticleSpeciesAssociated(-1),
148 fCheckMotherPDG(kTRUE),
149 fTrackletDphiCut(9999999.),
150 fSelectCharge(0),
151 fTriggerSelectCharge(0),
152 fAssociatedSelectCharge(0),
153 fTriggerRestrictEta(-1),
154 fEtaOrdering(kFALSE),
155 fCutConversions(kFALSE),
156 fCutResonances(kFALSE),
157 fRejectResonanceDaughters(-1),
158 fFillOnlyStep0(kFALSE),
159 fSkipStep6(kFALSE),
160 fRejectCentralityOutliers(kFALSE),
161 fRejectZeroTrackEvents(kFALSE),
162 fRemoveWeakDecays(kFALSE),
163 fRemoveDuplicates(kFALSE),
164 fSkipFastCluster(kFALSE),
165 fWeightPerEvent(kFALSE),
166 fCustomBinning(),
167 fPtOrder(kTRUE),
168 fTriggersFromDetector(0),
169 fAssociatedFromDetector(0),
170 fMCUseUncheckedCentrality(kFALSE),
171 fFillpT(kFALSE)
172 {
173   // Default constructor
174   // Define input and output slots here
175   // Input slot #0 works with a TChain
176   DefineInput(0, TChain::Class());
177   // Output slot #0 writes into a TList container
178   DefineOutput(0, TList::Class());
179 }
180
181 AliAnalysisTaskPhiCorrelations::~AliAnalysisTaskPhiCorrelations() 
182
183   // destructor
184   
185   if (fListOfHistos  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) 
186     delete fListOfHistos;
187 }
188
189 //____________________________________________________________________
190 void AliAnalysisTaskPhiCorrelations::ConnectInputData(Option_t* /*option*/)
191 {
192   
193   // Connect the input data  
194   if (fDebug > 1) AliInfo("ConnectInputData() ");
195   
196   // Since AODs can either be connected to the InputEventHandler
197   // or to the OutputEventHandler ( the AOD is created by a previus task in the train )
198   // we need to get the pointer to the AODEvent correctly.
199   
200   // Delta AODs are also accepted.
201   
202   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
203   
204   if( handler && handler->InheritsFrom("AliAODInputHandler") ) 
205   { // input AOD
206     fAOD = ((AliAODInputHandler*)handler)->GetEvent();
207     if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
208   } 
209   else 
210   {  //output AOD
211     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
212     if (handler && handler->InheritsFrom("AliAODHandler") ) 
213     {
214       fAOD = ((AliAODHandler*)handler)->GetAOD();
215       if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
216     } 
217     else 
218     {  // no AOD
219       AliWarning("I can't get any AOD Event Handler");
220     }
221   }
222   
223   if (handler && handler->InheritsFrom("AliESDInputHandler") ) 
224   { // input ESD
225     // pointer received per event in ::Exec
226     if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliESDInputHandler");
227   } 
228   
229   // Initialize common pointers
230   Initialize();
231 }
232
233 //____________________________________________________________________
234 void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
235 {
236   // Create the output container
237   
238   if (fDebug > 1) AliInfo("CreateOutputObjects()");
239    
240   // Initialize class with main algorithms, event and track selection. 
241   fAnalyseUE = new AliAnalyseLeadingTrackUE();
242   fAnalyseUE->SetParticleSelectionCriteria(fFilterBit, fUseChargeHadrons, fTrackEtaCut, fTrackEtaCutMin, fPtMin);
243   fAnalyseUE->SetDCAXYCut(fDCAXYCut);
244   fAnalyseUE->SetSharedClusterCut(fSharedClusterCut);
245   fAnalyseUE->SetCrossedRowsCut(fCrossedRowsCut);
246   fAnalyseUE->SetFoundFractionCut(fFoundFractionCut);
247   fAnalyseUE->SetTrackStatus(fTrackStatus);
248   fAnalyseUE->SetCheckMotherPDG(fCheckMotherPDG);
249   fAnalyseUE->SetDebug(fDebug); 
250   fAnalyseUE->DefineESDCuts(fFilterBit);
251   fAnalyseUE->SetEventSelection(fSelectBit);
252   fAnalyseUE->SetHelperPID(fHelperPID);
253   if(fTrackPhiCutEvPlMax!=0.)
254     fAnalyseUE->SetParticlePhiCutEventPlane(fTrackPhiCutEvPlMin,fTrackPhiCutEvPlMax);
255   if ((fParticleSpeciesTrigger != -1 || fParticleSpeciesAssociated != -1) && !fHelperPID)
256     AliFatal("HelperPID object should be set in the steering macro");
257
258   // Initialize output list of containers
259   if (fListOfHistos != NULL){
260         delete fListOfHistos;
261         fListOfHistos = NULL;
262         }
263   if (!fListOfHistos){
264         fListOfHistos = new TList();
265         fListOfHistos->SetOwner(kTRUE); 
266         }
267
268   // Initialize class to handle histograms 
269   TString histType = "4R";
270   if (fUseVtxAxis == 1)
271     histType = "5R";
272   else if (fUseVtxAxis == 2)
273     histType = "6R";
274   if (fCourseCentralityBinning)
275     histType += "C";
276   fHistos = new AliUEHistograms("AliUEHistogramsSame", histType, fCustomBinning);
277   fHistosMixed = new AliUEHistograms("AliUEHistogramsMixed", histType, fCustomBinning);
278   
279   fHistos->SetSelectCharge(fSelectCharge);
280   fHistosMixed->SetSelectCharge(fSelectCharge);
281   
282   fHistos->SetSelectTriggerCharge(fTriggerSelectCharge);
283   fHistosMixed->SetSelectTriggerCharge(fTriggerSelectCharge);
284   
285   fHistos->SetSelectAssociatedCharge(fAssociatedSelectCharge);
286   fHistosMixed->SetSelectAssociatedCharge(fAssociatedSelectCharge);
287
288   fHistos->SetTriggerRestrictEta(fTriggerRestrictEta);
289   fHistosMixed->SetTriggerRestrictEta(fTriggerRestrictEta);
290   
291   fHistos->SetOnlyOneEtaSide(fOnlyOneEtaSide);
292   fHistosMixed->SetOnlyOneEtaSide(fOnlyOneEtaSide);
293   
294   fHistos->SetEtaOrdering(fEtaOrdering);
295   fHistosMixed->SetEtaOrdering(fEtaOrdering);
296
297   fHistos->SetPairCuts(fCutConversions, fCutResonances);
298   fHistosMixed->SetPairCuts(fCutConversions, fCutResonances);
299   
300   fHistos->SetRejectResonanceDaughters(fRejectResonanceDaughters);
301   fHistosMixed->SetRejectResonanceDaughters(fRejectResonanceDaughters);
302   
303   fHistos->SetTrackEtaCut(fTrackEtaCut);
304   fHistosMixed->SetTrackEtaCut(fTrackEtaCut);
305   
306   fHistos->SetWeightPerEvent(fWeightPerEvent);
307   fHistosMixed->SetWeightPerEvent(fWeightPerEvent);
308   
309   fHistos->SetPtOrder(fPtOrder);
310   fHistosMixed->SetPtOrder(fPtOrder);
311   
312   fHistos->SetTwoTrackCutMinRadius(fTwoTrackCutMinRadius);
313   fHistosMixed->SetTwoTrackCutMinRadius(fTwoTrackCutMinRadius);
314   
315   if (fEfficiencyCorrectionTriggers)
316    {
317     fHistos->SetEfficiencyCorrectionTriggers(fEfficiencyCorrectionTriggers);
318     fHistosMixed->SetEfficiencyCorrectionTriggers((THnF*) fEfficiencyCorrectionTriggers->Clone());
319    }
320   if (fEfficiencyCorrectionAssociated)
321   {
322     fHistos->SetEfficiencyCorrectionAssociated(fEfficiencyCorrectionAssociated);
323     fHistosMixed->SetEfficiencyCorrectionAssociated((THnF*) fEfficiencyCorrectionAssociated->Clone());
324   }
325   
326   // add histograms to list
327   fListOfHistos->Add(fHistos);
328   fListOfHistos->Add(fHistosMixed);
329   // add HelperPID to list
330   if (fHelperPID)
331     fListOfHistos->Add(fHelperPID);
332   // add TMap to list
333   if (fMap)
334     fListOfHistos->Add(fMap);
335   
336   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));
337   fListOfHistos->Add(new TH1F("eventStat", ";;events", 4, -0.5, 3.5));
338   fListOfHistos->Add(new TH2F("mixedDist", ";centrality;tracks;events", 101, 0, 101, 200, 0, fMixingTracks * 1.5));
339   fListOfHistos->Add(new TH2F("referenceMultiplicity", ";centrality;tracks;events", 101, 0, 101, 200, 0, 200));
340   if (fCentralityMethod == "V0A_MANUAL")
341   {
342     fListOfHistos->Add(new TH2F("V0AMult", "V0A multiplicity;V0A multiplicity;V0A multiplicity (scaled)", 1000, -.5, 999.5, 1000, -.5, 999.5));
343     fListOfHistos->Add(new TH2F("V0AMultCorrelation", "V0A multiplicity;V0A multiplicity;SPD tracklets", 1000, -.5, 999.5, 1000, -.5, 999.5));
344   }
345   if (fTriggersFromDetector == 1 || fTriggersFromDetector == 2 || fAssociatedFromDetector == 1 || fAssociatedFromDetector == 2)
346     fListOfHistos->Add(new TH1F("V0SingleCells", "V0 single cell multiplicity;multiplicity;events", 100, -0.5, 99.5));
347   if (fTriggersFromDetector == 3 || fAssociatedFromDetector == 3)
348     fListOfHistos->Add(new TH1F("DphiTrklets", "tracklets Dphi;#Delta#phi,trklets (mrad);entries", 100, -100, 100));
349   
350   PostData(0,fListOfHistos);
351   
352   // Add task configuration to output list 
353   AddSettingsTree();
354
355   // event mixing
356   Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemention of AliEventPoolManager
357    
358   Int_t nCentralityBins  = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(1);
359   Double_t* centralityBins = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray();
360   
361   const Int_t kNZvtxBins  = 10+(1+10)*4;
362   // bins for further buffers are shifted by 100 cm
363   Double_t vertexBins[kNZvtxBins+1] = { -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10, 
364                                        90,  92,  94,  96,  98, 100, 102, 104, 106, 108, 110, 
365                                       190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210, 
366                                       290, 292, 294, 296, 298, 300, 302, 304, 306, 308, 310, 
367                                       390, 392, 394, 396, 398, 400, 402, 404, 406, 408, 410 };
368
369   Int_t nZvtxBins  = kNZvtxBins;
370   Double_t* zvtxbin = vertexBins;
371
372   if (fMode == 0 && fHistos->GetUEHist(2)->GetEventHist()->GetNVar() > 2)
373   {
374     nZvtxBins = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(2);
375     zvtxbin = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetXbins()->GetArray();
376   }
377
378   fPoolMgr = new AliEventPoolManager(poolsize, fMixingTracks, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
379   fPoolMgr->SetTargetValues(fMixingTracks, 0.1, 5);
380 }
381
382 //____________________________________________________________________
383 void  AliAnalysisTaskPhiCorrelations::Exec(Option_t */*option*/)
384 {
385   // exec (per event)
386   fAnalyseUE->NextEvent();
387   
388   // receive ESD pointer if we are not running AOD analysis
389   if (!fAOD)
390   {
391     AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
392     if (handler && handler->InheritsFrom("AliESDInputHandler"))
393       fESD = ((AliESDInputHandler*)handler)->GetEvent();
394   }
395
396   if (fMode)
397   {
398     // correction mode
399     
400     if (fMcHandler)
401       fMcEvent = fMcHandler->MCEvent();
402
403     if (fAOD)
404     {
405       // array of MC particles
406       fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
407       if (!fArrayMC)
408         AliFatal("No array of MC particles found !!!");
409     }
410     
411     AnalyseCorrectionMode();
412   }
413   else AnalyseDataMode();
414 }
415
416 /******************** ANALYSIS METHODS *****************************/
417
418 //____________________________________________________________________
419 void  AliAnalysisTaskPhiCorrelations::AddSettingsTree()
420 {
421   //Write settings to output list
422   TTree *settingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
423   settingsTree->Branch("fnTracksVertex", &fnTracksVertex,"nTracksVertex/I");
424   settingsTree->Branch("fZVertex", &fZVertex,"ZVertex/D");
425   settingsTree->Branch("fAcceptOnlyMuEvents", &fAcceptOnlyMuEvents,"AcceptOnlyMuEvents/O");
426   //settingsTree->Branch("fCentralityMethod", fCentralityMethod.Data(),"CentralityMethod/C");
427   settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
428   settingsTree->Branch("fTrackEtaCutMin", &fTrackEtaCutMin, "TrackEtaCutMin/D");
429   settingsTree->Branch("fOnlyOneEtaSide", &fOnlyOneEtaSide,"OnlyOneEtaSide/I");
430   settingsTree->Branch("fPtMin", &fPtMin, "PtMin/D");
431   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
432   settingsTree->Branch("fSharedClusterCut", &fSharedClusterCut,"SharedClusterCut/D");
433   settingsTree->Branch("fCrossedRowsCut", &fCrossedRowsCut,"CrossedRowsCut/I");
434   settingsTree->Branch("fFoundFractionCut", &fFoundFractionCut,"FoundFractionCut/D");
435   settingsTree->Branch("fTrackStatus", &fTrackStatus,"TrackStatus/I");
436   settingsTree->Branch("fSelectBit", &fSelectBit,"EventSelectionBit/I");
437   settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
438   settingsTree->Branch("fParticleSpeciesTrigger", &fParticleSpeciesTrigger,"ParticleSpeciesTrigger/I");
439   settingsTree->Branch("fParticleSpeciesAssociated", &fParticleSpeciesAssociated,"ParticleSpeciesAssociated/I");
440   settingsTree->Branch("fCheckMotherPDG", &fCheckMotherPDG,"CheckMotherPDG/I");
441   settingsTree->Branch("fSelectCharge", &fSelectCharge,"SelectCharge/I");
442   settingsTree->Branch("fTriggerSelectCharge", &fTriggerSelectCharge,"TriggerSelectCharge/I");
443   settingsTree->Branch("fAssociatedSelectCharge", &fAssociatedSelectCharge,"fAssociatedSelectCharge/I");
444   settingsTree->Branch("fTriggerRestrictEta", &fTriggerRestrictEta,"TriggerRestrictEta/D");
445   settingsTree->Branch("fEtaOrdering", &fEtaOrdering,"EtaOrdering/O");
446   settingsTree->Branch("fCutConversions", &fCutConversions,"CutConversions/O");
447   settingsTree->Branch("fCutResonances", &fCutResonances,"CutResonances/O");
448   settingsTree->Branch("fRejectResonanceDaughters", &fRejectResonanceDaughters,"RejectResonanceDaughters/I");
449   settingsTree->Branch("fFillpT", &fFillpT,"FillpT/O");
450   settingsTree->Branch("fMixingTracks", &fMixingTracks,"MixingTracks/I");
451   settingsTree->Branch("fSkipTrigger", &fSkipTrigger,"SkipTrigger/O");
452   settingsTree->Branch("fInjectedSignals", &fInjectedSignals,"InjectedSignals/O");
453   settingsTree->Branch("fRandomizeReactionPlane", &fRandomizeReactionPlane,"RandomizeReactionPlane/O");
454   settingsTree->Branch("fRejectCentralityOutliers", &fRejectCentralityOutliers,"RejectCentralityOutliers/O");
455   settingsTree->Branch("fRejectZeroTrackEvents", &fRejectZeroTrackEvents,"RejectZeroTrackEvents/O");
456   settingsTree->Branch("fRemoveWeakDecays", &fRemoveWeakDecays,"RemoveWeakDecays/O");
457   settingsTree->Branch("fRemoveDuplicates", &fRemoveDuplicates,"RemoveDuplicates/O");
458   settingsTree->Branch("fSkipFastCluster", &fSkipFastCluster,"SkipFastCluster/O");
459   settingsTree->Branch("fWeightPerEvent", &fWeightPerEvent,"WeightPerEvent/O");
460   settingsTree->Branch("fPtOrder", &fPtOrder,"PtOrder/O");
461   settingsTree->Branch("fTriggersFromDetector", &fTriggersFromDetector,"TriggersFromDetector/I");
462   settingsTree->Branch("fAssociatedFromDetector", &fAssociatedFromDetector,"AssociatedFromDetector/I");
463   settingsTree->Branch("fMCUseUncheckedCentrality", &fMCUseUncheckedCentrality,"MCUseUncheckedCentrality/O");
464   settingsTree->Branch("fTwoTrackEfficiencyCut", &fTwoTrackEfficiencyCut,"TwoTrackEfficiencyCut/D");
465   settingsTree->Branch("fTwoTrackCutMinRadius", &fTwoTrackCutMinRadius,"TwoTrackCutMinRadius/D");
466   
467   //fCustomBinning
468   
469   settingsTree->Fill();
470   fListOfHistos->Add(settingsTree);
471 }  
472
473 //____________________________________________________________________
474 void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
475 {
476   // Run the analysis on MC to get the correction maps
477   //
478  
479   if ( fDebug > 3 ) AliInfo( " Processing event in Corrections mode ..." );
480   
481   Double_t centrality = 0;
482   
483   if (fCentralityMethod.Length() > 0)
484   {
485     if (fCentralityMethod == "MC_b")
486     {
487       AliGenEventHeader* eventHeader = GetFirstHeader();
488       if (!eventHeader)
489       {
490         // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing 
491         // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
492         AliError("Event header not found. Skipping this event.");
493         fHistos->FillEvent(0, AliUEHist::kCFStepAnaTopology);
494         return;
495       }
496       
497       AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
498       if (!collGeometry)
499       {
500         eventHeader->Dump();
501         AliFatal("Asking for MC_b centrality, but event header has no collision geometry information");
502       }
503       
504       centrality = collGeometry->ImpactParameter();
505     }
506     else if (fCentralityMethod == "nano")
507     {
508       centrality = (Float_t) gROOT->ProcessLine(Form("100.0 + 100.0 * ((AliNanoAODHeader*) %p)->GetCentrality(\"%s\")", fAOD->GetHeader(), fCentralityMethod.Data())) / 100 - 1.0;
509     }
510     else
511     {
512       AliCentrality *centralityObj = 0;
513       if (fAOD)
514         centralityObj = fAOD->GetHeader()->GetCentralityP();
515       else if (fESD)
516         centralityObj = fESD->GetCentrality();
517       
518       if (centralityObj)
519       {
520         if (fMCUseUncheckedCentrality)
521           centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
522         else
523           centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
524       }
525       else
526       {
527         Printf("WARNING: Centrality object is 0");
528         centrality = -1;
529       }
530     }
531
532     AliInfo(Form("Centrality is %f", centrality));
533   }
534   
535   // Support for ESD and AOD based analysis
536   AliVEvent* inputEvent = fAOD;
537   if (!inputEvent)
538     inputEvent = fESD;
539
540   Float_t bSign = 0;
541   
542   if (inputEvent)
543   {
544     fHistos->SetRunNumber(inputEvent->GetRunNumber());
545     bSign = (inputEvent->GetMagneticField() > 0) ? 1 : -1;
546   }
547     
548   TObject* mc = fArrayMC;
549   if (!mc)
550     mc = fMcEvent;
551   
552   // count all events
553   fHistos->FillEvent(centrality, -1);
554   
555   if (centrality < 0)
556     return;
557
558   // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
559   TObject* vertexSupplier = fMcEvent;
560   if (fAOD) // AOD
561     vertexSupplier = fAOD->FindListObject(AliAODMCHeader::StdBranchName());
562     
563   if (!fAnalyseUE->VertexSelection(vertexSupplier, 0, fZVertex)) 
564     return;
565     
566     Float_t zVtx = 0;
567   if (fAOD)
568     zVtx = ((AliAODMCHeader*) vertexSupplier)->GetVtxZ();
569   else
570     zVtx = fMcEvent->GetPrimaryVertex()->GetZ();
571   Float_t weight = 1;
572   if (fFillpT)
573     weight = -1;
574   
575   // For productions with injected signals, figure out above which label to skip particles/tracks
576   Int_t skipParticlesAbove = 0;
577   if (fInjectedSignals)
578   {
579     AliGenEventHeader* eventHeader = 0;
580     Int_t headers = 0;
581
582     if (fMcEvent)
583     {
584       // ESD
585       AliHeader* header = (AliHeader*) fMcEvent->Header();
586       if (!header)
587         AliFatal("fInjectedSignals set but no MC header found");
588         
589       AliGenCocktailEventHeader* cocktailHeader = dynamic_cast<AliGenCocktailEventHeader*> (header->GenEventHeader());
590       if (!cocktailHeader)
591       {
592         header->Dump();
593         AliFatal("fInjectedSignals set but no MC cocktail header found");
594       }
595
596       headers = cocktailHeader->GetHeaders()->GetEntries();
597       eventHeader = dynamic_cast<AliGenEventHeader*> (cocktailHeader->GetHeaders()->First());
598       
599       if (fDebug > 4)
600       {
601         for (Int_t i=0; i<cocktailHeader->GetHeaders()->GetEntries(); i++)
602         {
603           AliGenEventHeader* headerTmp = dynamic_cast<AliGenEventHeader*> (cocktailHeader->GetHeaders()->At(i));
604           if (headerTmp)
605             Printf("%d particles in header:", headerTmp->NProduced());
606           cocktailHeader->GetHeaders()->At(i)->Dump();
607         }
608       }
609     }
610     else
611     {
612       // AOD
613       AliAODMCHeader* header = (AliAODMCHeader*) fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
614       if (!header)
615         AliFatal("fInjectedSignals set but no MC header found");
616       
617       headers = header->GetNCocktailHeaders();
618       eventHeader = header->GetCocktailHeader(0);
619     }
620     
621     if (!eventHeader)
622     {
623       // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing 
624       // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
625       AliError("First event header not found. Skipping this event.");
626       fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
627       return;
628     }
629     
630     skipParticlesAbove = eventHeader->NProduced();
631     AliInfo(Form("Injected signals in this event (%d headers). Keeping particles/tracks of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
632   }
633   
634   if (fCentralityWeights && !AcceptEventCentralityWeight(centrality))
635   {
636     AliInfo(Form("Rejecting event because of centrality weighting: %f", centrality));
637     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
638     return;
639   }
640   
641   // Get MC primaries
642   // triggers
643   TObjArray* tmpList = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
644   CleanUp(tmpList, mc, skipParticlesAbove);
645   TObjArray* tracksMC = CloneAndReduceTrackList(tmpList);
646   delete tmpList;
647   
648   // associated
649   TObjArray* tracksCorrelateMC = tracksMC;
650   if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
651   {
652     tmpList = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
653     CleanUp(tmpList, mc, skipParticlesAbove);
654     tracksCorrelateMC = CloneAndReduceTrackList(tmpList);
655     delete tmpList;
656   }
657   
658   if (fRandomizeReactionPlane)
659   {
660     Double_t centralityDigits = centrality*1000. - (Int_t)(centrality*1000.);
661     Double_t angle = TMath::TwoPi() * centralityDigits;
662     AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
663     ShiftTracks(tracksMC, angle);
664     if (tracksCorrelateMC != tracksMC)
665       ShiftTracks(tracksCorrelateMC, angle);
666   }
667   
668   if (fFillOnlyStep0)
669     zVtx = 0;
670   
671   // Event selection based on number of number of MC particles
672   if (fRejectZeroTrackEvents && tracksMC->GetEntriesFast() == 0)
673   {
674     AliInfo(Form("Rejecting event due to kinematic selection: %f %d", centrality, tracksMC->GetEntriesFast()));
675     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
676     if (tracksMC != tracksCorrelateMC)
677       delete tracksCorrelateMC;
678     delete tracksMC;
679     return;
680   }
681   
682   // (MC-true all particles)
683   // STEP 0
684   fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, tracksCorrelateMC, weight);
685   
686   // mixed event
687   if (fFillMixed)
688   {
689     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
690     if (fFillOnlyStep0)
691       ((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
692     if (pool->IsReady())
693       for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
694         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
695     pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateMC));
696   }
697   
698 //   Printf("trigger: %d", ((AliInputEventHandler*)fInputHandler)->IsEventSelected());
699   
700   // Trigger selection ************************************************
701   if (!fFillOnlyStep0 && (fSkipTrigger || fAnalyseUE->TriggerSelection(fInputHandler)))
702   {  
703     // (MC-true all particles)
704     // STEP 1
705     if (!fReduceMemoryFootprint)
706       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTriggered, tracksMC, tracksCorrelateMC, weight);
707     else
708       fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
709       
710     if (!inputEvent) {
711       AliFatal("UNEXPECTED: inputEvent is 0. Trigger selection should have failed");
712       return;
713     }
714     
715     // Vertex selection *************************************************
716     if (fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex))
717     {
718       // fill here for tracking efficiency
719       // loop over particle species
720       
721       for (Int_t particleSpecies = 0; particleSpecies < 4; particleSpecies++)
722       {
723         TObjArray* primMCParticles = fAnalyseUE->GetAcceptedParticles(mc, 0x0, kTRUE, particleSpecies, kTRUE);
724         TObjArray* primRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE, kFALSE);
725         TObjArray* allRecoTracksMatched  = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE, kFALSE);
726         TObjArray* primRecoTracksMatchedPID = 0;
727         TObjArray* allRecoTracksMatchedPID  = 0;
728         
729         if (fHelperPID)
730         {
731           primRecoTracksMatchedPID = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE, kTRUE);
732           allRecoTracksMatchedPID  = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE, kTRUE);
733         }
734         
735         CleanUp(primMCParticles, mc, skipParticlesAbove);
736         CleanUp(primRecoTracksMatched, mc, skipParticlesAbove);
737         CleanUp(allRecoTracksMatched, mc, skipParticlesAbove);
738         CleanUp(primRecoTracksMatchedPID, mc, skipParticlesAbove);
739         CleanUp(allRecoTracksMatchedPID, mc, skipParticlesAbove);
740         
741         // select charges
742         if (fTriggerSelectCharge != 0)
743         {
744           SelectCharge(primMCParticles);
745           SelectCharge(primRecoTracksMatched);
746           SelectCharge(allRecoTracksMatched);
747           SelectCharge(primRecoTracksMatchedPID);
748           SelectCharge(allRecoTracksMatchedPID);
749         }
750       
751         fHistos->FillTrackingEfficiency(primMCParticles, primRecoTracksMatched, allRecoTracksMatched, primRecoTracksMatchedPID, allRecoTracksMatchedPID, 0, particleSpecies, centrality, zVtx);
752         
753 //      Printf("%d --> %d %d %d", particleSpecies, primMCParticles->GetEntries(), primRecoTracksMatched->GetEntries(), allRecoTracksMatched->GetEntries());
754
755         delete primMCParticles;
756         delete primRecoTracksMatched;
757         delete allRecoTracksMatched;
758       }
759       
760       TObjArray* fakeParticles = fAnalyseUE->GetFakeParticles(inputEvent, mc, kFALSE, -1, kTRUE);
761       CleanUp((TObjArray*) fakeParticles->At(0), mc, skipParticlesAbove);
762       CleanUp((TObjArray*) fakeParticles->At(1), mc, skipParticlesAbove);
763
764       fHistos->FillTrackingEfficiency(0, 0, 0, 0, 0, (TObjArray*) fakeParticles->At(2), 0, centrality, zVtx);
765       fHistos->FillFakePt(fakeParticles, centrality);
766 //       Printf(">>>>> %d %d %d fakes", ((TObjArray*) fakeParticles->At(0))->GetEntriesFast(), ((TObjArray*) fakeParticles->At(1))->GetEntriesFast(), ((TObjArray*) fakeParticles->At(2))->GetEntriesFast());
767       delete fakeParticles;
768     
769       // (MC-true all particles)
770       // STEP 2
771       if (!fReduceMemoryFootprint)
772         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepVertex, tracksMC, tracksCorrelateMC, weight);
773       else
774         fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
775       
776       // Get MC primaries that match reconstructed track
777       // triggers
778       tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, fParticleSpeciesTrigger, kTRUE);
779       CleanUp(tmpList, mc, skipParticlesAbove);
780       TObjArray* tracksRecoMatchedPrim = CloneAndReduceTrackList(tmpList);
781       delete tmpList;
782       
783       // associated
784       TObjArray* tracksCorrelateRecoMatchedPrim = tracksRecoMatchedPrim;
785       if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
786       {
787         tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, fParticleSpeciesAssociated, kTRUE);
788         CleanUp(tmpList, mc, skipParticlesAbove);
789         tracksCorrelateRecoMatchedPrim = CloneAndReduceTrackList(tmpList);
790         delete tmpList;
791       }
792
793       // (RECO-matched (quantities from MC particle) primary particles)
794       // STEP 4
795       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, tracksCorrelateRecoMatchedPrim, weight);
796
797       // mixed event
798       if (fFillMixed)
799       {
800         AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx + 200);
801         if (pool->IsReady())
802           for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
803             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
804         pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateRecoMatchedPrim));
805       }
806       
807       // Get MC primaries + secondaries that match reconstructed track
808       // triggers
809       tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, fParticleSpeciesTrigger, kTRUE);
810       CleanUp(tmpList, mc, skipParticlesAbove);
811       TObjArray* tracksRecoMatchedAll = CloneAndReduceTrackList(tmpList);
812       delete tmpList;
813       
814       // associated
815       TObjArray* tracksCorrelateRecoMatchedAll = tracksRecoMatchedAll;
816       if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
817       {
818         tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, fParticleSpeciesAssociated, kTRUE);
819         CleanUp(tmpList, mc, skipParticlesAbove);
820         tracksCorrelateRecoMatchedAll = CloneAndReduceTrackList(tmpList);
821         delete tmpList;
822       }
823       
824       // (RECO-matched (quantities from MC particle) all particles)
825       // STEP 5
826       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, tracksCorrelateRecoMatchedAll, weight);
827       
828       // mixed event
829       if (fFillMixed)
830       {
831         AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx + 300);
832         if (pool->IsReady())
833           for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
834             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
835         pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateRecoMatchedAll));
836       }
837       
838       // Get RECO tracks
839       // triggers
840       tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
841       CleanUp(tmpList, mc, skipParticlesAbove);
842       TObjArray* tracks = CloneAndReduceTrackList(tmpList);
843       delete tmpList;
844       
845       // associated
846       TObjArray* tracksCorrelate = tracks;
847       if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
848       {
849         tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
850         CleanUp(tmpList, mc, skipParticlesAbove);
851         tracksCorrelate = CloneAndReduceTrackList(tmpList);
852         delete tmpList;
853       }
854      
855       // (RECO all tracks)
856       // STEP 6
857       if (!fSkipStep6)
858         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, tracksCorrelate, weight);
859       
860       // two track cut, STEP 8
861       if (fTwoTrackEfficiencyCut > 0)
862         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, tracksCorrelate, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut);
863
864       // apply correction efficiency, STEP 10
865       if (fEfficiencyCorrectionTriggers || fEfficiencyCorrectionAssociated)
866       {
867           // with or without two track efficiency depending on if fTwoTrackEfficiencyCut is set
868         Bool_t twoTrackCut = (fTwoTrackEfficiencyCut > 0);
869         
870         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepCorrected, tracks, tracksCorrelate, weight, kTRUE, twoTrackCut, bSign, fTwoTrackEfficiencyCut, kTRUE);
871       }
872       
873       // mixed event
874       if (fFillMixed)
875       {
876         AliEventPool* pool2 = fPoolMgr->GetEventPool(centrality, zVtx + 100);
877         ((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool2->NTracksInPool());
878         if (pool2->IsReady())
879         {
880           for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
881           {
882             // STEP 6
883             if (!fSkipStep6)
884               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0));
885             
886             // two track cut, STEP 8
887             if (fTwoTrackEfficiencyCut > 0)
888               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut);
889             
890             // apply correction efficiency, STEP 10
891             if (fEfficiencyCorrectionTriggers || fEfficiencyCorrectionAssociated)
892             {
893               // with or without two track efficiency depending on if fTwoTrackEfficiencyCut is set
894               Bool_t twoTrackCut = (fTwoTrackEfficiencyCut > 0);
895               
896               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepCorrected, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0), twoTrackCut, bSign, fTwoTrackEfficiencyCut, kTRUE);
897             }
898           }
899         }
900         pool2->UpdatePool(CloneAndReduceTrackList(tracksCorrelate));
901       }
902       
903       if (0 && !fReduceMemoryFootprint)
904       {
905         // make list of secondaries (matched with MC)
906         TObjArray* tracksRecoMatchedSecondaries = new TObjArray;
907         for (Int_t i=0; i<tracksRecoMatchedAll->GetEntriesFast(); i++)
908           if (((AliAODMCParticle*)tracksRecoMatchedAll->At(i))->IsPhysicalPrimary() == kFALSE)
909             tracksRecoMatchedSecondaries->Add(tracksRecoMatchedAll->At(i));
910       
911         // Study: Use only secondaries as trigger particles and plot the correlation vs. all particles; store in step 9
912         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy2, tracksRecoMatchedSecondaries, tracksRecoMatchedAll, weight);
913         
914         // Study: Use only primaries as trigger particles and plot the correlation vs. secondaries; store in step 8
915         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksRecoMatchedPrim, tracksRecoMatchedSecondaries, weight);
916       
917         // plot delta phi vs process id of secondaries
918         // trigger particles: primaries in 4 < pT < 10
919         // associated particles: secondaries in 1 < pT < 10
920         
921         for (Int_t i=0; i<tracksRecoMatchedPrim->GetEntriesFast(); i++)
922         {
923           AliVParticle* triggerParticle = (AliVParticle*) tracksRecoMatchedPrim->At(i);
924           
925           if (triggerParticle->Pt() < 4 || triggerParticle->Pt() > 10)
926             continue;
927           
928           for (Int_t j=0; j<tracksRecoMatchedSecondaries->GetEntriesFast(); j++)
929           {
930             AliAODMCParticle* particle = (AliAODMCParticle*) tracksRecoMatchedSecondaries->At(j);
931             
932             if (particle->Pt() < 1 || particle->Pt() > 10)
933               continue;
934             
935             if (particle->Pt() > triggerParticle->Pt())
936               continue;
937               
938             Double_t deltaPhi = triggerParticle->Phi() - particle->Phi();
939             if (deltaPhi > 1.5 * TMath::Pi()) 
940               deltaPhi -= TMath::TwoPi();
941             if (deltaPhi < -0.5 * TMath::Pi())
942               deltaPhi += TMath::TwoPi();
943               
944             Int_t processID = fMcEvent->Stack()->Particle(particle->GetLabel())->GetUniqueID();
945               
946             ((TH2F*) fListOfHistos->FindObject("processIDs"))->Fill(deltaPhi, processID);
947           }
948         }
949         
950         delete tracksRecoMatchedSecondaries;
951       }
952   
953       if (tracksCorrelateRecoMatchedPrim != tracksRecoMatchedPrim)
954         delete tracksCorrelateRecoMatchedPrim;
955       delete tracksRecoMatchedPrim;
956
957       if (tracksCorrelateRecoMatchedAll != tracksRecoMatchedAll)
958         delete tracksCorrelateRecoMatchedAll;
959       delete tracksRecoMatchedAll;
960       
961       if (tracksCorrelate != tracks)
962         delete tracksCorrelate;
963       delete tracks;
964     }
965   }
966   
967   if (tracksMC != tracksCorrelateMC)
968     delete tracksCorrelateMC;
969   delete tracksMC;
970 }
971
972 //____________________________________________________________________
973 AliGenEventHeader* AliAnalysisTaskPhiCorrelations::GetFirstHeader()
974 {
975   // get first MC header from either ESD/AOD (including cocktail header if available)
976   
977   if (fMcEvent)
978   {
979     // ESD
980     AliHeader* header = (AliHeader*) fMcEvent->Header();
981     if (!header)
982       return 0;
983       
984     AliGenCocktailEventHeader* cocktailHeader = dynamic_cast<AliGenCocktailEventHeader*> (header->GenEventHeader());
985     if (cocktailHeader)
986       return dynamic_cast<AliGenEventHeader*> (cocktailHeader->GetHeaders()->First());
987
988     return dynamic_cast<AliGenEventHeader*> (header->GenEventHeader());
989   }
990   else
991   {
992     // AOD
993     AliAODMCHeader* header = (AliAODMCHeader*) fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
994     if (!header)
995       return 0;
996     
997     return header->GetCocktailHeader(0);
998   }
999 }
1000
1001 //____________________________________________________________________
1002 void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
1003 {
1004
1005   // Run the analysis on DATA or MC to get raw distributions
1006  
1007   if ( fDebug > 3 ) AliInfo( " Processing event in Data mode ..." );
1008
1009   ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(0);
1010
1011   if (!fInputHandler)
1012     return;
1013     
1014   // skip not selected events here (the AOD is not updated for those)
1015   if (!fSkipTrigger && !(fInputHandler->IsEventSelected() & fSelectBit))
1016     return;
1017   
1018   // skip fast cluster events here if requested
1019   if (fSkipFastCluster && (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly))
1020     return;
1021  
1022   // Support for ESD and AOD based analysis
1023   AliVEvent* inputEvent = fAOD;
1024   if (!inputEvent)
1025     inputEvent = fESD;
1026
1027   Double_t centrality = 0;
1028   
1029   AliCentrality *centralityObj = 0;
1030   if (fCentralityMethod.Length() > 0)
1031   {
1032     if (fCentralityMethod == "ZNA_MANUAL")
1033     {
1034       Bool_t zna = kFALSE;
1035       for(Int_t j = 0; j < 4; ++j) {
1036         if (fESD->GetZDCData()->GetZDCTDCData(12,j) != 0) {
1037           zna = kTRUE;
1038         }
1039       }
1040
1041 //       Printf("%d %f", zna, fZNAtower[0]);
1042       if (zna)
1043       {
1044         // code from Chiara O (23.10.12)
1045         const Double_t *fZNAtower = fESD->GetZDCData()->GetZN2TowerEnergy();
1046         Float_t znacut[4] = {681., 563., 413., 191.};
1047         
1048         if(fZNAtower[0]>znacut[0]) centrality = 1;
1049         else if(fZNAtower[0]>znacut[1]) centrality = 21;
1050         else if(fZNAtower[0]>znacut[2]) centrality = 41;
1051         else if(fZNAtower[0]>znacut[3]) centrality = 61;
1052         else centrality = 81;
1053       }
1054       else
1055         centrality = -1;
1056     }
1057     else if (fCentralityMethod == "TRACKS_MANUAL")
1058     {
1059       // for pp
1060       TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
1061       centrality = tracks->GetEntriesFast();
1062       if (centrality > 40)
1063         centrality = 41;
1064 //       Printf("%d %f", tracks->GetEntriesFast(), centrality);
1065       delete tracks;
1066     }
1067     else if (fCentralityMethod == "V0A_MANUAL")
1068     {
1069       // for pp
1070       
1071       //Total multiplicity in the VZERO A detector
1072       Float_t MV0A=inputEvent->GetVZEROData()->GetMTotV0A();
1073       Float_t MV0AScaled=0.;
1074       if (fMap){
1075         TParameter<float>* sf=(TParameter<float>*)fMap->GetValue(Form("%d",inputEvent->GetRunNumber()));
1076         if(sf)MV0AScaled=MV0A*sf->GetVal();
1077       }
1078       
1079       if (MV0AScaled > 0)
1080         centrality = MV0AScaled;
1081       else
1082         centrality = -1;
1083     }
1084     else if (fCentralityMethod == "nano")
1085     {
1086 //       fAOD->GetHeader()->Dump();
1087 //       Printf("%p %p %d", dynamic_cast<AliNanoAODHeader*> (fAOD->GetHeader()), dynamic_cast<AliNanoAODHeader*> ((TObject*) (fAOD->GetHeader())), fAOD->GetHeader()->InheritsFrom("AliNanoAODHeader"));
1088
1089       Int_t error = 0;
1090       centrality = (Float_t) gROOT->ProcessLine(Form("100.0 + 100.0 * ((AliNanoAODHeader*) %p)->GetCentrality(\"%s\")", fAOD->GetHeader(), fCentralityMethod.Data()), &error) / 100 - 1.0;
1091       if (error != TInterpreter::kNoError)
1092         centrality = -1;
1093     }
1094     else
1095     {
1096       if (fAOD)
1097         centralityObj = fAOD->GetHeader()->GetCentralityP();
1098       else if (fESD)
1099         centralityObj = fESD->GetCentrality();
1100       
1101       if (centralityObj)
1102         centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
1103         //centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
1104       else
1105         centrality = -1;
1106       
1107       if (fAOD)
1108       {
1109         // remove outliers
1110         if (centrality == 0)
1111         {
1112           if (fAOD->GetVZEROData())
1113           {
1114             Float_t multV0 = 0;
1115             for (Int_t i=0; i<64; i++)
1116               multV0 += fAOD->GetVZEROData()->GetMultiplicity(i);
1117             if (multV0 < 19500)
1118             {
1119               centrality = -1;
1120               AliInfo("Rejecting event due to too small V0 multiplicity");
1121             }
1122           }
1123         }
1124       }
1125     }
1126     
1127     AliInfo(Form("Centrality is %f", centrality));
1128   }
1129   
1130   Float_t bSign = (inputEvent->GetMagneticField() > 0) ? 1 : -1;
1131
1132   fHistos->SetRunNumber(inputEvent->GetRunNumber());
1133   
1134   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
1135   fHistos->FillEvent(centrality, AliUEHist::kCFStepAll);
1136   
1137   // Trigger selection ************************************************
1138   if (!fSkipTrigger && !fAnalyseUE->TriggerSelection(fInputHandler)) return;
1139
1140   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
1141   fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
1142   
1143   // Pileup selection ************************************************
1144   if (fAnalysisUtils && fAnalysisUtils->IsPileUpEvent(inputEvent)) 
1145   {
1146     // count the removed events
1147     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1148
1149     return;
1150   }
1151   
1152   // Reject events without a muon in the muon arm ************************************************
1153   if(fAcceptOnlyMuEvents && !IsMuEvent())return;
1154   
1155   // Vertex selection *************************************************
1156   if(!fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex)) return;
1157   
1158   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
1159   fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
1160  
1161   // fill V0 control histograms
1162   if (fCentralityMethod == "V0A_MANUAL")
1163   {
1164     ((TH2F*) fListOfHistos->FindObject("V0AMult"))->Fill(inputEvent->GetVZEROData()->GetMTotV0A(), centrality);
1165     if (fAOD)
1166       ((TH2F*) fListOfHistos->FindObject("V0AMultCorrelation"))->Fill(inputEvent->GetVZEROData()->GetMTotV0A(), fAOD->GetTracklets()->GetNumberOfTracklets());
1167   }
1168     
1169   // optimization
1170   if (centrality < 0)
1171     return;
1172
1173   TObjArray* tracks = 0;
1174  
1175   Double_t evtPlanePhi = -999.; //A value outside [-pi/2,pi/2] will be ignored
1176   if(fTrackPhiCutEvPlMax!=0.) {
1177     AliEventplane* evtPlane = inputEvent->GetEventplane();
1178     Double_t qx = 0; Double_t qy = 0;
1179     if(evtPlane) evtPlanePhi = evtPlane->CalculateVZEROEventPlane(inputEvent, 10, 2, qx, qy);
1180     //Reject event if the plane is not available
1181     else return; 
1182   }
1183  
1184   if (fTriggersFromDetector == 0)
1185     tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesTrigger, kTRUE, kTRUE, evtPlanePhi);
1186   else if (fTriggersFromDetector <= 4)
1187     tracks=GetParticlesFromDetector(inputEvent,fTriggersFromDetector);
1188   else
1189     AliFatal(Form("Invalid setting for fTriggersFromDetector: %d", fTriggersFromDetector));
1190   
1191   //Printf("Accepted %d tracks", tracks->GetEntries());
1192   
1193   // check for outlier in centrality vs number of tracks (rough constants extracted from correlation histgram)
1194   Bool_t reject = kFALSE;
1195   if (fRejectCentralityOutliers)
1196   {
1197     if (centrality > 40 && centrality <= 50 && tracks->GetEntriesFast() > 1160)
1198       reject = kTRUE;
1199     if (centrality > 50 && centrality <= 60 && tracks->GetEntriesFast() > 650)
1200       reject = kTRUE;
1201     if (centrality > 60 && centrality <= 70 && tracks->GetEntriesFast() > 370)
1202       reject = kTRUE;
1203     if (centrality > 70 && centrality <= 80 && tracks->GetEntriesFast() > 220)
1204       reject = kTRUE;
1205     if (centrality > 80 && centrality <= 90 && tracks->GetEntriesFast() > 130)
1206       reject = kTRUE;
1207     if (centrality > 90 && tracks->GetEntriesFast() > 75)
1208       reject = kTRUE;
1209   }
1210   
1211   if (reject)
1212   {
1213     AliInfo(Form("Rejecting event due to centrality vs tracks correlation: %f %d", centrality, tracks->GetEntriesFast()));
1214     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1215     delete tracks;
1216     return;
1217   }
1218   
1219   if (fRejectZeroTrackEvents && tracks->GetEntriesFast() == 0)
1220   {
1221     AliInfo(Form("Rejecting event because it has no tracks: %f %d", centrality, tracks->GetEntriesFast()));
1222     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1223     delete tracks;
1224     return;
1225   }
1226   
1227   if (fCentralityWeights && !AcceptEventCentralityWeight(centrality))
1228   {
1229     AliInfo(Form("Rejecting event because of centrality weighting: %f", centrality));
1230     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1231     delete tracks;
1232     return;
1233   }
1234   
1235   // correlate particles with...
1236   TObjArray* tracksCorrelate = 0;
1237   if(fAssociatedFromDetector==0){
1238     if (fParticleSpeciesAssociated != fParticleSpeciesTrigger || fTriggersFromDetector > 0 || fTrackPhiCutEvPlMax != 0.)
1239       tracksCorrelate = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
1240   }
1241   else if (fAssociatedFromDetector <= 4){
1242     if(fAssociatedFromDetector != fTriggersFromDetector)
1243       tracksCorrelate = GetParticlesFromDetector(inputEvent,fAssociatedFromDetector);
1244   }
1245   else
1246     AliFatal(Form("Invalid setting for fAssociatedFromDetector: %d", fAssociatedFromDetector));
1247   
1248   // reference multiplicity
1249   Int_t referenceMultiplicity = -1;
1250   if (fESD)
1251     referenceMultiplicity = AliESDtrackCuts::GetReferenceMultiplicity(fESD);
1252   else if (fAOD)
1253     referenceMultiplicity = tracks->GetEntriesFast(); // TODO to be replaced by the estimator once available in the AOD
1254 //    referenceMultiplicity = fAOD->GetHeader()->GetRefMultiplicityComb05();
1255
1256   ((TH2F*) fListOfHistos->FindObject("referenceMultiplicity"))->Fill(centrality, referenceMultiplicity);
1257   
1258   const AliVVertex* vertex = inputEvent->GetPrimaryVertex();
1259   Double_t zVtx = vertex->GetZ();
1260   
1261   Float_t weight = 1;
1262   if (fFillpT)
1263     weight = -1;
1264
1265   // Fill containers at STEP 6 (reconstructed)
1266   if (centrality >= 0)
1267   {
1268     if (!fSkipStep6)
1269       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, tracksCorrelate, weight, kTRUE, kFALSE, 0, 0.02, kTRUE);
1270     
1271     ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(1);
1272     
1273     if (fTwoTrackEfficiencyCut > 0)
1274       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, tracksCorrelate, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut, kTRUE);
1275   }
1276
1277   // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
1278   TObjArray* tracksClone = CloneAndReduceTrackList(tracks);
1279   delete tracks;
1280   
1281   if (fFillMixed)
1282   {
1283     // event mixing
1284     
1285     // 1. First get an event pool corresponding in mult (cent) and
1286     //    zvertex to the current event. Once initialized, the pool
1287     //    should contain nMix (reduced) events. This routine does not
1288     //    pre-scan the chain. The first several events of every chain
1289     //    will be skipped until the needed pools are filled to the
1290     //    specified depth. If the pool categories are not too rare, this
1291     //    should not be a problem. If they are rare, you could lose
1292     //    statistics.
1293
1294     // 2. Collect the whole pool's content of tracks into one TObjArray
1295     //    (bgTracks), which is effectively a single background super-event.
1296
1297     // 3. The reduced and bgTracks arrays must both be passed into
1298     //    FillCorrelations(). Also nMix should be passed in, so a weight
1299     //    of 1./nMix can be applied.
1300
1301     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
1302     
1303     if (!pool)
1304       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, zVtx));
1305     
1306 //     pool->SetDebug(1);
1307      
1308     if (pool->IsReady()) 
1309     {
1310       
1311       Int_t nMix = pool->GetCurrentNEvents();
1312 //       cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
1313       
1314       ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
1315       ((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
1316       if (pool->IsReady())
1317         ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
1318     
1319       // Fill mixed-event histos here  
1320       for (Int_t jMix=0; jMix<nMix; jMix++) 
1321       {
1322         TObjArray* bgTracks = pool->GetEvent(jMix);
1323         
1324         if (!fSkipStep6)
1325           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kFALSE, 0, 0.02, kTRUE);
1326
1327         if (fTwoTrackEfficiencyCut > 0)
1328           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut, kTRUE);
1329       }
1330     }
1331     
1332     // ownership is with the pool now
1333     if (tracksCorrelate)
1334     {
1335       pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelate));
1336       delete tracksClone;
1337     }
1338     else
1339       pool->UpdatePool(tracksClone);
1340     //pool->PrintInfo();
1341   }
1342   else
1343   {
1344     delete tracksClone;
1345     if (tracksCorrelate)
1346       delete tracksCorrelate;
1347   }
1348 }
1349
1350 TObjArray* AliAnalysisTaskPhiCorrelations::CloneAndReduceTrackList(TObjArray* tracks)
1351 {
1352   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
1353   
1354   TObjArray* tracksClone = new TObjArray;
1355   tracksClone->SetOwner(kTRUE);
1356   
1357   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
1358   {
1359     AliVParticle* particle = (AliVParticle*) tracks->UncheckedAt(i);
1360     AliDPhiBasicParticle* copy = new AliDPhiBasicParticle(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge());
1361     copy->SetUniqueID(particle->GetUniqueID());
1362     tracksClone->Add(copy);
1363   }
1364   
1365   return tracksClone;
1366 }
1367
1368 //____________________________________________________________________
1369 void  AliAnalysisTaskPhiCorrelations::Initialize()
1370 {
1371   // input handler
1372   fInputHandler = (AliInputEventHandler*)
1373          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1374   // MC handler
1375   fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1376 }
1377
1378 //____________________________________________________________________
1379 void AliAnalysisTaskPhiCorrelations::RemoveDuplicates(TObjArray* tracks)
1380 {
1381   // remove particles with the same label
1382   
1383   Int_t before = tracks->GetEntriesFast();
1384
1385   for (Int_t i=0; i<before; ++i) 
1386   {
1387     AliVParticle* part = (AliVParticle*) tracks->At(i);
1388     
1389     for (Int_t j=i+1; j<before; ++j) 
1390     {
1391       AliVParticle* part2 = (AliVParticle*) tracks->At(j);
1392       
1393       if (part->GetLabel() == part2->GetLabel())
1394       {
1395         Printf("Removing %d with label %d (duplicated in %d)", i, part->GetLabel(), j); part->Dump(); part2->Dump();
1396         TObject* object = tracks->RemoveAt(i);
1397         if (tracks->IsOwner())
1398           delete object;
1399         break;
1400       }
1401     }
1402   }
1403  
1404   tracks->Compress();
1405   
1406   if (before > tracks->GetEntriesFast())
1407     AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast())); 
1408 }
1409
1410 void AliAnalysisTaskPhiCorrelations::CleanUp(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
1411 {
1412   // calls RemoveInjectedSignals, RemoveWeakDecays and RemoveDuplicates
1413   
1414   if (!tracks)
1415     return;
1416   
1417   if (fInjectedSignals)
1418     fAnalyseUE->RemoveInjectedSignals(tracks, mcObj, maxLabel);
1419   if (fRemoveWeakDecays)
1420     fAnalyseUE->RemoveWeakDecays(tracks, mcObj);
1421   if (fRemoveDuplicates)
1422     RemoveDuplicates(tracks);
1423 }
1424
1425 //____________________________________________________________________
1426 void AliAnalysisTaskPhiCorrelations::SelectCharge(TObjArray* tracks)
1427 {
1428   // remove particles with charge not selected (depending on fTriggerSelectCharge)
1429   
1430   if (!tracks)
1431     return;
1432   
1433   Int_t before = tracks->GetEntriesFast();
1434
1435   for (Int_t i=0; i<before; ++i) 
1436   {
1437     AliVParticle* part = (AliVParticle*) tracks->At(i);
1438     
1439     if (part->Charge() * fTriggerSelectCharge < -1)
1440     {
1441 //       Printf("Removing %d with charge %d", i, part->Charge());
1442       TObject* object = tracks->RemoveAt(i);
1443       if (tracks->IsOwner())
1444         delete object;
1445     }
1446   }
1447  
1448   tracks->Compress();
1449   
1450   if (before > tracks->GetEntriesFast())
1451     AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast())); 
1452 }
1453
1454 //____________________________________________________________________
1455 Bool_t AliAnalysisTaskPhiCorrelations::AcceptEventCentralityWeight(Double_t centrality)
1456 {
1457   // rejects "randomly" events such that the centrality gets flat
1458   // uses fCentralityWeights histogram
1459
1460   // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
1461   
1462   Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
1463   Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
1464   
1465   Bool_t result = kFALSE;
1466   if (centralityDigits < weight) 
1467     result = kTRUE;
1468   
1469   AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
1470   
1471   return result;
1472 }
1473
1474 //____________________________________________________________________
1475 void AliAnalysisTaskPhiCorrelations::ShiftTracks(TObjArray* tracks, Double_t angle)
1476 {
1477   // shifts the phi angle of all tracks by angle
1478   // 0 <= angle <= 2pi
1479   
1480   for (Int_t i=0; i<tracks->GetEntriesFast(); ++i) 
1481   {
1482     AliDPhiBasicParticle* part = (AliDPhiBasicParticle*) tracks->At(i);
1483     Double_t newAngle = part->Phi() + angle; 
1484     if (newAngle >= TMath::TwoPi())
1485       newAngle -= TMath::TwoPi();
1486     
1487     part->SetPhi(newAngle);
1488   }
1489 }
1490
1491 //____________________________________________________________________
1492 TObjArray* AliAnalysisTaskPhiCorrelations::GetParticlesFromDetector(AliVEvent* inputEvent, Int_t idet)
1493 {
1494   //1 = VZERO_A; 2 = VZERO_C; 3 = SPD tracklets
1495   TObjArray* obj = new TObjArray;
1496   obj->SetOwner(kTRUE);
1497   
1498   if (idet == 1 || idet == 2)
1499   {
1500     AliVVZERO* vZero = inputEvent->GetVZEROData();
1501     
1502     const Int_t vZeroStart = (idet == 1) ? 32 : 0;
1503     
1504     TH1F* singleCells = (TH1F*) fListOfHistos->FindObject("V0SingleCells");
1505     for (Int_t i=vZeroStart; i<vZeroStart+32; i++)
1506       {
1507         Float_t weight = vZero->GetMultiplicity(i);
1508         singleCells->Fill(weight);
1509         
1510         // rough estimate of multiplicity
1511         for (Int_t j=0; j<TMath::Nint(weight); j++)
1512           {
1513             AliDPhiBasicParticle* particle = new AliDPhiBasicParticle((AliVVZERO::GetVZEROEtaMax(i) + AliVVZERO::GetVZEROEtaMin(i)) / 2, AliVVZERO::GetVZEROAvgPhi(i), 1.1, 0); // fit pT = 1.1 and charge = 0
1514             particle->SetUniqueID(fAnalyseUE->GetEventCounter()* 100000 + j + i * 1000);
1515             
1516             obj->Add(particle);
1517           }
1518       }    
1519   }
1520   else if (idet == 3)
1521     {
1522       if(!fAOD)
1523         AliFatal("Tracklets only available on AOD");
1524       AliAODTracklets* trklets=(AliAODTracklets*)fAOD->GetTracklets();
1525       if(!trklets)
1526         AliFatal("AliAODTracklets not found");
1527       for(Int_t itrklets=0;itrklets<trklets->GetNumberOfTracklets();itrklets++)
1528         {
1529           Float_t eta=-TMath::Log(TMath::Tan(trklets->GetTheta(itrklets)/2));
1530           if(TMath::Abs(eta)>fTrackEtaCut)continue;
1531           Float_t pT=1000*TMath::Abs(trklets->GetDeltaPhi(itrklets));//in mrad
1532           if(pT>fTrackletDphiCut)continue;
1533           TH1F* DphiTrklets = (TH1F*)fListOfHistos->FindObject("DphiTrklets");
1534           DphiTrklets->Fill(1000*trklets->GetDeltaPhi(itrklets)); //in mrad
1535           Float_t phi=trklets->GetPhi(itrklets);
1536           phi+=trklets->GetDeltaPhi(itrklets)*39./34.; //correction dphi*39./34. (Dphi in rad)
1537           if (phi<0) phi+=TMath::TwoPi();
1538           if (phi>TMath::TwoPi()) phi-=TMath::TwoPi();
1539           
1540           AliDPhiBasicParticle* particle = new AliDPhiBasicParticle(eta,phi, pT, 0); // pT = TMath::Abs(trklets->GetDeltaPhi(itrklets)) in mrad and charge = 0
1541           particle->SetUniqueID(fAnalyseUE->GetEventCounter()* 100000 + itrklets);
1542           
1543           obj->Add(particle);
1544         }      
1545     }
1546   else if (idet == 4)
1547     {
1548       if(!fAOD)
1549         AliFatal("Muon selection only implemented on AOD");//FIXME to be implemented also for ESDs as in AliAnalyseLeadingTrackUE::GetAcceptedPArticles
1550       for (Int_t iTrack = 0; iTrack < fAOD->GetNTracks(); iTrack++) {
1551         AliAODTrack* track = fAOD->GetTrack(iTrack);
1552         if (!track->IsMuonTrack()) continue;
1553         //Float_t dca    = track->DCA();
1554         //Float_t chi2   = track->Chi2perNDF();
1555         Float_t rabs   = track->GetRAtAbsorberEnd();
1556         Float_t eta    = track->Eta();
1557         Int_t   matching   = track->GetMatchTrigger();
1558         if (rabs < 17.6 || rabs > 89.5) continue;
1559         if (eta < -4 || eta > -2.5) continue;
1560         if (matching < 2) continue;
1561         
1562         AliDPhiBasicParticle* particle = new AliDPhiBasicParticle(eta,track->Phi(), track->Pt(), track->Charge()); 
1563         particle->SetUniqueID(fAnalyseUE->GetEventCounter() * 100000 + iTrack);
1564         
1565         obj->Add(particle);
1566       }
1567     }      
1568   else
1569     AliFatal(Form("GetParticlesFromDetector: Invalid idet value: %d", idet));
1570   
1571   return obj;  
1572 }
1573
1574 //____________________________________________________________________
1575 Bool_t AliAnalysisTaskPhiCorrelations::IsMuEvent(){
1576   
1577   if(!fAOD)
1578     AliFatal("Muon selection only implemented on AOD");//FIXME to be implemented also for ESDs as in AliAnalyseLeadingTrackUE::GetAcceptedPArticles
1579   for (Int_t iTrack = 0; iTrack < fAOD->GetNTracks(); iTrack++) {
1580     AliAODTrack* track = fAOD->GetTrack(iTrack);
1581     if (!track->IsMuonTrack()) continue;
1582     //Float_t dca    = track->DCA();
1583     //Float_t chi2   = track->Chi2perNDF();
1584     Float_t rabs   = track->GetRAtAbsorberEnd();
1585     Float_t eta    = track->Eta();
1586     Int_t   matching   = track->GetMatchTrigger();
1587     if (rabs < 17.6 || rabs > 89.5) continue;
1588     if (eta < -4 || eta > -2.5) continue;
1589     if (matching < 2) continue;
1590     return kTRUE;
1591   }
1592   return kFALSE;
1593   
1594 }