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