Implementing destructor
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetServices.cxx
CommitLineData
bf7b8731 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
63ClassImp(AliAnalysisTaskJetServices)
64
65AliAnalysisTaskJetServices::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
82AliAnalysisTaskJetServices::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
102Bool_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
133void 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
203void 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
214void 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
334void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
335{
336 // Terminate analysis
337 //
338 if (fDebug > 1) printf("AnalysisJetServices: Terminate() \n");
339}