]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGCF/Correlations/DPhi/AliAnalysisTaskPhiCorrelations.cxx
optimization and bug fix in two track cut
[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 <TRandom.h>
26
27 #include "AliAnalysisTaskPhiCorrelations.h"
28 #include "AliAnalyseLeadingTrackUE.h"
29 #include "AliUEHistograms.h"
30 #include "AliUEHist.h"
31
32 #include "AliAnalysisManager.h"
33 #include "AliAODHandler.h"
34 #include "AliAODInputHandler.h"
35 #include "AliAODMCParticle.h"
36 #include "AliInputEventHandler.h"
37 #include "AliLog.h"
38 #include "AliMCEventHandler.h"
39 #include "AliVParticle.h"
40 #include "AliCFContainer.h"
41
42 #include "AliESDEvent.h"
43 #include "AliESDInputHandler.h"
44 #include "AliMultiplicity.h"
45 #include "AliCentrality.h"
46 #include "AliStack.h"
47 #include "AliAODMCHeader.h"
48
49 #include "AliEventPoolManager.h"
50
51
52 ////////////////////////////////////////////////////////////////////////
53 //
54 // Analysis class for azimuthal correlation studies
55 // Based on the UE task from Sara Vallero and Jan Fiete
56 //
57 // This class needs input AODs.
58 // The output is a list of analysis-specific containers.
59 //
60 // The AOD can be either connected to the InputEventHandler  
61 // for a chain of AOD files 
62 // or 
63 // to the OutputEventHandler
64 // for a chain of ESD files,
65 // in this case the class should be in the train after the jet-finder
66 //
67 //    Authors:
68 //    Jan Fiete Grosse-Oetringhaus
69 // 
70 ////////////////////////////////////////////////////////////////////////
71
72
73 ClassImp( AliAnalysisTaskPhiCorrelations )
74 ClassImp( AliDPhiBasicParticle )
75
76 //____________________________________________________________________
77 AliAnalysisTaskPhiCorrelations:: AliAnalysisTaskPhiCorrelations(const char* name):
78 AliAnalysisTask(name,""),
79 // general configuration
80 fDebug(0),
81 fMode(0),
82 fReduceMemoryFootprint(kFALSE),
83 fFillMixed(kTRUE),
84 fMixingTracks(50000),
85 fCompareCentralities(kFALSE),
86 fTwoTrackEfficiencyStudy(kFALSE),
87 fTwoTrackEfficiencyCut(kFALSE),
88 fUseVtxAxis(kFALSE),
89 fSkipTrigger(kFALSE),
90 // pointers to UE classes
91 fAnalyseUE(0x0),
92 fHistos(0x0),
93 fHistosMixed(0),
94 fkTrackingEfficiency(0x0),
95 // handlers and events
96 fAOD(0x0),
97 fESD(0x0),
98 fArrayMC(0x0),
99 fInputHandler(0x0),
100 fMcEvent(0x0),
101 fMcHandler(0x0),
102 fPoolMgr(0x0),
103 // histogram settings
104 fListOfHistos(0x0), 
105 // event QA
106 fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default) 
107 fZVertex(7.),
108 fCentralityMethod("V0M"),
109 // track cuts
110 fTrackEtaCut(0.8),
111 fPtMin(0.5),
112 fFilterBit(0xFF),
113 fSelectBit(0),
114 fUseChargeHadrons(kFALSE),
115 fSelectCharge(0),
116 fTriggerRestrictEta(-1),
117 fEtaOrdering(kFALSE),
118 fCutConversions(kFALSE),
119 fCutResonances(kFALSE),
120 fFillpT(kFALSE)
121 {
122   // Default constructor
123   // Define input and output slots here
124   // Input slot #0 works with a TChain
125   DefineInput(0, TChain::Class());
126   // Output slot #0 writes into a TList container
127   DefineOutput(0, TList::Class());
128 }
129
130 AliAnalysisTaskPhiCorrelations::~AliAnalysisTaskPhiCorrelations() 
131
132   // destructor
133   
134   if (fListOfHistos  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) 
135     delete fListOfHistos;
136 }
137
138 //____________________________________________________________________
139 void AliAnalysisTaskPhiCorrelations::ConnectInputData(Option_t* /*option*/)
140 {
141   
142   // Connect the input data  
143   if (fDebug > 1) AliInfo("ConnectInputData() ");
144   
145   // Since AODs can either be connected to the InputEventHandler
146   // or to the OutputEventHandler ( the AOD is created by a previus task in the train )
147   // we need to get the pointer to the AODEvent correctly.
148   
149   // Delta AODs are also accepted.
150   
151   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
152   
153   if( handler && handler->InheritsFrom("AliAODInputHandler") ) 
154   { // input AOD
155     fAOD = ((AliAODInputHandler*)handler)->GetEvent();
156     if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
157   } 
158   else 
159   {  //output AOD
160     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
161     if (handler && handler->InheritsFrom("AliAODHandler") ) 
162     {
163       fAOD = ((AliAODHandler*)handler)->GetAOD();
164       if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
165     } 
166     else 
167     {  // no AOD
168       AliWarning("I can't get any AOD Event Handler");
169     }
170   }
171   
172   if (handler && handler->InheritsFrom("AliESDInputHandler") ) 
173   { // input ESD
174     // pointer received per event in ::Exec
175     if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliESDInputHandler");
176   } 
177   
178   // Initialize common pointers
179   Initialize();
180 }
181
182 //____________________________________________________________________
183 void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
184 {
185   // Create the output container
186   
187   if (fDebug > 1) AliInfo("CreateOutputObjects()");
188    
189   // Initialize class with main algorithms, event and track selection. 
190   fAnalyseUE = new AliAnalyseLeadingTrackUE();
191   fAnalyseUE->SetParticleSelectionCriteria(fFilterBit, fUseChargeHadrons, fTrackEtaCut, fPtMin);
192   fAnalyseUE->SetDebug(fDebug); 
193   fAnalyseUE->DefineESDCuts(fFilterBit);
194
195   // Initialize output list of containers
196   if (fListOfHistos != NULL){
197         delete fListOfHistos;
198         fListOfHistos = NULL;
199         }
200   if (!fListOfHistos){
201         fListOfHistos = new TList();
202         fListOfHistos->SetOwner(kTRUE); 
203         }
204
205   // Initialize class to handle histograms 
206   const char* histType = "4";
207   if (fUseVtxAxis)
208     histType = "5";
209   fHistos = new AliUEHistograms("AliUEHistogramsSame", histType);
210   fHistosMixed = new AliUEHistograms("AliUEHistogramsMixed", histType);
211   
212   fHistos->SetSelectCharge(fSelectCharge);
213   fHistosMixed->SetSelectCharge(fSelectCharge);
214   
215   fHistos->SetTriggerRestrictEta(fTriggerRestrictEta);
216   fHistosMixed->SetTriggerRestrictEta(fTriggerRestrictEta);
217   
218   fHistos->SetEtaOrdering(fEtaOrdering);
219   fHistosMixed->SetEtaOrdering(fEtaOrdering);
220
221   fHistos->SetPairCuts(fCutConversions, fCutResonances);
222   fHistosMixed->SetPairCuts(fCutConversions, fCutResonances);
223   
224   // add histograms to list
225   fListOfHistos->Add(fHistos);
226   fListOfHistos->Add(fHistosMixed);
227   
228   fListOfHistos->Add(new TH2F("trackletsVsV0Cent", ";L1 clusters;v0 centrality", 100, -0.5, 9999.5, 101, 0, 101));
229   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));
230   fListOfHistos->Add(new TH1F("eventStat", ";;events", 4, -0.5, 3.5));
231   fListOfHistos->Add(new TH2F("mixedDist", ";centrality;tracks;events", 101, 0, 101, 200, 0, fMixingTracks * 1.5));
232   fListOfHistos->Add(new TH1F("pids", ";pdg;tracks", 2001, -1000.5, 1000.5));
233   
234   PostData(0,fListOfHistos);
235   
236   // Add task configuration to output list 
237   AddSettingsTree();
238
239   // event mixing
240   Int_t trackDepth = fMixingTracks; 
241   Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPool
242    
243   Int_t nCentralityBins  = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(1);
244   Double_t* centralityBins = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray();
245   
246   Int_t nZvtxBins  = 7+1+7;
247   // bins for second buffer are shifted by 100 cm
248   Double_t vertexBins[] = { -7, -5, -3, -1, 1, 3, 5, 7, 93, 95, 97, 99, 101, 103, 105, 107 };
249   Double_t* zvtxbin = vertexBins;
250
251   if (fHistos->GetUEHist(2)->GetEventHist()->GetNVar() > 2)
252   {
253     nZvtxBins = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(2);
254     zvtxbin = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetXbins()->GetArray();
255   }
256
257   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
258 }
259
260 //____________________________________________________________________
261 void  AliAnalysisTaskPhiCorrelations::Exec(Option_t */*option*/)
262 {
263   // receive ESD pointer if we are not running AOD analysis
264   if (!fAOD)
265   {
266     AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
267     if (handler && handler->InheritsFrom("AliESDInputHandler"))
268       fESD = ((AliESDInputHandler*)handler)->GetEvent();
269   }
270
271   if (fMode)
272   {
273     // correction mode
274     
275     if (fMcHandler)
276       fMcEvent = fMcHandler->MCEvent();
277
278     if (fAOD)
279     {
280       // array of MC particles
281       fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
282       if (!fArrayMC)
283         AliFatal("No array of MC particles found !!!");
284     }
285     
286     AnalyseCorrectionMode();
287   }
288   else AnalyseDataMode();
289 }
290
291 /******************** ANALYSIS METHODS *****************************/
292
293 //____________________________________________________________________
294 void  AliAnalysisTaskPhiCorrelations::AddSettingsTree()
295 {
296   //Write settings to output list
297   TTree *settingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
298   settingsTree->Branch("fnTracksVertex", &fnTracksVertex,"nTracksVertex/I");
299   settingsTree->Branch("fZVertex", &fZVertex,"ZVertex/D");
300   //settingsTree->Branch("fCentralityMethod", fCentralityMethod.Data(),"CentralityMethod/C");
301   settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
302   settingsTree->Branch("fPtMin", &fPtMin, "PtMin/D");
303   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
304   settingsTree->Branch("fSelectBit", &fSelectBit,"EventSelectionBit/I");
305   settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
306   settingsTree->Branch("fSelectCharge", &fSelectCharge,"SelectCharge/I");
307   settingsTree->Branch("fTriggerRestrictEta", &fTriggerRestrictEta,"TriggerRestrictEta/D");
308   settingsTree->Branch("fEtaOrdering", &fEtaOrdering,"EtaOrdering/O");
309   settingsTree->Branch("fCutConversions", &fCutConversions,"CutConversions/O");
310   settingsTree->Branch("fCutResonances", &fCutResonances,"CutResonances/O");
311   settingsTree->Branch("fFillpT", &fFillpT,"FillpT/O");
312   settingsTree->Branch("fkTrackingEfficiency", "TH1D", &fkTrackingEfficiency);
313   settingsTree->Branch("fMixingTracks", &fMixingTracks,"MixingTracks/I");
314   settingsTree->Branch("fSkipTrigger", &fSkipTrigger,"SkipTrigger/O");
315   
316   settingsTree->Fill();
317   fListOfHistos->Add(settingsTree);
318 }  
319
320 //____________________________________________________________________
321 void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
322 {
323   // Run the analysis on MC to get the correction maps
324   //
325  
326   if ( fDebug > 3 ) AliInfo( " Processing event in Corrections mode ..." );
327   
328   Double_t centrality = 0;
329   
330   if (fCentralityMethod.Length() > 0)
331   {
332     AliCentrality *centralityObj = 0;
333     if (fAOD)
334       centralityObj = fAOD->GetHeader()->GetCentralityP();
335     else if (fESD)
336       centralityObj = fESD->GetCentrality();
337     
338     if (centralityObj)
339     {
340       centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
341       AliInfo(Form("Centrality is %f", centrality));
342     }
343     else
344     {
345       Printf("WARNING: Centrality object is 0");
346       centrality = -1;
347      }
348   }
349   
350   // Support for ESD and AOD based analysis
351   AliVEvent* inputEvent = fAOD;
352   if (!inputEvent)
353     inputEvent = fESD;
354   
355   TObject* mc = fArrayMC;
356   if (!mc)
357     mc = fMcEvent;
358   
359   // count all events
360   fHistos->FillEvent(centrality, -1);
361   
362   // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
363   TObject* vertexSupplier = fMcEvent;
364   if (fAOD) // AOD
365     vertexSupplier = fAOD->FindListObject(AliAODMCHeader::StdBranchName());
366     
367   if (!fAnalyseUE->VertexSelection(vertexSupplier, 0, fZVertex)) 
368     return;
369   
370   Float_t zVtx = 0;
371   if (fAOD)
372     zVtx = ((AliAODMCHeader*) vertexSupplier)->GetVtxZ();
373   else
374     zVtx = fMcEvent->GetPrimaryVertex()->GetZ();
375   Float_t weight = 1;
376   if (fFillpT)
377     weight = -1;
378     
379   // Get MC primaries
380   TObjArray* tracksMC = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, -1, kTRUE);
381   
382   if (fAOD)
383   {
384     for (Int_t i=0; i<fArrayMC->GetEntriesFast(); i++)
385       ((TH1F*) fListOfHistos->FindObject("pids"))->Fill(((AliAODMCParticle*) fArrayMC->At(i))->PdgCode());
386   }
387   else
388   {
389     for (Int_t i=0; i<fMcEvent->GetNumberOfTracks(); i++)
390       ((TH1F*) fListOfHistos->FindObject("pids"))->Fill(fMcEvent->GetTrack(i)->PdgCode());
391   }
392   
393   // (MC-true all particles)
394   // STEP 0
395   fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, 0, weight);
396   
397   // mixed event
398   if (fFillMixed)
399   {
400     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
401     //pool->PrintInfo();
402     if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5) 
403       for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++) 
404         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
405     pool->UpdatePool(CloneAndReduceTrackList(tracksMC));
406   }
407   
408 //   Printf("trigger: %d", ((AliInputEventHandler*)fInputHandler)->IsEventSelected());
409   
410   // Trigger selection ************************************************
411   if (fSkipTrigger || fAnalyseUE->TriggerSelection(fInputHandler))
412   {  
413     // (MC-true all particles)
414     // STEP 1
415     if (!fReduceMemoryFootprint)
416       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTriggered, tracksMC, 0, weight);
417     else
418       fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
419       
420     if (!inputEvent) {
421       AliFatal("UNEXPECTED: inputEvent is 0. Trigger selection should have failed");
422       return;
423     }
424     
425     // Vertex selection *************************************************
426     if (fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex))
427     {
428       // fill here for tracking efficiency
429       // loop over particle species
430       
431       for (Int_t particleSpecies = 0; particleSpecies < 4; particleSpecies++)
432       {
433         TObjArray* primMCParticles = fAnalyseUE->GetAcceptedParticles(mc, 0x0, kTRUE, particleSpecies, kTRUE);
434         TObjArray* primRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE);
435         TObjArray* allRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE);
436       
437         fHistos->FillTrackingEfficiency(primMCParticles, primRecoTracksMatched, allRecoTracksMatched, particleSpecies, centrality);
438         
439 //      Printf("%d --> %d %d %d", particleSpecies, primMCParticles->GetEntries(), primRecoTracksMatched->GetEntries(), allRecoTracksMatched->GetEntries());
440
441         delete primMCParticles;
442         delete primRecoTracksMatched;
443         delete allRecoTracksMatched;
444       }
445     
446       // (MC-true all particles)
447       // STEP 2
448       if (!fReduceMemoryFootprint)
449         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepVertex, tracksMC, 0, weight);
450       else
451         fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
452       
453       // Get MC primaries that match reconstructed track
454       TObjArray* tracksRecoMatchedPrim = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, -1, kTRUE);
455       
456       // (RECO-matched (quantities from MC particle) primary particles)
457       // STEP 4
458       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, 0, weight);
459       
460       // Get MC primaries + secondaries that match reconstructed track
461       TObjArray* tracksRecoMatchedAll = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, -1, kTRUE);
462       
463       // (RECO-matched (quantities from MC particle) all particles)
464       // STEP 5
465       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, 0, weight);
466       
467       // Get RECO tracks
468       TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
469       
470       // (RECO all tracks)
471       // STEP 6
472       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, 0, weight);
473       
474       // mixed event
475       if (fFillMixed)
476       {
477         AliEventPool* pool2 = fPoolMgr->GetEventPool(centrality, zVtx + 100);
478         //pool2->PrintInfo();
479         if (pool2->IsReady() || pool2->NTracksInPool() > fMixingTracks / 10 || pool2->GetCurrentNEvents() >= 5) 
480           for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++) 
481             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0));
482         pool2->UpdatePool(CloneAndReduceTrackList(tracks));
483       }
484       
485       if (0 && !fReduceMemoryFootprint)
486       {
487         // make list of secondaries (matched with MC)
488         TObjArray* tracksRecoMatchedSecondaries = new TObjArray;
489         for (Int_t i=0; i<tracksRecoMatchedAll->GetEntriesFast(); i++)
490           if (((AliAODMCParticle*)tracksRecoMatchedAll->At(i))->IsPhysicalPrimary() == kFALSE)
491             tracksRecoMatchedSecondaries->Add(tracksRecoMatchedAll->At(i));
492       
493         // Study: Use only secondaries as trigger particles and plot the correlation vs. all particles; store in step 9
494         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy2, tracksRecoMatchedSecondaries, tracksRecoMatchedAll, weight);
495         
496         // Study: Use only primaries as trigger particles and plot the correlation vs. secondaries; store in step 8
497         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksRecoMatchedPrim, tracksRecoMatchedSecondaries, weight);
498       
499         // plot delta phi vs process id of secondaries
500         // trigger particles: primaries in 4 < pT < 10
501         // associated particles: secondaries in 1 < pT < 10
502         
503         for (Int_t i=0; i<tracksRecoMatchedPrim->GetEntriesFast(); i++)
504         {
505           AliVParticle* triggerParticle = (AliVParticle*) tracksRecoMatchedPrim->At(i);
506           
507           if (triggerParticle->Pt() < 4 || triggerParticle->Pt() > 10)
508             continue;
509           
510           for (Int_t j=0; j<tracksRecoMatchedSecondaries->GetEntriesFast(); j++)
511           {
512             AliAODMCParticle* particle = (AliAODMCParticle*) tracksRecoMatchedSecondaries->At(j);
513             
514             if (particle->Pt() < 1 || particle->Pt() > 10)
515               continue;
516             
517             if (particle->Pt() > triggerParticle->Pt())
518               continue;
519               
520             Double_t deltaPhi = triggerParticle->Phi() - particle->Phi();
521             if (deltaPhi > 1.5 * TMath::Pi()) 
522               deltaPhi -= TMath::TwoPi();
523             if (deltaPhi < -0.5 * TMath::Pi())
524               deltaPhi += TMath::TwoPi();
525               
526             Int_t processID = fMcEvent->Stack()->Particle(particle->GetLabel())->GetUniqueID();
527               
528             ((TH2F*) fListOfHistos->FindObject("processIDs"))->Fill(deltaPhi, processID);
529           }
530         }
531         
532         delete tracksRecoMatchedSecondaries;
533       }
534   
535       delete tracksRecoMatchedPrim;
536       delete tracksRecoMatchedAll;
537       delete tracks;
538     }
539   }
540   
541   delete tracksMC;
542 }
543
544 //____________________________________________________________________
545 void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
546 {
547
548   // Run the analysis on DATA or MC to get raw distributions
549  
550   if ( fDebug > 3 ) AliInfo( " Processing event in Data mode ..." );
551
552   ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(0);
553
554   // skip not selected events here (the AOD is not updated for those)
555   if (!fInputHandler)
556     return;
557     
558   if (!(fInputHandler->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kUserDefined)))
559     return;
560
561   Double_t centrality = 0;
562   
563   AliCentrality *centralityObj = 0;
564   if (fCentralityMethod.Length() > 0)
565   {
566     if (fAOD)
567       centralityObj = fAOD->GetHeader()->GetCentralityP();
568     else if (fESD)
569       centralityObj = fESD->GetCentrality();
570     
571     if (centralityObj)
572       centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
573       //centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
574     else
575       centrality = -1;
576
577     if (fAOD)
578     {
579       // remove outliers
580       if (centrality == 0)
581       {
582         if (fAOD->GetVZEROData())
583         {
584           Float_t multV0 = 0;
585           for (Int_t i=0; i<64; i++)
586             multV0 += fAOD->GetVZEROData()->GetMultiplicity(i);
587           if (multV0 < 19500)
588           {
589             centrality = -1;
590             AliInfo("Rejecting event due to too small V0 multiplicity");
591           }
592         }
593       }
594     }
595     
596     AliInfo(Form("Centrality is %f", centrality));
597   }
598   
599   // Support for ESD and AOD based analysis
600   AliVEvent* inputEvent = fAOD;
601   if (!inputEvent)
602     inputEvent = fESD;
603
604   Float_t bSign = (inputEvent->GetMagneticField() > 0) ? 1 : -1;
605
606   fHistos->SetRunNumber(inputEvent->GetRunNumber());
607   
608   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
609   fHistos->FillEvent(centrality, AliUEHist::kCFStepAll);
610   
611   // Trigger selection ************************************************
612   if (!fAnalyseUE->TriggerSelection(fInputHandler)) return;
613   
614   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
615   fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
616   
617   // Vertex selection *************************************************
618   if(!fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex)) return;
619   
620   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
621   fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
622  
623   // optimization
624   if (centrality < 0 && !fCompareCentralities)
625     return;
626
627   TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
628   // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
629   TObjArray* tracksClone = CloneAndReduceTrackList(tracks);
630   delete tracks;
631   
632   //Printf("Accepted %d tracks", tracks->GetEntries());
633   
634   const AliVVertex* vertex = inputEvent->GetPrimaryVertex();
635   Double_t zVtx = vertex->GetZ();
636   
637   Float_t weight = 1;
638   if (fFillpT)
639     weight = -1;
640
641   // fill two different centralities (syst study)
642   // the zvtx axis is used to distinguish the results of both centralities: configured centrality in zvtx = 0, SPD in zvtx = 2
643   if (fCompareCentralities)
644     zVtx = 0;
645   
646   // Fill containers at STEP 6 (reconstructed)
647   if (centrality >= 0)
648   {
649     fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, 0, weight);
650     ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(1);
651     
652     if (fTwoTrackEfficiencyCut)
653       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, 0, weight, kTRUE, kTRUE, bSign);
654   }
655
656   // fill second time with SPD centrality
657   if (fCompareCentralities && centralityObj)
658   {
659     centrality = centralityObj->GetCentralityPercentile("CL1");
660     if (centrality >= 0)
661       fHistos->FillCorrelations(centrality, 2, AliUEHist::kCFStepReconstructed, tracksClone, 0, weight);
662   }
663   
664   if (fFillMixed)
665   {
666     // event mixing
667     
668     // 1. First get an event pool corresponding in mult (cent) and
669     //    zvertex to the current event. Once initialized, the pool
670     //    should contain nMix (reduced) events. This routine does not
671     //    pre-scan the chain. The first several events of every chain
672     //    will be skipped until the needed pools are filled to the
673     //    specified depth. If the pool categories are not too rare, this
674     //    should not be a problem. If they are rare, you could lose
675     //    statistics.
676
677     // 2. Collect the whole pool's content of tracks into one TObjArray
678     //    (bgTracks), which is effectively a single background super-event.
679
680     // 3. The reduced and bgTracks arrays must both be passed into
681     //    FillCorrelations(). Also nMix should be passed in, so a weight
682     //    of 1./nMix can be applied.
683
684     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
685     
686     if (!pool)
687       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, zVtx));
688     
689     //pool->SetDebug(1);
690      
691     if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5) 
692     {
693       
694       Int_t nMix = pool->GetCurrentNEvents();
695 //       cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
696       
697       ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
698       ((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
699       if (pool->IsReady())
700         ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
701     
702       // Fill mixed-event histos here  
703       for (Int_t jMix=0; jMix<nMix; jMix++) 
704       {
705         TObjArray* bgTracks = pool->GetEvent(jMix);
706         
707         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0));
708
709         if (fTwoTrackEfficiencyCut)
710           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kTRUE, bSign);
711       }
712     }
713     
714     // ownership is with the pool now
715     pool->UpdatePool(tracksClone);
716     //pool->PrintInfo();
717   }
718   else
719     delete tracksClone;
720 }
721
722 TObjArray* AliAnalysisTaskPhiCorrelations::CloneAndReduceTrackList(TObjArray* tracks)
723 {
724   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
725   
726   TObjArray* tracksClone = new TObjArray;
727   tracksClone->SetOwner(kTRUE);
728   
729   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
730   {
731     AliVParticle* particle = (AliVParticle*) tracks->At(i);
732     tracksClone->Add(new AliDPhiBasicParticle(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge()));
733   }
734   
735   return tracksClone;
736 }
737
738 //____________________________________________________________________
739 void  AliAnalysisTaskPhiCorrelations::Initialize()
740 {
741   // input handler
742   fInputHandler = (AliInputEventHandler*)
743          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
744   // MC handler
745   fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
746 }