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