a469338197be74dbd9cb807c29aba39d9bca3646
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskUE.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 //class AliAODInputHandler;
19 //class AliAODHandler;
20 //class AliLog;
21 //class TROOT;
22
23 #include <TROOT.h>
24 #include <TChain.h>
25 #include <TFile.h>
26 #include <TList.h>
27 #include <TMath.h>
28 #include <TProfile.h>
29 #include <TTree.h>
30 #include <TVector3.h>
31 #include <TH1F.h>
32
33 #include "AliAnalyseUE.h"
34 #include "AliHistogramsUE.h"
35 #include "AliAnalysisTaskUE.h"
36 #include "AliHistogramsUE.h"
37
38 #include "AliAnalysisManager.h"
39 #include "AliMCEventHandler.h"
40
41 #include "AliAnalysisHelperJetTasks.h"
42 #include "AliAODHandler.h"
43 #include "AliAODInputHandler.h"
44 #include "AliGenPythiaEventHeader.h"
45 #include "AliLog.h"
46 #include "AliInputEventHandler.h"
47
48
49 ////////////////////////////////////////////////////////////////////////
50 //
51 // Analysis class for Underlying Event studies
52 //
53 // Look for correlations on the tranverse regions to 
54 // the leading charged jet
55 //
56 // This class needs as input AOD with track and Jets.
57 // The output is a list of histograms
58 //
59 // AOD can be either connected to the InputEventHandler  
60 // for a chain of AOD files 
61 // or 
62 // to the OutputEventHandler
63 // for a chain of ESD files, so this case class should be 
64 // in the train after the Jet finder
65 //
66 //    Arian.Abrahantes.Quintana@cern.ch 
67 //    Ernesto.Lopez.Torres@cern.ch
68 //    vallero@physi.uni-heidelberg.de
69 // 
70 ////////////////////////////////////////////////////////////////////////
71
72 ClassImp( AliAnalysisTaskUE)
73
74 // Define global pointer
75 AliAnalysisTaskUE* AliAnalysisTaskUE::fgTaskUE=NULL;
76
77 //____________________________________________________________________
78 AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name):
79 AliAnalysisTask(name,""),
80 fHistosUE(0x0),
81 fAnaUE(0x0),
82 fAOD(0x0),            
83 fAODBranch("jets"),
84 fDebug(0),
85 fListOfHistos(0x0), 
86 //Configuration
87 fBinsPtInHist(30),     
88 fIsNorm2Area(kTRUE),
89 fMaxJetPtInHist(300.), 
90 fMinJetPtInHist(0.),
91 fConstrainDistance(kTRUE),
92 fMinDistance(0.2),
93 fSimulateChJetPt(kFALSE),
94 fUseAliStack(kTRUE),
95 fUseMCParticleBranch(kFALSE),
96 fnTracksVertex(3),  // QA tracks pointing to principal vertex (= 3 default) 
97 fZVertex(5.),
98 fAnaType(1),         
99 fConePosition(1),
100 fConeRadius(0.7),
101 fFilterBit(0xFF),
102 fJetsOnFly(kFALSE),
103 fRegionType(1),
104 fUseChargeHadrons(kFALSE),
105 fUseChPartJet(kFALSE),
106 fUsePositiveCharge(kTRUE),
107 fUseSingleCharge(kFALSE),
108 fOrdering(1),
109 fChJetPtMin(5.0),
110 fJet1EtaCut(0.2),
111 fJet2DeltaPhiCut(2.616),    // 150 degrees
112 fJet2RatioPtCut(0.8),
113 fJet3PtCut(15.),
114 fTrackEtaCut(0.9),
115 fTrackPtCut(0.),
116 //For MC
117 fAvgTrials(1)
118 /*//Histograms
119 fhNJets(0x0),
120 fhEleadingPt(0x0),
121 fhMinRegPtDist(0x0),
122 fhRegionMultMin(0x0),
123 fhMinRegAvePt(0x0),
124 fhMinRegSumPt(0x0),
125 fhMinRegMaxPtPart(0x0),
126 fhMinRegSumPtvsMult(0x0),
127 fhdNdEtaPhiDist(0x0),
128 fhFullRegPartPtDistVsEt(0x0),
129 fhTransRegPartPtDistVsEt(0x0),
130 fhRegionSumPtMaxVsEt(0x0),
131 fhRegionMultMax(0x0),
132 fhRegionMultMaxVsEt(0x0),
133 fhRegionSumPtMinVsEt(0x0),
134 fhRegionMultMinVsEt(0x0),
135 fhRegionAveSumPtVsEt(0x0),
136 fhRegionDiffSumPtVsEt(0x0),
137 fhRegionAvePartPtMaxVsEt(0x0),
138 fhRegionAvePartPtMinVsEt(0x0),
139 fhRegionMaxPartPtMaxVsEt(0x0),    
140 fhRegForwardMult(0x0),
141 fhRegForwardSumPtvsMult(0x0),
142 fhRegBackwardMult(0x0),
143 fhRegBackwardSumPtvsMult(0x0),
144 fhRegForwardPartPtDistVsEt(0x0),
145 fhRegBackwardPartPtDistVsEt(0x0),
146 fhRegTransMult(0x0),
147 fhRegTransSumPtVsMult(0x0),
148 fhMinRegSumPtJetPtBin(0x0),
149 fhMaxRegSumPtJetPtBin(0x0),
150 fhVertexMult(0x0),
151 fh1Xsec(0x0),   
152 fh1Trials(0x0)
153 */
154 {
155   // Default constructor
156   // Define input and output slots here
157   // Input slot #0 works with a TChain
158   DefineInput(0, TChain::Class());
159   // Output slot #0 writes into a TList container
160   DefineOutput(0, TList::Class());
161
162 }
163
164 //____________________________________________________________________
165 AliAnalysisTaskUE:: AliAnalysisTaskUE(const AliAnalysisTaskUE & original):
166 AliAnalysisTask(),
167 fHistosUE(original.fHistosUE),
168 fAnaUE(original.fAnaUE),
169 fAOD(original.fAOD),            
170 fAODBranch(original.fAODBranch),
171 fDebug(original.fDebug),
172 fListOfHistos(original.fListOfHistos),  
173 //Configuration
174 fBinsPtInHist(original.fBinsPtInHist),     
175 fIsNorm2Area(original.fIsNorm2Area),
176 fMaxJetPtInHist(original.fMaxJetPtInHist), 
177 fMinJetPtInHist(original.fMinJetPtInHist),
178 fConstrainDistance(original.fConstrainDistance),
179 fMinDistance(original.fMinDistance),
180 fSimulateChJetPt(original.fSimulateChJetPt),
181 fUseAliStack(original.fUseAliStack),
182 fUseMCParticleBranch(original.fUseMCParticleBranch),
183 fnTracksVertex(original.fnTracksVertex),  // QA tracks pointing to principal vertex (= 3 default) 
184 fZVertex(original.fZVertex),
185 fAnaType(original.fAnaType),         
186 fConePosition(original.fConePosition),
187 fConeRadius(original.fConeRadius),
188 fFilterBit(original.fFilterBit),
189 fJetsOnFly(original.fJetsOnFly),
190 fRegionType(original.fRegionType),
191 fUseChargeHadrons(original.fUseChargeHadrons),
192 fUseChPartJet(original.fUseChPartJet),
193 fUsePositiveCharge(original.fUsePositiveCharge),
194 fUseSingleCharge(original.fUseSingleCharge),
195 fOrdering(original.fOrdering),
196 fChJetPtMin(original.fChJetPtMin),
197 fJet1EtaCut(original.fJet1EtaCut),
198 fJet2DeltaPhiCut(original.fJet2DeltaPhiCut),    // 150 degrees
199 fJet2RatioPtCut(original.fJet2RatioPtCut),
200 fJet3PtCut(original.fJet3PtCut),
201 fTrackEtaCut(original.fTrackEtaCut),
202 fTrackPtCut(original.fTrackPtCut),
203 //For MC
204 fAvgTrials(original.fAvgTrials)
205 /*//Histograms
206 fhNJets(original.fhNJets),
207 fhEleadingPt(original.fhEleadingPt),
208 fhMinRegPtDist(original.fhMinRegPtDist),
209 fhRegionMultMin(original.fhRegionMultMin),
210 fhMinRegAvePt(original.fhMinRegAvePt),
211 fhMinRegSumPt(original.fhMinRegSumPt),
212 fhMinRegMaxPtPart(original.fhMinRegMaxPtPart),
213 fhMinRegSumPtvsMult(original.fhMinRegSumPtvsMult),
214 fhdNdEtaPhiDist(original.fhdNdEtaPhiDist),
215 fhFullRegPartPtDistVsEt(original.fhFullRegPartPtDistVsEt),
216 fhTransRegPartPtDistVsEt(original.fhTransRegPartPtDistVsEt),
217 fhRegionSumPtMaxVsEt(original.fhRegionSumPtMaxVsEt),
218 fhRegionMultMax(original.fhRegionMultMax),
219 fhRegionMultMaxVsEt(original.fhRegionMultMaxVsEt),
220 fhRegionSumPtMinVsEt(original.fhRegionSumPtMinVsEt),
221 fhRegionMultMinVsEt(original.fhRegionMultMinVsEt),
222 fhRegionAveSumPtVsEt(original.fhRegionAveSumPtVsEt),
223 fhRegionDiffSumPtVsEt(original.fhRegionDiffSumPtVsEt),
224 fhRegionAvePartPtMaxVsEt(original.fhRegionAvePartPtMaxVsEt),
225 fhRegionAvePartPtMinVsEt(original.fhRegionAvePartPtMinVsEt),
226 fhRegionMaxPartPtMaxVsEt(original.fhRegionMaxPartPtMaxVsEt),
227 fhRegForwardMult(original.fhRegForwardMult),
228 fhRegForwardSumPtvsMult(original.fhRegForwardSumPtvsMult),
229 fhRegBackwardMult(original.fhRegBackwardMult),
230 fhRegBackwardSumPtvsMult(original.fhRegBackwardSumPtvsMult),
231 fhRegForwardPartPtDistVsEt(original.fhRegForwardPartPtDistVsEt),
232 fhRegBackwardPartPtDistVsEt(original.fhRegBackwardPartPtDistVsEt),
233 fhRegTransMult(original.fhRegTransMult),
234 fhRegTransSumPtVsMult(original.fhRegTransSumPtVsMult),
235 fhMinRegSumPtJetPtBin(original.fhMinRegSumPtJetPtBin),
236 fhMaxRegSumPtJetPtBin(original.fhMaxRegSumPtJetPtBin),
237 fhVertexMult(original.fhVertexMult),
238 fh1Xsec(original.fh1Xsec),      
239 fh1Trials(original.fh1Trials)
240 */
241 {
242   // Copy constructor
243 }
244
245
246 //______________________________________________________________
247 AliAnalysisTaskUE & AliAnalysisTaskUE::operator = (const AliAnalysisTaskUE & /*source*/)
248 {
249   // assignment operator
250   return *this;
251 }
252
253 //______________________________________________________________
254 Bool_t AliAnalysisTaskUE::Notify()
255 {
256   //
257   // Implemented Notify() to read the cross sections
258   // and number of trials from pyxsec.root
259   // Copy from AliAnalysisTaskJFSystematics
260   fAvgTrials = 1;
261   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
262   Float_t xsection = 0;
263   Float_t trials  = 1;
264   if(tree){
265         TFile *curfile = tree->GetCurrentFile();
266         if (!curfile) {
267                 Error("Notify","No current file");
268                 return kFALSE;
269                 }
270         
271         AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials); 
272         
273         fHistosUE->GetXsec()->Fill("<#sigma>",xsection);
274
275         // construct average trials
276         Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
277         if(trials>=nEntries && nEntries>0.)fAvgTrials = trials/nEntries;
278         }
279   
280   return kTRUE;
281 }
282
283 //____________________________________________________________________
284 void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
285 {
286   // Connect the input data  
287   
288   // We need AOD with tracks and jets.
289   // Since AOD can be either connected to the InputEventHandler (input chain fron AOD files)
290   // or to the OutputEventHandler ( AOD is create by a previus task in the train )
291   // we need to check where it is and get the pointer to AODEvent in the right way
292   
293   // Delta AODs are also accepted
294   
295  
296   if (fDebug > 1) AliInfo("ConnectInputData() ");
297   
298   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
299   
300   if( handler && handler->InheritsFrom("AliAODInputHandler") ) { // input AOD
301         fAOD = ((AliAODInputHandler*)handler)->GetEvent();
302         if(!fJetsOnFly){
303                 if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
304                 }else{
305                 if (fDebug > 1) AliInfo(" ==== Tracks from AliAODInputHandler / Jets on-the-fly");
306                 }
307         } else {  //output AOD
308         handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
309         if( handler && handler->InheritsFrom("AliAODHandler") ) {
310                 fAOD = ((AliAODHandler*)handler)->GetAOD();
311                 if (!fJetsOnFly){
312                         if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
313                         } else {
314                         if (fDebug > 1) AliInfo(" ==== Tracks from AliAODHandler / Jets on-the-fly");
315                         }
316                 }else {
317                 AliFatal("I can't get any AOD Event Handler");
318                 return;
319                 }
320         }       
321
322    fAnaUE->Initialize( *this );
323    
324 }
325
326 //____________________________________________________________________
327 void  AliAnalysisTaskUE::CreateOutputObjects()
328 {
329   // Create the output container
330   
331   if (fDebug > 1) AliInfo("CreateOutPutData()");
332    
333   //Initialize AliAnalysisUE, a class implementing the main analysis algorithms
334   fAnaUE = new AliAnalyseUE();
335   fHistosUE = new AliHistogramsUE();
336  
337   if (fListOfHistos != NULL){
338         delete fListOfHistos;
339         fListOfHistos = NULL;
340         }
341   if (!fListOfHistos){
342         fListOfHistos = new TList();
343         fListOfHistos->SetOwner(kTRUE); 
344         }
345   
346   //Initialize output histograms
347   fHistosUE->CreateHistograms(fListOfHistos,fBinsPtInHist, fMinJetPtInHist, fMaxJetPtInHist, fTrackEtaCut);
348   AddSettingsTree();
349
350   PostData(0,fListOfHistos);
351 }
352
353 //____________________________________________________________________
354 void  AliAnalysisTaskUE::Exec(Option_t */*option*/)
355 {
356   // Trigger selection ************************************************
357   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
358          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
359   if( !inputHandler->InheritsFrom("AliAODInputHandler") ) { // input AOD
360         if (inputHandler->IsEventSelected()) {
361                 if (fDebug > 1) AliInfo(" Trigger Selection: event ACCEPTED ... ");
362         } else {
363                 if (fDebug > 1) AliInfo(" Trigger Selection: event REJECTED ... ");
364                 return;
365                 }
366         }                                
367   // Event selection (vertex) *****************************************
368    
369   //if(!fAnaUE->VertexSelection(fAOD,fnTracksVertex,fZVertex)) return;
370   if(!fAnaUE->VertexSelectionOld(fAOD)) return; // temporary to compare with old task and to have same cuts for MC !!!
371   
372   // Execute analysis for current event ******************************
373   
374   if ( fDebug > 3 ) AliInfo( " Processing event..." );
375
376   // fetch the pythia header info and get the trials
377   AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
378   Float_t nTrials = 1;
379   if (mcHandler) {
380         AliMCEvent* mcEvent = mcHandler->MCEvent();
381         if (mcEvent) {
382                 AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
383                 if(pythiaGenHeader){
384                         nTrials = pythiaGenHeader->Trials();
385                         }
386                 }
387         }
388   fHistosUE->GetTrials()->Fill("#sum{ntrials}",fAvgTrials);
389   
390   //analyse the event
391   AnalyseUE();
392  
393  PostData(0,fListOfHistos);
394 }
395
396 //____________________________________________________________________
397 void  AliAnalysisTaskUE::AddSettingsTree()
398 {
399   //Write settings to output list
400   TTree *settingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
401   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
402   settingsTree->Branch("fConeRadius", &fConeRadius,"Rad/D");
403   settingsTree->Branch("fJet1EtaCut", &fJet1EtaCut, "LeadJetEtaCut/D");
404   settingsTree->Branch("fJet2DeltaPhiCut", &fJet2DeltaPhiCut, "DeltaPhi/D");
405   settingsTree->Branch("fJet2RatioPtCut", &fJet2RatioPtCut, "Jet2Ratio/D");
406   settingsTree->Branch("fJet3PtCut", &fJet3PtCut, "Jet3PtCut/D");
407   settingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D");
408   settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
409   settingsTree->Branch("fAnaType", &fAnaType, "Ana/I");        
410   settingsTree->Branch("fRegionType", &fRegionType,"Reg/I");
411   settingsTree->Branch("fOrdering", &fOrdering,"OrderMeth/I");
412   settingsTree->Branch("fUseChPartJet", &fUseChPartJet,"UseChPart/O");
413   settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
414   settingsTree->Branch("fUseSingleCharge", &fUseSingleCharge,"UseSingleCh/O");
415   settingsTree->Branch("fUsePositiveCharge", &fUsePositiveCharge,"UsePositiveCh/O");
416   settingsTree->Fill();
417   fListOfHistos->Add(settingsTree);
418 }  
419
420 //____________________________________________________________________
421 void  AliAnalysisTaskUE::AnalyseUE()
422 {
423
424   
425    // Get jets and order by pT
426    TVector3 jetVect[3];
427    *jetVect = fAnaUE->GetOrderedClusters(fAODBranch, fUseChPartJet, fChJetPtMin );
428   
429    if( jetVect[0].Pt() < 0. ) {
430         if( fDebug > 1 ) AliInfo("\n   Skipping Event, not jet found...");
431         return;
432         } else {
433                 if (fDebug >1 ) AliInfo(Form("\n   Pt Leading Jet = %6.1f eta=%5.3f ",  jetVect[0].Pt(), jetVect[0].Eta() ));
434                 }
435
436    // Select events according to analysis type
437    if ( ! (fAnaUE->AnaTypeSelection(jetVect))) return;
438
439   // Find max and min regions with real tracks
440   if (!fUseMCParticleBranch){
441         // Primary vertex distribution
442         AliAODVertex* vertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
443         fHistosUE->FillHistogram("hVertexMult",vertex->GetNContributors());
444         
445         // Fill leading "jet" histogram
446         fHistosUE->FillHistogram("hEleadingPt",jetVect[0].Pt());
447
448         fAnaUE->FindMaxMinRegions( jetVect,  fConePosition );
449
450   }else { 
451     
452     // this is the part we only use when we have MC information
453     // More than a test for values of it also resumes the reconstruction efficiency of jets
454     // As commented bellow if available for the data, we try to pair reconstructed jets with simulated ones
455     // afterwards we kept angular variables of MC jet to perform UE analysis over MC particles
456     // TODO: Handle Multiple jet environment. 06/2009 just suited for inclusive jet condition ( fAnaType = 1 ) 
457     
458       AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
459       if (!mcHandler) {
460         Printf("ERROR: Could not retrieve MC event handler");
461         return;
462       }
463       
464       AliMCEvent* mcEvent = mcHandler->MCEvent();
465       if (!mcEvent) {
466         Printf("ERROR: Could not retrieve MC event");
467         return;
468       }
469     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
470     if(!pythiaGenHeader){
471       return;
472     }
473     
474     fAnaUE->AnalyseMC(jetVect,mcEvent,pythiaGenHeader, fConePosition, fUseAliStack, fConstrainDistance, fMinDistance);
475
476   }
477
478   fAnaUE->FillRegions(fIsNorm2Area, jetVect);
479   
480 }
481
482 //____________________________________________________________________
483 AliAnalysisTaskUE* AliAnalysisTaskUE::Instance()
484
485   
486   //Create instance
487   if (fgTaskUE) {
488         return fgTaskUE;
489   } else {
490         fgTaskUE = new AliAnalysisTaskUE();
491         return fgTaskUE;
492         }
493 }
494
495
496 //____________________________________________________________________
497 void  AliAnalysisTaskUE::Terminate(Option_t */*option*/)
498 {
499   
500   // Terminate analysis
501
502   fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
503   if (!fListOfHistos){
504         AliError("Histogram List is not available");
505         return;
506   }
507   if (!gROOT->IsBatch()){
508   //call class AliHistogramsUE
509   AliHistogramsUE *histos=new AliHistogramsUE(fListOfHistos);
510   histos->DrawUE(0);
511   } else {
512         AliInfo(" Batch mode, not histograms will be shown...");
513   }
514
515   if( fDebug > 1 ) 
516     AliInfo("End analysis");
517  
518 }
519
520 void  AliAnalysisTaskUE::WriteSettings()
521
522
523   //Print analysis settings on screen
524   if (fDebug > 5){
525     AliInfo(" All Analysis Settings in Saved Tree");
526     fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
527     if (!fListOfHistos){
528         AliError("Histogram List is not available");
529         return;
530         }
531      TTree *tree = (TTree*)(fListOfHistos->FindObject("UEAnalysisSettings"));
532      tree->Scan();
533      }
534 }