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