Moving AliTRDCalPID to libTRDbase.pkg since it is used by AliTRDcalibDB (loadlibs.C)
[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 }
cf049e94 303 Float_t r2 = xvtx *xvtx + yvtx *yvtx;
304 if(TMath::Abs(zvtx)<fZVtxCut&&r2<1)fh2ESDTriggerRun->Fill(run,iTrig+3);
b5a3f310 305 }
306 else{
307 fh2ESDTriggerRun->Fill(run,0);
308 }
309
310
311 }
bf7b8731 312
313 // loop over all possible trigger and
314 for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
315 Bool_t esdTrig = kFALSE;
316 if(esd){
317 esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
318 if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
319 }
320 Bool_t aodTrig = kFALSE;
321 if(aod){
322 aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
323 if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
324 }
325
326 // Check wether we have also an SPD vertex
327
328 if(aod){
329 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
330 // Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
331
332 if(vtxAOD->GetNContributors()>0){
333 if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
334 Float_t zvtx = vtxAOD->GetZ();
335 Float_t yvtx = vtxAOD->GetY();
336 Float_t xvtx = vtxAOD->GetX();
337 fh2TriggerVtx->Fill(it,zvtx);
cf049e94 338 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
339 if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&r2<1){
bf7b8731 340 fh2TriggerCount->Fill(it,kTriggeredVertexIn);
341 }
342 }
343 }
344 if(esd){
345 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
346 // Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
216601f0 347 Bool_t cand = fInputHandler->IsEventSelected();
348 if(cand){
349 fh2ESDTriggerCount->Fill(it,kSelectedALICE);
216601f0 350 }
cc0649e4 351 if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
352
bf7b8731 353 if(vtxESD->GetNContributors()>0){
354 if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
355 Float_t zvtx = vtxESD->GetZ();
356 Float_t yvtx = vtxESD->GetY();
357 Float_t xvtx = vtxESD->GetX();
a923bd34 358 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
bf7b8731 359 fh2ESDTriggerVtx->Fill(it,zvtx);
a923bd34 360 if(TMath::Abs(zvtx)<fZVtxCut&&esdTrig&&r2<1){
bf7b8731 361 fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
a923bd34 362 if(cand){
363 fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
364 fh2ESDTriggerCount->Fill(it,kSelected);
365 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
366 }
bf7b8731 367 // here we select based on ESD info...
bf7b8731 368 }
369 }
bf7b8731 370 }
bf7b8731 371 }
372
373
374
375 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
376
377
378 Double_t ptHard = 0;
379 Double_t nTrials = 1; // Trials for MC trigger
380
381 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
382 AliMCEvent* mcEvent = MCEvent();
383 // AliStack *pStack = 0;
384 if(mcEvent){
385 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
386 if(pythiaGenHeader){
387 nTrials = pythiaGenHeader->Trials();
388 ptHard = pythiaGenHeader->GetPtHard();
389 int iProcessType = pythiaGenHeader->ProcessType();
390 // 11 f+f -> f+f
391 // 12 f+barf -> f+barf
392 // 13 f+barf -> g+g
393 // 28 f+g -> f+g
394 // 53 g+g -> f+barf
395 // 68 g+g -> g+g
396 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
397 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
398 fh1PtHard->Fill(ptHard);
399 fh1PtHardTrials->Fill(ptHard,nTrials);
400
401 }// if pythia gen header
402 }
403
404 // trigger selection
405
406
407 PostData(1, fHistList);
408}
409
410
411void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
412{
413 // Terminate analysis
414 //
bf7b8731 415}