]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG4/JetTasks/AliAnalysisTaskPhiCorrelations.cxx
coverity fixes
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / 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 #include <TMCProcess.h>
27
28 #include "AliAnalysisTaskPhiCorrelations.h"
29 #include "AliAnalyseLeadingTrackUE.h"
30 #include "AliUEHistograms.h"
31 #include "AliUEHist.h"
32
33 #include "AliAnalysisHelperJetTasks.h"
34 #include "AliAnalysisManager.h"
35 #include "AliAODHandler.h"
36 #include "AliAODInputHandler.h"
37 #include "AliAODMCParticle.h"
38 #include "AliGenPythiaEventHeader.h"
39 #include "AliInputEventHandler.h"
40 #include "AliLog.h"
41 #include "AliMCEventHandler.h"
42 #include "AliVParticle.h"
43 #include "AliCFContainer.h"
44
45 #include "AliESDEvent.h"
46 #include "AliESDInputHandler.h"
47 #include "AliMultiplicity.h"
48 #include "AliCentrality.h"
49 #include "AliStack.h"
50
51 #include "AliEventPoolManager.h"
52
53
54 ////////////////////////////////////////////////////////////////////////
55 //
56 // Analysis class for azimuthal correlation studies
57 // Based on the UE task from Sara Vallero and Jan Fiete
58 //
59 // This class needs input AODs.
60 // The output is a list of analysis-specific containers.
61 //
62 // The AOD can be either connected to the InputEventHandler  
63 // for a chain of AOD files 
64 // or 
65 // to the OutputEventHandler
66 // for a chain of ESD files,
67 // in this case the class should be in the train after the jet-finder
68 //
69 //    Authors:
70 //    Jan Fiete Grosse-Oetringhaus
71 // 
72 ////////////////////////////////////////////////////////////////////////
73
74
75 ClassImp( AliAnalysisTaskPhiCorrelations )
76
77 //____________________________________________________________________
78 AliAnalysisTaskPhiCorrelations:: AliAnalysisTaskPhiCorrelations(const char* name):
79 AliAnalysisTask(name,""),
80 // general configuration
81 fDebug(0),
82 fMode(0),
83 fReduceMemoryFootprint(kFALSE),
84 // pointers to UE classes
85 fAnalyseUE(0x0),
86 fHistos(0x0),
87 fHistosMixed(0),
88 fkTrackingEfficiency(0x0),
89 // handlers and events
90 fAOD(0x0),
91 fESD(0x0),
92 fArrayMC(0x0),
93 fInputHandler(0x0),
94 fMcEvent(0x0),
95 fMcHandler(0x0),
96 fPoolMgr(0x0),
97 // histogram settings
98 fListOfHistos(0x0), 
99 // event QA
100 fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default) 
101 fZVertex(7.),
102 fCentralityMethod("V0M"),
103 // track cuts
104 fTrackEtaCut(0.8),
105 fPtMin(0.5),
106 fFilterBit(0xFF),
107 fSelectBit(0),
108 fUseChargeHadrons(kFALSE),
109 fSelectCharge(0)
110 {
111   // Default constructor
112   // Define input and output slots here
113   // Input slot #0 works with a TChain
114   DefineInput(0, TChain::Class());
115   // Output slot #0 writes into a TList container
116   DefineOutput(0, TList::Class());
117 }
118
119 AliAnalysisTaskPhiCorrelations::~AliAnalysisTaskPhiCorrelations() 
120
121   // destructor
122   
123   if (fListOfHistos  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) 
124     delete fListOfHistos;
125 }
126
127 //____________________________________________________________________
128 void AliAnalysisTaskPhiCorrelations::ConnectInputData(Option_t* /*option*/)
129 {
130   
131   // Connect the input data  
132   if (fDebug > 1) AliInfo("ConnectInputData() ");
133   
134   // Since AODs can either be connected to the InputEventHandler
135   // or to the OutputEventHandler ( the AOD is created by a previus task in the train )
136   // we need to get the pointer to the AODEvent correctly.
137   
138   // Delta AODs are also accepted.
139   
140   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
141   
142   if( handler && handler->InheritsFrom("AliAODInputHandler") ) 
143   { // input AOD
144     fAOD = ((AliAODInputHandler*)handler)->GetEvent();
145     if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
146   } 
147   else 
148   {  //output AOD
149     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
150     if (handler && handler->InheritsFrom("AliAODHandler") ) 
151     {
152       fAOD = ((AliAODHandler*)handler)->GetAOD();
153       if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
154     } 
155     else 
156     {  // no AOD
157       AliWarning("I can't get any AOD Event Handler");
158     }
159   }
160   
161   if (handler && handler->InheritsFrom("AliESDInputHandler") ) 
162   { // input ESD
163     // pointer received per event in ::Exec
164     if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliESDInputHandler");
165   } 
166   
167   // Initialize common pointers
168   Initialize();
169 }
170
171 //____________________________________________________________________
172 void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
173 {
174   // Create the output container
175   
176   if (fDebug > 1) AliInfo("CreateOutputObjects()");
177    
178   // Initialize class with main algorithms, event and track selection. 
179   fAnalyseUE = new AliAnalyseLeadingTrackUE();
180   fAnalyseUE->SetParticleSelectionCriteria(fFilterBit, fUseChargeHadrons, fTrackEtaCut, fPtMin);
181   fAnalyseUE->SetDebug(fDebug); 
182   fAnalyseUE->DefineESDCuts(fFilterBit);
183
184   // Initialize output list of containers
185   if (fListOfHistos != NULL){
186         delete fListOfHistos;
187         fListOfHistos = NULL;
188         }
189   if (!fListOfHistos){
190         fListOfHistos = new TList();
191         fListOfHistos->SetOwner(kTRUE); 
192         }
193
194   // Initialize class to handle histograms 
195   fHistos = new AliUEHistograms("AliUEHistogramsSame", "4");
196   fHistosMixed = new AliUEHistograms("AliUEHistogramsMixed", "4");
197   
198   fHistos->SetSelectCharge(fSelectCharge);
199   fHistosMixed->SetSelectCharge(fSelectCharge);
200   
201   // add histograms to list
202   fListOfHistos->Add(fHistos);
203   fListOfHistos->Add(fHistosMixed);
204   
205   fListOfHistos->Add(new TH2F("trackletsVsV0Cent", ";L1 clusters;v0 centrality", 100, -0.5, 9999.5, 101, 0, 101));
206   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));
207   
208   // Add task configuration to output list 
209   AddSettingsTree();
210
211   // event mixing
212   //  Int_t trackDepth = 100; // Require e.g. 20 5-track events, or 2 50-track events
213   Int_t trackDepth = 5000; 
214   Int_t poolsize   = 100;  // Maximum number of events
215   
216   Int_t nCentralityBins  = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(1);
217   Double_t* centralityBins = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray();
218   
219   Int_t nZvtxBins  = 7;
220   Double_t vertexBins[] = { -7, -5, -3, -1, 1, 3, 5, 7 };
221   Double_t* zvtxbin = vertexBins;
222   
223   if (fHistos->GetUEHist(2)->GetEventHist()->GetNVar() > 2)
224   {
225     nZvtxBins = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(2);
226     zvtxbin = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetXbins()->GetArray();
227   }
228
229   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
230 }
231
232 //____________________________________________________________________
233 void  AliAnalysisTaskPhiCorrelations::Exec(Option_t */*option*/)
234 {
235   // array of MC particles
236   if (fMcHandler) {
237     if (fAOD)
238     {
239       fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
240       if (!fArrayMC)
241         AliFatal("No array of MC particles found !!!");
242     }
243     fMcEvent = fMcHandler->MCEvent();
244   }
245
246   // receive ESD pointer if we are not running AOD analysis
247   if (!fAOD)
248   {
249     AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
250     if (handler && handler->InheritsFrom("AliESDInputHandler"))
251       fESD = ((AliESDInputHandler*)handler)->GetEvent();
252   }
253
254   // Analyse the event
255   if (fMode) AnalyseCorrectionMode();
256   else AnalyseDataMode();
257
258   PostData(0,fListOfHistos);
259 }
260
261 /******************** ANALYSIS METHODS *****************************/
262
263 //____________________________________________________________________
264 void  AliAnalysisTaskPhiCorrelations::AddSettingsTree()
265 {
266   //Write settings to output list
267   TTree *settingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
268   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
269   settingsTree->Branch("fSelectBit", &fSelectBit,"EventSelectionBit/I");
270   settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
271   settingsTree->Branch("fPtMin", &fPtMin, "PtMin/D");
272   settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
273   settingsTree->Branch("fkTrackingEfficiency", "TH1D", &fkTrackingEfficiency);
274   settingsTree->Fill();
275   fListOfHistos->Add(settingsTree);
276 }  
277
278 //____________________________________________________________________
279 void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
280 {
281   // Run the analysis on MC to get the correction maps
282   //
283  
284   if ( fDebug > 3 ) AliInfo( " Processing event in Corrections mode ..." );
285   
286   Double_t centrality = 0;
287   
288   if (fCentralityMethod.Length() > 0)
289   {
290     AliCentrality *centralityObj = 0;
291     if (fAOD)
292       centralityObj = fAOD->GetHeader()->GetCentralityP();
293     else if (fESD)
294       centralityObj = fESD->GetCentrality();
295     
296     if (centralityObj)
297     {
298       centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
299       AliInfo(Form("Centrality is %f", centrality));
300     }
301     else
302     {
303       Printf("WARNING: Centrality object is 0");
304       centrality = -1;
305      }
306   }
307   
308   // Support for ESD and AOD based analysis
309   AliVEvent* inputEvent = fAOD;
310   if (!inputEvent)
311     inputEvent = fESD;
312   
313   TObject* mc = fArrayMC;
314   if (!mc)
315     mc = fMcEvent;
316   
317   // count all events
318   fHistos->FillEvent(centrality, -1);
319   
320   // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
321   if (!fAnalyseUE->VertexSelection(fMcEvent, 0, fZVertex)) 
322     return;
323   
324   Float_t zVtx = fMcEvent->GetPrimaryVertex()->GetZ();
325     
326   // Get MC primaries
327   TObjArray* tracksMC = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, -1, kTRUE);
328   
329   // (MC-true all particles)
330   // STEP 0
331   fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC);
332   
333   // Trigger selection ************************************************
334   if (fAnalyseUE->TriggerSelection(fInputHandler))
335   {  
336     // (MC-true all particles)
337     // STEP 1
338     if (!fReduceMemoryFootprint)
339       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTriggered, tracksMC);
340     else
341       fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
342       
343     if (!inputEvent)
344       AliFatal("UNEXPECTED: inputEvent is 0. Trigger selection should have failed");
345     
346     // Vertex selection *************************************************
347     if (fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex))
348     {
349       // fill here for tracking efficiency
350       // loop over particle species
351       
352       for (Int_t particleSpecies = 0; particleSpecies < 4; particleSpecies++)
353       {
354         TObjArray* primMCParticles = fAnalyseUE->GetAcceptedParticles(mc, 0x0, kTRUE, particleSpecies, kTRUE);
355         TObjArray* primRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE);
356         TObjArray* allRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE);
357       
358         fHistos->FillTrackingEfficiency(primMCParticles, primRecoTracksMatched, allRecoTracksMatched, particleSpecies, centrality);
359         
360         delete primMCParticles;
361         delete primRecoTracksMatched;
362         delete allRecoTracksMatched;
363       }
364     
365       // (MC-true all particles)
366       // STEP 2
367       if (!fReduceMemoryFootprint)
368         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepVertex, tracksMC);
369       else
370         fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
371       
372       // Get MC primaries that match reconstructed track
373       TObjArray* tracksRecoMatchedPrim = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, -1, kTRUE);
374       
375       // (RECO-matched (quantities from MC particle) primary particles)
376       // STEP 4
377       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim);
378       
379       // Get MC primaries + secondaries that match reconstructed track
380       TObjArray* tracksRecoMatchedAll = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, -1, kTRUE);
381       
382       // (RECO-matched (quantities from MC particle) all particles)
383       // STEP 5
384       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll);
385       
386       // Get RECO tracks
387       TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
388       
389       // (RECO all tracks)
390       // STEP 6
391       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks);
392       
393       if (0 && !fReduceMemoryFootprint)
394       {
395         // make list of secondaries (matched with MC)
396         TObjArray* tracksRecoMatchedSecondaries = new TObjArray;
397         for (Int_t i=0; i<tracksRecoMatchedAll->GetEntries(); i++)
398           if (((AliAODMCParticle*)tracksRecoMatchedAll->At(i))->IsPhysicalPrimary() == kFALSE)
399             tracksRecoMatchedSecondaries->Add(tracksRecoMatchedAll->At(i));
400       
401         // Study: Use only secondaries as trigger particles and plot the correlation vs. all particles; store in step 9
402         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy2, tracksRecoMatchedSecondaries, tracksRecoMatchedAll);
403         
404         // Study: Use only primaries as trigger particles and plot the correlation vs. secondaries; store in step 8
405         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksRecoMatchedPrim, tracksRecoMatchedSecondaries);
406       
407         // plot delta phi vs process id of secondaries
408         // trigger particles: primaries in 4 < pT < 10
409         // associated particles: secondaries in 1 < pT < 10
410         
411         for (Int_t i=0; i<tracksRecoMatchedPrim->GetEntries(); i++)
412         {
413           AliVParticle* triggerParticle = (AliVParticle*) tracksRecoMatchedPrim->At(i);
414           
415           if (triggerParticle->Pt() < 4 || triggerParticle->Pt() > 10)
416             continue;
417           
418           for (Int_t j=0; j<tracksRecoMatchedSecondaries->GetEntries(); j++)
419           {
420             AliAODMCParticle* particle = (AliAODMCParticle*) tracksRecoMatchedSecondaries->At(j);
421             
422             if (particle->Pt() < 1 || particle->Pt() > 10)
423               continue;
424             
425             if (particle->Pt() > triggerParticle->Pt())
426               continue;
427               
428             Double_t deltaPhi = triggerParticle->Phi() - particle->Phi();
429             if (deltaPhi > 1.5 * TMath::Pi()) 
430               deltaPhi -= TMath::TwoPi();
431             if (deltaPhi < -0.5 * TMath::Pi())
432               deltaPhi += TMath::TwoPi();
433               
434             Int_t processID = fMcEvent->Stack()->Particle(particle->GetLabel())->GetUniqueID();
435               
436             ((TH2F*) fListOfHistos->FindObject("processIDs"))->Fill(deltaPhi, processID);
437           }
438         }
439         
440         delete tracksRecoMatchedSecondaries;
441       }
442   
443       delete tracksRecoMatchedPrim;
444       delete tracksRecoMatchedAll;
445       delete tracks;
446     }
447   }
448   
449   delete tracksMC;
450 }
451
452 //____________________________________________________________________
453 void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
454 {
455
456   // Run the analysis on DATA or MC to get raw distributions
457  
458   if ( fDebug > 3 ) AliInfo( " Processing event in Data mode ..." );
459
460   // skip not selected events here (the AOD is not updated for those)
461   if (!fInputHandler)
462     return;
463     
464   if (!(fInputHandler->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kUserDefined)))
465     return;
466
467   Double_t centrality = 0;
468   
469   if (fCentralityMethod.Length() > 0)
470   {
471     AliCentrality *centralityObj = 0;
472     if (fAOD)
473       centralityObj = fAOD->GetHeader()->GetCentralityP();
474     else if (fESD)
475       centralityObj = fESD->GetCentrality();
476     
477     if (centralityObj)
478       centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
479       //centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
480     else
481       centrality = -1;
482     AliInfo(Form("Centrality is %f", centrality));
483   }
484   
485   // Support for ESD and AOD based analysis
486   AliVEvent* inputEvent = fAOD;
487   if (!inputEvent)
488     inputEvent = fESD;
489
490   fHistos->SetRunNumber(inputEvent->GetRunNumber());
491   
492   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
493   fHistos->FillEvent(centrality, AliUEHist::kCFStepAll);
494   
495   // Trigger selection ************************************************
496   if (!fAnalyseUE->TriggerSelection(fInputHandler)) return;
497   
498   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
499   fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
500   
501   // Vertex selection *************************************************
502   if(!fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex)) return;
503   
504   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
505   fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
506  
507   if (centrality < 0)
508     return;
509
510   TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
511   //Printf("Accepted %d tracks", tracks->GetEntries());
512   
513   // event mixing
514   
515   // 1. First get an event pool corresponding in mult (cent) and
516   //    zvertex to the current event. Once initialized, the pool
517   //    should contain nMix (reduced) events. This routine does not
518   //    pre-scan the chain. The first several events of every chain
519   //    will be skipped until the needed pools are filled to the
520   //    specified depth. If the pool categories are not too rare, this
521   //    should not be a problem. If they are rare, you could lose
522   //    statistics.
523
524   // 2. Collect the whole pool's content of tracks into one TObjArray
525   //    (bgTracks), which is effectively a single background super-event.
526
527   // 3. The reduced and bgTracks arrays must both be passed into
528   //    FillCorrelations(). Also nMix should be passed in, so a weight
529   //    of 1./nMix can be applied.
530
531   const AliVVertex* vertex = inputEvent->GetPrimaryVertex();
532   Double_t zVtx = vertex->GetZ();
533   
534   AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
535   
536   if (!pool)
537     AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, zVtx));
538   
539   //pool->SetDebug(1);
540     
541   // Fill containers at STEP 6 (reconstructed)
542   fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks);
543   
544   if (pool->IsReady()) 
545   {
546     
547     Int_t nMix = pool->GetCurrentNEvents();
548     //cout << "nMix = " << nMix << endl;
549   
550     // Fill mixed-event histos here  
551     for (Int_t jMix=0; jMix<nMix; jMix++) {
552       TObjArray* bgTracks = pool->GetEvent(jMix);
553       fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, bgTracks, 1.0 / nMix, (jMix == 0));
554     }
555   }
556   
557   // clone and give ownership to event pool
558   TObjArray* tracksClone = (TObjArray*) tracks->Clone();
559   tracksClone->SetOwner(kTRUE);
560   
561   pool->UpdatePool(tracksClone);
562   //pool->PrintInfo();
563
564   delete tracks;
565 }
566
567 //____________________________________________________________________
568 void  AliAnalysisTaskPhiCorrelations::Initialize()
569 {
570   // input handler
571   fInputHandler = (AliInputEventHandler*)
572          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
573   // MC handler
574   fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
575 }