New book-keeping class
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum2.cxx
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
63 ClassImp(AliAnalysisTaskJetSpectrum2)
64
65 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
66   fJetHeaderRec(0x0),
67   fJetHeaderGen(0x0),
68   fAOD(0x0),
69   fhnCorrelation(0x0),
70   fBranchRec("jets"),
71   fBranchGen(""),
72   fUseAODInput(kFALSE),
73   fUseGlobalSelection(kFALSE),
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
102 AliAnalysisTaskJetSpectrum2::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),
111   fUseGlobalSelection(kFALSE),
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
143 Bool_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 && nEntries>0.)fAvgTrials = ftrials/nEntries;
170   }  
171   return kTRUE;
172 }
173
174 void 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   
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
288 void 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
299 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
300 {
301
302   if(!AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
303     // no selection by the service task, we continue
304     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
305     PostData(1, fHistList);
306     return;
307   }
308
309   //
310   // Execute analysis for current event
311   //
312   AliESDEvent *fESD = 0;
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     }
328     if(fDebug>0){
329       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
330     }
331
332
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   ///////////////////////////
364
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
555   
556   
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
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       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
573       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
574         AliAODTrack *tr = fAOD->GetTrack(it);
575         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
576         tr->Print();
577         tr->Dump();
578         if(fESD){
579           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
580           esdTr->Print("");
581           esdTr->Dump();
582         }
583       }
584     }
585   
586
587     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
588     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
589     if(TMath::Abs(etaRec)<fRecEtaWindow){
590       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
591       fh1PtRecIn[ir]->Fill(ptRec,eventW);
592       // fill the fragmentation function
593       for(int it = 0;it<recParticles.GetEntries();++it){
594         AliVParticle *part = (AliVParticle*)recParticles.At(it);
595         if(recJets[ir].DeltaR(part)<radiusRec){
596           Float_t z = part->Pt()/ptRec;
597           Float_t lnz =  -1.*TMath::Log(z);
598           fh2FragRec[ir]->Fill(z,ptRec,eventW);
599           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
600         }
601       }
602     }
603
604
605     // Fill Correlation
606     Int_t ig = iGenIndex[ir];
607     if(ig>=0 && ig<nGenJets){
608       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
609       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
610       Double_t ptGen  = genJets[ig].Pt();
611       Double_t phiGen = genJets[ig].Phi();
612       if(phiGen<0)phiGen+=TMath::Pi()*2.; 
613       Double_t etaGen = genJets[ig].Eta();
614
615       container[3] = ptGen;
616       container[4] = etaGen;
617       container[5] = phiGen;
618
619       // 
620       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
621       // 
622
623       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
624       if(TMath::Abs(etaRec)<fRecEtaWindow){
625         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
626         fhnCorrelation->Fill(container);
627       }// if etarec in window
628
629     } 
630     ////////////////////////////////////////////////////  
631     else{
632
633     }
634   }// loop over reconstructed jets
635
636
637   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
638   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
639   PostData(1, fHistList);
640 }
641
642
643 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
644   //
645   // Create the particle container for the correction framework manager and 
646   // link it
647   //
648   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
649   const Double_t kPtmin = 10.0, kPtmax = 260.; // we do not want to have empty bins at the beginning...
650   const Double_t kEtamin = -3.0, kEtamax = 3.0;
651   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
652
653   // can we neglect migration in eta and phi?
654   // phi should be no problem since we cover full phi and are phi symmetric
655   // eta migration is more difficult  due to needed acceptance correction
656   // in limited eta range
657
658
659   //arrays for the number of bins in each dimension
660   Int_t iBin[kNvar];
661   iBin[0] = 100; //bins in pt
662   iBin[1] =  1; //bins in eta 
663   iBin[2] = 1; // bins in phi
664
665   //arrays for lower bounds :
666   Double_t* binEdges[kNvar];
667   for(Int_t ivar = 0; ivar < kNvar; ivar++)
668     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
669
670   //values for bin lower bounds
671   //  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);  
672   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/iBin[0]*(Double_t)i;
673   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
674   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
675
676
677   for(int i = 0;i < kMaxStep*2;++i){
678     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
679     for (int k=0; k<kNvar; k++) {
680       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
681     }
682   }
683   //create correlation matrix for unfolding
684   Int_t thnDim[2*kNvar];
685   for (int k=0; k<kNvar; k++) {
686     //first half  : reconstructed 
687     //second half : MC
688     thnDim[k]      = iBin[k];
689     thnDim[k+kNvar] = iBin[k];
690   }
691
692   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
693   for (int k=0; k<kNvar; k++) {
694     fhnCorrelation->SetBinEdges(k,binEdges[k]);
695     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
696   }
697   fhnCorrelation->Sumw2();
698
699   // Add a histogram for Fake jets
700   //  thnDim[3] = AliPID::kSPECIES;
701   //  fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
702   // for(Int_t idim = 0; idim < kNvar; idim++)
703   //  fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
704 }
705
706 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
707 {
708 // Terminate analysis
709 //
710     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
711 }
712
713
714 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
715
716
717   Int_t iCount = 0;
718   if(type==kTrackAODIn||type==kTrackAODOut){
719     AliAODEvent *aod = 0;
720     if(type==kTrackAODIn)aod = dynamic_cast<AliAODEvent*>(InputEvent());
721     else if (type==kTrackAODOut) aod = AODEvent();
722     if(!aod){
723       return iCount;
724     }
725     for(int it = 0;it < aod->GetNumberOfTracks();++it){
726       AliAODTrack *tr = aod->GetTrack(it);
727       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
728       list->Add(tr);
729       iCount++;
730     }
731   }
732   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
733     AliMCEvent* mcEvent = MCEvent();
734     if(!mcEvent)return iCount;
735     // we want to have alivpartilces so use get track
736     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
737       if(!mcEvent->IsPhysicalPrimary(it))continue;
738       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
739       if(type == kTrackKineAll){
740         list->Add(part);
741         iCount++;
742       }
743       else if(type == kTrackKineCharged){
744         if(part->Particle()->GetPDG()->Charge()==0)continue;
745         list->Add(part);
746         iCount++;
747       }
748     }
749   }
750   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll ) {
751     AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
752     if(!aod) aod = AODEvent();
753     if(!aod)return iCount;
754     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
755     if(!tca)return iCount;
756     for(int it = 0;it < tca->GetEntriesFast();++it){
757       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
758       if(!part->IsPhysicalPrimary())continue;
759       if(type == kTrackAODMCAll){
760         list->Add(part);
761         iCount++;
762       }
763       else if (type == kTrackAODMCCharged){
764         if(part->Charge()==0)continue;
765         list->Add(part);
766         iCount++;
767       }
768     }
769   }// AODMCparticle
770   return iCount;
771
772 }