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