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