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