878ed669ebfb661f40244c9e2cd9547ea4c5ad8d
[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   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
101 AliAnalysisTaskJetSpectrum2::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
141 Bool_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
172 void 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   
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
286 void 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
297 void 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   ///////////////////////////
349
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
607 void 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
670 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
671 {
672 // Terminate analysis
673 //
674     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
675 }
676
677
678 Int_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 == kTrackAODMCAll ) {
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 == kTrackAODMCAll){
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   return iCount;
736
737 }