]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added fit macro from M. Putis
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum.cxx
CommitLineData
77bbf113 1// **************************************
2// Task used for the correction of determiantion of reconstructed jet spectra
3// Compares input (gen) and output (rec) jets
4// *******************************************
5
6
f3050824 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
f3050824 22
23#include <TROOT.h>
1d676c86 24#include <TRandom.h>
f3050824 25#include <TSystem.h>
26#include <TInterpreter.h>
27#include <TChain.h>
28#include <TFile.h>
60c2f4e5 29#include <TKey.h>
f3050824 30#include <TH1F.h>
31#include <TH2F.h>
32#include <TH3F.h>
4dbfdecc 33#include <TProfile.h>
f3050824 34#include <TList.h>
35#include <TLorentzVector.h>
36#include <TClonesArray.h>
37#include "TDatabasePDG.h"
38
39#include "AliAnalysisTaskJetSpectrum.h"
40#include "AliAnalysisManager.h"
41#include "AliJetFinder.h"
42#include "AliJetHeader.h"
43#include "AliJetReader.h"
44#include "AliJetReaderHeader.h"
45#include "AliUA1JetHeaderV1.h"
f3050824 46#include "AliESDEvent.h"
47#include "AliAODEvent.h"
48#include "AliAODHandler.h"
49#include "AliAODTrack.h"
50#include "AliAODJet.h"
51#include "AliMCEventHandler.h"
52#include "AliMCEvent.h"
53#include "AliStack.h"
54#include "AliGenPythiaEventHeader.h"
55#include "AliJetKineReaderHeader.h"
56#include "AliGenCocktailEventHeader.h"
57#include "AliInputEventHandler.h"
58
db6bcb0e 59
f3050824 60#include "AliAnalysisHelperJetTasks.h"
61
62ClassImp(AliAnalysisTaskJetSpectrum)
63
df65bddb 64const Float_t AliAnalysisTaskJetSpectrum::fgkJetNpartCut[AliAnalysisTaskJetSpectrum::kMaxCorrelation] = {5,10,1E+09};
65
f3050824 66AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
67 fJetHeaderRec(0x0),
68 fJetHeaderGen(0x0),
69 fAOD(0x0),
70 fBranchRec("jets"),
f3050824 71 fBranchGen(""),
f3050824 72 fUseAODInput(kFALSE),
73 fUseExternalWeightOnly(kFALSE),
74 fLimitGenJetEta(kFALSE),
75 fAnalysisType(0),
76 fExternalWeight(1),
d40dc1ba 77 fRecEtaWindow(0.5),
4dbfdecc 78 fh1Xsec(0x0),
79 fh1Trials(0x0),
f3050824 80 fh1PtHard(0x0),
77bbf113 81 fh1PtHardNoW(0x0),
82 fh1PtHardTrials(0x0),
f3050824 83 fh1NGenJets(0x0),
84 fh1NRecJets(0x0),
db6bcb0e 85 fHistList(0x0) ,
df65bddb 86 ////////////////
87 fh1JetMultiplicity(0x0) ,
88 fh2ERecZRec(0x0) ,
89 fh2EGenZGen(0x0) ,
90 fh2Efficiency(0x0) ,
91 fh3EGenERecN(0x0)
92 ////////////////
f3050824 93{
94 // Default constructor
95 for(int ij = 0;ij<kMaxJets;++ij){
df65bddb 96 fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
d40dc1ba 97 fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] = fh2PtGenDeltaPhi[ij] = fh2PtGenDeltaEta[ij] = 0;
77bbf113 98 fh3PtRecGenHard[ij] = fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
df65bddb 99 }
100 for(int ic = 0;ic < kMaxCorrelation;ic++){
101 fhnCorrelation[ic] = 0;
102 }
f3050824 103
104}
105
106AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
107 AliAnalysisTaskSE(name),
108 fJetHeaderRec(0x0),
109 fJetHeaderGen(0x0),
110 fAOD(0x0),
111 fBranchRec("jets"),
f3050824 112 fBranchGen(""),
f3050824 113 fUseAODInput(kFALSE),
114 fUseExternalWeightOnly(kFALSE),
115 fLimitGenJetEta(kFALSE),
116 fAnalysisType(0),
117 fExternalWeight(1),
d40dc1ba 118 fRecEtaWindow(0.5),
4dbfdecc 119 fh1Xsec(0x0),
120 fh1Trials(0x0),
f3050824 121 fh1PtHard(0x0),
77bbf113 122 fh1PtHardNoW(0x0),
123 fh1PtHardTrials(0x0),
f3050824 124 fh1NGenJets(0x0),
125 fh1NRecJets(0x0),
df65bddb 126 fHistList(0x0) ,
127 ////////////////
128 fh1JetMultiplicity(0x0) ,
129 fh2ERecZRec(0x0) ,
130 fh2EGenZGen(0x0) ,
131 fh2Efficiency(0x0) ,
132 fh3EGenERecN(0x0)
133 ////////////////
f3050824 134{
135 // Default constructor
136 for(int ij = 0;ij<kMaxJets;++ij){
137 fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
d40dc1ba 138 fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtGenDeltaPhi[ij] = fh2PtGenDeltaEta[ij] = fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] =0;
f3050824 139
77bbf113 140 fh3PtRecGenHard[ij] = fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
f3050824 141 }
142
df65bddb 143 for(int ic = 0;ic < kMaxCorrelation;ic++){
144 fhnCorrelation[ic] = 0;
145 }
146
f3050824 147 DefineOutput(1, TList::Class());
148}
149
150
151
4dbfdecc 152Bool_t AliAnalysisTaskJetSpectrum::Notify()
153{
77bbf113 154 //
155 // Implemented Notify() to read the cross sections
156 // and number of trials from pyxsec.root
157 //
4dbfdecc 158 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
159 Double_t xsection = 0;
160 UInt_t ntrials = 0;
60c2f4e5 161 Float_t ftrials = 0;
4dbfdecc 162 if(tree){
163 TFile *curfile = tree->GetCurrentFile();
164 if (!curfile) {
165 Error("Notify","No current file");
166 return kFALSE;
167 }
168 if(!fh1Xsec||!fh1Trials){
169 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
170 return kFALSE;
171 }
172
173 TString fileName(curfile->GetName());
174 if(fileName.Contains("AliESDs.root")){
60c2f4e5 175 fileName.ReplaceAll("AliESDs.root", "");
5385a44b 176 }
177 else if(fileName.Contains("AliAOD.root")){
178 fileName.ReplaceAll("AliAOD.root", "");
4dbfdecc 179 }
60c2f4e5 180 else if(fileName.Contains("AliAODs.root")){
181 fileName.ReplaceAll("AliAODs.root", "");
4dbfdecc 182 }
183 else if(fileName.Contains("galice.root")){
184 // for running with galice and kinematics alone...
60c2f4e5 185 fileName.ReplaceAll("galice.root", "");
4dbfdecc 186 }
60c2f4e5 187 TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"));
4dbfdecc 188 if(!fxsec){
60c2f4e5 189 if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec.root"));
190 // next trial fetch the histgram file
191 fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
192 if(!fxsec){
193 // not a severe condition
194 if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
195 return kTRUE;
196 }
197 else{
198 // find the tlist we want to be independtent of the name so use the Tkey
199 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
200 if(!key){
201 if(fDebug>0)Printf("%s:%d key not found in the file",(char*)__FILE__,__LINE__);
202 return kTRUE;
203 }
204 TList *list = dynamic_cast<TList*>(key->ReadObj());
205 if(!list){
206 if(fDebug>0)Printf("%s:%d key is not a tlist",(char*)__FILE__,__LINE__);
207 return kTRUE;
208 }
209 xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
210 ftrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
211 }
4dbfdecc 212 }
60c2f4e5 213 else{
214 TTree *xtree = (TTree*)fxsec->Get("Xsection");
215 if(!xtree){
216 Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
217 }
e855f5c5 218 else{
219 xtree->SetBranchAddress("xsection",&xsection);
220 xtree->SetBranchAddress("ntrials",&ntrials);
221 ftrials = ntrials;
222 xtree->GetEntry(0);
223 }
4dbfdecc 224 }
4dbfdecc 225 fh1Xsec->Fill("<#sigma>",xsection);
60c2f4e5 226 fh1Trials->Fill("#sum{ntrials}",ftrials);
4dbfdecc 227 }
228
229 return kTRUE;
230}
f3050824 231
232void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
233{
234
235 //
236 // Create the output container
237 //
238
239
240 // Connect the AOD
241
242 if(fUseAODInput){
243 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
244 if(!fAOD){
245 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
246 return;
247 }
248 // fethc the header
249 fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
250 if(!fJetHeaderRec){
251 Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
252 }
253 }
254 else{
255 // assume that the AOD is in the general output...
256 fAOD = AODEvent();
257 if(!fAOD){
258 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
259 return;
260 }
261 fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
262 if(!fJetHeaderRec){
263 Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
264 }
265 else{
266 if(fDebug>10)fJetHeaderRec->Dump();
267 }
268 }
269
270
271
272 if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
273
274 OpenFile(1);
275 if(!fHistList)fHistList = new TList();
276
277 Bool_t oldStatus = TH1::AddDirectoryStatus();
278 TH1::AddDirectory(kFALSE);
279
280 //
281 // Histogram
282
283 const Int_t nBinPt = 100;
284 Double_t binLimitsPt[nBinPt+1];
285 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
286 if(iPt == 0){
287 binLimitsPt[iPt] = 0.0;
288 }
289 else {// 1.0
de6f8090 290 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
f3050824 291 }
292 }
d40dc1ba 293
294 const Int_t nBinEta = 26;
295 Double_t binLimitsEta[nBinEta+1] = {
296 -1.6,-1.4,-1.2,-1.0,
297 -0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
298 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
299 1.0, 1.2, 1.4, 1.6
300 };
301
f3050824 302
22ceb537 303 const Int_t nBinPhi = 30;
f3050824 304 Double_t binLimitsPhi[nBinPhi+1];
305 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
306 if(iPhi==0){
307 binLimitsPhi[iPhi] = 0;
308 }
309 else{
310 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
311 }
312 }
313
314 const Int_t nBinFrag = 25;
315
316
4dbfdecc 317 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
318 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
319
320 fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
321 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
322
f3050824 323 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
324
77bbf113 325 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
f3050824 326
77bbf113 327 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
f3050824 328
329 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
330
331 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
332
333
334 for(int ij = 0;ij<kMaxJets;++ij){
335 fh1E[ij] = new TH1F(Form("fh1E_j%d",ij),"Jet Energy;E_{jet} (GeV);N",nBinPt,binLimitsPt);
336 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
337 fh1PtRecOut[ij] = new TH1F(Form("fh1PtRecOut_j%d",ij),"rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
338 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
339 fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
340
341
d40dc1ba 342
f3050824 343 fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
344 nBinPt,binLimitsPt,nBinPt,binLimitsPt);
345
de6f8090 346 fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};#phi_{gen}",
22ceb537 347 nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
f3050824 348
de6f8090 349 fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};#eta_{gen}",
22ceb537 350 nBinEta,binLimitsEta,nBinEta,binLimitsEta);
f3050824 351
22ceb537 352
de6f8090 353 fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaPhi_j%d",ij),"delta phi vs. P_{T,gen};p_{T,gen} (GeV/c);#phi_{gen}-#phi_{rec}",
77bbf113 354 nBinPt,binLimitsPt,100,-1.0,1.0);
de6f8090 355 fh2PtGenDeltaEta[ij] = new TH2F(Form("fh2PtGenDeltaEta_j%d",ij),"delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-#eta_{rec}",
77bbf113 356 nBinPt,binLimitsPt,100,-1.0,1.0);
22ceb537 357
358
de6f8090 359 fh2PtRecDeltaR[ij] = new TH2F(Form("fh2PtRecDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,rec,j};#Delta R",
d40dc1ba 360 nBinPt,binLimitsPt,60,0,6.0);
de6f8090 361 fh2PtGenDeltaR[ij] = new TH2F(Form("fh2PtGenDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,gen,j};#Delta R",
d40dc1ba 362 nBinPt,binLimitsPt,60,0,6.0);
f3050824 363
364
365
366 fh3PtRecGenHard[ij] = new TH3F(Form("fh3PtRecGenHard_j%d",ij), "Pt hard vs. pt gen vs. pt rec;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
367
368
369
77bbf113 370 fh3PtRecGenHardNoW[ij] = new TH3F(Form("fh3PtRecGenHardNoW_j%d",ij), "Pt hard vs. pt gen vs. pt rec no weight;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
f3050824 371
372
373 fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
374 nBinFrag,0.,1.,nBinPt,binLimitsPt);
375
376 fh2FragLn[ij] = new TH2F(Form("fh2FragLn_j%d",ij),"Jet Fragmentation Ln;#xi=ln(E_{jet}/E_{i});E_{jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
377 nBinFrag,0.,10.,nBinPt,binLimitsPt);
378
379 fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
380 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
381
382
383
77bbf113 384 fh3RecEtaPhiPtNoGen[ij] = new TH3F(Form("fh3RecEtaPhiPtNoGen_j%d",ij),"No generated for found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
f3050824 385 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
386
387
77bbf113 388 fh3GenEtaPhiPtNoFound[ij] = new TH3F(Form("fh3GenEtaPhiPtNoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
f3050824 389 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
390
391
392
de6f8090 393 fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)",
f3050824 394 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
395
396 }
df65bddb 397
1d676c86 398
399 // tmp histos do not add to the header
400 TH2F *hCorrPt = new TH2F("fh2PtRecPhiCorrPt","#Delta#phi correlation pt weighted",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
401 fHistList->Add(hCorrPt);
402 TH2F *hCorrRanPt = new TH2F("fh2PtRecPhiCorrPtRan","#Delta#phi Random correlation pt weighted",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
403 fHistList->Add(hCorrRanPt);
404
405 TH2F *hCorr = new TH2F("fh2PtRecPhiCorr","#Delta#phi correlation",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
406 fHistList->Add(hCorr);
407 TH2F *hCorrRan = new TH2F("fh2PtRecPhiCorrRan","#Delta#phi Random correlation",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
408 fHistList->Add(hCorrRan);
409
410
df65bddb 411 /////////////////////////////////////////////////////////////////
412 fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
f3050824 413
df65bddb 414 fh2ERecZRec = new TH2F("fh2ERecZRec", " ; E^{jet}_{rec} [GeV]; z^{lp}_{rec}", 100, 0., 250., 100, 0., 2.);
415 fh2EGenZGen = new TH2F("fh2EGenZGen", " ; E^{jet}_{gen} [GeV]; z^{lp}_{gen}", 100, 0., 250., 100, 0., 2.);
416 fh2Efficiency = new TH2F("fh2Efficiency", "ERec/EGen;E^{jet}_{gen} [GeV];E^{jet}_{rec}/E^{jet}_{gen}", 100, 0., 250., 100, 0., 1.5);
f3050824 417
22ceb537 418 fh3EGenERecN = new TH3F("fh3EGenERecN", "Efficiency vs. Jet Multiplicity", 100, 0., 250., 100, 0., 250., 51, 0., 50.);
df65bddb 419
420 // Response map
421 //arrays for bin limits
422 const Int_t nbin[4] = {100, 100, 100, 100};
77bbf113 423 Double_t vLowEdge[4] = {0.,0.,0.,0.};
424 Double_t vUpEdge[4] = {250., 250., 1., 1.};
df65bddb 425
426 for(int ic = 0;ic < kMaxCorrelation;ic++){
77bbf113 427 fhnCorrelation[ic] = new THnSparseF(Form("fhnCorrelation_%d",ic), "Response Map", 4, nbin, vLowEdge, vUpEdge);
df65bddb 428 if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
429 else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
430 }
22ceb537 431 const Int_t saveLevel = 3; // large save level more histos
f3050824 432 if(saveLevel>0){
4dbfdecc 433 fHistList->Add(fh1Xsec);
434 fHistList->Add(fh1Trials);
f3050824 435 fHistList->Add(fh1PtHard);
77bbf113 436 fHistList->Add(fh1PtHardNoW);
437 fHistList->Add(fh1PtHardTrials);
f3050824 438 fHistList->Add(fh1NGenJets);
439 fHistList->Add(fh1NRecJets);
df65bddb 440 ////////////////////////
441 fHistList->Add(fh1JetMultiplicity);
442 fHistList->Add(fh2ERecZRec);
443 fHistList->Add(fh2EGenZGen);
444 fHistList->Add(fh2Efficiency);
445 fHistList->Add(fh3EGenERecN);
446
447 for(int ic = 0;ic < kMaxCorrelation;++ic){
448 fHistList->Add(fhnCorrelation[ic]);
449 }
450 ////////////////////////
f3050824 451 for(int ij = 0;ij<kMaxJets;++ij){
452 fHistList->Add(fh1E[ij]);
453 fHistList->Add(fh1PtRecIn[ij]);
454 fHistList->Add(fh1PtRecOut[ij]);
455 fHistList->Add(fh1PtGenIn[ij]);
456 fHistList->Add(fh1PtGenOut[ij]);
457 fHistList->Add(fh2PtFGen[ij]);
22ceb537 458 fHistList->Add(fh2PhiFGen[ij]);
459 fHistList->Add(fh2EtaFGen[ij]);
460 fHistList->Add(fh2PtGenDeltaEta[ij]);
461 fHistList->Add(fh2PtGenDeltaPhi[ij]);
d40dc1ba 462 fHistList->Add(fh2PtRecDeltaR[ij]);
463 fHistList->Add(fh2PtGenDeltaR[ij]);
22ceb537 464 fHistList->Add(fh3RecEtaPhiPt[ij]);
465 fHistList->Add(fh3GenEtaPhiPt[ij]);
f3050824 466 if(saveLevel>2){
77bbf113 467 fHistList->Add(fh3RecEtaPhiPtNoGen[ij]);
468 fHistList->Add(fh3GenEtaPhiPtNoFound[ij]);
f3050824 469 }
470 }
471 }
472
473 // =========== Switch on Sumw2 for all histos ===========
474 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
475 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
476 if (h1){
477 // Printf("%s ",h1->GetName());
478 h1->Sumw2();
4dbfdecc 479 continue;
f3050824 480 }
4dbfdecc 481 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
482 if(hn)hn->Sumw2();
f3050824 483 }
484
485 TH1::AddDirectory(oldStatus);
486
487}
488
489void AliAnalysisTaskJetSpectrum::Init()
490{
491 //
492 // Initialization
493 //
494
495 Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
496 if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
497
498}
499
500void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
501{
502 //
503 // Execute analysis for current event
504 //
505
506
507
508 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
509
510
511 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
512
513 if(!aodH){
514 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
515 return;
516 }
517
f3050824 518 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
c4396fd4 519
f3050824 520 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
521 if(!aodRecJets){
522 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
523 return;
524 }
525
526 // ==== General variables needed
527
528
529 // We use statice array, not to fragment the memory
530 AliAODJet genJets[kMaxJets];
531 Int_t nGenJets = 0;
532 AliAODJet recJets[kMaxJets];
533 Int_t nRecJets = 0;
df65bddb 534 ///////////////////////////
535 Int_t nTracks = 0;
536 //////////////////////////
f3050824 537
538 Double_t eventW = 1;
539 Double_t ptHard = 0;
540 Double_t nTrials = 1; // Trials for MC trigger weigth for real data
541
542 if(fUseExternalWeightOnly){
543 eventW = fExternalWeight;
544 }
545
546
547 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
548 if((fAnalysisType&kAnaMC)==kAnaMC){
549 // this is the part we only use when we have MC information
550 AliMCEvent* mcEvent = MCEvent();
551 // AliStack *pStack = 0;
552 if(!mcEvent){
553 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
554 return;
555 }
556 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
557 if(!pythiaGenHeader){
558 return;
559 }
560
561 nTrials = pythiaGenHeader->Trials();
562 ptHard = pythiaGenHeader->GetPtHard();
de6f8090 563 int iProcessType = pythiaGenHeader->ProcessType();
564 // 11 f+f -> f+f
565 // 12 f+barf -> f+barf
566 // 13 f+barf -> g+g
567 // 28 f+g -> f+g
568 // 53 g+g -> f+barf
569 // 68 g+g -> g+g
570 /*
571 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
572 // if(iProcessType != 13 && iProcessType != 68){ // allow only glue
573 if(iProcessType != 11 && iProcessType != 12 && iProcessType != 53){ // allow only quark
574 // if(iProcessType != 28){ // allow only -> f+g
575 PostData(1, fHistList);
576 return;
577 }
578 */
579 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
f3050824 580
581 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
582
583 if(!fUseExternalWeightOnly){
584 // case were we combine more than one p_T hard bin...
585 }
586
587 // fetch the pythia generated jets only to be used here
588 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
589 AliAODJet pythiaGenJets[kMaxJets];
590 Int_t iCount = 0;
591 for(int ip = 0;ip < nPythiaGenJets;++ip){
592 if(iCount>=kMaxJets)continue;
593 Float_t p[4];
594 pythiaGenHeader->TriggerJet(ip,p);
595 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
596
597 if(fLimitGenJetEta){
598 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
599 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
600 }
601
602
603 if(fBranchGen.Length()==0){
604 // if we have MC particles and we do not read from the aod branch
605 // use the pythia jets
606 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
607 }
608 iCount++;
609 }
610 if(fBranchGen.Length()==0)nGenJets = iCount;
611
612 }// (fAnalysisType&kMC)==kMC)
613
614 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
615 fh1PtHard->Fill(ptHard,eventW);
77bbf113 616 fh1PtHardNoW->Fill(ptHard,1);
617 fh1PtHardTrials->Fill(ptHard,nTrials);
f3050824 618
619 // If we set a second branch for the input jets fetch this
620 if(fBranchGen.Length()>0){
621 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
622 if(aodGenJets){
623 Int_t iCount = 0;
624 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
625 if(iCount>=kMaxJets)continue;
626 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
627 if(!tmp)continue;
628 if(fLimitGenJetEta){
629 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
630 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
631 }
632 genJets[iCount] = *tmp;
633 iCount++;
634 }
635 nGenJets = iCount;
636 }
637 else{
638 Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
639 }
640 }
641
642 fh1NGenJets->Fill(nGenJets);
643 // We do not want to exceed the maximum jet number
644 nGenJets = TMath::Min(nGenJets,kMaxJets);
645
646 // Fetch the reconstructed jets...
647
648
649 nRecJets = aodRecJets->GetEntries();
650 fh1NRecJets->Fill(nRecJets);
651 nRecJets = TMath::Min(nRecJets,kMaxJets);
df65bddb 652 //////////////////////////////////////////
653 nTracks = fAOD->GetNumberOfTracks();
654 ///////////////////////////////////////////
f3050824 655
656 for(int ir = 0;ir < nRecJets;++ir){
657 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
658 if(!tmp)continue;
659 recJets[ir] = *tmp;
660 }
661
662
663 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
664 // Relate the jets
665 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
666 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
667
668 for(int i = 0;i<kMaxJets;++i){
669 iGenIndex[i] = iRecIndex[i] = -1;
670 }
671
672
db6bcb0e 673 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
f3050824 674 iGenIndex,iRecIndex,fDebug);
675 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
676
677 if(fDebug){
678 for(int i = 0;i<kMaxJets;++i){
679 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
680 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
681 }
682 }
683
684 // loop over reconstructed jets
685 for(int ir = 0;ir < nRecJets;++ir){
686 Double_t ptRec = recJets[ir].Pt();
687 Double_t phiRec = recJets[ir].Phi();
1d676c86 688 Double_t phiRecRan = TMath::Pi()*gRandom->Rndm(); // better take real jet axis from previous events (TPC acceptance in phi)
f3050824 689 if(phiRec<0)phiRec+=TMath::Pi()*2.;
690 Double_t etaRec = recJets[ir].Eta();
77bbf113 691 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
f3050824 692 fh1E[ir]->Fill(recJets[ir].E(),eventW);
693 fh1PtRecIn[ir]->Fill(ptRec,eventW);
694 fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
d40dc1ba 695 for(int irr = ir+1;irr<nRecJets;irr++){
696 fh2PtRecDeltaR[ir]->Fill(recJets[irr].Pt(),recJets[ir].DeltaR(&recJets[irr]));
697 }
f3050824 698 // Fill Correlation
699 Int_t ig = iGenIndex[ir];
df65bddb 700 if(ig>=0 && ig<nGenJets){
77bbf113 701 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
702 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
f3050824 703 fh1PtRecOut[ir]->Fill(ptRec,eventW);
df65bddb 704 Double_t ptGen = genJets[ig].Pt();
705 Double_t phiGen = genJets[ig].Phi();
706 if(phiGen<0)phiGen+=TMath::Pi()*2.;
707 Double_t etaGen = genJets[ig].Eta();
d40dc1ba 708
709 //
710 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
711 //
712
713 if(TMath::Abs(etaRec)<fRecEtaWindow){
714
f3050824 715 fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
22ceb537 716 fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
717 fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
77bbf113 718 fh2PtGenDeltaEta[ir]->Fill(ptGen,etaGen-etaRec,eventW);
77bbf113 719 fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
f3050824 720 fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
77bbf113 721 fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
df65bddb 722 /////////////////////////////////////////////////////
de6f8090 723
724 // Double_t eRec = recJets[ir].E();
725 // Double_t eGen = genJets[ig].E();
726 // CKB use p_T not Energy
727 // TODO recname variabeles and histos
77bbf113 728 Double_t eRec = recJets[ir].E();
729 Double_t eGen = genJets[ig].E();
d40dc1ba 730
77bbf113 731 fh2Efficiency->Fill(eGen, eRec/eGen);
77bbf113 732
733 if (eGen>=0. && eGen<=250.){
734 Double_t eLeading = -1;
df65bddb 735 Double_t ptleading = -1;
736 Int_t nPart=0;
737 // loop over tracks
738 for (Int_t it = 0; it< nTracks; it++){
d40dc1ba 739 // if (fAOD->GetTrack(it)->E() > eGen) continue; // CKB. Not allowed! cannot cut on gen properties in real events!
df65bddb 740 // find leading particle
de6f8090 741 // if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){
742 // TODO implement esd filter flag to be the same as in the jet finder
743 // allow also for MC particles...
744 Float_t r = recJets[ir].DeltaR(fAOD->GetTrack(it));
745 if (r<0.4 && fAOD->GetTrack(it)->Pt()>ptleading){
77bbf113 746 eLeading = fAOD->GetTrack(it)->E();
df65bddb 747 ptleading = fAOD->GetTrack(it)->Pt();
748 }
d40dc1ba 749 // if (fAOD->GetTrack(it)->Pt()>0.03*eGen && fAOD->GetTrack(it)->E()<=eGen && r<0.7) // CKB cannot cut on gen properties
de6f8090 750 if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->Pt()<=eRec && r<0.7)
df65bddb 751 nPart++;
1d676c86 752
753
754 // correlate jet axis of leading jet with particles
755 if(ir==0){
756 Float_t phi = fAOD->GetTrack(it)->Phi();
757 Float_t dPhi = phi - phiRec;
758 if(dPhi>TMath::Pi()/1.5)dPhi = dPhi - 2.*TMath::Pi();
759 if(dPhi<(-0.5*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
760 Float_t dPhiRan = phi - phiRecRan;
761 if(dPhiRan>TMath::Pi()/1.5)dPhiRan = dPhiRan - 2.*TMath::Pi();
762 if(dPhiRan<(-0.5*TMath::Pi()))dPhiRan = dPhiRan + 2.*TMath::Pi();
763 ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorr"))->Fill(ptRec,dPhi);
764 ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrRan"))->Fill(ptRec,dPhiRan);
765 ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPt"))->Fill(ptRec,dPhi,fAOD->GetTrack(it)->Pt());
766 ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPtRan"))->Fill(ptRec,dPhiRan,fAOD->GetTrack(it)->Pt());
767
768 }
df65bddb 769 }
77bbf113 770 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
771
df65bddb 772 // fill Response Map (4D histogram) and Energy vs z distributions
77bbf113 773 Double_t var[4] = {eGen, eRec, ptleading/eGen, ptleading/eRec};
df65bddb 774 fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
775 fh2EGenZGen->Fill(var[0],var[2]);
776 fh1JetMultiplicity->Fill(nPart);
77bbf113 777 fh3EGenERecN->Fill(eGen, eRec, nPart);
df65bddb 778 for(int ic = 0;ic <kMaxCorrelation;ic++){
d40dc1ba 779 if (nPart<=fgkJetNpartCut[ic]){ // is this corrected for CKB
77bbf113 780 fhnCorrelation[ic]->Fill(var);
781 break;
782 }
df65bddb 783 }
784 }
d40dc1ba 785
786 }// if etarec in window
787
788 }
df65bddb 789 ////////////////////////////////////////////////////
f3050824 790 else{
77bbf113 791 fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
f3050824 792 }
793 }// loop over reconstructed jets
794
795
796 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
797 for(int ig = 0;ig < nGenJets;++ig){
798 Double_t ptGen = genJets[ig].Pt();
799 // Fill Correlation
800 Double_t phiGen = genJets[ig].Phi();
801 if(phiGen<0)phiGen+=TMath::Pi()*2.;
802 Double_t etaGen = genJets[ig].Eta();
803 fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
804 fh1PtGenIn[ig]->Fill(ptGen,eventW);
d40dc1ba 805 for(int igg = ig+1;igg<nGenJets;igg++){
806 fh2PtGenDeltaR[ig]->Fill(genJets[igg].Pt(),genJets[ig].DeltaR(&genJets[igg]));
807 }
f3050824 808 Int_t ir = iRecIndex[ig];
809 if(ir>=0&&ir<nRecJets){
810 fh1PtGenOut[ig]->Fill(ptGen,eventW);
811 }
812 else{
77bbf113 813 fh3GenEtaPhiPtNoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
f3050824 814 }
815 }// loop over reconstructed jets
816
817 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
818 PostData(1, fHistList);
819}
820
821void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
822{
823// Terminate analysis
824//
825 if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
826}