]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG4/JetTasks/AliAnalysisTaskPhiCorrelations.cxx
updated dphi correlation analysis
[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; // Require e.g. 20 5-track events, or 2 50-track events
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  = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(2);
220   Double_t* zvtxbin = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetXbins()->GetArray();
221
222   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
223   
224   delete[] zvtxbin;
225 }
226
227 //____________________________________________________________________
228 void  AliAnalysisTaskPhiCorrelations::Exec(Option_t */*option*/)
229 {
230   // array of MC particles
231   if (fMcHandler) {
232     if (fAOD)
233     {
234       fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
235       if (!fArrayMC)
236         AliFatal("No array of MC particles found !!!");
237     }
238     fMcEvent = fMcHandler->MCEvent();
239   }
240
241   // receive ESD pointer if we are not running AOD analysis
242   if (!fAOD)
243   {
244     AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
245     if (handler && handler->InheritsFrom("AliESDInputHandler"))
246       fESD = ((AliESDInputHandler*)handler)->GetEvent();
247   }
248
249   // Analyse the event
250   if (fMode) AnalyseCorrectionMode();
251   else AnalyseDataMode();
252
253   PostData(0,fListOfHistos);
254 }
255
256 /******************** ANALYSIS METHODS *****************************/
257
258 //____________________________________________________________________
259 void  AliAnalysisTaskPhiCorrelations::AddSettingsTree()
260 {
261   //Write settings to output list
262   TTree *settingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
263   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
264   settingsTree->Branch("fSelectBit", &fSelectBit,"EventSelectionBit/I");
265   settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
266   settingsTree->Branch("fPtMin", &fPtMin, "PtMin/D");
267   settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
268   settingsTree->Branch("fkTrackingEfficiency", "TH1D", &fkTrackingEfficiency);
269   settingsTree->Fill();
270   fListOfHistos->Add(settingsTree);
271 }  
272
273 //____________________________________________________________________
274 void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
275 {
276   // Run the analysis on MC to get the correction maps
277   //
278  
279   if ( fDebug > 3 ) AliInfo( " Processing event in Corrections mode ..." );
280   
281   Double_t centrality = 0;
282   
283   if (fCentralityMethod.Length() > 0)
284   {
285     AliCentrality *centralityObj = 0;
286     if (fAOD)
287       centralityObj = fAOD->GetHeader()->GetCentralityP();
288     else if (fESD)
289       centralityObj = fESD->GetCentrality();
290     
291     if (centralityObj)
292     {
293       centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
294       Printf("Centrality is %f", centrality);  
295     }
296     else
297     {
298       Printf("WARNING: Centrality object is 0");
299       centrality = -1;
300      }
301   }
302   
303   // Support for ESD and AOD based analysis
304   AliVEvent* inputEvent = fAOD;
305   if (!inputEvent)
306     inputEvent = fESD;
307   
308   TObject* mc = fArrayMC;
309   if (!mc)
310     mc = fMcEvent;
311   
312   // count all events
313   fHistos->FillEvent(centrality, -1);
314   
315   // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
316   if (!fAnalyseUE->VertexSelection(fMcEvent, 0, fZVertex)) 
317     return;
318   
319   Float_t zVtx = fMcEvent->GetPrimaryVertex()->GetZ();
320     
321   // Get MC primaries
322   TObjArray* tracksMC = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, -1, kTRUE);
323   
324   // (MC-true all particles)
325   // STEP 0
326   fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC);
327   
328   // Trigger selection ************************************************
329   if (fAnalyseUE->TriggerSelection(fInputHandler))
330   {  
331     fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
332     
333     // Vertex selection *************************************************
334     if (fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex))
335     {
336       // fill here for tracking efficiency
337       // loop over particle species
338       
339       for (Int_t particleSpecies = 0; particleSpecies < 4; particleSpecies++)
340       {
341         TObjArray* primMCParticles = fAnalyseUE->GetAcceptedParticles(mc, 0x0, kTRUE, particleSpecies, kTRUE);
342         TObjArray* primRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE);
343         TObjArray* allRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE);
344       
345         fHistos->FillTrackingEfficiency(primMCParticles, primRecoTracksMatched, allRecoTracksMatched, particleSpecies, centrality);
346         
347         delete primMCParticles;
348         delete primRecoTracksMatched;
349         delete allRecoTracksMatched;
350       }
351     
352       // (MC-true all particles)
353       // STEP 2
354       if (!fReduceMemoryFootprint)
355         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepVertex, tracksMC);
356       
357       // Get MC primaries that match reconstructed track
358       TObjArray* tracksRecoMatchedPrim = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, -1, kTRUE);
359       
360       // (RECO-matched (quantities from MC particle) primary particles)
361       // STEP 4
362       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim);
363       
364       // Get MC primaries + secondaries that match reconstructed track
365       TObjArray* tracksRecoMatchedAll = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, -1, kTRUE);
366       
367       // (RECO-matched (quantities from MC particle) all particles)
368       // STEP 5
369       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll);
370       
371       // Get RECO tracks
372       TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
373       
374       // (RECO all tracks)
375       // STEP 6
376       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks);
377       
378       if (0 && !fReduceMemoryFootprint)
379       {
380         // make list of secondaries (matched with MC)
381         TObjArray* tracksRecoMatchedSecondaries = new TObjArray;
382         for (Int_t i=0; i<tracksRecoMatchedAll->GetEntries(); i++)
383           if (((AliAODMCParticle*)tracksRecoMatchedAll->At(i))->IsPhysicalPrimary() == kFALSE)
384             tracksRecoMatchedSecondaries->Add(tracksRecoMatchedAll->At(i));
385       
386         // Study: Use only secondaries as trigger particles and plot the correlation vs. all particles; store in step 9
387         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy2, tracksRecoMatchedSecondaries, tracksRecoMatchedAll);
388         
389         // Study: Use only primaries as trigger particles and plot the correlation vs. secondaries; store in step 8
390         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksRecoMatchedPrim, tracksRecoMatchedSecondaries);
391       
392         // plot delta phi vs process id of secondaries
393         // trigger particles: primaries in 4 < pT < 10
394         // associated particles: secondaries in 1 < pT < 10
395         
396         for (Int_t i=0; i<tracksRecoMatchedPrim->GetEntries(); i++)
397         {
398           AliVParticle* triggerParticle = (AliVParticle*) tracksRecoMatchedPrim->At(i);
399           
400           if (triggerParticle->Pt() < 4 || triggerParticle->Pt() > 10)
401             continue;
402           
403           for (Int_t j=0; j<tracksRecoMatchedSecondaries->GetEntries(); j++)
404           {
405             AliAODMCParticle* particle = (AliAODMCParticle*) tracksRecoMatchedSecondaries->At(j);
406             
407             if (particle->Pt() < 1 || particle->Pt() > 10)
408               continue;
409             
410             if (particle->Pt() > triggerParticle->Pt())
411               continue;
412               
413             Double_t deltaPhi = triggerParticle->Phi() - particle->Phi();
414             if (deltaPhi > 1.5 * TMath::Pi()) 
415               deltaPhi -= TMath::TwoPi();
416             if (deltaPhi < -0.5 * TMath::Pi())
417               deltaPhi += TMath::TwoPi();
418               
419             Int_t processID = fMcEvent->Stack()->Particle(particle->GetLabel())->GetUniqueID();
420               
421             ((TH2F*) fListOfHistos->FindObject("processIDs"))->Fill(deltaPhi, processID);
422           }
423         }
424         
425         delete tracksRecoMatchedSecondaries;
426       }
427   
428       delete tracksRecoMatchedPrim;
429       delete tracksRecoMatchedAll;
430       delete tracks;
431     }
432   }
433   
434   delete tracksMC;
435 }
436
437 //____________________________________________________________________
438 void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
439 {
440
441   // Run the analysis on DATA or MC to get raw distributions
442  
443   if ( fDebug > 3 ) AliInfo( " Processing event in Data mode ..." );
444
445   // skip not selected events here (the AOD is not updated for those)
446   if (!fInputHandler)
447     return;
448     
449   if (!(fInputHandler->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kUserDefined)))
450     return;
451
452   Double_t centrality = 0;
453   
454   if (fCentralityMethod.Length() > 0)
455   {
456     AliCentrality *centralityObj = 0;
457     if (fAOD)
458       centralityObj = fAOD->GetHeader()->GetCentralityP();
459     else if (fESD)
460       centralityObj = fESD->GetCentrality();
461     
462     if (centralityObj)
463       centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
464       //centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
465     else
466       centrality = -1;
467     Printf("Centrality is %f", centrality);  
468   }
469   
470   // Support for ESD and AOD based analysis
471   AliVEvent* inputEvent = fAOD;
472   if (!inputEvent)
473     inputEvent = fESD;
474
475   fHistos->SetRunNumber(inputEvent->GetRunNumber());
476   
477   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
478   fHistos->FillEvent(centrality, AliUEHist::kCFStepAll);
479   
480   // Trigger selection ************************************************
481   if (!fAnalyseUE->TriggerSelection(fInputHandler)) return;
482   
483   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
484   fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
485   
486   // Vertex selection *************************************************
487   if(!fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex)) return;
488   
489   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
490   fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
491  
492   if (centrality < 0)
493     return;
494
495   TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
496   //Printf("Accepted %d tracks", tracks->GetEntries());
497   
498   // event mixing
499   
500   // 1. First get an event pool corresponding in mult (cent) and
501   //    zvertex to the current event. Once initialized, the pool
502   //    should contain nMix (reduced) events. This routine does not
503   //    pre-scan the chain. The first several events of every chain
504   //    will be skipped until the needed pools are filled to the
505   //    specified depth. If the pool categories are not too rare, this
506   //    should not be a problem. If they are rare, you could lose
507   //    statistics.
508
509   // 2. Collect the whole pool's content of tracks into one TObjArray
510   //    (bgTracks), which is effectively a single background super-event.
511
512   // 3. The reduced and bgTracks arrays must both be passed into
513   //    FillCorrelations(). Also nMix should be passed in, so a weight
514   //    of 1./nMix can be applied.
515
516   const AliVVertex* vertex = inputEvent->GetPrimaryVertex();
517   Double_t zVtx = vertex->GetZ();
518   
519   AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
520   
521   if (!pool)
522     AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, zVtx));
523   
524   //pool->SetDebug(1);
525     
526   // Fill containers at STEP 6 (reconstructed)
527   fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks);
528   
529   if (pool->IsReady()) 
530   {
531     
532     Int_t nMix = pool->GetCurrentNEvents();
533     //cout << "nMix = " << nMix << endl;
534   
535     // Fill mixed-event histos here  
536     for (Int_t jMix=0; jMix<nMix; jMix++) {
537       TObjArray* bgTracks = pool->GetEvent(jMix);
538       fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, bgTracks, 1.0 / nMix, (jMix == 0));
539     }
540   }
541   
542   // clone and give ownership to event pool
543   TObjArray* tracksClone = (TObjArray*) tracks->Clone();
544   tracksClone->SetOwner(kTRUE);
545   
546   pool->UpdatePool(tracksClone);
547   //pool->PrintInfo();
548
549   delete tracks;
550 }
551
552 //____________________________________________________________________
553 void  AliAnalysisTaskPhiCorrelations::Initialize()
554 {
555   // input handler
556   fInputHandler = (AliInputEventHandler*)
557          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
558   // MC handler
559   fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
560 }