Adding some histograms for jet shape
[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   fhnCorrelationPhiZRec(0x0),
71   f1PtScale(0x0),
72                                                             fBranchRec("jets"),
73                                                             fBranchGen(""),
74                                                             fUseAODJetInput(kFALSE),
75                                                             fUseAODTrackInput(kFALSE),
76                                                             fUseAODMCInput(kFALSE),
77                                                             fUseGlobalSelection(kFALSE),
78                                                             fUseExternalWeightOnly(kFALSE),
79                                                             fLimitGenJetEta(kFALSE),
80                                                             fFilterMask(0),
81                                                             fAnalysisType(0),
82   fTrackTypeRec(kTrackUndef),
83   fTrackTypeGen(kTrackUndef),
84   fAvgTrials(1),
85   fExternalWeight(1),    
86                                                             fRecEtaWindow(0.5),
87   fMinJetPt(0),
88                                                             fDeltaPhiWindow(20./180.*TMath::Pi()),
89   fh1Xsec(0x0),
90   fh1Trials(0x0),
91   fh1PtHard(0x0),
92   fh1PtHardNoW(0x0),  
93   fh1PtHardTrials(0x0),
94   fh1NGenJets(0x0),
95   fh1NRecJets(0x0),
96   fh1PtTrackRec(0x0),   
97   fh1SumPtTrackRec(0x0),  
98   fh1SumPtTrackAreaRec(0x0),  
99   fh1TmpRho(0x0),
100   fh1PtJetsRecIn(0x0),
101   fh1PtJetsLeadingRecIn(0x0),
102   fh1PtTracksRecIn(0x0),
103   fh1PtTracksLeadingRecIn(0x0),
104   fh1PtTracksGenIn(0x0),
105   fh2NRecJetsPt(0x0),
106   fh2NRecTracksPt(0x0),
107   fh2JetsLeadingPhiEta(0x0),
108   fh2JetsLeadingPhiPt(0x0),
109   fh2TracksLeadingPhiEta(0x0),
110   fh2TracksLeadingPhiPt(0x0),
111   fh2TracksLeadingJetPhiPt(0x0),
112   fh2JetPtJetPhi(0x0),
113   fh2TrackPtTrackPhi(0x0),
114   fh2DijetDeltaPhiPt(0x0),      
115   fh2DijetAsymPt(0x0),          
116   fh2DijetAsymPtCut(0x0),       
117   fh2DijetDeltaPhiDeltaEta(0x0),
118   fh2DijetPt2vsPt1(0x0),        
119   fh2DijetDifvsSum(0x0),        
120   fh1DijetMinv(0x0),            
121   fh1DijetMinvCut(0x0),         
122   fHistList(0x0)  
123 {
124   for(int i = 0;i < kMaxStep*2;++i){
125     fhnJetContainer[i] = 0;
126   }
127   for(int i = 0;i < kMaxJets;++i){
128     fh2PhiPt[i] = fh2PhiEta[i] = fh2RhoPtRec[i] = fh2PsiPtRec[i] = fh2FragRec[i] = fh2PsiPtGen[i] = fh2FragGen[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
129     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
130   }  
131
132 }
133
134 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
135   AliAnalysisTaskSE(name),
136   fJetHeaderRec(0x0),
137   fJetHeaderGen(0x0),
138   fAOD(0x0),
139   fhnCorrelation(0x0),
140   fhnCorrelationPhiZRec(0x0),
141   f1PtScale(0x0),
142   fBranchRec("jets"),
143   fBranchGen(""),
144   fUseAODJetInput(kFALSE),
145   fUseAODTrackInput(kFALSE),
146   fUseAODMCInput(kFALSE),
147   fUseGlobalSelection(kFALSE),
148   fUseExternalWeightOnly(kFALSE),
149   fLimitGenJetEta(kFALSE),
150   fFilterMask(0),
151   fAnalysisType(0),
152   fTrackTypeRec(kTrackUndef),
153   fTrackTypeGen(kTrackUndef),
154   fAvgTrials(1),
155   fExternalWeight(1),    
156   fRecEtaWindow(0.5),
157   fMinJetPt(0),
158   fDeltaPhiWindow(20./180.*TMath::Pi()),
159   fh1Xsec(0x0),
160   fh1Trials(0x0),
161   fh1PtHard(0x0),
162   fh1PtHardNoW(0x0),  
163   fh1PtHardTrials(0x0),
164   fh1NGenJets(0x0),
165   fh1NRecJets(0x0),
166   fh1PtTrackRec(0x0),   
167   fh1SumPtTrackRec(0x0),  
168   fh1SumPtTrackAreaRec(0x0),  
169   fh1TmpRho(0x0),
170   fh1PtJetsRecIn(0x0),
171   fh1PtJetsLeadingRecIn(0x0),
172   fh1PtTracksRecIn(0x0),
173   fh1PtTracksLeadingRecIn(0x0),
174   fh1PtTracksGenIn(0x0),
175   fh2NRecJetsPt(0x0),
176   fh2NRecTracksPt(0x0),
177   fh2JetsLeadingPhiEta(0x0),
178   fh2JetsLeadingPhiPt(0x0),
179   fh2TracksLeadingPhiEta(0x0),
180   fh2TracksLeadingPhiPt(0x0),
181   fh2TracksLeadingJetPhiPt(0x0),
182   fh2JetPtJetPhi(0x0),
183   fh2TrackPtTrackPhi(0x0),
184   fh2DijetDeltaPhiPt(0x0),      
185   fh2DijetAsymPt(0x0),          
186   fh2DijetAsymPtCut(0x0),       
187   fh2DijetDeltaPhiDeltaEta(0x0),
188   fh2DijetPt2vsPt1(0x0),        
189   fh2DijetDifvsSum(0x0),        
190   fh1DijetMinv(0x0),            
191   fh1DijetMinvCut(0x0),         
192   fHistList(0x0)
193 {
194
195   for(int i = 0;i < kMaxStep*2;++i){
196     fhnJetContainer[i] = 0;
197   }  
198   for(int i = 0;i < kMaxJets;++i){
199     fh2PhiPt[i] = fh2PhiEta[i] = fh2RhoPtRec[i] = fh2PsiPtRec[i] = fh2FragRec[i] = fh2PsiPtGen[i] =  fh2RhoPtGen[i] = fh2FragGen[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
200     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
201   }
202   DefineOutput(1, TList::Class());  
203 }
204
205
206
207 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
208 {
209   //
210   // Implemented Notify() to read the cross sections
211   // and number of trials from pyxsec.root
212   // 
213
214   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
215   Float_t xsection = 0;
216   Float_t ftrials  = 1;
217
218   fAvgTrials = 1;
219   if(tree){
220     TFile *curfile = tree->GetCurrentFile();
221     if (!curfile) {
222       Error("Notify","No current file");
223       return kFALSE;
224     }
225     if(!fh1Xsec||!fh1Trials){
226       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
227       return kFALSE;
228     }
229     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
230     fh1Xsec->Fill("<#sigma>",xsection);
231     // construct a poor man average trials 
232     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
233     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
234   }  
235   return kTRUE;
236 }
237
238 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
239 {
240
241   //
242   // Create the output container
243   //
244
245
246   // Connect the AOD
247
248
249   MakeJetContainer();
250
251
252   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
253
254   OpenFile(1);
255   if(!fHistList)fHistList = new TList();
256   fHistList->SetOwner(kTRUE);
257   Bool_t oldStatus = TH1::AddDirectoryStatus();
258   TH1::AddDirectory(kFALSE);
259
260   //
261   //  Histogram
262     
263   const Int_t nBinPt = 240;
264   Double_t binLimitsPt[nBinPt+1];
265   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
266     if(iPt == 0){
267       binLimitsPt[iPt] = 0.0;
268     }
269     else {// 1.0
270       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
271     }
272   }
273   
274   const Int_t nBinPhi = 90;
275   Double_t binLimitsPhi[nBinPhi+1];
276   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
277     if(iPhi==0){
278       binLimitsPhi[iPhi] = -1.*TMath::Pi();
279     }
280     else{
281       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
282     }
283   }
284
285
286   const Int_t nBinPhi2 = 360;
287   Double_t binLimitsPhi2[nBinPhi2+1];
288   for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
289     if(iPhi2==0){
290       binLimitsPhi2[iPhi2] = 0.;
291     }
292     else{
293       binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
294     }
295   }
296
297
298
299   const Int_t nBinEta = 40;
300   Double_t binLimitsEta[nBinEta+1];
301   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
302     if(iEta==0){
303       binLimitsEta[iEta] = -2.0;
304     }
305     else{
306       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
307     }
308   }
309
310   const Int_t nBinFrag = 25;
311
312
313   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
314   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
315
316   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
317   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
318
319   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
320   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
321   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
322
323   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
324   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
325
326   fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
327   fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
328   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);
329   
330   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
331   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
332   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
333   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
334   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
335
336   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
337   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
338   // 
339
340   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
341                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
342   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
343                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
344   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
345                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
346   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
347                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
348
349   fh2JetPtJetPhi = new TH2F("fh2JetPtJetPhi","Reconstructed jet phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
350   fh2TrackPtTrackPhi = new TH2F("fh2TrackPtTrackPhi","Reconstructed track phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
351
352
353   for(int ij = 0;ij < kMaxJets;++ij){
354     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
355     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
356
357     fh2PhiPt[ij] =  new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
358                              nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
359
360     fh2PhiEta[ij] =  new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
361                               nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
362
363     fh2RhoPtRec[ij] =  new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
364                                 20,0.,1.,nBinPt,binLimitsPt);
365     fh2PsiPtRec[ij] =  new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
366                                 20,0.,1.,nBinPt,binLimitsPt);
367
368     fh2RhoPtGen[ij] =  new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
369                                 20,0.,1.,nBinPt,binLimitsPt);
370     fh2PsiPtGen[ij] =  new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
371                                 20,0.,1.,nBinPt,binLimitsPt);
372     if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",20,0.,1);
373
374
375     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",
376                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
377     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",
378                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
379
380     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",
381                               nBinFrag,0.,1.,nBinPt,binLimitsPt);
382     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",
383                                 nBinFrag,0.,10.,nBinPt,binLimitsPt);
384   }
385
386   // Dijet histograms
387
388   fh2DijetDeltaPhiPt       = new TH2F("fh2DeltaPhiPt","Difference in the azimuthal angle;#Delta#phi;p_{T,1};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
389   fh2DijetAsymPt            = new TH2F("fh2DijetAsym","Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
390   fh2DijetAsymPtCut         = new TH2F("fh2DijetAsymCut","Pt asymmetry after delta phi cut;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
391   fh2DijetDeltaPhiDeltaEta = new TH2F("fh2DijetDeltaPhiDeltaEta","Difference in the azimuthal angle;#Delta#phi;Entries",180,0.,TMath::Pi(),20,-2.,2.);
392   fh2DijetPt2vsPt1          = new TH2F("fh2DijetPt2vsPt1","Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
393   fh2DijetDifvsSum          = new TH2F("fh2DijetDifvsSum","Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
394   fh1DijetMinv               = new TH1F("fh1DijetMinv","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
395   fh1DijetMinvCut           = new TH1F("fh1DijetMinvCut","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
396
397
398
399
400   const Int_t saveLevel = 3; // large save level more histos
401   if(saveLevel>0){
402     fHistList->Add(fh1Xsec);
403     fHistList->Add(fh1Trials);
404     fHistList->Add(fh1PtHard);
405     fHistList->Add(fh1PtHardNoW);
406     fHistList->Add(fh1PtHardTrials);
407     if(fBranchGen.Length()>0){
408       fHistList->Add(fh1NGenJets);
409       fHistList->Add(fh1PtTracksGenIn);
410     }
411     fHistList->Add(fh1PtJetsRecIn);
412     fHistList->Add(fh1PtJetsLeadingRecIn);
413     fHistList->Add(fh1PtTracksRecIn);
414     fHistList->Add(fh1PtTracksLeadingRecIn);
415     fHistList->Add(fh1NRecJets);
416     fHistList->Add(fh1PtTrackRec);
417     fHistList->Add(fh1SumPtTrackRec);
418     fHistList->Add(fh1SumPtTrackAreaRec);
419     fHistList->Add(fh2NRecJetsPt);
420     fHistList->Add(fh2NRecTracksPt);
421     fHistList->Add(fh2JetsLeadingPhiEta );
422     fHistList->Add(fh2JetsLeadingPhiPt );
423     fHistList->Add(fh2TracksLeadingPhiEta);
424     fHistList->Add(fh2TracksLeadingPhiPt);
425     for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
426     for(int ij = 0;ij<kMaxJets;++ij){
427       fHistList->Add( fh1PtRecIn[ij]);
428
429       if(fBranchGen.Length()>0){
430         fHistList->Add( fh1PtGenIn[ij]);
431         fHistList->Add( fh2FragGen[ij]);
432         fHistList->Add( fh2FragLnGen[ij]);
433         fHistList->Add(fh2RhoPtGen[ij]);
434         fHistList->Add(fh2PsiPtGen[ij]);
435         fHistList->Add( fh2FragGen[ij]);
436       }
437       fHistList->Add( fh2PhiPt[ij]);
438       fHistList->Add( fh2PhiEta[ij]);
439       fHistList->Add(fh2RhoPtRec[ij]);
440       fHistList->Add(fh2PsiPtRec[ij]);
441       fHistList->Add( fh2FragRec[ij]);
442       fHistList->Add( fh2FragLnRec[ij]);
443     }
444     fHistList->Add(fhnCorrelation);
445     fHistList->Add(fhnCorrelationPhiZRec);
446     fHistList->Add(fh2JetPtJetPhi);
447     fHistList->Add(fh2TrackPtTrackPhi);
448
449     fHistList->Add(fh2DijetDeltaPhiPt);       
450     fHistList->Add(fh2DijetAsymPt);       
451     fHistList->Add(fh2DijetAsymPtCut);               
452     fHistList->Add(fh2DijetDeltaPhiDeltaEta);        
453     fHistList->Add(fh2DijetPt2vsPt1);                
454     fHistList->Add(fh2DijetDifvsSum);                
455     fHistList->Add(fh1DijetMinv);                    
456     fHistList->Add(fh1DijetMinvCut);                 
457   }
458
459   // =========== Switch on Sumw2 for all histos ===========
460   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
461     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
462     if (h1){
463       h1->Sumw2();
464       continue;
465     }
466     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
467     if(hn)hn->Sumw2();
468   }
469   TH1::AddDirectory(oldStatus);
470 }
471
472 void AliAnalysisTaskJetSpectrum2::Init()
473 {
474   //
475   // Initialization
476   //
477
478   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
479
480 }
481
482 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
483 {
484
485   if(!AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
486     // no selection by the service task, we continue
487     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
488     PostData(1, fHistList);
489     return;
490   }
491
492   //
493   // Execute analysis for current event
494   //
495   AliESDEvent *fESD = 0;
496   if(fUseAODJetInput){    
497     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
498     if(!fAOD){
499       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
500       return;
501     }
502     // fethc the header
503   }
504   else{
505     //  assume that the AOD is in the general output...
506     fAOD  = AODEvent();
507     if(!fAOD){
508       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
509       return;
510     }
511     if(fDebug>0){
512       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
513     }
514   }
515   
516
517
518
519   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
520
521   
522   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
523
524   if(!aodH){
525     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
526     return;
527   }
528
529   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
530   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
531   if(!aodRecJets){
532     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
533     return;
534   }
535
536   // ==== General variables needed
537
538
539   // We use statice array, not to fragment the memory
540   AliAODJet genJets[kMaxJets];
541   Int_t nGenJets = 0;
542   AliAODJet recJets[kMaxJets];
543   Int_t nRecJets = 0;
544   ///////////////////////////
545
546
547   Double_t eventW = 1;
548   Double_t ptHard = 0; 
549   Double_t nTrials = 1; // Trials for MC trigger 
550
551   if(fUseExternalWeightOnly){
552     eventW = fExternalWeight;
553   }
554
555   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
556   //  if(fDebug>0)aodH->SetFillAOD(kFALSE);
557   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
558   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
559     // this is the part we only use when we have MC information
560     AliMCEvent* mcEvent = MCEvent();
561     //    AliStack *pStack = 0; 
562     if(!mcEvent){
563       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
564       return;
565     }
566     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
567     Int_t iCount = 0;  
568     if(pythiaGenHeader){
569       nTrials = pythiaGenHeader->Trials();
570       ptHard  = pythiaGenHeader->GetPtHard();
571       int iProcessType = pythiaGenHeader->ProcessType();
572       // 11 f+f -> f+f
573       // 12 f+barf -> f+barf
574       // 13 f+barf -> g+g
575       // 28 f+g -> f+g
576       // 53 g+g -> f+barf
577       // 68 g+g -> g+g
578       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
579       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
580       
581       // fetch the pythia generated jets only to be used here
582       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
583       AliAODJet pythiaGenJets[kMaxJets];
584       for(int ip = 0;ip < nPythiaGenJets;++ip){
585         if(iCount>=kMaxJets)continue;
586         Float_t p[4];
587         pythiaGenHeader->TriggerJet(ip,p);
588         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
589         
590         if(fBranchGen.Length()==0){
591           /*    
592           if(fLimitGenJetEta){
593             if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
594                pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
595           }
596           */
597           if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
598           // if we have MC particles and we do not read from the aod branch
599           // use the pythia jets
600           genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
601           iCount++;
602         }
603       }
604     }
605     if(fBranchGen.Length()==0)nGenJets = iCount;    
606   }// (fAnalysisType&kMCESD)==kMCESD)
607
608
609   // we simply fetch the tracks/mc particles as a list of AliVParticles
610
611   TList recParticles;
612   TList genParticles;
613
614
615
616
617   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
618   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
619   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
620   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
621
622
623   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
624   fh1PtHard->Fill(ptHard,eventW);
625   fh1PtHardNoW->Fill(ptHard,1);
626   fh1PtHardTrials->Fill(ptHard,nTrials);
627
628   // If we set a second branch for the input jets fetch this 
629   if(fBranchGen.Length()>0){
630     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
631     if(aodGenJets){
632       Int_t iCount = 0;
633       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
634         if(iCount>=kMaxJets)continue;
635         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
636         if(!tmp)continue;
637         /*
638         if(fLimitGenJetEta){
639           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
640              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
641         }
642         */
643         if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
644         genJets[iCount] = *tmp;
645         iCount++;
646       }
647       nGenJets = iCount;
648     }
649     else{
650       if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
651       if(fDebug>2)fAOD->Print();
652     }
653   }
654
655   fh1NGenJets->Fill(nGenJets);
656   // We do not want to exceed the maximum jet number
657   nGenJets = TMath::Min(nGenJets,kMaxJets);
658
659   // Fetch the reconstructed jets...
660
661   nRecJets = aodRecJets->GetEntries();
662
663   nRecJets = aodRecJets->GetEntries();
664   fh1NRecJets->Fill(nRecJets);
665
666   // Do something with the all rec jets
667   Int_t nRecOver = nRecJets;
668
669   // check that the jets are sorted
670   Float_t ptOld = 999999;
671   for(int ir = 0;ir < nRecJets;ir++){
672     AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
673     Float_t tmpPt = tmp->Pt();
674     if(tmpPt>ptOld){
675       Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
676     }
677     ptOld = tmpPt;
678   }
679
680
681   if(nRecOver>0){
682     TIterator *recIter = aodRecJets->MakeIterator();
683     AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());  
684     Float_t pt = tmpRec->Pt();
685     if(tmpRec){
686       for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
687         Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
688         while(pt<ptCut&&tmpRec){
689           nRecOver--;
690           tmpRec = (AliAODJet*)(recIter->Next()); 
691           if(tmpRec){
692             pt = tmpRec->Pt();
693           }
694         }
695         if(nRecOver<=0)break;
696         fh2NRecJetsPt->Fill(ptCut,nRecOver);
697       }
698     }
699     recIter->Reset();
700     AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
701     Float_t phi = leading->Phi();
702     if(phi<0)phi+=TMath::Pi()*2.;    
703     Float_t eta = leading->Eta();
704     pt = leading->Pt();
705      while((tmpRec = (AliAODJet*)(recIter->Next()))){
706       Float_t tmpPt = tmpRec->Pt();
707       fh1PtJetsRecIn->Fill(tmpPt);
708       if(tmpRec==leading){
709         fh1PtJetsLeadingRecIn->Fill(tmpPt);
710         continue;
711       }
712       // correlation
713       Float_t tmpPhi =  tmpRec->Phi();
714
715       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
716       Float_t dPhi = phi - tmpRec->Phi();
717       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
718       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
719       Float_t dEta = eta - tmpRec->Eta();
720       fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
721       fh2JetsLeadingPhiPt->Fill(dPhi,pt);
722     }  
723     delete recIter;
724   }
725
726   Int_t nTrackOver = recParticles.GetSize();
727   // do the same for tracks and jets
728   if(nTrackOver>0){
729     TIterator *recIter = recParticles.MakeIterator();
730     AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());  
731     Float_t pt = tmpRec->Pt();
732     //    Printf("Leading track p_t %3.3E",pt);
733     for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
734       Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
735       while(pt<ptCut&&tmpRec){
736         nTrackOver--;
737         tmpRec = (AliVParticle*)(recIter->Next()); 
738         if(tmpRec){
739           pt = tmpRec->Pt();
740         }
741       }
742       if(nTrackOver<=0)break;
743       fh2NRecTracksPt->Fill(ptCut,nTrackOver);
744     }
745     
746     recIter->Reset();
747     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
748     Float_t phi = leading->Phi();
749     if(phi<0)phi+=TMath::Pi()*2.;    
750     Float_t eta = leading->Eta();
751     pt  = leading->Pt();
752     while((tmpRec = (AliVParticle*)(recIter->Next()))){
753       Float_t tmpPt = tmpRec->Pt();
754       fh1PtTracksRecIn->Fill(tmpPt);
755       if(tmpRec==leading){
756         fh1PtTracksLeadingRecIn->Fill(tmpPt);
757         continue;
758       }
759       // correlation
760       Float_t tmpPhi =  tmpRec->Phi();
761
762       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
763       Float_t dPhi = phi - tmpRec->Phi();
764       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
765       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
766       Float_t dEta = eta - tmpRec->Eta();
767       fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
768       fh2TracksLeadingPhiPt->Fill(dPhi,pt);
769     }  
770     delete recIter;
771   }
772   
773   if(genParticles.GetSize()){
774     TIterator *genIter = genParticles.MakeIterator();
775     AliVParticle *tmpGen = 0;
776     while((tmpGen = (AliVParticle*)(genIter->Next()))){
777       if(TMath::Abs(tmpGen->Eta())<0.9){
778         Float_t tmpPt = tmpGen->Pt();
779         fh1PtTracksGenIn->Fill(tmpPt);
780       }  
781     }
782     delete genIter;
783   }
784   
785   nRecJets = TMath::Min(nRecJets,kMaxJets);
786   
787   Int_t iCountRec = 0;
788   for(int ir = 0;ir < nRecJets;++ir){
789     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
790     if(!tmp)continue;
791     if(tmp->Pt()<fMinJetPt)continue;
792     recJets[ir] = *tmp;
793     iCountRec++;
794   }
795   nRecJets = iCountRec;
796
797
798   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
799   // Relate the jets
800   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
801   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
802   
803
804   for(int i = 0;i<kMaxJets;++i){    
805     iGenIndex[i] = iRecIndex[i] = -1;
806   }
807
808   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
809                                             iGenIndex,iRecIndex,fDebug);
810   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
811
812   if(fDebug){
813     for(int i = 0;i<kMaxJets;++i){
814       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]); 
815       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]); 
816     }
817   }
818
819
820
821
822   Double_t container[6];
823   Double_t containerPhiZ[6];
824
825   // loop over generated jets
826
827   // radius; tmp, get from the jet header itself
828   // at some point, how todeal woht FastJet on MC?
829   Float_t radiusGen =0.4;
830   Float_t radiusRec =0.4;
831
832   for(int ig = 0;ig < nGenJets;++ig){
833     Double_t ptGen = genJets[ig].Pt();
834     Double_t phiGen = genJets[ig].Phi();
835     if(phiGen<0)phiGen+=TMath::Pi()*2.;    
836     Double_t etaGen = genJets[ig].Eta();
837     
838     container[3] = ptGen;
839     container[4] = etaGen;
840     container[5] = phiGen;
841     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
842     Int_t ir = iRecIndex[ig];
843
844     if(TMath::Abs(etaGen)<fRecEtaWindow){
845       fh1TmpRho->Reset();
846
847       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
848       fh1PtGenIn[ig]->Fill(ptGen,eventW);
849       // fill the fragmentation function
850       for(int it = 0;it<genParticles.GetEntries();++it){
851         AliVParticle *part = (AliVParticle*)genParticles.At(it);
852         Float_t deltaR = genJets[ig].DeltaR(part);
853         fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
854         if(deltaR<radiusGen){
855           Float_t z = part->Pt()/ptGen;
856           Float_t lnz =  -1.*TMath::Log(z);
857           fh2FragGen[ig]->Fill(z,ptGen,eventW);
858           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
859         }
860
861       }
862       Float_t rhoSum = 0;
863       for(int ibx = 1;ibx<fh2RhoPtGen[ir]->GetNbinsX();ibx++){
864         Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
865         Float_t rho = fh1TmpRho->GetBinContent(ibx);
866         rhoSum += rho;
867         fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
868         fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
869       }
870     }
871     if(ir>=0&&ir<nRecJets){   
872       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
873       Double_t etaRec = recJets[ir].Eta();
874       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
875     }
876   }// loop over generated jets
877   
878   Float_t sumPt = 0;
879   for(int it = 0;it<recParticles.GetEntries();++it){
880     AliVParticle *part = (AliVParticle*)recParticles.At(it);
881     // fill sum pt and P_t of all paritles
882     if(TMath::Abs(part->Eta())<0.9){
883       Float_t pt = part->Pt();
884       fh1PtTrackRec->Fill(pt,eventW);
885       fh2TrackPtTrackPhi->Fill(pt,part->Phi());
886       sumPt += pt;
887     }
888   }
889   fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
890   fh1SumPtTrackRec->Fill(sumPt,eventW);
891
892   
893   // loop over reconstructed jets
894   for(int ir = 0;ir < nRecJets;++ir){
895     Double_t etaRec = recJets[ir].Eta();
896     Double_t ptRec = recJets[ir].Pt();
897     Double_t phiRec = recJets[ir].Phi();
898     if(phiRec<0)phiRec+=TMath::Pi()*2.;    
899
900
901   // do something with dijets...
902   if(ir==1){
903     Double_t etaRec1 = recJets[0].Eta();
904     Double_t ptRec1 = recJets[0].Pt();
905     Double_t phiRec1 = recJets[0].Phi();
906     if(phiRec1<0)phiRec1+=TMath::Pi()*2.;    
907     
908   
909     if(TMath::Abs(etaRec1)<fRecEtaWindow
910       &&TMath::Abs(etaRec)<fRecEtaWindow){
911   
912       Float_t deltaPhi = phiRec1 - phiRec;
913       
914       if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
915       if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();      
916       deltaPhi = TMath::Abs(deltaPhi);
917       fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);      
918       Float_t asym = (ptRec1-ptRec)/(ptRec1+ptRec);
919       fh2DijetAsymPt->Fill(asym,ptRec1);
920       fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
921       fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);        
922       fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);        
923       Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
924                          recJets[0].Px()*recJets[1].Px()- 
925                          recJets[0].Py()*recJets[1].Py()- 
926                          recJets[0].Pz()*recJets[1].Py());
927       minv = TMath::Sqrt(minv);
928       // with mass == 0;
929       
930       fh1DijetMinv->Fill(minv);            
931       if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
932         fh1DijetMinvCut->Fill(minv);         
933         fh2DijetAsymPtCut->Fill(asym,ptRec1);      
934       }
935     }
936   }
937
938
939     container[0] = ptRec;
940     container[1] = etaRec;
941     container[2] = phiRec;
942     containerPhiZ[0] = ptRec;
943     containerPhiZ[1] = phiRec;
944     if(ptRec>30.&&fDebug>0){
945       // need to cast to int, otherwise the printf overwrites
946       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
947       Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
948       if(fESD)Printf("ESDEvent  GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
949       //  aodH->SetFillAOD(kTRUE);
950       fAOD->GetHeader()->Print();
951       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
952       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
953         AliAODTrack *tr = fAOD->GetTrack(it);
954         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
955         tr->Print();
956         //      tr->Dump();
957         if(fESD){
958           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
959           esdTr->Print("");
960           // esdTr->Dump();
961         }
962       }
963     }
964   
965
966     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
967     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
968
969     Float_t zLeading = -1;
970     if(TMath::Abs(etaRec)<fRecEtaWindow){
971       fh2JetPtJetPhi->Fill(ptRec,phiRec);
972       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
973       fh1PtRecIn[ir]->Fill(ptRec,eventW);
974       // fill the fragmentation function
975       
976       fh1TmpRho->Reset();
977
978       for(int it = 0;it<recParticles.GetEntries();++it){
979         AliVParticle *part = (AliVParticle*)recParticles.At(it);
980         Float_t eta = part->Eta();
981         if(TMath::Abs(eta)<0.9){
982           Float_t phi = part->Phi();
983           if(phi<0)phi+=TMath::Pi()*2.;    
984           Float_t dPhi = phi - phiRec;
985           Float_t dEta = eta - etaRec;
986           if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
987           if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
988           fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
989           fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
990         }
991
992         Float_t deltaR = recJets[ir].DeltaR(part);
993         fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
994
995
996         if(deltaR<radiusRec){
997           Float_t z = part->Pt()/ptRec;
998           if(z>zLeading)zLeading=z;
999           Float_t lnz =  -1.*TMath::Log(z);
1000           fh2FragRec[ir]->Fill(z,ptRec,eventW);
1001           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1002         }
1003       }
1004       // fill the jet shapes
1005       Float_t rhoSum = 0;
1006       for(int ibx = 1;ibx<fh2RhoPtRec[ir]->GetNbinsX();ibx++){
1007         Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1008         Float_t rho = fh1TmpRho->GetBinContent(ibx);
1009         rhoSum += rho;
1010         fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1011         fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1012       }
1013     }
1014
1015
1016     containerPhiZ[2] = zLeading;
1017
1018     // Fill Correlation
1019     Int_t ig = iGenIndex[ir];
1020     if(ig>=0 && ig<nGenJets){
1021       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1022       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1023       Double_t ptGen  = genJets[ig].Pt();
1024       Double_t phiGen = genJets[ig].Phi();
1025       if(phiGen<0)phiGen+=TMath::Pi()*2.; 
1026       Double_t etaGen = genJets[ig].Eta();
1027
1028       container[3] = ptGen;
1029       container[4] = etaGen;
1030       container[5] = phiGen;
1031       containerPhiZ[3] = ptGen;
1032       // 
1033       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1034       // 
1035
1036       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1037       if(TMath::Abs(etaRec)<fRecEtaWindow){
1038         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1039         fhnCorrelation->Fill(container);
1040         if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1041
1042       }// if etarec in window
1043
1044     } 
1045     else{
1046       containerPhiZ[3] = 0;
1047       if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1048     }
1049   }// loop over reconstructed jets
1050
1051
1052   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1053   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1054   PostData(1, fHistList);
1055 }
1056
1057
1058 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1059   //
1060   // Create the particle container for the correction framework manager and 
1061   // link it
1062   //
1063   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1064   const Double_t kPtmin = 0.0, kPtmax = 160.; // we do not want to have empty bins at the beginning...
1065   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1066   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1067   const Double_t kZmin = 0., kZmax = 1;
1068
1069   // can we neglect migration in eta and phi?
1070   // phi should be no problem since we cover full phi and are phi symmetric
1071   // eta migration is more difficult  due to needed acceptance correction
1072   // in limited eta range
1073
1074   //arrays for the number of bins in each dimension
1075   Int_t iBin[kNvar];
1076   iBin[0] = 160; //bins in pt
1077   iBin[1] =  1; //bins in eta 
1078   iBin[2] = 1; // bins in phi
1079
1080   //arrays for lower bounds :
1081   Double_t* binEdges[kNvar];
1082   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1083     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1084
1085   //values for bin lower bounds
1086   //  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);  
1087   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1088   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1089   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1090
1091
1092   for(int i = 0;i < kMaxStep*2;++i){
1093     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
1094     for (int k=0; k<kNvar; k++) {
1095       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1096     }
1097   }
1098   //create correlation matrix for unfolding
1099   Int_t thnDim[2*kNvar];
1100   for (int k=0; k<kNvar; k++) {
1101     //first half  : reconstructed 
1102     //second half : MC
1103     thnDim[k]      = iBin[k];
1104     thnDim[k+kNvar] = iBin[k];
1105   }
1106
1107   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1108   for (int k=0; k<kNvar; k++) {
1109     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1110     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1111   }
1112   fhnCorrelation->Sumw2();
1113
1114   // for second correlation histogram
1115
1116
1117   const Int_t kNvarPhiZ   = 4; 
1118   //arrays for the number of bins in each dimension
1119   Int_t iBinPhiZ[kNvarPhiZ];
1120   iBinPhiZ[0] = 80; //bins in pt
1121   iBinPhiZ[1] = 72; //bins in phi 
1122   iBinPhiZ[2] = 20; // bins in Z
1123   iBinPhiZ[3] = 80; //bins in ptgen
1124
1125   //arrays for lower bounds :
1126   Double_t* binEdgesPhiZ[kNvarPhiZ];
1127   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1128     binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1129
1130   for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1131   for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1132   for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin  + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1133   for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1134
1135   fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1136   for (int k=0; k<kNvarPhiZ; k++) {
1137     fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1138   }
1139   fhnCorrelationPhiZRec->Sumw2();
1140
1141
1142   // Add a histogram for Fake jets
1143
1144   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1145     delete [] binEdges[ivar];
1146
1147   for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1148     delete [] binEdgesPhiZ[ivar];
1149
1150 }
1151
1152 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1153 {
1154 // Terminate analysis
1155 //
1156     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1157 }
1158
1159
1160 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1161
1162   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1163
1164   Int_t iCount = 0;
1165   if(type==kTrackAOD){
1166     AliAODEvent *aod = 0;
1167     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1168     else aod = AODEvent();
1169     if(!aod){
1170       return iCount;
1171     }
1172     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1173       AliAODTrack *tr = aod->GetTrack(it);
1174       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1175       if(TMath::Abs(tr->Eta())>0.9)continue;
1176       if(fDebug>0){
1177         if(tr->Pt()>20){
1178           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1179           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
1180           tr->Print();
1181           //    tr->Dump();
1182           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1183           if(fESD){
1184             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1185             esdTr->Print("");
1186             // esdTr->Dump();
1187           }
1188         }
1189       }
1190       list->Add(tr);
1191       iCount++;
1192     }
1193   }
1194   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1195     AliMCEvent* mcEvent = MCEvent();
1196     if(!mcEvent)return iCount;
1197     // we want to have alivpartilces so use get track
1198     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1199       if(!mcEvent->IsPhysicalPrimary(it))continue;
1200       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1201       if(type == kTrackKineAll){
1202         list->Add(part);
1203         iCount++;
1204       }
1205       else if(type == kTrackKineCharged){
1206         if(part->Particle()->GetPDG()->Charge()==0)continue;
1207         list->Add(part);
1208         iCount++;
1209       }
1210     }
1211   }
1212   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1213     AliAODEvent *aod = 0;
1214     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1215     else aod = AODEvent();
1216     if(!aod)return iCount;
1217     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1218     if(!tca)return iCount;
1219     for(int it = 0;it < tca->GetEntriesFast();++it){
1220       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1221       if(!part->IsPhysicalPrimary())continue;
1222       if(type == kTrackAODMCAll){
1223         list->Add(part);
1224         iCount++;
1225       }
1226       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1227         if(part->Charge()==0)continue;
1228         if(kTrackAODMCCharged){
1229           list->Add(part);
1230         }
1231         else {
1232           if(TMath::Abs(part->Eta())>0.9)continue;
1233           list->Add(part);
1234         }
1235         iCount++;
1236       }
1237     }
1238   }// AODMCparticle
1239   list->Sort();
1240   return iCount;
1241
1242 }