Adding AliPhysicsSelectionTask, some clean up wiht output file managing on the grid...
[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),
70 fAvgTrials(1),
fe669ac6 71 fZVtxCut(8.),
72 fRealData(kFALSE),
6d5f6d70 73 fPhysicsSelection(0),
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),
92 fAvgTrials(1),
fe669ac6 93 fZVtxCut(8.),
94 fRealData(kFALSE),
6d5f6d70 95 fPhysicsSelection(0),
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);
160
fe669ac6 161 if(!fPhysicsSelection)
162 fPhysicsSelection = new AliPhysicsSelection();
163 fPhysicsSelection->SetName("AliPhysicsSelection_outputlist"); // to prevent conflict with object that is automatically streamed back
164 //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
165 fHistList->Add(fPhysicsSelection);
166
bf7b8731 167 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
168 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
169 fHistList->Add(fh1Xsec);
170
171 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
172 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
173 fHistList->Add(fh1Trials);
174
175 const Int_t nBinPt = 100;
176 Double_t binLimitsPt[nBinPt+1];
177 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
178 if(iPt == 0){
179 binLimitsPt[iPt] = 0.0;
180 }
181 else {// 1.0
182 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
183 }
184 }
185
186 fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
187 fHistList->Add(fh2TriggerCount);
188
189 fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
190 fHistList->Add(fh2ESDTriggerCount);
191
192 fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
193 fHistList->Add(fh2TriggerVtx);
194
195 fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
196 fHistList->Add(fh2ESDTriggerVtx);
197
198
199 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
200 fHistList->Add(fh1PtHard);
201 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
202 fHistList->Add(fh1PtHardTrials);
203
b5a3f310 204 // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
205 // 3 triggers BB BE/EB EE
206
207 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);
208 fHistList->Add(fh2ESDTriggerRun);
209
210 fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
211 fHistList->Add(fh2VtxXY);
bf7b8731 212 // =========== Switch on Sumw2 for all histos ===========
213 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
214 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
215 if (h1){
216 h1->Sumw2();
217 continue;
218 }
219 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
220 if(hn)hn->Sumw2();
221 }
222
223
224 TH1::AddDirectory(oldStatus);
225}
226
227void AliAnalysisTaskJetServices::Init()
228{
229 //
230 // Initialization
231 //
232
233 Printf(">>> AnalysisTaskJetServices::Init() debug level %d\n",fDebug);
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
275 if(vtxESD->GetNContributors()>0){
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);
306 if(vtxESD->GetNContributors()>0){
307 fh2ESDTriggerRun->Fill(run,iTrig+2);
308 fh2VtxXY->Fill(xvtx,yvtx);
309 }
310 if(TMath::Abs(zvtx)<fZVtxCut&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5)fh2ESDTriggerRun->Fill(run,iTrig+3);
311 }
312 else{
313 fh2ESDTriggerRun->Fill(run,0);
314 }
315
316
317 }
bf7b8731 318
319 // loop over all possible trigger and
320 for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
321 Bool_t esdTrig = kFALSE;
322 if(esd){
323 esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
324 if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
325 }
326 Bool_t aodTrig = kFALSE;
327 if(aod){
328 aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
329 if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
330 }
331
332 // Check wether we have also an SPD vertex
333
334 if(aod){
335 const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
336 // Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
337
338 if(vtxAOD->GetNContributors()>0){
339 if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
340 Float_t zvtx = vtxAOD->GetZ();
341 Float_t yvtx = vtxAOD->GetY();
342 Float_t xvtx = vtxAOD->GetX();
343 fh2TriggerVtx->Fill(it,zvtx);
344 if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
345 fh2TriggerCount->Fill(it,kTriggeredVertexIn);
346 }
347 }
348 }
349 if(esd){
350 const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
351 // Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
6d5f6d70 352 Bool_t cand = true;
353 if(fRealData)cand = fPhysicsSelection->IsCollisionCandidate(esd);
354 if(cand) fh2ESDTriggerCount->Fill(it,kSelectedALICE);
bf7b8731 355 if(vtxESD->GetNContributors()>0){
356 if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
357 Float_t zvtx = vtxESD->GetZ();
358 Float_t yvtx = vtxESD->GetY();
359 Float_t xvtx = vtxESD->GetX();
360 fh2ESDTriggerVtx->Fill(it,zvtx);
361 if(TMath::Abs(zvtx)<fZVtxCut&&esdTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
362 fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
363 // here we select based on ESD info...
6d5f6d70 364 if(esdTrig&&(it==AliAnalysisHelperJetTasks::kMB1)){
fe669ac6 365 fh2ESDTriggerCount->Fill(it,kSelected);
366 AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
367 }
bf7b8731 368 }
fe669ac6 369
370
bf7b8731 371 }
372
373 }
374
375 }
376
377
378
379 if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
380
381
382 Double_t ptHard = 0;
383 Double_t nTrials = 1; // Trials for MC trigger
384
385 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
386 AliMCEvent* mcEvent = MCEvent();
387 // AliStack *pStack = 0;
388 if(mcEvent){
389 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
390 if(pythiaGenHeader){
391 nTrials = pythiaGenHeader->Trials();
392 ptHard = pythiaGenHeader->GetPtHard();
393 int iProcessType = pythiaGenHeader->ProcessType();
394 // 11 f+f -> f+f
395 // 12 f+barf -> f+barf
396 // 13 f+barf -> g+g
397 // 28 f+g -> f+g
398 // 53 g+g -> f+barf
399 // 68 g+g -> g+g
400 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
401 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
402 fh1PtHard->Fill(ptHard);
403 fh1PtHardTrials->Fill(ptHard,nTrials);
404
405 }// if pythia gen header
406 }
407
408 // trigger selection
409
410
411 PostData(1, fHistList);
412}
413
414
415void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
416{
417 // Terminate analysis
418 //
fe669ac6 419
420 TDirectory* owd = gDirectory;
421
bf7b8731 422 if (fDebug > 1) printf("AnalysisJetServices: Terminate() \n");
fe669ac6 423
424 fHistList = dynamic_cast<TList*> (GetOutputData(1));
425 if (!fHistList)
e1a13e49 426 Printf("ERROR: fHistList not available");
fe669ac6 427
428
429
430 AliAnalysisDataContainer *cont = GetOutputSlot(1)->GetContainer();
431 TString filename = cont->GetFileName();
432 TFile *f = NULL;
433 // Check first if the file is already opened
434 f = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
435 if (f) {
436 // Cd to file
437 f->cd();
438 // Check for a folder request
439 TString dir = cont->GetFolderName();
440 if (!dir.IsNull()) {
441 if (!f->GetDirectory(dir)) f->mkdir(dir);
442 f->cd(dir);
443 }
444 }
445
446 if (fHistList)
447 {
448 fPhysicsSelection = dynamic_cast<AliPhysicsSelection*> (fHistList->FindObject("AliPhysicsSelection_outputlist"));
449 }
450
451 if (fPhysicsSelection)
452 {
453 fPhysicsSelection->SaveHistograms("physics_selection");
454 fPhysicsSelection->Print();
455 }
456
457 owd->cd();
458
bf7b8731 459}