Switched on vertex constrained again for Jet SPectrum, limit some printout in three...
[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
3e3987f1 613 if(ptRec>10.&&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);
a923bd34 616 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
7fa8b2da 617 fAOD->GetHeader()->Print();
a221560c 618 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
7fa8b2da 619 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
620 AliAODTrack *tr = fAOD->GetTrack(it);
621 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
622 tr->Print();
623 tr->Dump();
624 if(fESD){
625 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
626 esdTr->Print("");
627 esdTr->Dump();
628 }
629 }
630 }
631
632
3b7ffecf 633 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
634 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
635 if(TMath::Abs(etaRec)<fRecEtaWindow){
636 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
637 fh1PtRecIn[ir]->Fill(ptRec,eventW);
638 // fill the fragmentation function
639 for(int it = 0;it<recParticles.GetEntries();++it){
640 AliVParticle *part = (AliVParticle*)recParticles.At(it);
edfbe476 641 Float_t eta = part->Eta();
642 if(TMath::Abs(eta)<0.9){
643 Float_t phi = part->Phi();
644 if(phi<0)phi+=TMath::Pi()*2.;
645 Float_t dPhi = phi - phiRec;
646 Float_t dEta = eta - etaRec;
647 if(dPhi<(-1.*TMath::Pi()/2))phiRec+=TMath::Pi()*2.;
648 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
649 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
650 }
3b7ffecf 651 if(recJets[ir].DeltaR(part)<radiusRec){
652 Float_t z = part->Pt()/ptRec;
653 Float_t lnz = -1.*TMath::Log(z);
654 fh2FragRec[ir]->Fill(z,ptRec,eventW);
655 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
656 }
657 }
658 }
659
660
661 // Fill Correlation
662 Int_t ig = iGenIndex[ir];
663 if(ig>=0 && ig<nGenJets){
664 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
665 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
666 Double_t ptGen = genJets[ig].Pt();
667 Double_t phiGen = genJets[ig].Phi();
668 if(phiGen<0)phiGen+=TMath::Pi()*2.;
669 Double_t etaGen = genJets[ig].Eta();
670
671 container[3] = ptGen;
672 container[4] = etaGen;
673 container[5] = phiGen;
674
675 //
676 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
677 //
678
679 if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
680 if(TMath::Abs(etaRec)<fRecEtaWindow){
681 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
682 fhnCorrelation->Fill(container);
683 }// if etarec in window
684
685 }
686 ////////////////////////////////////////////////////
687 else{
688
689 }
690 }// loop over reconstructed jets
691
692
693 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
694 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
695 PostData(1, fHistList);
696}
697
698
699void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
700 //
701 // Create the particle container for the correction framework manager and
702 // link it
703 //
704 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
a923bd34 705 const Double_t kPtmin = 5.0, kPtmax = 105.; // we do not want to have empty bins at the beginning...
3b7ffecf 706 const Double_t kEtamin = -3.0, kEtamax = 3.0;
707 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
708
709 // can we neglect migration in eta and phi?
710 // phi should be no problem since we cover full phi and are phi symmetric
711 // eta migration is more difficult due to needed acceptance correction
712 // in limited eta range
713
714
715 //arrays for the number of bins in each dimension
716 Int_t iBin[kNvar];
717 iBin[0] = 100; //bins in pt
718 iBin[1] = 1; //bins in eta
719 iBin[2] = 1; // bins in phi
720
721 //arrays for lower bounds :
722 Double_t* binEdges[kNvar];
723 for(Int_t ivar = 0; ivar < kNvar; ivar++)
724 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
725
726 //values for bin lower bounds
727 // 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);
a923bd34 728 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
3b7ffecf 729 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
730 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
731
732
733 for(int i = 0;i < kMaxStep*2;++i){
734 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
735 for (int k=0; k<kNvar; k++) {
736 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
737 }
738 }
739 //create correlation matrix for unfolding
740 Int_t thnDim[2*kNvar];
741 for (int k=0; k<kNvar; k++) {
742 //first half : reconstructed
743 //second half : MC
744 thnDim[k] = iBin[k];
745 thnDim[k+kNvar] = iBin[k];
746 }
747
748 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
749 for (int k=0; k<kNvar; k++) {
750 fhnCorrelation->SetBinEdges(k,binEdges[k]);
751 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
752 }
753 fhnCorrelation->Sumw2();
754
755 // Add a histogram for Fake jets
756 // thnDim[3] = AliPID::kSPECIES;
757 // fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
758 // for(Int_t idim = 0; idim < kNvar; idim++)
759 // fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
760}
761
762void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
763{
764// Terminate analysis
765//
766 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
767}
768
769
770Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
771
772
773 Int_t iCount = 0;
774 if(type==kTrackAODIn||type==kTrackAODOut){
775 AliAODEvent *aod = 0;
776 if(type==kTrackAODIn)aod = dynamic_cast<AliAODEvent*>(InputEvent());
777 else if (type==kTrackAODOut) aod = AODEvent();
778 if(!aod){
779 return iCount;
780 }
781 for(int it = 0;it < aod->GetNumberOfTracks();++it){
782 AliAODTrack *tr = aod->GetTrack(it);
783 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
784 list->Add(tr);
785 iCount++;
786 }
787 }
788 else if (type == kTrackKineAll||type == kTrackKineCharged){
789 AliMCEvent* mcEvent = MCEvent();
790 if(!mcEvent)return iCount;
791 // we want to have alivpartilces so use get track
792 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
793 if(!mcEvent->IsPhysicalPrimary(it))continue;
794 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
795 if(type == kTrackKineAll){
796 list->Add(part);
797 iCount++;
798 }
799 else if(type == kTrackKineCharged){
800 if(part->Particle()->GetPDG()->Charge()==0)continue;
801 list->Add(part);
802 iCount++;
803 }
804 }
805 }
c4631cdb 806 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll ) {
3b7ffecf 807 AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
5010a3f7 808 if(!aod) aod = AODEvent();
809 if(!aod)return iCount;
3b7ffecf 810 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
811 if(!tca)return iCount;
812 for(int it = 0;it < tca->GetEntriesFast();++it){
813 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
814 if(!part->IsPhysicalPrimary())continue;
c4631cdb 815 if(type == kTrackAODMCAll){
3b7ffecf 816 list->Add(part);
817 iCount++;
818 }
819 else if (type == kTrackAODMCCharged){
820 if(part->Charge()==0)continue;
821 list->Add(part);
822 iCount++;
823 }
824 }
825 }// AODMCparticle
c4631cdb 826 return iCount;
3b7ffecf 827
828}