]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Make tasks aware of AliPhysicsSlection, fixes to Three Jets Task to fill output in...
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum2.cxx
CommitLineData
3b7ffecf 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 "AliAnalysisTaskJetSpectrum2.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"
46#include "AliESDEvent.h"
47#include "AliAODEvent.h"
48#include "AliAODHandler.h"
49#include "AliAODTrack.h"
50#include "AliAODJet.h"
51#include "AliAODMCParticle.h"
52#include "AliMCEventHandler.h"
53#include "AliMCEvent.h"
54#include "AliStack.h"
55#include "AliGenPythiaEventHeader.h"
56#include "AliJetKineReaderHeader.h"
57#include "AliGenCocktailEventHeader.h"
58#include "AliInputEventHandler.h"
59
60
61#include "AliAnalysisHelperJetTasks.h"
62
63ClassImp(AliAnalysisTaskJetSpectrum2)
64
65AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
66 fJetHeaderRec(0x0),
67 fJetHeaderGen(0x0),
68 fAOD(0x0),
69 fhnCorrelation(0x0),
70 fBranchRec("jets"),
71 fBranchGen(""),
72 fUseAODInput(kFALSE),
b5cc0c6d 73 fUseGlobalSelection(kFALSE),
3b7ffecf 74 fUseExternalWeightOnly(kFALSE),
75 fLimitGenJetEta(kFALSE),
76 fFilterMask(0),
77 fAnalysisType(0),
78 fTrackTypeRec(kTrackUndef),
79 fTrackTypeGen(kTrackUndef),
80 fAvgTrials(1),
81 fExternalWeight(1),
82 fRecEtaWindow(0.5),
83 fh1Xsec(0x0),
84 fh1Trials(0x0),
85 fh1PtHard(0x0),
86 fh1PtHardNoW(0x0),
87 fh1PtHardTrials(0x0),
88 fh1NGenJets(0x0),
89 fh1NRecJets(0x0),
edfbe476 90 fh1PtTrackRec(0x0),
91 fh1SumPtTrackRec(0x0),
92 fh1SumPtTrackAreaRec(0x0),
3b7ffecf 93 fHistList(0x0)
94{
95 for(int i = 0;i < kMaxStep*2;++i){
96 fhnJetContainer[i] = 0;
97 }
98 for(int i = 0;i < kMaxJets;++i){
edfbe476 99 fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
3b7ffecf 100 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
101 }
102
103}
104
105AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
106 AliAnalysisTaskSE(name),
107 fJetHeaderRec(0x0),
108 fJetHeaderGen(0x0),
109 fAOD(0x0),
110 fhnCorrelation(0x0),
111 fBranchRec("jets"),
112 fBranchGen(""),
113 fUseAODInput(kFALSE),
b5cc0c6d 114 fUseGlobalSelection(kFALSE),
3b7ffecf 115 fUseExternalWeightOnly(kFALSE),
116 fLimitGenJetEta(kFALSE),
117 fFilterMask(0),
118 fAnalysisType(0),
119 fTrackTypeRec(kTrackUndef),
120 fTrackTypeGen(kTrackUndef),
121 fAvgTrials(1),
122 fExternalWeight(1),
123 fRecEtaWindow(0.5),
124 fh1Xsec(0x0),
125 fh1Trials(0x0),
126 fh1PtHard(0x0),
127 fh1PtHardNoW(0x0),
128 fh1PtHardTrials(0x0),
129 fh1NGenJets(0x0),
130 fh1NRecJets(0x0),
edfbe476 131 fh1PtTrackRec(0x0),
132 fh1SumPtTrackRec(0x0),
133 fh1SumPtTrackAreaRec(0x0),
3b7ffecf 134 fHistList(0x0)
135{
136
137 for(int i = 0;i < kMaxStep*2;++i){
138 fhnJetContainer[i] = 0;
139 }
140 for(int i = 0;i < kMaxJets;++i){
edfbe476 141 fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
3b7ffecf 142 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
143 }
144 DefineOutput(1, TList::Class());
145}
146
147
148
149Bool_t AliAnalysisTaskJetSpectrum2::Notify()
150{
151 //
152 // Implemented Notify() to read the cross sections
153 // and number of trials from pyxsec.root
154 //
155
156 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
157 Float_t xsection = 0;
158 Float_t ftrials = 1;
159
160 fAvgTrials = 1;
161 if(tree){
162 TFile *curfile = tree->GetCurrentFile();
163 if (!curfile) {
164 Error("Notify","No current file");
165 return kFALSE;
166 }
167 if(!fh1Xsec||!fh1Trials){
168 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
169 return kFALSE;
170 }
171 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
172 fh1Xsec->Fill("<#sigma>",xsection);
173 // construct a poor man average trials
174 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
a221560c 175 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
3b7ffecf 176 }
177 return kTRUE;
178}
179
180void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
181{
182
183 //
184 // Create the output container
185 //
186
187
188 // Connect the AOD
189
190
191 MakeJetContainer();
192
193
194 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
195
196 OpenFile(1);
197 if(!fHistList)fHistList = new TList();
198
199 Bool_t oldStatus = TH1::AddDirectoryStatus();
200 TH1::AddDirectory(kFALSE);
201
202 //
203 // Histogram
204
205 const Int_t nBinPt = 100;
206 Double_t binLimitsPt[nBinPt+1];
207 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
208 if(iPt == 0){
209 binLimitsPt[iPt] = 0.0;
210 }
211 else {// 1.0
edfbe476 212 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
3b7ffecf 213 }
214 }
215
edfbe476 216 const Int_t nBinPhi = 90;
3b7ffecf 217 Double_t binLimitsPhi[nBinPhi+1];
218 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
219 if(iPhi==0){
edfbe476 220 binLimitsPhi[iPhi] = -0.5*TMath::Pi();
3b7ffecf 221 }
222 else{
223 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
224 }
225 }
226
edfbe476 227
228
229 const Int_t nBinEta = 40;
230 Double_t binLimitsEta[nBinEta+1];
231 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
232 if(iEta==0){
233 binLimitsEta[iEta] = -2.0;
234 }
235 else{
236 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 1/(Float_t)nBinEta + 0.1;
237 }
238 }
239
3b7ffecf 240 const Int_t nBinFrag = 25;
241
242
243 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
244 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
245
246 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
247 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
248
249 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
250 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
251 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
252
253 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
254 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
255
edfbe476 256 fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
257 fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
258 fh1SumPtTrackAreaRec = new TH1F("fh1SumPtTrackAreaRec","Sum Rec track P_T #eta <0.9 / 1.8 * 2 * 0.4*0.4;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
259
3b7ffecf 260 //
261 for(int ij = 0;ij < kMaxJets;++ij){
262 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
263 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
264
edfbe476 265 fh2PhiPt[ij] = new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
266 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
267
268 fh2PhiEta[ij] = new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;p_{T,jet}",
269 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
270
271
3b7ffecf 272 fh2FragRec[ij] = new TH2F(Form("fh2FragRec_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
273 nBinFrag,0.,1.,nBinPt,binLimitsPt);
274 fh2FragLnRec[ij] = new TH2F(Form("fh2FragLnRec_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
275 nBinFrag,0.,10.,nBinPt,binLimitsPt);
276
277 fh2FragGen[ij] = new TH2F(Form("fh2FragGen_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
278 nBinFrag,0.,1.,nBinPt,binLimitsPt);
279 fh2FragLnGen[ij] = new TH2F(Form("fh2FragLnGen_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
280 nBinFrag,0.,10.,nBinPt,binLimitsPt);
281 }
282
283
284 const Int_t saveLevel = 3; // large save level more histos
285 if(saveLevel>0){
286 fHistList->Add(fh1Xsec);
287 fHistList->Add(fh1Trials);
288 fHistList->Add(fh1PtHard);
289 fHistList->Add(fh1PtHardNoW);
290 fHistList->Add(fh1PtHardTrials);
291 fHistList->Add(fh1NGenJets);
292 fHistList->Add(fh1NRecJets);
edfbe476 293 fHistList->Add(fh1PtTrackRec);
294 fHistList->Add(fh1SumPtTrackRec);
295 fHistList->Add(fh1SumPtTrackAreaRec);
3b7ffecf 296 for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
297 for(int ij = 0;ij<kMaxJets;++ij){
298 fHistList->Add( fh1PtRecIn[ij]);
299 fHistList->Add( fh1PtGenIn[ij]);
edfbe476 300 fHistList->Add( fh2PhiPt[ij]);
301 fHistList->Add( fh2PhiEta[ij]);
3b7ffecf 302 fHistList->Add( fh2FragRec[ij]);
303 fHistList->Add( fh2FragLnRec[ij]);
304 fHistList->Add( fh2FragGen[ij]);
305 fHistList->Add( fh2FragLnGen[ij]);
306 }
307 fHistList->Add(fhnCorrelation);
308 }
309
310 // =========== Switch on Sumw2 for all histos ===========
311 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
312 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
313 if (h1){
314 h1->Sumw2();
315 continue;
316 }
317 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
318 if(hn)hn->Sumw2();
319 }
320 TH1::AddDirectory(oldStatus);
321}
322
323void AliAnalysisTaskJetSpectrum2::Init()
324{
325 //
326 // Initialization
327 //
328
329 Printf(">>> AnalysisTaskJetSpectrum2::Init() debug level %d\n",fDebug);
330 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
331
332}
333
334void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
335{
b5cc0c6d 336
a221560c 337 if(!AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
b5cc0c6d 338 // no selection by the service task, we continue
a221560c 339 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
b5cc0c6d 340 PostData(1, fHistList);
341 return;
342 }
343
3b7ffecf 344 //
345 // Execute analysis for current event
346 //
7fa8b2da 347 AliESDEvent *fESD = 0;
3b7ffecf 348 if(fUseAODInput){
349 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
350 if(!fAOD){
351 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
352 return;
353 }
354 // fethc the header
355 }
356 else{
357 // assume that the AOD is in the general output...
358 fAOD = AODEvent();
359 if(!fAOD){
360 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
361 return;
362 }
7fa8b2da 363 if(fDebug>0){
364 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
365 }
3b7ffecf 366 }
367
368
369
370
371 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
372
373
374 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
375
376 if(!aodH){
377 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
378 return;
379 }
380
381 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
382 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
383 if(!aodRecJets){
384 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
385 return;
386 }
387
388 // ==== General variables needed
389
390
391 // We use statice array, not to fragment the memory
392 AliAODJet genJets[kMaxJets];
393 Int_t nGenJets = 0;
394 AliAODJet recJets[kMaxJets];
395 Int_t nRecJets = 0;
396 ///////////////////////////
599396fa 397
3b7ffecf 398
399 Double_t eventW = 1;
400 Double_t ptHard = 0;
401 Double_t nTrials = 1; // Trials for MC trigger
402
403 if(fUseExternalWeightOnly){
404 eventW = fExternalWeight;
405 }
406
407 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
408
409 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
410 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
411 // this is the part we only use when we have MC information
412 AliMCEvent* mcEvent = MCEvent();
413 // AliStack *pStack = 0;
414 if(!mcEvent){
415 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
416 return;
417 }
418 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
419 Int_t iCount = 0;
420 if(pythiaGenHeader){
421 nTrials = pythiaGenHeader->Trials();
422 ptHard = pythiaGenHeader->GetPtHard();
423 int iProcessType = pythiaGenHeader->ProcessType();
424 // 11 f+f -> f+f
425 // 12 f+barf -> f+barf
426 // 13 f+barf -> g+g
427 // 28 f+g -> f+g
428 // 53 g+g -> f+barf
429 // 68 g+g -> g+g
430 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
431 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
432
433 // fetch the pythia generated jets only to be used here
434 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
435 AliAODJet pythiaGenJets[kMaxJets];
436 for(int ip = 0;ip < nPythiaGenJets;++ip){
437 if(iCount>=kMaxJets)continue;
438 Float_t p[4];
439 pythiaGenHeader->TriggerJet(ip,p);
440 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
441
442 if(fLimitGenJetEta){
443 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
444 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
445 }
446
447
448 if(fBranchGen.Length()==0){
449 // if we have MC particles and we do not read from the aod branch
450 // use the pythia jets
451 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
452 }
453 iCount++;
454 }
455 }
456 if(fBranchGen.Length()==0)nGenJets = iCount;
457 }// (fAnalysisType&kMCESD)==kMCESD)
458
459
460 // we simply fetch the tracks/mc particles as a list of AliVParticles
461
462 TList recParticles;
463 TList genParticles;
464
465 GetListOfTracks(&recParticles,fTrackTypeRec);
466 GetListOfTracks(&genParticles,fTrackTypeGen);
467
468 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
469 fh1PtHard->Fill(ptHard,eventW);
470 fh1PtHardNoW->Fill(ptHard,1);
471 fh1PtHardTrials->Fill(ptHard,nTrials);
472
473 // If we set a second branch for the input jets fetch this
474 if(fBranchGen.Length()>0){
475 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
476 if(aodGenJets){
477 Int_t iCount = 0;
478 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
479 if(iCount>=kMaxJets)continue;
480 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
481 if(!tmp)continue;
482 if(fLimitGenJetEta){
483 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
484 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
485 }
486 genJets[iCount] = *tmp;
487 iCount++;
488 }
489 nGenJets = iCount;
490 }
491 else{
492 Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
493 fAOD->Print();
494 }
495 }
496
497 fh1NGenJets->Fill(nGenJets);
498 // We do not want to exceed the maximum jet number
499 nGenJets = TMath::Min(nGenJets,kMaxJets);
500
501 // Fetch the reconstructed jets...
502
503 nRecJets = aodRecJets->GetEntries();
504
505
506
507
508 fh1NRecJets->Fill(nRecJets);
509 nRecJets = TMath::Min(nRecJets,kMaxJets);
510
511 for(int ir = 0;ir < nRecJets;++ir){
512 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
513 if(!tmp)continue;
514 recJets[ir] = *tmp;
515 }
516
517
518 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
519 // Relate the jets
520 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
521 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
522
523 for(int i = 0;i<kMaxJets;++i){
524 iGenIndex[i] = iRecIndex[i] = -1;
525 }
526
527
528 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
529 iGenIndex,iRecIndex,fDebug);
530 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
531
532 if(fDebug){
533 for(int i = 0;i<kMaxJets;++i){
534 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
535 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
536 }
537 }
538
539
540
541
542 Double_t container[6];
543
544 // loop over generated jets
545
546 // radius; tmp, get from the jet header itself
547 // at some point, how todeal woht FastJet on MC?
548 Float_t radiusGen =0.4;
549 Float_t radiusRec =0.4;
550
551 for(int ig = 0;ig < nGenJets;++ig){
552 Double_t ptGen = genJets[ig].Pt();
553 Double_t phiGen = genJets[ig].Phi();
554 if(phiGen<0)phiGen+=TMath::Pi()*2.;
555 Double_t etaGen = genJets[ig].Eta();
556
557 container[3] = ptGen;
558 container[4] = etaGen;
559 container[5] = phiGen;
560 fhnJetContainer[kStep0]->Fill(&container[3],eventW);
561 Int_t ir = iRecIndex[ig];
562
563 if(TMath::Abs(etaGen)<fRecEtaWindow){
564 fhnJetContainer[kStep1]->Fill(&container[3],eventW);
565 fh1PtGenIn[ig]->Fill(ptGen,eventW);
566 // fill the fragmentation function
567 for(int it = 0;it<genParticles.GetEntries();++it){
568 AliVParticle *part = (AliVParticle*)genParticles.At(it);
569 if(genJets[ig].DeltaR(part)<radiusGen){
570 Float_t z = part->Pt()/ptGen;
571 Float_t lnz = -1.*TMath::Log(z);
572 fh2FragGen[ig]->Fill(z,ptGen,eventW);
573 fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
574 }
575 }
576 if(ir>=0&&ir<nRecJets){
577 fhnJetContainer[kStep3]->Fill(&container[3],eventW);
578 }
579 }
580
581 if(ir>=0&&ir<nRecJets){
582 fhnJetContainer[kStep2]->Fill(&container[3],eventW);
583 Double_t etaRec = recJets[ir].Eta();
584 if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
585 }
586 }// loop over generated jets
587
7fa8b2da 588
edfbe476 589 Float_t sumPt = 0;
590 for(int it = 0;it<recParticles.GetEntries();++it){
591 AliVParticle *part = (AliVParticle*)recParticles.At(it);
592 // fill sum pt and P_t of all paritles
593 if(TMath::Abs(part->Eta())<0.9){
594 Float_t pt = part->Pt();
595 fh1PtTrackRec->Fill(pt,eventW);
596 sumPt += pt;
597 }
598 }
599 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
600 fh1SumPtTrackRec->Fill(sumPt,eventW);
601
3b7ffecf 602
603 // loop over reconstructed jets
604 for(int ir = 0;ir < nRecJets;++ir){
605 Double_t etaRec = recJets[ir].Eta();
606 Double_t ptRec = recJets[ir].Pt();
607 Double_t phiRec = recJets[ir].Phi();
608 if(phiRec<0)phiRec+=TMath::Pi()*2.;
609 container[0] = ptRec;
610 container[1] = etaRec;
611 container[2] = phiRec;
612
edfbe476 613 if(ptRec>15.&&fDebug>0){
7fa8b2da 614 // need to cast to int, otherwise the printf overwrites
615 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
616 fAOD->GetHeader()->Print();
a221560c 617 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
7fa8b2da 618 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
619 AliAODTrack *tr = fAOD->GetTrack(it);
620 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
621 tr->Print();
622 tr->Dump();
623 if(fESD){
624 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
625 esdTr->Print("");
626 esdTr->Dump();
627 }
628 }
629 }
630
631
3b7ffecf 632 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
633 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
634 if(TMath::Abs(etaRec)<fRecEtaWindow){
635 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
636 fh1PtRecIn[ir]->Fill(ptRec,eventW);
637 // fill the fragmentation function
638 for(int it = 0;it<recParticles.GetEntries();++it){
639 AliVParticle *part = (AliVParticle*)recParticles.At(it);
edfbe476 640 Float_t eta = part->Eta();
641 if(TMath::Abs(eta)<0.9){
642 Float_t phi = part->Phi();
643 if(phi<0)phi+=TMath::Pi()*2.;
644 Float_t dPhi = phi - phiRec;
645 Float_t dEta = eta - etaRec;
646 if(dPhi<(-1.*TMath::Pi()/2))phiRec+=TMath::Pi()*2.;
647 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
648 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
649 }
3b7ffecf 650 if(recJets[ir].DeltaR(part)<radiusRec){
651 Float_t z = part->Pt()/ptRec;
652 Float_t lnz = -1.*TMath::Log(z);
653 fh2FragRec[ir]->Fill(z,ptRec,eventW);
654 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
655 }
656 }
657 }
658
659
660 // Fill Correlation
661 Int_t ig = iGenIndex[ir];
662 if(ig>=0 && ig<nGenJets){
663 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
664 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
665 Double_t ptGen = genJets[ig].Pt();
666 Double_t phiGen = genJets[ig].Phi();
667 if(phiGen<0)phiGen+=TMath::Pi()*2.;
668 Double_t etaGen = genJets[ig].Eta();
669
670 container[3] = ptGen;
671 container[4] = etaGen;
672 container[5] = phiGen;
673
674 //
675 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
676 //
677
678 if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
679 if(TMath::Abs(etaRec)<fRecEtaWindow){
680 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
681 fhnCorrelation->Fill(container);
682 }// if etarec in window
683
684 }
685 ////////////////////////////////////////////////////
686 else{
687
688 }
689 }// loop over reconstructed jets
690
691
692 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
693 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
694 PostData(1, fHistList);
695}
696
697
698void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
699 //
700 // Create the particle container for the correction framework manager and
701 // link it
702 //
703 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
704 const Double_t kPtmin = 10.0, kPtmax = 260.; // we do not want to have empty bins at the beginning...
705 const Double_t kEtamin = -3.0, kEtamax = 3.0;
706 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
707
708 // can we neglect migration in eta and phi?
709 // phi should be no problem since we cover full phi and are phi symmetric
710 // eta migration is more difficult due to needed acceptance correction
711 // in limited eta range
712
713
714 //arrays for the number of bins in each dimension
715 Int_t iBin[kNvar];
716 iBin[0] = 100; //bins in pt
717 iBin[1] = 1; //bins in eta
718 iBin[2] = 1; // bins in phi
719
720 //arrays for lower bounds :
721 Double_t* binEdges[kNvar];
722 for(Int_t ivar = 0; ivar < kNvar; ivar++)
723 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
724
725 //values for bin lower bounds
726 // for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
727 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/iBin[0]*(Double_t)i;
728 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
729 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
730
731
732 for(int i = 0;i < kMaxStep*2;++i){
733 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
734 for (int k=0; k<kNvar; k++) {
735 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
736 }
737 }
738 //create correlation matrix for unfolding
739 Int_t thnDim[2*kNvar];
740 for (int k=0; k<kNvar; k++) {
741 //first half : reconstructed
742 //second half : MC
743 thnDim[k] = iBin[k];
744 thnDim[k+kNvar] = iBin[k];
745 }
746
747 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
748 for (int k=0; k<kNvar; k++) {
749 fhnCorrelation->SetBinEdges(k,binEdges[k]);
750 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
751 }
752 fhnCorrelation->Sumw2();
753
754 // Add a histogram for Fake jets
755 // thnDim[3] = AliPID::kSPECIES;
756 // fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
757 // for(Int_t idim = 0; idim < kNvar; idim++)
758 // fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
759}
760
761void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
762{
763// Terminate analysis
764//
765 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
766}
767
768
769Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
770
771
772 Int_t iCount = 0;
773 if(type==kTrackAODIn||type==kTrackAODOut){
774 AliAODEvent *aod = 0;
775 if(type==kTrackAODIn)aod = dynamic_cast<AliAODEvent*>(InputEvent());
776 else if (type==kTrackAODOut) aod = AODEvent();
777 if(!aod){
778 return iCount;
779 }
780 for(int it = 0;it < aod->GetNumberOfTracks();++it){
781 AliAODTrack *tr = aod->GetTrack(it);
782 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
783 list->Add(tr);
784 iCount++;
785 }
786 }
787 else if (type == kTrackKineAll||type == kTrackKineCharged){
788 AliMCEvent* mcEvent = MCEvent();
789 if(!mcEvent)return iCount;
790 // we want to have alivpartilces so use get track
791 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
792 if(!mcEvent->IsPhysicalPrimary(it))continue;
793 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
794 if(type == kTrackKineAll){
795 list->Add(part);
796 iCount++;
797 }
798 else if(type == kTrackKineCharged){
799 if(part->Particle()->GetPDG()->Charge()==0)continue;
800 list->Add(part);
801 iCount++;
802 }
803 }
804 }
c4631cdb 805 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll ) {
3b7ffecf 806 AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
5010a3f7 807 if(!aod) aod = AODEvent();
808 if(!aod)return iCount;
3b7ffecf 809 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
810 if(!tca)return iCount;
811 for(int it = 0;it < tca->GetEntriesFast();++it){
812 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
813 if(!part->IsPhysicalPrimary())continue;
c4631cdb 814 if(type == kTrackAODMCAll){
3b7ffecf 815 list->Add(part);
816 iCount++;
817 }
818 else if (type == kTrackAODMCCharged){
819 if(part->Charge()==0)continue;
820 list->Add(part);
821 iCount++;
822 }
823 }
824 }// AODMCparticle
c4631cdb 825 return iCount;
3b7ffecf 826
827}