New book-keeping class
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetServices.cxx
1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets   
4 // *******************************************
5
6
7 /**************************************************************************
8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9  *                                                                        *
10  * Author: The ALICE Off-line Project.                                    *
11  * Contributors are mentioned in the code where appropriate.              *
12  *                                                                        *
13  * Permission to use, copy, modify and distribute this software and its   *
14  * documentation strictly for non-commercial purposes is hereby granted   *
15  * without fee, provided that the above copyright notice appears in all   *
16  * copies and that both the copyright notice and this permission notice   *
17  * appear in the supporting documentation. The authors make no claims     *
18  * about the suitability of this software for any purpose. It is          *
19  * provided "as is" without express or implied warranty.                  *
20  **************************************************************************/
21
22  
23 #include <TROOT.h>
24 #include <TRandom.h>
25 #include <TSystem.h>
26 #include <TInterpreter.h>
27 #include <TChain.h>
28 #include <TFile.h>
29 #include <TKey.h>
30 #include <TH1F.h>
31 #include <TH2F.h>
32 #include <TH3F.h>
33 #include <TProfile.h>
34 #include <TList.h>
35 #include <TLorentzVector.h>
36 #include <TClonesArray.h>
37 #include  "TDatabasePDG.h"
38
39 #include "AliAnalysisTaskJetServices.h"
40 #include "AliAnalysisManager.h"
41 #include "AliJetFinder.h"
42 #include "AliJetHeader.h"
43 #include "AliJetReader.h"
44 #include "AliJetReaderHeader.h"
45 #include "AliUA1JetHeaderV1.h"
46 #include "AliESDEvent.h"
47 #include "AliAODEvent.h"
48 #include "AliAODHandler.h"
49 #include "AliAODTrack.h"
50 #include "AliAODJet.h"
51 #include "AliAODMCParticle.h"
52 #include "AliMCEventHandler.h"
53 #include "AliMCEvent.h"
54 #include "AliStack.h"
55 #include "AliGenPythiaEventHeader.h"
56 #include "AliJetKineReaderHeader.h"
57 #include "AliGenCocktailEventHeader.h"
58 #include "AliInputEventHandler.h"
59
60
61 #include "AliAnalysisHelperJetTasks.h"
62
63 ClassImp(AliAnalysisTaskJetServices)
64
65 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
66   fUseAODInput(kFALSE),
67   fAvgTrials(1),
68   fZVtxCut(10.),
69   fh1Xsec(0x0),
70   fh1Trials(0x0),
71   fh1PtHard(0x0),
72   fh1PtHardTrials(0x0),
73   fh2TriggerCount(0x0),
74   fh2ESDTriggerCount(0x0),
75   fh2TriggerVtx(0x0),
76   fh2ESDTriggerVtx(0x0),
77   fHistList(0x0)  
78 {
79
80 }
81
82 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
83   AliAnalysisTaskSE(name),
84   fUseAODInput(kFALSE),
85   fAvgTrials(1),
86   fZVtxCut(8),
87   fh1Xsec(0x0),
88   fh1Trials(0x0),
89   fh1PtHard(0x0),
90   fh1PtHardTrials(0x0),
91   fh2TriggerCount(0x0),
92   fh2ESDTriggerCount(0x0),
93   fh2TriggerVtx(0x0),
94   fh2ESDTriggerVtx(0x0),
95   fHistList(0x0)  
96 {
97   DefineOutput(1,TList::Class());
98 }
99
100
101
102 Bool_t AliAnalysisTaskJetServices::Notify()
103 {
104   //
105   // Implemented Notify() to read the cross sections
106   // and number of trials from pyxsec.root
107   // 
108
109   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
110   Float_t xsection = 0;
111   Float_t ftrials  = 1;
112
113   fAvgTrials = 1;
114   if(tree){
115     TFile *curfile = tree->GetCurrentFile();
116     if (!curfile) {
117       Error("Notify","No current file");
118       return kFALSE;
119     }
120     if(!fh1Xsec||!fh1Trials){
121       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
122       return kFALSE;
123     }
124     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
125     fh1Xsec->Fill("<#sigma>",xsection);
126     // construct a poor man average trials 
127     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
128     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
129   }  
130   return kTRUE;
131 }
132
133 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
134 {
135
136   //
137   // Create the output container
138   //
139
140
141   // Connect the AOD
142
143   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
144
145   OpenFile(1);
146   if(!fHistList)fHistList = new TList();
147
148   Bool_t oldStatus = TH1::AddDirectoryStatus();
149   TH1::AddDirectory(kFALSE);
150
151   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
152   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
153   fHistList->Add(fh1Xsec);
154
155   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
156   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
157   fHistList->Add(fh1Trials);
158
159   const Int_t nBinPt = 100;
160   Double_t binLimitsPt[nBinPt+1];
161   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
162     if(iPt == 0){
163       binLimitsPt[iPt] = 0.0;
164     }
165     else {// 1.0
166       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
167     }
168   }
169   
170   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
171   fHistList->Add(fh2TriggerCount);
172
173   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
174   fHistList->Add(fh2ESDTriggerCount);
175
176   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); 
177   fHistList->Add(fh2TriggerVtx);
178
179   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); 
180   fHistList->Add(fh2ESDTriggerVtx);
181   
182
183   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
184   fHistList->Add(fh1PtHard);
185   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
186   fHistList->Add(fh1PtHardTrials);
187
188   // =========== Switch on Sumw2 for all histos ===========
189   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
190     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
191     if (h1){
192       h1->Sumw2();
193       continue;
194     }
195     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
196     if(hn)hn->Sumw2();
197   }
198
199
200   TH1::AddDirectory(oldStatus);
201 }
202
203 void AliAnalysisTaskJetServices::Init()
204 {
205   //
206   // Initialization
207   //
208
209   Printf(">>> AnalysisTaskJetServices::Init() debug level %d\n",fDebug);
210   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
211
212 }
213
214 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
215 {
216
217   //
218   // Execute analysis for current event
219   //
220  
221   AliAODEvent *aod = 0;
222   AliESDEvent *esd = 0;
223
224   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
225
226   if(fUseAODInput){    
227     aod = dynamic_cast<AliAODEvent*>(InputEvent());
228     if(!aod){
229       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
230       return;
231     }
232     // fethc the header
233   }
234   else{
235     //  assume that the AOD is in the general output...
236     aod  = AODEvent();
237     if(!aod){
238       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
239       return;
240     }
241     esd = dynamic_cast<AliESDEvent*>(InputEvent());
242   }
243   
244   
245   // loop over all possible trigger and 
246   for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
247     Bool_t esdTrig = kFALSE;
248     if(esd){
249       esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
250       if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
251     }
252     Bool_t aodTrig = kFALSE;
253     if(aod){
254       aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
255       if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
256     }
257
258     // Check wether we have also an SPD vertex
259     
260     if(aod){
261       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
262       //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
263       
264       if(vtxAOD->GetNContributors()>0){
265         if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
266         Float_t zvtx = vtxAOD->GetZ();
267         Float_t yvtx = vtxAOD->GetY();
268         Float_t xvtx = vtxAOD->GetX();
269         fh2TriggerVtx->Fill(it,zvtx);
270         if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
271           fh2TriggerCount->Fill(it,kTriggeredVertexIn);
272         }
273       }
274     }
275     if(esd){
276       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
277       //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
278       if(vtxESD->GetNContributors()>0){
279         if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
280         Float_t zvtx = vtxESD->GetZ();
281         Float_t yvtx = vtxESD->GetY();
282         Float_t xvtx = vtxESD->GetX();
283         fh2ESDTriggerVtx->Fill(it,zvtx);
284         if(TMath::Abs(zvtx)<fZVtxCut&&esdTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
285           fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
286           // here we select based on ESD info...
287           fh2ESDTriggerCount->Fill(it,kSelected);
288           AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
289         }
290       }
291
292     }
293
294   }
295
296
297
298   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
299
300   
301   Double_t ptHard = 0; 
302   Double_t nTrials = 1; // Trials for MC trigger 
303
304   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
305   AliMCEvent* mcEvent = MCEvent();
306   //    AliStack *pStack = 0; 
307   if(mcEvent){
308     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
309     if(pythiaGenHeader){
310       nTrials = pythiaGenHeader->Trials();
311       ptHard  = pythiaGenHeader->GetPtHard();
312       int iProcessType = pythiaGenHeader->ProcessType();
313       // 11 f+f -> f+f
314       // 12 f+barf -> f+barf
315       // 13 f+barf -> g+g
316       // 28 f+g -> f+g
317       // 53 g+g -> f+barf
318       // 68 g+g -> g+g
319       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
320       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
321       fh1PtHard->Fill(ptHard);
322       fh1PtHardTrials->Fill(ptHard,nTrials);
323
324     }// if pythia gen header
325   }
326
327   // trigger selection
328   
329
330   PostData(1, fHistList);
331 }
332
333
334 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
335 {
336   // Terminate analysis
337   //
338   if (fDebug > 1) printf("AnalysisJetServices: Terminate() \n");
339 }