Fixed compilation warning
[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>
b1cd90b5 25#include <TString.h>
bf7b8731 26#include <TSystem.h>
27#include <TInterpreter.h>
28#include <TChain.h>
29#include <TFile.h>
30#include <TKey.h>
31#include <TH1F.h>
32#include <TH2F.h>
33#include <TH3F.h>
34#include <TProfile.h>
35#include <TList.h>
36#include <TLorentzVector.h>
37#include <TClonesArray.h>
38#include "TDatabasePDG.h"
39
40#include "AliAnalysisTaskJetServices.h"
fe669ac6 41#include "AliAnalysisDataContainer.h"
42#include "AliAnalysisDataSlot.h"
bf7b8731 43#include "AliAnalysisManager.h"
44#include "AliJetFinder.h"
45#include "AliJetHeader.h"
46#include "AliJetReader.h"
47#include "AliJetReaderHeader.h"
48#include "AliUA1JetHeaderV1.h"
49#include "AliESDEvent.h"
50#include "AliAODEvent.h"
51#include "AliAODHandler.h"
52#include "AliAODTrack.h"
53#include "AliAODJet.h"
54#include "AliAODMCParticle.h"
55#include "AliMCEventHandler.h"
56#include "AliMCEvent.h"
57#include "AliStack.h"
58#include "AliGenPythiaEventHeader.h"
59#include "AliJetKineReaderHeader.h"
60#include "AliGenCocktailEventHeader.h"
61#include "AliInputEventHandler.h"
fe669ac6 62#include "AliPhysicsSelection.h"
bf7b8731 63
64
65#include "AliAnalysisHelperJetTasks.h"
66
67ClassImp(AliAnalysisTaskJetServices)
68
69AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
70 fUseAODInput(kFALSE),
cc0649e4 71 fUsePhysicsSelection(kFALSE),
72 fRealData(kFALSE),
bf7b8731 73 fAvgTrials(1),
fe669ac6 74 fZVtxCut(8.),
bf7b8731 75 fh1Xsec(0x0),
76 fh1Trials(0x0),
77 fh1PtHard(0x0),
78 fh1PtHardTrials(0x0),
79 fh2TriggerCount(0x0),
80 fh2ESDTriggerCount(0x0),
81 fh2TriggerVtx(0x0),
82 fh2ESDTriggerVtx(0x0),
b5a3f310 83 fh2ESDTriggerRun(0x0),
84 fh2VtxXY(0x0),
bf7b8731 85 fHistList(0x0)
86{
b5a3f310 87 fRunRange[0] = fRunRange[1] = 0;
bf7b8731 88}
89
90AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
91 AliAnalysisTaskSE(name),
92 fUseAODInput(kFALSE),
cc0649e4 93 fUsePhysicsSelection(kFALSE),
94 fRealData(kFALSE),
bf7b8731 95 fAvgTrials(1),
fe669ac6 96 fZVtxCut(8.),
bf7b8731 97 fh1Xsec(0x0),
98 fh1Trials(0x0),
99 fh1PtHard(0x0),
100 fh1PtHardTrials(0x0),
101 fh2TriggerCount(0x0),
102 fh2ESDTriggerCount(0x0),
103 fh2TriggerVtx(0x0),
104 fh2ESDTriggerVtx(0x0),
b5a3f310 105 fh2ESDTriggerRun(0x0),
106 fh2VtxXY(0x0),
bf7b8731 107 fHistList(0x0)
108{
b5a3f310 109 fRunRange[0] = fRunRange[1] = 0;
bf7b8731 110 DefineOutput(1,TList::Class());
111}
112
113
114
115Bool_t AliAnalysisTaskJetServices::Notify()
116{
117 //
118 // Implemented Notify() to read the cross sections
119 // and number of trials from pyxsec.root
120 //
121
122 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
123 Float_t xsection = 0;
124 Float_t ftrials = 1;
125
126 fAvgTrials = 1;
127 if(tree){
128 TFile *curfile = tree->GetCurrentFile();
129 if (!curfile) {
130 Error("Notify","No current file");
131 return kFALSE;
132 }
133 if(!fh1Xsec||!fh1Trials){
134 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
135 return kFALSE;
136 }
137 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
138 fh1Xsec->Fill("<#sigma>",xsection);
139 // construct a poor man average trials
140 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
141 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
142 }
143 return kTRUE;
144}
145
146void AliAnalysisTaskJetServices::UserCreateOutputObjects()
147{
148
149 //
150 // Create the output container
151 //
152
153
bf7b8731 154 if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
155
156 OpenFile(1);
157 if(!fHistList)fHistList = new TList();
158
159 Bool_t oldStatus = TH1::AddDirectoryStatus();
160 TH1::AddDirectory(kFALSE);
bf7b8731 161 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
162 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
163 fHistList->Add(fh1Xsec);
164
165 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
166 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
167 fHistList->Add(fh1Trials);
168
169 const Int_t nBinPt = 100;
170 Double_t binLimitsPt[nBinPt+1];
171 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
172 if(iPt == 0){
173 binLimitsPt[iPt] = 0.0;
174 }
175 else {// 1.0
176 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
177 }
178 }
179
180 fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
181 fHistList->Add(fh2TriggerCount);
182
183 fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
184 fHistList->Add(fh2ESDTriggerCount);
185
186 fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
187 fHistList->Add(fh2TriggerVtx);
188
189 fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
190 fHistList->Add(fh2ESDTriggerVtx);
191
192
193 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
194 fHistList->Add(fh1PtHard);
195 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
196 fHistList->Add(fh1PtHardTrials);
197
b5a3f310 198 // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
199 // 3 triggers BB BE/EB EE
200
201 fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Trigger vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
202 fHistList->Add(fh2ESDTriggerRun);
203
204 fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
205 fHistList->Add(fh2VtxXY);
bf7b8731 206 // =========== Switch on Sumw2 for all histos ===========
207 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
208 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
209 if (h1){
210 h1->Sumw2();
211 continue;
212 }
213 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
214 if(hn)hn->Sumw2();
215 }
216
217
218 TH1::AddDirectory(oldStatus);
219}
220
221void AliAnalysisTaskJetServices::Init()
222{
223 //
224 // Initialization
225 //
bf7b8731 226 if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
227
228}
229
230void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
231{
232
233 //
234 // Execute analysis for current event
235 //
236
237 AliAODEvent *aod = 0;
238 AliESDEvent *esd = 0;
239
240 AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
241
242 if(fUseAODInput){
243 aod = dynamic_cast<AliAODEvent*>(InputEvent());
244 if(!aod){
245 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
246 return;
247 }
248 // fethc the header
249 }
250 else{
251 // assume that the AOD is in the general output...
252 aod = AODEvent();
253 if(!aod){
254 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
255 return;
256 }
257 esd = dynamic_cast<AliESDEvent*>(InputEvent());
258 }
259
b5a3f310 260 if(esd){
261 Float_t run = (Float_t)esd->GetRunNumber();
262 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
263 Float_t zvtx = -999;
264 Float_t xvtx = -999;
265 Float_t yvtx = -999;
266
267 if(vtxESD->GetNContributors()>0){
268 zvtx = vtxESD->GetZ();
269 yvtx = vtxESD->GetY();
270 xvtx = vtxESD->GetX();
271 }
272
273 Int_t iTrig = -1;
274 if(esd->GetFiredTriggerClasses().Contains("CINT1B")
275 ||esd->GetFiredTriggerClasses().Contains("CSMBB")
276 ||esd->GetFiredTriggerClasses().Contains("MB1")
277 ||esd->GetFiredTriggerClasses().Contains("CINT6B")){
278 iTrig = 0;
279 }
280 else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
281 ||esd->GetFiredTriggerClasses().Contains("CSMBA")
282 ||esd->GetFiredTriggerClasses().Contains("CINT6A")
283 ||esd->GetFiredTriggerClasses().Contains("CINT1C")
284 ||esd->GetFiredTriggerClasses().Contains("CSMBC")
285 ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
286 // empty bunch or bunch empty
287 iTrig = 1;
288 }
289 if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
290 ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
291 iTrig = 2;
292 }
293
294
295 if(iTrig>=0){
296 iTrig *= 3;
297 fh2ESDTriggerRun->Fill(run,iTrig+1);
298 if(vtxESD->GetNContributors()>0){
299 fh2ESDTriggerRun->Fill(run,iTrig+2);
300 fh2VtxXY->Fill(xvtx,yvtx);
301 }
cf049e94 302 Float_t r2 = xvtx *xvtx + yvtx *yvtx;
303 if(TMath::Abs(zvtx)<fZVtxCut&&r2<1)fh2ESDTriggerRun->Fill(run,iTrig+3);
b5a3f310 304 }
305 else{
306 fh2ESDTriggerRun->Fill(run,0);
307 }
308
309
310 }
bf7b8731 311
312 // loop over all possible trigger and
313 for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
314 Bool_t esdTrig = kFALSE;
315 if(esd){
316 esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
317 if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
318 }
319 Bool_t aodTrig = kFALSE;
320 if(aod){
321 aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
322 if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
323 }
324
325 // Check wether we have also an SPD vertex
326
327 if(aod){
328 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
329 // Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
330
331 if(vtxAOD->GetNContributors()>0){
332 if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
333 Float_t zvtx = vtxAOD->GetZ();
334 Float_t yvtx = vtxAOD->GetY();
335 Float_t xvtx = vtxAOD->GetX();
336 fh2TriggerVtx->Fill(it,zvtx);
cf049e94 337 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
338 if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&r2<1){
bf7b8731 339 fh2TriggerCount->Fill(it,kTriggeredVertexIn);
340 }
341 }
342 }
343 if(esd){
344 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
345 // Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
216601f0 346 Bool_t cand = fInputHandler->IsEventSelected();
347 if(cand){
348 fh2ESDTriggerCount->Fill(it,kSelectedALICE);
216601f0 349 }
cc0649e4 350 if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
351
bf7b8731 352 if(vtxESD->GetNContributors()>0){
353 if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
354 Float_t zvtx = vtxESD->GetZ();
355 Float_t yvtx = vtxESD->GetY();
356 Float_t xvtx = vtxESD->GetX();
a923bd34 357 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
bf7b8731 358 fh2ESDTriggerVtx->Fill(it,zvtx);
5cb0f438 359 Bool_t vtxIn = TMath::Abs(zvtx)<fZVtxCut&&r2<1;
360
361 if(cand&&vtxIn){
362 fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
363 fh2ESDTriggerCount->Fill(it,kSelected);
364 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
365 }
366 if(vtxIn&&esdTrig){
bf7b8731 367 fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
368 // here we select based on ESD info...
bf7b8731 369 }
370 }
bf7b8731 371 }
bf7b8731 372 }
373
374
375
376 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
377
378
379 Double_t ptHard = 0;
380 Double_t nTrials = 1; // Trials for MC trigger
381
382 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
383 AliMCEvent* mcEvent = MCEvent();
384 // AliStack *pStack = 0;
385 if(mcEvent){
386 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
387 if(pythiaGenHeader){
388 nTrials = pythiaGenHeader->Trials();
389 ptHard = pythiaGenHeader->GetPtHard();
390 int iProcessType = pythiaGenHeader->ProcessType();
391 // 11 f+f -> f+f
392 // 12 f+barf -> f+barf
393 // 13 f+barf -> g+g
394 // 28 f+g -> f+g
395 // 53 g+g -> f+barf
396 // 68 g+g -> g+g
397 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
398 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
399 fh1PtHard->Fill(ptHard);
400 fh1PtHardTrials->Fill(ptHard,nTrials);
401
402 }// if pythia gen header
403 }
404
405 // trigger selection
406
407
408 PostData(1, fHistList);
409}
410
411
412void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
413{
414 // Terminate analysis
415 //
bf7b8731 416}