- moving macro to programs folder as it is more a tool then an example
[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>
c0997643 38#include "TDatabasePDG.h"
bf7b8731 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),
c0997643 74 fVtxXMean(0),
75 fVtxYMean(0),
76 fVtxZMean(0),
77 fVtxRCut(1.),
78 fVtxZCut(8.),
bf7b8731 79 fh1Xsec(0x0),
80 fh1Trials(0x0),
81 fh1PtHard(0x0),
82 fh1PtHardTrials(0x0),
83 fh2TriggerCount(0x0),
84 fh2ESDTriggerCount(0x0),
85 fh2TriggerVtx(0x0),
86 fh2ESDTriggerVtx(0x0),
b5a3f310 87 fh2ESDTriggerRun(0x0),
88 fh2VtxXY(0x0),
bf7b8731 89 fHistList(0x0)
90{
b5a3f310 91 fRunRange[0] = fRunRange[1] = 0;
bf7b8731 92}
93
94AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
95 AliAnalysisTaskSE(name),
96 fUseAODInput(kFALSE),
cc0649e4 97 fUsePhysicsSelection(kFALSE),
98 fRealData(kFALSE),
bf7b8731 99 fAvgTrials(1),
c0997643 100 fVtxXMean(0),
101 fVtxYMean(0),
102 fVtxZMean(0),
103 fVtxRCut(1.),
104 fVtxZCut(8.),
bf7b8731 105 fh1Xsec(0x0),
106 fh1Trials(0x0),
107 fh1PtHard(0x0),
108 fh1PtHardTrials(0x0),
109 fh2TriggerCount(0x0),
110 fh2ESDTriggerCount(0x0),
111 fh2TriggerVtx(0x0),
112 fh2ESDTriggerVtx(0x0),
b5a3f310 113 fh2ESDTriggerRun(0x0),
114 fh2VtxXY(0x0),
bf7b8731 115 fHistList(0x0)
116{
b5a3f310 117 fRunRange[0] = fRunRange[1] = 0;
bf7b8731 118 DefineOutput(1,TList::Class());
119}
120
121
122
123Bool_t AliAnalysisTaskJetServices::Notify()
124{
125 //
126 // Implemented Notify() to read the cross sections
127 // and number of trials from pyxsec.root
128 //
129
130 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
131 Float_t xsection = 0;
132 Float_t ftrials = 1;
133
134 fAvgTrials = 1;
135 if(tree){
136 TFile *curfile = tree->GetCurrentFile();
137 if (!curfile) {
138 Error("Notify","No current file");
139 return kFALSE;
140 }
141 if(!fh1Xsec||!fh1Trials){
142 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
143 return kFALSE;
144 }
145 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
146 fh1Xsec->Fill("<#sigma>",xsection);
147 // construct a poor man average trials
148 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
149 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
150 }
151 return kTRUE;
152}
153
154void AliAnalysisTaskJetServices::UserCreateOutputObjects()
155{
156
157 //
158 // Create the output container
159 //
160
161
bf7b8731 162 if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
163
164 OpenFile(1);
165 if(!fHistList)fHistList = new TList();
166
167 Bool_t oldStatus = TH1::AddDirectoryStatus();
168 TH1::AddDirectory(kFALSE);
bf7b8731 169 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
170 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
171 fHistList->Add(fh1Xsec);
172
173 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
174 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
175 fHistList->Add(fh1Trials);
176
177 const Int_t nBinPt = 100;
178 Double_t binLimitsPt[nBinPt+1];
179 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
180 if(iPt == 0){
181 binLimitsPt[iPt] = 0.0;
182 }
183 else {// 1.0
184 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
185 }
186 }
187
188 fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
189 fHistList->Add(fh2TriggerCount);
190
191 fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
192 fHistList->Add(fh2ESDTriggerCount);
193
194 fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
195 fHistList->Add(fh2TriggerVtx);
196
197 fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
198 fHistList->Add(fh2ESDTriggerVtx);
199
200
201 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
202 fHistList->Add(fh1PtHard);
203 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
204 fHistList->Add(fh1PtHardTrials);
205
b5a3f310 206 // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
207 // 3 triggers BB BE/EB EE
208
209 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);
210 fHistList->Add(fh2ESDTriggerRun);
211
212 fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
213 fHistList->Add(fh2VtxXY);
bf7b8731 214 // =========== Switch on Sumw2 for all histos ===========
215 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
216 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
217 if (h1){
218 h1->Sumw2();
219 continue;
220 }
221 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
222 if(hn)hn->Sumw2();
223 }
224
225
226 TH1::AddDirectory(oldStatus);
227}
228
229void AliAnalysisTaskJetServices::Init()
230{
231 //
232 // Initialization
233 //
bf7b8731 234 if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
235
236}
237
238void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
239{
240
241 //
242 // Execute analysis for current event
243 //
244
245 AliAODEvent *aod = 0;
246 AliESDEvent *esd = 0;
247
248 AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
249
250 if(fUseAODInput){
251 aod = dynamic_cast<AliAODEvent*>(InputEvent());
252 if(!aod){
253 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
254 return;
255 }
256 // fethc the header
257 }
258 else{
259 // assume that the AOD is in the general output...
260 aod = AODEvent();
261 if(!aod){
262 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
263 return;
264 }
265 esd = dynamic_cast<AliESDEvent*>(InputEvent());
266 }
267
b5a3f310 268 if(esd){
269 Float_t run = (Float_t)esd->GetRunNumber();
270 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
271 Float_t zvtx = -999;
272 Float_t xvtx = -999;
273 Float_t yvtx = -999;
274
3568c3d6 275 if(vtxESD->GetNContributors()>2){
b5a3f310 276 zvtx = vtxESD->GetZ();
277 yvtx = vtxESD->GetY();
278 xvtx = vtxESD->GetX();
279 }
280
281 Int_t iTrig = -1;
282 if(esd->GetFiredTriggerClasses().Contains("CINT1B")
283 ||esd->GetFiredTriggerClasses().Contains("CSMBB")
284 ||esd->GetFiredTriggerClasses().Contains("MB1")
285 ||esd->GetFiredTriggerClasses().Contains("CINT6B")){
286 iTrig = 0;
287 }
288 else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
289 ||esd->GetFiredTriggerClasses().Contains("CSMBA")
290 ||esd->GetFiredTriggerClasses().Contains("CINT6A")
291 ||esd->GetFiredTriggerClasses().Contains("CINT1C")
292 ||esd->GetFiredTriggerClasses().Contains("CSMBC")
293 ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
294 // empty bunch or bunch empty
295 iTrig = 1;
296 }
297 if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
298 ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
299 iTrig = 2;
300 }
301
302
303 if(iTrig>=0){
304 iTrig *= 3;
305 fh2ESDTriggerRun->Fill(run,iTrig+1);
3568c3d6 306 if(vtxESD->GetNContributors()>2){
b5a3f310 307 fh2ESDTriggerRun->Fill(run,iTrig+2);
308 fh2VtxXY->Fill(xvtx,yvtx);
309 }
c0997643 310 xvtx -= fVtxXMean;
311 yvtx -= fVtxYMean;
312 zvtx -= fVtxZMean;
cf049e94 313 Float_t r2 = xvtx *xvtx + yvtx *yvtx;
c0997643 314 if(TMath::Abs(zvtx)<fVtxZCut&&r2<fVtxRCut)fh2ESDTriggerRun->Fill(run,iTrig+3);
b5a3f310 315 }
316 else{
317 fh2ESDTriggerRun->Fill(run,0);
318 }
319
320
321 }
bf7b8731 322
bf7b8731 323
c0997643 324 // Apply additional constraints
325 Bool_t esdEventSelected = IsEventSelectedESD(esd);
326 Bool_t aodEventSelected = IsEventSelectedAOD(aod);
327
328
329 if(esd&&aod&&fDebug){
330 if(esdEventSelected&&!aodEventSelected){
331 Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
332 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
bf7b8731 333 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
c0997643 334 Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
335 vtxESD->Print();
336 Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
337 vtxAOD->Print();
bf7b8731 338 }
c0997643 339 }
340
341 // loop over all possible triggers for esd
342
343 if(esd){
344 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
345 // Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
346 TString vtxName(vtxESD->GetName());
347 for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
348 Bool_t esdTrig = kFALSE;
349 esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
350 if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
216601f0 351 Bool_t cand = fInputHandler->IsEventSelected();
352 if(cand){
353 fh2ESDTriggerCount->Fill(it,kSelectedALICE);
216601f0 354 }
cc0649e4 355 if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
3568c3d6 356 if(vtxESD->GetNContributors()>2&&!vtxName.Contains("TPCVertex")){
c0997643 357 if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredVertex);
bf7b8731 358 Float_t zvtx = vtxESD->GetZ();
c0997643 359 if(esdEventSelected&&esdTrig){
bf7b8731 360 fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
cc89bb69 361 fh2ESDTriggerVtx->Fill(it,zvtx);
bf7b8731 362 }
363 }
c0997643 364 if(cand&&esdEventSelected){
365 fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
366 fh2ESDTriggerCount->Fill(it,kSelected);
367 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
368 }
bf7b8731 369 }
bf7b8731 370 }
371
c0997643 372 if(aod){
373 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
374 // Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
375 TString vtxTitle(vtxAOD->GetTitle());
376 for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
377 Bool_t aodTrig = kFALSE;
378 aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
379 if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
3568c3d6 380 if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
c0997643 381 if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex);
382 Float_t zvtx = vtxAOD->GetZ();
c0997643 383 if(aodTrig&&aodEventSelected){
cc89bb69 384 fh2TriggerVtx->Fill(it,zvtx);
c0997643 385 fh2TriggerCount->Fill(it,kTriggeredVertexIn);
386 }
387 }
388 }
389 }
bf7b8731 390
391 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
392
393
394 Double_t ptHard = 0;
395 Double_t nTrials = 1; // Trials for MC trigger
396
397 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
398 AliMCEvent* mcEvent = MCEvent();
399 // AliStack *pStack = 0;
400 if(mcEvent){
401 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
402 if(pythiaGenHeader){
403 nTrials = pythiaGenHeader->Trials();
404 ptHard = pythiaGenHeader->GetPtHard();
405 int iProcessType = pythiaGenHeader->ProcessType();
406 // 11 f+f -> f+f
407 // 12 f+barf -> f+barf
408 // 13 f+barf -> g+g
409 // 28 f+g -> f+g
410 // 53 g+g -> f+barf
411 // 68 g+g -> g+g
412 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
413 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
414 fh1PtHard->Fill(ptHard);
415 fh1PtHardTrials->Fill(ptHard,nTrials);
416
417 }// if pythia gen header
418 }
419
420 // trigger selection
421
422
423 PostData(1, fHistList);
424}
425
c0997643 426Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){
427 if(!esd)return kFALSE;
428 const AliESDVertex *vtx = esd->GetPrimaryVertex();
429 // We can check the type of the vertex by its name and title
430 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
431 // vtx name title
432 // Tracks PrimaryVertex VertexerTracksNoConstraint
433 // TPC TPCVertex VertexerTracksNoConstraint
434 // SPD SPDVertex vertexer: 3D
435 // SPD SPDVertex vertexer: Z
436
437
438
3568c3d6 439 if(vtx->GetNContributors()<3)return kFALSE;
c0997643 440
441 // do not want tpc only primary vertex
442 TString vtxName(vtx->GetName());
443 if(vtxName.Contains("TPCVertex"))return kFALSE;
444 Float_t zvtx = vtx->GetZ();
445 Float_t yvtx = vtx->GetY();
446 Float_t xvtx = vtx->GetX();
447
448 // here we may fill the histos for run dependence....
449
450 xvtx -= fVtxXMean;
451 yvtx -= fVtxYMean;
452 zvtx -= fVtxZMean;
453
454 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
455
456 Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
457 return eventSel;
458}
459
460Bool_t AliAnalysisTaskJetServices::IsEventSelectedAOD(AliAODEvent* aod){
461 if(!aod)return kFALSE;
462 const AliAODVertex *vtx = aod->GetPrimaryVertex();
463 // We can check the type of the vertex by its name and title
464 // if vertexer failed title is empty (default c'tor) and ncontributors is 0
465 // vtx name title
466 // Tracks PrimaryVertex VertexerTracksNoConstraint
467 // TPC TPCVertex VertexerTracksNoConstraint
468 // SPD SPDVertex vertexer: 3D
469 // SPD SPDVertex vertexer: Z
470
471
472
3568c3d6 473 if(vtx->GetNContributors()<3)return kFALSE;
c0997643 474
475 // do not want tpc only primary vertex
476 TString vtxName(vtx->GetName());
477 if(vtxName.Contains("TPCVertex"))return kFALSE;
478
479 Float_t zvtx = vtx->GetZ();
480 Float_t yvtx = vtx->GetY();
481 Float_t xvtx = vtx->GetX();
482
483 // here we may fill the histos for run dependence....
484
485 xvtx -= fVtxXMean;
486 yvtx -= fVtxYMean;
487 zvtx -= fVtxZMean;
488
489 Float_t r2 = yvtx*yvtx+xvtx*xvtx;
490 Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
491 return eventSel;
492}
493
494
bf7b8731 495
496void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
497{
498 // Terminate analysis
499 //
bf7b8731 500}