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