]> 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
Adding some histograms for jet shape
[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 }
218 xtree->SetBranchAddress("xsection",&xsection);
219 xtree->SetBranchAddress("ntrials",&ntrials);
220 ftrials = ntrials;
221 xtree->GetEntry(0);
4dbfdecc 222 }
4dbfdecc 223 fh1Xsec->Fill("<#sigma>",xsection);
60c2f4e5 224 fh1Trials->Fill("#sum{ntrials}",ftrials);
4dbfdecc 225 }
226
227 return kTRUE;
228}
f3050824 229
230void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
231{
232
233 //
234 // Create the output container
235 //
236
237
238 // Connect the AOD
239
240 if(fUseAODInput){
241 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
242 if(!fAOD){
243 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
244 return;
245 }
246 // fethc the header
247 fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
248 if(!fJetHeaderRec){
249 Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
250 }
251 }
252 else{
253 // assume that the AOD is in the general output...
254 fAOD = AODEvent();
255 if(!fAOD){
256 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
257 return;
258 }
259 fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
260 if(!fJetHeaderRec){
261 Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
262 }
263 else{
264 if(fDebug>10)fJetHeaderRec->Dump();
265 }
266 }
267
268
269
270 if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
271
272 OpenFile(1);
273 if(!fHistList)fHistList = new TList();
274
275 Bool_t oldStatus = TH1::AddDirectoryStatus();
276 TH1::AddDirectory(kFALSE);
277
278 //
279 // Histogram
280
281 const Int_t nBinPt = 100;
282 Double_t binLimitsPt[nBinPt+1];
283 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
284 if(iPt == 0){
285 binLimitsPt[iPt] = 0.0;
286 }
287 else {// 1.0
de6f8090 288 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
f3050824 289 }
290 }
d40dc1ba 291
292 const Int_t nBinEta = 26;
293 Double_t binLimitsEta[nBinEta+1] = {
294 -1.6,-1.4,-1.2,-1.0,
295 -0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
296 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
297 1.0, 1.2, 1.4, 1.6
298 };
299
f3050824 300
22ceb537 301 const Int_t nBinPhi = 30;
f3050824 302 Double_t binLimitsPhi[nBinPhi+1];
303 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
304 if(iPhi==0){
305 binLimitsPhi[iPhi] = 0;
306 }
307 else{
308 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
309 }
310 }
311
312 const Int_t nBinFrag = 25;
313
314
4dbfdecc 315 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
316 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
317
318 fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
319 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
320
f3050824 321 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
322
77bbf113 323 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
f3050824 324
77bbf113 325 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
f3050824 326
327 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
328
329 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
330
331
332 for(int ij = 0;ij<kMaxJets;++ij){
333 fh1E[ij] = new TH1F(Form("fh1E_j%d",ij),"Jet Energy;E_{jet} (GeV);N",nBinPt,binLimitsPt);
334 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
335 fh1PtRecOut[ij] = new TH1F(Form("fh1PtRecOut_j%d",ij),"rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
336 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
337 fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
338
339
d40dc1ba 340
f3050824 341 fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
342 nBinPt,binLimitsPt,nBinPt,binLimitsPt);
343
de6f8090 344 fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};#phi_{gen}",
22ceb537 345 nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
f3050824 346
de6f8090 347 fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};#eta_{gen}",
22ceb537 348 nBinEta,binLimitsEta,nBinEta,binLimitsEta);
f3050824 349
22ceb537 350
de6f8090 351 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 352 nBinPt,binLimitsPt,100,-1.0,1.0);
de6f8090 353 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 354 nBinPt,binLimitsPt,100,-1.0,1.0);
22ceb537 355
356
de6f8090 357 fh2PtRecDeltaR[ij] = new TH2F(Form("fh2PtRecDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,rec,j};#Delta R",
d40dc1ba 358 nBinPt,binLimitsPt,60,0,6.0);
de6f8090 359 fh2PtGenDeltaR[ij] = new TH2F(Form("fh2PtGenDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,gen,j};#Delta R",
d40dc1ba 360 nBinPt,binLimitsPt,60,0,6.0);
f3050824 361
362
363
364 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);
365
366
367
77bbf113 368 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 369
370
371 fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
372 nBinFrag,0.,1.,nBinPt,binLimitsPt);
373
374 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",
375 nBinFrag,0.,10.,nBinPt,binLimitsPt);
376
377 fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
378 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
379
380
381
77bbf113 382 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 383 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
384
385
77bbf113 386 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 387 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
388
389
390
de6f8090 391 fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)",
f3050824 392 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
393
394 }
df65bddb 395
1d676c86 396
397 // tmp histos do not add to the header
398 TH2F *hCorrPt = new TH2F("fh2PtRecPhiCorrPt","#Delta#phi correlation pt weighted",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
399 fHistList->Add(hCorrPt);
400 TH2F *hCorrRanPt = new TH2F("fh2PtRecPhiCorrPtRan","#Delta#phi Random correlation pt weighted",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
401 fHistList->Add(hCorrRanPt);
402
403 TH2F *hCorr = new TH2F("fh2PtRecPhiCorr","#Delta#phi correlation",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
404 fHistList->Add(hCorr);
405 TH2F *hCorrRan = new TH2F("fh2PtRecPhiCorrRan","#Delta#phi Random correlation",nBinPt,binLimitsPt,180,TMath::Pi()/-2,1.5*TMath::Pi());
406 fHistList->Add(hCorrRan);
407
408
df65bddb 409 /////////////////////////////////////////////////////////////////
410 fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
f3050824 411
df65bddb 412 fh2ERecZRec = new TH2F("fh2ERecZRec", " ; E^{jet}_{rec} [GeV]; z^{lp}_{rec}", 100, 0., 250., 100, 0., 2.);
413 fh2EGenZGen = new TH2F("fh2EGenZGen", " ; E^{jet}_{gen} [GeV]; z^{lp}_{gen}", 100, 0., 250., 100, 0., 2.);
414 fh2Efficiency = new TH2F("fh2Efficiency", "ERec/EGen;E^{jet}_{gen} [GeV];E^{jet}_{rec}/E^{jet}_{gen}", 100, 0., 250., 100, 0., 1.5);
f3050824 415
22ceb537 416 fh3EGenERecN = new TH3F("fh3EGenERecN", "Efficiency vs. Jet Multiplicity", 100, 0., 250., 100, 0., 250., 51, 0., 50.);
df65bddb 417
418 // Response map
419 //arrays for bin limits
420 const Int_t nbin[4] = {100, 100, 100, 100};
77bbf113 421 Double_t vLowEdge[4] = {0.,0.,0.,0.};
422 Double_t vUpEdge[4] = {250., 250., 1., 1.};
df65bddb 423
424 for(int ic = 0;ic < kMaxCorrelation;ic++){
77bbf113 425 fhnCorrelation[ic] = new THnSparseF(Form("fhnCorrelation_%d",ic), "Response Map", 4, nbin, vLowEdge, vUpEdge);
df65bddb 426 if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
427 else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
428 }
22ceb537 429 const Int_t saveLevel = 3; // large save level more histos
f3050824 430 if(saveLevel>0){
4dbfdecc 431 fHistList->Add(fh1Xsec);
432 fHistList->Add(fh1Trials);
f3050824 433 fHistList->Add(fh1PtHard);
77bbf113 434 fHistList->Add(fh1PtHardNoW);
435 fHistList->Add(fh1PtHardTrials);
f3050824 436 fHistList->Add(fh1NGenJets);
437 fHistList->Add(fh1NRecJets);
df65bddb 438 ////////////////////////
439 fHistList->Add(fh1JetMultiplicity);
440 fHistList->Add(fh2ERecZRec);
441 fHistList->Add(fh2EGenZGen);
442 fHistList->Add(fh2Efficiency);
443 fHistList->Add(fh3EGenERecN);
444
445 for(int ic = 0;ic < kMaxCorrelation;++ic){
446 fHistList->Add(fhnCorrelation[ic]);
447 }
448 ////////////////////////
f3050824 449 for(int ij = 0;ij<kMaxJets;++ij){
450 fHistList->Add(fh1E[ij]);
451 fHistList->Add(fh1PtRecIn[ij]);
452 fHistList->Add(fh1PtRecOut[ij]);
453 fHistList->Add(fh1PtGenIn[ij]);
454 fHistList->Add(fh1PtGenOut[ij]);
455 fHistList->Add(fh2PtFGen[ij]);
22ceb537 456 fHistList->Add(fh2PhiFGen[ij]);
457 fHistList->Add(fh2EtaFGen[ij]);
458 fHistList->Add(fh2PtGenDeltaEta[ij]);
459 fHistList->Add(fh2PtGenDeltaPhi[ij]);
d40dc1ba 460 fHistList->Add(fh2PtRecDeltaR[ij]);
461 fHistList->Add(fh2PtGenDeltaR[ij]);
22ceb537 462 fHistList->Add(fh3RecEtaPhiPt[ij]);
463 fHistList->Add(fh3GenEtaPhiPt[ij]);
f3050824 464 if(saveLevel>2){
77bbf113 465 fHistList->Add(fh3RecEtaPhiPtNoGen[ij]);
466 fHistList->Add(fh3GenEtaPhiPtNoFound[ij]);
f3050824 467 }
468 }
469 }
470
471 // =========== Switch on Sumw2 for all histos ===========
472 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
473 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
474 if (h1){
475 // Printf("%s ",h1->GetName());
476 h1->Sumw2();
4dbfdecc 477 continue;
f3050824 478 }
4dbfdecc 479 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
480 if(hn)hn->Sumw2();
f3050824 481 }
482
483 TH1::AddDirectory(oldStatus);
484
485}
486
487void AliAnalysisTaskJetSpectrum::Init()
488{
489 //
490 // Initialization
491 //
492
493 Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
494 if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
495
496}
497
498void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
499{
500 //
501 // Execute analysis for current event
502 //
503
504
505
506 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
507
508
509 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
510
511 if(!aodH){
512 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
513 return;
514 }
515
f3050824 516 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
c4396fd4 517
f3050824 518 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
519 if(!aodRecJets){
520 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
521 return;
522 }
523
524 // ==== General variables needed
525
526
527 // We use statice array, not to fragment the memory
528 AliAODJet genJets[kMaxJets];
529 Int_t nGenJets = 0;
530 AliAODJet recJets[kMaxJets];
531 Int_t nRecJets = 0;
df65bddb 532 ///////////////////////////
533 Int_t nTracks = 0;
534 //////////////////////////
f3050824 535
536 Double_t eventW = 1;
537 Double_t ptHard = 0;
538 Double_t nTrials = 1; // Trials for MC trigger weigth for real data
539
540 if(fUseExternalWeightOnly){
541 eventW = fExternalWeight;
542 }
543
544
545 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
546 if((fAnalysisType&kAnaMC)==kAnaMC){
547 // this is the part we only use when we have MC information
548 AliMCEvent* mcEvent = MCEvent();
549 // AliStack *pStack = 0;
550 if(!mcEvent){
551 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
552 return;
553 }
554 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
555 if(!pythiaGenHeader){
556 return;
557 }
558
559 nTrials = pythiaGenHeader->Trials();
560 ptHard = pythiaGenHeader->GetPtHard();
de6f8090 561 int iProcessType = pythiaGenHeader->ProcessType();
562 // 11 f+f -> f+f
563 // 12 f+barf -> f+barf
564 // 13 f+barf -> g+g
565 // 28 f+g -> f+g
566 // 53 g+g -> f+barf
567 // 68 g+g -> g+g
568 /*
569 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
570 // if(iProcessType != 13 && iProcessType != 68){ // allow only glue
571 if(iProcessType != 11 && iProcessType != 12 && iProcessType != 53){ // allow only quark
572 // if(iProcessType != 28){ // allow only -> f+g
573 PostData(1, fHistList);
574 return;
575 }
576 */
577 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
f3050824 578
579 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
580
581 if(!fUseExternalWeightOnly){
582 // case were we combine more than one p_T hard bin...
583 }
584
585 // fetch the pythia generated jets only to be used here
586 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
587 AliAODJet pythiaGenJets[kMaxJets];
588 Int_t iCount = 0;
589 for(int ip = 0;ip < nPythiaGenJets;++ip){
590 if(iCount>=kMaxJets)continue;
591 Float_t p[4];
592 pythiaGenHeader->TriggerJet(ip,p);
593 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
594
595 if(fLimitGenJetEta){
596 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
597 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
598 }
599
600
601 if(fBranchGen.Length()==0){
602 // if we have MC particles and we do not read from the aod branch
603 // use the pythia jets
604 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
605 }
606 iCount++;
607 }
608 if(fBranchGen.Length()==0)nGenJets = iCount;
609
610 }// (fAnalysisType&kMC)==kMC)
611
612 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
613 fh1PtHard->Fill(ptHard,eventW);
77bbf113 614 fh1PtHardNoW->Fill(ptHard,1);
615 fh1PtHardTrials->Fill(ptHard,nTrials);
f3050824 616
617 // If we set a second branch for the input jets fetch this
618 if(fBranchGen.Length()>0){
619 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
620 if(aodGenJets){
621 Int_t iCount = 0;
622 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
623 if(iCount>=kMaxJets)continue;
624 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
625 if(!tmp)continue;
626 if(fLimitGenJetEta){
627 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
628 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
629 }
630 genJets[iCount] = *tmp;
631 iCount++;
632 }
633 nGenJets = iCount;
634 }
635 else{
636 Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
637 }
638 }
639
640 fh1NGenJets->Fill(nGenJets);
641 // We do not want to exceed the maximum jet number
642 nGenJets = TMath::Min(nGenJets,kMaxJets);
643
644 // Fetch the reconstructed jets...
645
646
647 nRecJets = aodRecJets->GetEntries();
648 fh1NRecJets->Fill(nRecJets);
649 nRecJets = TMath::Min(nRecJets,kMaxJets);
df65bddb 650 //////////////////////////////////////////
651 nTracks = fAOD->GetNumberOfTracks();
652 ///////////////////////////////////////////
f3050824 653
654 for(int ir = 0;ir < nRecJets;++ir){
655 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
656 if(!tmp)continue;
657 recJets[ir] = *tmp;
658 }
659
660
661 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
662 // Relate the jets
663 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
664 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
665
666 for(int i = 0;i<kMaxJets;++i){
667 iGenIndex[i] = iRecIndex[i] = -1;
668 }
669
670
db6bcb0e 671 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
f3050824 672 iGenIndex,iRecIndex,fDebug);
673 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
674
675 if(fDebug){
676 for(int i = 0;i<kMaxJets;++i){
677 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
678 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
679 }
680 }
681
682 // loop over reconstructed jets
683 for(int ir = 0;ir < nRecJets;++ir){
684 Double_t ptRec = recJets[ir].Pt();
685 Double_t phiRec = recJets[ir].Phi();
1d676c86 686 Double_t phiRecRan = TMath::Pi()*gRandom->Rndm(); // better take real jet axis from previous events (TPC acceptance in phi)
f3050824 687 if(phiRec<0)phiRec+=TMath::Pi()*2.;
688 Double_t etaRec = recJets[ir].Eta();
77bbf113 689 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
f3050824 690 fh1E[ir]->Fill(recJets[ir].E(),eventW);
691 fh1PtRecIn[ir]->Fill(ptRec,eventW);
692 fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
d40dc1ba 693 for(int irr = ir+1;irr<nRecJets;irr++){
694 fh2PtRecDeltaR[ir]->Fill(recJets[irr].Pt(),recJets[ir].DeltaR(&recJets[irr]));
695 }
f3050824 696 // Fill Correlation
697 Int_t ig = iGenIndex[ir];
df65bddb 698 if(ig>=0 && ig<nGenJets){
77bbf113 699 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
700 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
f3050824 701 fh1PtRecOut[ir]->Fill(ptRec,eventW);
df65bddb 702 Double_t ptGen = genJets[ig].Pt();
703 Double_t phiGen = genJets[ig].Phi();
704 if(phiGen<0)phiGen+=TMath::Pi()*2.;
705 Double_t etaGen = genJets[ig].Eta();
d40dc1ba 706
707 //
708 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
709 //
710
711 if(TMath::Abs(etaRec)<fRecEtaWindow){
712
f3050824 713 fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
22ceb537 714 fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
715 fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
77bbf113 716 fh2PtGenDeltaEta[ir]->Fill(ptGen,etaGen-etaRec,eventW);
77bbf113 717 fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
f3050824 718 fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
77bbf113 719 fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
df65bddb 720 /////////////////////////////////////////////////////
de6f8090 721
722 // Double_t eRec = recJets[ir].E();
723 // Double_t eGen = genJets[ig].E();
724 // CKB use p_T not Energy
725 // TODO recname variabeles and histos
77bbf113 726 Double_t eRec = recJets[ir].E();
727 Double_t eGen = genJets[ig].E();
d40dc1ba 728
77bbf113 729 fh2Efficiency->Fill(eGen, eRec/eGen);
77bbf113 730
731 if (eGen>=0. && eGen<=250.){
732 Double_t eLeading = -1;
df65bddb 733 Double_t ptleading = -1;
734 Int_t nPart=0;
735 // loop over tracks
736 for (Int_t it = 0; it< nTracks; it++){
d40dc1ba 737 // if (fAOD->GetTrack(it)->E() > eGen) continue; // CKB. Not allowed! cannot cut on gen properties in real events!
df65bddb 738 // find leading particle
de6f8090 739 // if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){
740 // TODO implement esd filter flag to be the same as in the jet finder
741 // allow also for MC particles...
742 Float_t r = recJets[ir].DeltaR(fAOD->GetTrack(it));
743 if (r<0.4 && fAOD->GetTrack(it)->Pt()>ptleading){
77bbf113 744 eLeading = fAOD->GetTrack(it)->E();
df65bddb 745 ptleading = fAOD->GetTrack(it)->Pt();
746 }
d40dc1ba 747 // if (fAOD->GetTrack(it)->Pt()>0.03*eGen && fAOD->GetTrack(it)->E()<=eGen && r<0.7) // CKB cannot cut on gen properties
de6f8090 748 if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->Pt()<=eRec && r<0.7)
df65bddb 749 nPart++;
1d676c86 750
751
752 // correlate jet axis of leading jet with particles
753 if(ir==0){
754 Float_t phi = fAOD->GetTrack(it)->Phi();
755 Float_t dPhi = phi - phiRec;
756 if(dPhi>TMath::Pi()/1.5)dPhi = dPhi - 2.*TMath::Pi();
757 if(dPhi<(-0.5*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
758 Float_t dPhiRan = phi - phiRecRan;
759 if(dPhiRan>TMath::Pi()/1.5)dPhiRan = dPhiRan - 2.*TMath::Pi();
760 if(dPhiRan<(-0.5*TMath::Pi()))dPhiRan = dPhiRan + 2.*TMath::Pi();
761 ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorr"))->Fill(ptRec,dPhi);
762 ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrRan"))->Fill(ptRec,dPhiRan);
763 ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPt"))->Fill(ptRec,dPhi,fAOD->GetTrack(it)->Pt());
764 ((TH2F*)fHistList->FindObject("fh2PtRecPhiCorrPtRan"))->Fill(ptRec,dPhiRan,fAOD->GetTrack(it)->Pt());
765
766 }
df65bddb 767 }
77bbf113 768 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
769
df65bddb 770 // fill Response Map (4D histogram) and Energy vs z distributions
77bbf113 771 Double_t var[4] = {eGen, eRec, ptleading/eGen, ptleading/eRec};
df65bddb 772 fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
773 fh2EGenZGen->Fill(var[0],var[2]);
774 fh1JetMultiplicity->Fill(nPart);
77bbf113 775 fh3EGenERecN->Fill(eGen, eRec, nPart);
df65bddb 776 for(int ic = 0;ic <kMaxCorrelation;ic++){
d40dc1ba 777 if (nPart<=fgkJetNpartCut[ic]){ // is this corrected for CKB
77bbf113 778 fhnCorrelation[ic]->Fill(var);
779 break;
780 }
df65bddb 781 }
782 }
d40dc1ba 783
784 }// if etarec in window
785
786 }
df65bddb 787 ////////////////////////////////////////////////////
f3050824 788 else{
77bbf113 789 fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
f3050824 790 }
791 }// loop over reconstructed jets
792
793
794 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
795 for(int ig = 0;ig < nGenJets;++ig){
796 Double_t ptGen = genJets[ig].Pt();
797 // Fill Correlation
798 Double_t phiGen = genJets[ig].Phi();
799 if(phiGen<0)phiGen+=TMath::Pi()*2.;
800 Double_t etaGen = genJets[ig].Eta();
801 fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
802 fh1PtGenIn[ig]->Fill(ptGen,eventW);
d40dc1ba 803 for(int igg = ig+1;igg<nGenJets;igg++){
804 fh2PtGenDeltaR[ig]->Fill(genJets[igg].Pt(),genJets[ig].DeltaR(&genJets[igg]));
805 }
f3050824 806 Int_t ir = iRecIndex[ig];
807 if(ir>=0&&ir<nRecJets){
808 fh1PtGenOut[ig]->Fill(ptGen,eventW);
809 }
810 else{
77bbf113 811 fh3GenEtaPhiPtNoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
f3050824 812 }
813 }// loop over reconstructed jets
814
815 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
816 PostData(1, fHistList);
817}
818
819void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
820{
821// Terminate analysis
822//
823 if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
824}