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