flag jets with too high p_T track
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum2.cxx
1 // **************************************
2 // 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 <TRandom3.h>
26 #include <TSystem.h>
27 #include <TInterpreter.h>
28 #include <TChain.h>
29 #include <TFile.h>
30 #include <TKey.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TH3F.h>
34 #include <TProfile.h>
35 #include <TProfile2D.h>
36 #include <TList.h>
37 #include <TLorentzVector.h>
38 #include <TClonesArray.h>
39 #include <TRefArray.h>
40 #include  "TDatabasePDG.h"
41
42 #include "AliAnalysisTaskJetSpectrum2.h"
43 #include "AliAnalysisManager.h"
44 #include "AliJetFinder.h"
45 #include "AliTHn.h"
46 #include "AliJetHeader.h"
47 #include "AliJetReader.h"
48 #include "AliJetReaderHeader.h"
49 #include "AliUA1JetHeaderV1.h"
50 #include "AliESDEvent.h"
51 #include "AliAODEvent.h"
52 #include "AliAODHandler.h"
53 #include "AliAODTrack.h"
54 #include "AliAODJet.h"
55 #include "AliAODJetEventBackground.h"
56 #include "AliAODMCParticle.h"
57 #include "AliMCEventHandler.h"
58 #include "AliMCEvent.h"
59 #include "AliStack.h"
60 #include "AliGenPythiaEventHeader.h"
61 #include "AliJetKineReaderHeader.h"
62 #include "AliGenCocktailEventHeader.h"
63 #include "AliInputEventHandler.h"
64
65
66 #include "AliAnalysisHelperJetTasks.h"
67
68 ClassImp(AliAnalysisTaskJetSpectrum2)
69
70 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): 
71   AliAnalysisTaskSE(),
72   fJetHeaderRec(0x0),
73   fJetHeaderGen(0x0),
74   fAODIn(0x0),
75   fAODOut(0x0),
76   fAODExtension(0x0),
77   fhnJetContainer(0x0),
78   fhnCorrelation(0x0),
79   fhnEvent(0x0),
80   f1PtScale(0x0),
81   fBranchRec("jets"),
82   fBranchGen(""),
83   fBranchBkgRec(""), 
84   fBranchBkgGen(""), 
85   fNonStdFile(""),
86   fRandomizer(0x0),
87   fUseAODJetInput(kFALSE),
88   fUseAODTrackInput(kFALSE),
89   fUseAODMCInput(kFALSE),
90   fUseGlobalSelection(kFALSE),
91   fUseExternalWeightOnly(kFALSE),
92   fLimitGenJetEta(kFALSE),
93   fDoMatching(kFALSE),
94   fNMatchJets(5),
95   fNRPBins(3),
96   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
97   fFilterMask(0),
98   fEventSelectionMask(0),
99   fAnalysisType(0),
100   fTrackTypeRec(kTrackUndef),
101   fTrackTypeGen(kTrackUndef),
102   fEventClass(0),
103   fRPMethod(0),
104   fAvgTrials(1),
105   fExternalWeight(1),    
106   fJetRecEtaWindow(0.5),
107   fTrackRecEtaWindow(0.5),
108   fMinJetPt(0),
109   fMinTrackPt(0.15),
110   fDeltaPhiWindow(90./180.*TMath::Pi()),
111   fCentrality(100),
112   fRPAngle(0),
113   fMultRec(0),
114   fMultGen(0),
115   fh1Xsec(0x0),
116   fh1Trials(0x0),
117   fh1PtHard(0x0),
118   fh1PtHardNoW(0x0),  
119   fh1PtHardTrials(0x0),
120   fh1ZVtx(0x0),
121   fh1RP(0x0),
122   fh1Centrality(0x0),
123   fh1TmpRho(0x0),
124   fh2MultRec(0x0),
125   fh2MultGen(0x0),
126   fh2RPCentrality(0x0),
127   fh2PtFGen(0x0),
128   fh2RelPtFGen(0x0),
129   fHistList(0x0)  
130 {
131
132   for(int ij = 0;ij <kJetTypes;++ij){    
133     fFlagJetType[ij] = 1; // default = on
134     fh1NJets[ij] = 0;
135     fh1SumPtTrack[ij] = 0;
136     fh1PtJetsIn[ij] = 0;
137     fh1PtJetsInRej[ij] = 0;
138     fh1PtTracksIn[ij] = 0;
139     fh1PtTracksInLow[ij] = 0;
140     fh2NJetsPt[ij]  = 0;
141     fh2NTracksPt[ij]  = 0;
142     fp2MultRPPhiTrackPt[ij] = 0;
143     fp2CentRPPhiTrackPt[ij] = 0;
144     fhnJetPt[ij] = 0;
145     fhnJetPtQA[ij] = 0;
146     fhnTrackPt[ij] = 0;
147     fhnTrackPtQA[ij] = 0;
148     for(int i = 0;i <= kMaxJets;++i){
149       fh2LTrackPtJetPt[ij][i] = 0;
150       fh1PtIn[ij][i] = 0;
151     }
152
153     fh1DijetMinv[ij] = 0;      
154     fh2DijetDeltaPhiPt[ij] = 0;  
155     fh2DijetAsymPt[ij] = 0; 
156     fh2DijetPt2vsPt1[ij] = 0;
157     fh2DijetDifvsSum[ij] = 0;
158
159   }
160 }
161
162 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
163   AliAnalysisTaskSE(name),
164   fJetHeaderRec(0x0),
165   fJetHeaderGen(0x0),
166   fAODIn(0x0),
167   fAODOut(0x0),
168   fAODExtension(0x0),
169   fhnJetContainer(0x0),
170   fhnCorrelation(0x0),
171   fhnEvent(0x0),
172   f1PtScale(0x0),
173   fBranchRec("jets"),
174   fBranchGen(""),
175   fBranchBkgRec(""),
176   fBranchBkgGen(""),
177   fNonStdFile(""),
178   fRandomizer(0x0),
179   fUseAODJetInput(kFALSE),
180   fUseAODTrackInput(kFALSE),
181   fUseAODMCInput(kFALSE),
182   fUseGlobalSelection(kFALSE),
183   fUseExternalWeightOnly(kFALSE),
184   fLimitGenJetEta(kFALSE),
185   fDoMatching(kFALSE),
186   fNMatchJets(5),
187   fNRPBins(3),
188   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
189   fFilterMask(0),
190   fEventSelectionMask(0),
191   fAnalysisType(0),
192   fTrackTypeRec(kTrackUndef),
193   fTrackTypeGen(kTrackUndef),
194   fEventClass(0),
195   fRPMethod(0),
196   fAvgTrials(1),
197   fExternalWeight(1),    
198   fJetRecEtaWindow(0.5),
199   fTrackRecEtaWindow(0.5),
200   fMinJetPt(0),
201   fMinTrackPt(0.15),
202   fDeltaPhiWindow(90./180.*TMath::Pi()),
203   fCentrality(100),
204   fRPAngle(0),
205   fMultRec(0),
206   fMultGen(0),
207   fh1Xsec(0x0),
208   fh1Trials(0x0),
209   fh1PtHard(0x0),
210   fh1PtHardNoW(0x0),  
211   fh1PtHardTrials(0x0),
212   fh1ZVtx(0x0),
213   fh1RP(0x0),
214   fh1Centrality(0x0),
215   fh1TmpRho(0x0),
216   fh2MultRec(0x0),
217   fh2MultGen(0x0),
218   fh2RPCentrality(0x0),
219   fh2PtFGen(0x0),
220   fh2RelPtFGen(0x0),
221   fHistList(0x0)
222 {
223
224   for(int ij = 0;ij <kJetTypes;++ij){    
225     fFlagJetType[ij] = 1; // default = on
226     fh1NJets[ij] = 0;
227     fh1SumPtTrack[ij] = 0;
228     fh1PtJetsIn[ij] = 0;
229     fh1PtJetsInRej[ij] = 0;
230     fh1PtTracksIn[ij] = 0;
231     fh1PtTracksInLow[ij] = 0;
232     fh2NJetsPt[ij]  = 0;
233     fh2NTracksPt[ij]  = 0;
234     fp2MultRPPhiTrackPt[ij] = 0;
235     fp2CentRPPhiTrackPt[ij] = 0;
236     fhnJetPt[ij] = 0;
237     fhnJetPtQA[ij] = 0;
238     fhnTrackPt[ij] = 0;
239     fhnTrackPtQA[ij] = 0;
240     for(int i = 0;i <= kMaxJets;++i){
241       fh2LTrackPtJetPt[ij][i] = 0;
242       fh1PtIn[ij][i] = 0;
243     }
244
245     fh1DijetMinv[ij] = 0;      
246     fh2DijetDeltaPhiPt[ij] = 0;  
247     fh2DijetAsymPt[ij] = 0; 
248     fh2DijetPt2vsPt1[ij] = 0;
249     fh2DijetDifvsSum[ij] = 0;
250   } 
251
252   DefineOutput(1, TList::Class());  
253 }
254
255
256
257 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
258 {
259
260
261
262   //
263   // Implemented Notify() to read the cross sections
264   // and number of trials from pyxsec.root
265   // 
266   
267   // Fetch the aod also from the input in,
268   // have todo it in notify
269   
270   
271   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
272   //  assume that the AOD is in the general output...
273   fAODOut  = AODEvent();
274   
275   if(fNonStdFile.Length()!=0){
276     // case that we have an AOD extension we need can fetch the jets from the extended output
277     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
278     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);    
279     if(!fAODExtension){
280       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
281     }
282   }
283
284
285   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
286   Float_t xsection = 0;
287   Float_t ftrials  = 1;
288
289   fAvgTrials = 1;
290   if(tree){
291     TFile *curfile = tree->GetCurrentFile();
292     if (!curfile) {
293       Error("Notify","No current file");
294       return kFALSE;
295     }
296     if(!fh1Xsec||!fh1Trials){
297       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
298       return kFALSE;
299     }
300     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
301     fh1Xsec->Fill("<#sigma>",xsection);
302     // construct a poor man average trials 
303     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
304     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
305   }  
306
307   if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
308
309   return kTRUE;
310 }
311
312 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
313 {
314
315   
316   // Connect the AOD
317
318   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
319   OpenFile(1);
320   if(!fHistList)fHistList = new TList(); 
321   PostData(1, fHistList); // post data in any case once
322
323   if(!fRandomizer)fRandomizer = new TRandom3(0);
324
325   fHistList->SetOwner(kTRUE);
326   Bool_t oldStatus = TH1::AddDirectoryStatus(); 
327   TH1::AddDirectory(kFALSE);
328
329
330   
331   // event npsparse cent, mult
332   const Int_t nBinsSparse0 = 2;
333   const Int_t nBins0[nBinsSparse0] = {     100, 500};
334   const Double_t xmin0[nBinsSparse0]  = {    0,   0};
335   const Double_t xmax0[nBinsSparse0]  = {  100,5000};
336       
337
338   fhnEvent = new THnSparseF("fhnEvent",";cent;mult",nBinsSparse0,
339                             nBins0,xmin0,xmax0);
340   fHistList->Add(fhnEvent);
341
342   if(fDoMatching){
343     MakeJetContainer();
344     fHistList->Add(fhnCorrelation);
345     fHistList->Add(fhnJetContainer);
346   }
347   //
348   //  Histogram
349     
350
351
352   const Int_t nBinPt = 100;
353   Double_t binLimitsPt[nBinPt+1];
354   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
355     if(iPt == 0){
356       binLimitsPt[iPt] = 0.0;
357     }
358     else {// 1.0
359       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
360     }
361   }
362   const Int_t nBinPhi = 90;
363   Double_t binLimitsPhi[nBinPhi+1];
364   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
365     if(iPhi==0){
366       binLimitsPhi[iPhi] = 0;
367     }
368     else{
369       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
370     }
371   }
372
373   const Int_t nBinEta = 40;
374   Double_t binLimitsEta[nBinEta+1];
375   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
376     if(iEta==0){
377       binLimitsEta[iEta] = -2.0;
378     }
379     else{
380       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
381     }
382   }
383
384
385   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
386   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
387   fHistList->Add(fh1Xsec);
388   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
389   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
390   fHistList->Add(fh1Trials);
391   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
392   fHistList->Add(fh1PtHard);
393   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
394   fHistList->Add(fh1PtHardNoW);
395   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
396   fHistList->Add(fh1PtHardTrials);
397   
398   fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
399   fHistList->Add(fh1ZVtx);
400
401
402   fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
403   fHistList->Add(fh1RP);
404
405   fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",103,-1,102);
406   fHistList->Add(fh1Centrality);
407
408   fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",500,0,5000,500,0.,5000);
409   fHistList->Add(fh2MultRec);
410   fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,0,5000,500,0.,5000);
411   fHistList->Add(fh2MultGen);
412
413
414   fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
415   fHistList->Add(fh2RPCentrality);
416
417   fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
418   fHistList->Add(fh2PtFGen);
419
420   fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
421   fHistList->Add(fh2RelPtFGen);
422   
423   for(int ij = 0;ij <kJetTypes;++ij){    
424     TString cAdd = "";
425     TString cJetBranch = "";
426     if(ij==kJetRec){
427       cAdd = "Rec";
428       cJetBranch = fBranchRec.Data();
429     }
430     else if (ij==kJetGen){
431       cAdd = "Gen";
432       cJetBranch = fBranchGen.Data();
433     }
434     else if (ij==kJetRecFull){
435       cAdd = "RecFull";
436       cJetBranch = fBranchRec.Data();
437     }
438     else if (ij==kJetGenFull){
439       cAdd = "GenFull";
440       cJetBranch = fBranchGen.Data();
441     }
442
443     if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
444     if(!fFlagJetType[ij])continue;
445
446     fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
447     fHistList->Add(fh1NJets[ij]);
448     
449     fh1PtJetsIn[ij]  = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
450     fHistList->Add(fh1PtJetsIn[ij]);
451
452     fh1PtJetsInRej[ij]  = new TH1F(Form("fh1PtJets%sInRej",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
453     fHistList->Add(fh1PtJetsInRej[ij]);
454     
455     fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
456     fHistList->Add(fh1PtTracksIn[ij]);
457     
458     fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),100,0.,5.);
459     fHistList->Add(fh1PtTracksInLow[ij]);
460     
461     fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),1000,0.,3000.);
462     fHistList->Add(fh1SumPtTrack[ij]);
463     
464     fh2NJetsPt[ij]  = new TH2F(Form("fh2N%sJetsPt",cAdd.Data()),Form("Number of %s jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",cAdd.Data()),nBinPt,binLimitsPt,50,-0.5,49.5);
465     fHistList->Add(fh2NJetsPt[ij]);
466     
467     fh2NTracksPt[ij]  = new TH2F(Form("fh2N%sTracksPt",cAdd.Data()),Form("Number of %s tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",cAdd.Data()),nBinPt,binLimitsPt,1000,0.,5000);
468     fHistList->Add(fh2NTracksPt[ij]);
469
470     
471     fp2MultRPPhiTrackPt[ij] = new TProfile2D(Form("fp2MultRPPhiTrackPt%s",cAdd.Data()),"RP phi vs Number of tracks;# tracks;#Delta#phi_{RP}; <p_{T}>",20,0,4000,181,-1./180.*TMath::Pi(),TMath::Pi(),"S");
472     fHistList->Add(fp2MultRPPhiTrackPt[ij]);
473     fp2CentRPPhiTrackPt[ij] = new TProfile2D(Form("fp2CentRPPhiTrackPt%s",cAdd.Data()),"RP phi vs cent;# cent;#Delta#phi_{RP}; <p_{T}>",10,0,100,181,-1./180.*TMath::Pi(),TMath::Pi(),"S");
474     fHistList->Add(fp2CentRPPhiTrackPt[ij]);    
475
476     // Bins:  Jet number: pTJet, cent, mult, RP, Area.   total bins = 4.5M
477     const Int_t nBinsSparse1 = 6;
478     const Int_t nBins1[nBinsSparse1] = {     kMaxJets+1,100, 10,  25,    fNRPBins, 10};
479     const Double_t xmin1[nBinsSparse1]  = {        -0.5,  0,  0,   0,        -0.5, 0.};
480     const Double_t xmax1[nBinsSparse1]  = {kMaxJets+0.5,250,100,5000,fNRPBins-0.5,1.0};
481     
482     fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area",nBinsSparse1,nBins1,xmin1,xmax1);
483     fHistList->Add(fhnJetPt[ij]);
484     
485     // Bins:  Jet number: pTJet, cent, eta, phi, Area.   total bins = 9.72 M
486     const Int_t nBinsSparse2 = 6;
487       const Int_t nBins2[nBinsSparse2] = {     kMaxJets+1, 25,   5,  18,             360, 10};
488       const Double_t xmin2[nBinsSparse2]  = {        -0.5,  0,   0,-0.9,              0,  0.};
489       const Double_t xmax2[nBinsSparse2]  = {kMaxJets+0.5,250, 100, 0.9, 2.*TMath::Pi(),1.0};
490       fhnJetPtQA[ij] = new THnSparseF(Form("fhnJetPtQA%s",cAdd.Data()),";jet number;p_{T,jet};cent;#eta;#phi;area",nBinsSparse2,nBins2,xmin2,xmax2);
491       fHistList->Add(fhnJetPtQA[ij]);
492
493       // Bins:track number  pTtrack, cent, mult, RP.   total bins = 224 k
494       const Int_t nBinsSparse3 = 5;
495       const Int_t nBins3[nBinsSparse3] = {       2,    100,     10,  20,    fNRPBins};
496       const Double_t xmin3[nBinsSparse3]  = { -0.5,     0,   0,      0,        -0.5};
497       const Double_t xmax3[nBinsSparse3]  = { 1.5,    200, 100,   4000,fNRPBins-0.5};  
498
499       // change the binning ot the pT axis:
500       Double_t *xPt3 = new Double_t[nBins3[1]+1];
501       xPt3[0] = 0.;
502       for(int i = 1; i<=nBins3[1];i++){
503         if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.05; // 1 - 40
504         else if(xPt3[i-1]<4)xPt3[i] = xPt3[i-1] + 0.2; // 41 - 50
505         else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 0.5; // 50 - 62
506         else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] +  1.; // 62 - 72
507         else if(xPt3[i-1]<30)xPt3[i] = xPt3[i-1] + 2.5; // 74 - 78
508         else xPt3[i] = xPt3[i-1] + 5.; // 78 - 100 = 140 
509       }
510       
511       fhnTrackPt[ij] = new THnSparseF(Form("fhnTrackPt%s",cAdd.Data()),";track number;p_{T};cent;#tracks;RP",nBinsSparse3,nBins3,xmin3,xmax3);
512       fhnTrackPt[ij]->SetBinEdges(1,xPt3);
513       fHistList->Add(fhnTrackPt[ij]);
514       delete [] xPt3;
515
516       // Track QA bins track nr, pTrack, cent, eta, phi bins 5.4 M
517       const Int_t nBinsSparse4 = 5;
518       const Int_t nBins4[nBinsSparse4] =    {    2, 50,  10,  20, 360};
519       const Double_t xmin4[nBinsSparse4]  = { -0.5,  0,   0, -1.0,   0.};
520       const Double_t xmax4[nBinsSparse4]  = {  1.5,150, 100,  1.0,2.*TMath::Pi()};  
521
522       // change the binning ot the pT axis:
523       Double_t *xPt4 = new Double_t[nBins4[1]+1];
524       xPt4[0] = 0.;
525       for(int i = 1; i<=nBins4[1];i++){
526         if(xPt4[i-1]<2)xPt4[i] = xPt4[i-1] + 0.1;
527         else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 0.5;
528         else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] +  1.;
529         else if(xPt4[i-1]<30)xPt4[i] = xPt4[i-1] +  2.5;
530         else xPt4[i] = xPt4[i-1] + 5.;
531       }
532       fhnTrackPtQA[ij] = new THnSparseF(Form("fhnTrackPtQA%s",cAdd.Data()),";track number;p_{T};cent;#eta;#phi",nBinsSparse4,nBins4,xmin4,xmax4);
533       fhnTrackPtQA[ij]->SetBinEdges(1,xPt4);
534       fHistList->Add(fhnTrackPtQA[ij]);
535       delete [] xPt4;
536
537       for(int i = 0;i <= kMaxJets;++i){
538         fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
539         fHistList->Add(fh1PtIn[ij][i]);
540
541
542         if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
543         fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
544                                            Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
545                                                 cAdd.Data()),
546                                            200,0.,200.,nBinPt,binLimitsPt);
547         fHistList->Add(fh2LTrackPtJetPt[ij][i]);
548       }
549
550
551       fh1DijetMinv[ij]                = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
552       fHistList->Add(fh1DijetMinv[ij]);
553
554       fh2DijetDeltaPhiPt[ij]       = new TH2F(Form("fh2Dijet%sDeltaPhiPt",cAdd.Data()),"Difference in the azimuthal angle;#Delta#phi;p_{T,2};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
555       fHistList->Add(fh2DijetDeltaPhiPt[ij]);
556
557       fh2DijetAsymPt[ij]            = new TH2F(Form("fh2Dijet%sAsym",cAdd.Data()),"Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
558       fHistList->Add(fh2DijetAsymPt[ij]);
559
560       fh2DijetPt2vsPt1[ij]          = new TH2F(Form("fh2Dijet%sPt2vsPt1",cAdd.Data()),"Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
561       fHistList->Add(fh2DijetPt2vsPt1[ij]);
562       fh2DijetDifvsSum[ij]         = new TH2F(Form("fh2Dijet%sDifvsSum",cAdd.Data()),"Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
563       fHistList->Add( fh2DijetDifvsSum[ij]);
564     }   
565   // =========== Switch on Sumw2 for all histos ===========
566   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
567     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
568     if (h1){
569       h1->Sumw2();
570       continue;
571     }
572     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
573     if(hn)hn->Sumw2();
574   }
575   TH1::AddDirectory(oldStatus);
576 }
577
578 void AliAnalysisTaskJetSpectrum2::Init()
579 {
580   //
581   // Initialization
582   //
583
584   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
585
586 }
587
588 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
589
590
591   Bool_t selected = kTRUE;
592
593   if(fUseGlobalSelection&&fEventSelectionMask==0){
594     selected = AliAnalysisHelperJetTasks::Selected();
595   }
596   else if(fUseGlobalSelection&&fEventSelectionMask>0){
597     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
598   }
599
600   if(fEventClass>0){
601     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
602   }
603
604   if(!selected){
605     // no selection by the service task, we continue
606     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
607     PostData(1, fHistList);
608     return;
609   }
610
611
612   static AliAODEvent* aod = 0;
613
614   // take all other information from the aod we take the tracks from
615   if(!aod){
616    if(fUseAODTrackInput)aod = fAODIn;
617    else aod = fAODOut;
618   }
619
620
621   //
622   // Execute analysis for current event
623   //
624   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);  
625   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
626   if(!aodH){
627     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
628     return;
629   }
630
631   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
632   TClonesArray *aodRecJets = 0;
633
634   if(fAODOut&&!aodRecJets){
635     aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
636   }
637   if(fAODExtension&&!aodRecJets){
638     aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
639   }
640   if(fAODIn&&!aodRecJets){
641     aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
642   }
643
644
645
646   if(!aodRecJets){
647     if(fDebug){
648
649         Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
650         if(fAODIn){
651           Printf("Input AOD >>>>");
652           fAODIn->Print();    
653         }
654         if(fAODExtension){
655           Printf("AOD Extension >>>>");
656           fAODExtension->Print();
657         }
658         if(fAODOut){
659           Printf("Output AOD >>>>");
660           fAODOut->Print();    
661         }
662     }
663         return;
664   }
665
666   TClonesArray *aodGenJets = 0;
667   if(fBranchGen.Length()>0){
668     if(fAODOut&&!aodGenJets){
669       aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
670     }
671     if(fAODExtension&&!aodGenJets){
672       aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
673     }
674     if(fAODIn&&!aodGenJets){
675       aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
676     }
677
678     if(!aodGenJets){
679       Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
680       return;
681     }
682   }
683
684   TClonesArray *aodBackRecJets = 0;
685   if(fBranchBkgRec.Length()>0){
686     if(fAODOut&&!aodBackRecJets){
687       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
688     }
689     if(fAODExtension&&!aodBackRecJets){
690       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
691     }
692     if(fAODIn&&!aodBackRecJets){
693       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
694     }
695
696     if(!aodBackRecJets){
697       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
698       return;
699     }
700   }
701
702
703   TClonesArray *aodBackGenJets = 0;
704
705   if(fBranchBkgGen.Length()>0){
706     if(fAODOut&&!aodBackGenJets){
707       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
708     }
709     if(fAODExtension&&!aodBackGenJets){
710       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
711     }
712     if(fAODIn&&!aodBackGenJets){
713       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
714     }
715
716     if(!aodBackGenJets){
717       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
718       return;
719     }
720   }
721
722  
723   // new Scheme
724   // first fill all the pure  histograms, i.e. full jets 
725   // in case of the correaltion limit to the n-leading jets
726
727   // reconstructed
728
729   
730   // generated
731
732
733   // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
734
735
736
737   TList genJetsList;         // full acceptance
738   TList genJetsListCut;      // acceptance cut
739   TList recJetsList;         // full acceptance
740   TList recJetsListCut;      // acceptance cut
741
742
743   GetListOfJets(&recJetsList,aodRecJets,0);
744   GetListOfJets(&recJetsListCut,aodRecJets,1);
745
746   if(fBranchGen.Length()>0){
747     GetListOfJets(&genJetsList,aodGenJets,0);
748     GetListOfJets(&genJetsListCut,aodGenJets,1);
749   }
750
751   Double_t eventW = 1;
752   Double_t ptHard = 0; 
753   Double_t nTrials = 1; // Trials for MC trigger 
754   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
755
756   // Getting some global properties
757   fCentrality = GetCentrality();
758   if(fCentrality<=0)fCentrality = 0;
759   fh1Centrality->Fill(fCentrality);
760
761
762   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
763     // this is the part we only use when we have MC information
764     AliMCEvent* mcEvent = MCEvent();
765     //    AliStack *pStack = 0; 
766     if(!mcEvent){
767       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
768       return;
769     }
770     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
771     if(pythiaGenHeader){
772       nTrials = pythiaGenHeader->Trials();
773       ptHard  = pythiaGenHeader->GetPtHard();
774       int iProcessType = pythiaGenHeader->ProcessType();
775       // 11 f+f -> f+f
776       // 12 f+barf -> f+barf
777       // 13 f+barf -> g+g
778       // 28 f+g -> f+g
779       // 53 g+g -> f+barf
780       // 68 g+g -> g+g
781       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
782       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
783     } 
784   }// (fAnalysisType&kMCESD)==kMCESD)
785
786
787   // we simply fetch the tracks/mc particles as a list of AliVParticles
788
789   TList recParticles;
790   TList genParticles;
791
792   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
793   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
794   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
795   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
796
797   //  CalculateReactionPlaneAngle(&recParticles);
798   fRPAngle = 0;
799   
800   if(fRPMethod==0)fRPAngle = aod->GetHeader()->GetEventplane();
801   else if(fRPMethod==1||fRPMethod==2){
802     fRPAngle = aod->GetHeader()->GetQTheta(fRPMethod);
803   }
804   fh1RP->Fill(fRPAngle);
805   fh2RPCentrality->Fill(fCentrality,fRPAngle);
806   // Event control and counting ...  
807   // MC
808   fh1PtHard->Fill(ptHard,eventW);
809   fh1PtHardNoW->Fill(ptHard,1);
810   fh1PtHardTrials->Fill(ptHard,nTrials);
811
812   // Real
813   if(aod->GetPrimaryVertex()){// No vtx for pure MC
814     fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
815   }
816
817
818   Int_t recMult1 = recParticles.GetEntries();
819   Int_t genMult1 = genParticles.GetEntries();
820
821   Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
822   Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
823
824   fh2MultRec->Fill(recMult1,recMult2);
825   fh2MultGen->Fill(genMult1,genMult2);
826   fMultRec = recMult1;
827   if(fMultRec<=0)fMultRec = recMult2;
828   fMultGen = genMult1;
829   if(fMultGen<=0)fMultGen = genMult2;
830
831   Double_t var0[2] = {0,};
832   var0[0] = fCentrality;
833   var0[1] = fMultRec;
834   fhnEvent->Fill(var0);
835
836   // the loops for rec and gen should be indentical... pass it to a separate
837   // function ...
838   // Jet Loop
839   // Track Loop
840   // Jet Jet Loop
841   // Jet Track loop
842
843   FillJetHistos(recJetsListCut,recParticles,kJetRec);
844   FillJetHistos(recJetsList,recParticles,kJetRecFull);
845   FillTrackHistos(recParticles,kJetRec);
846
847   FillJetHistos(genJetsListCut,genParticles,kJetGen);
848   FillJetHistos(genJetsList,genParticles,kJetGenFull);
849   FillTrackHistos(genParticles,kJetGen);
850
851   // Here follows the jet matching part
852   // delegate to a separated method?
853
854   if(fDoMatching){
855     FillMatchHistos(recJetsList,genJetsList);
856   }
857
858   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
859   PostData(1, fHistList);
860 }
861
862 void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
863
864   if(iType>=kJetTypes){
865     return;
866   }
867   if(!fFlagJetType[iType])return;
868
869   Int_t refMult = fMultRec;
870   if(iType==kJetGen||iType==kJetGenFull){
871     refMult = fMultGen;
872   }
873
874   Int_t nJets = jetsList.GetEntries(); 
875   fh1NJets[iType]->Fill(nJets);
876
877   if(nJets<=0)return;
878   
879   Float_t ptOld = 1.E+32;
880   Float_t pT0 = 0;
881   Float_t pT1 = 0;
882   Float_t phi0 = 0;
883   Float_t phi1 = 0;
884   Int_t ij0 = -1;
885   Int_t ij1 = -1;
886
887   Double_t var1[6] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area
888   var1[2] = fCentrality; 
889   var1[3] = refMult;
890
891   Double_t var2[6] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area
892   var2[2] = fCentrality;
893
894   for(int ij = 0;ij < nJets;ij++){
895     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
896     Float_t ptJet = jet->Pt();
897     if(jet->Trigger()&fJetTriggerExcludeMask){
898       fh1PtJetsInRej[iType]->Fill(ptJet);
899       continue;
900     }
901     fh1PtJetsIn[iType]->Fill(ptJet);
902     if(ptJet>ptOld){
903       Printf("%s:%d Jets Type %d Not Sorted !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,iType,ij,ptJet,ij-1,ptOld);
904     }
905     ptOld = ptJet;
906     
907     // find the dijets assume sorting and acceptance cut...
908     if(ij==0){
909       ij0 = ij;
910       pT0 = ptJet;
911       phi0 = jet->Phi();
912       if(phi0<0)phi0+=TMath::Pi()*2.;
913     }
914     else if(ptJet>pT1){
915       // take only the backward hemisphere??                                                        
916       phi1 = jet->Phi();
917       if(phi1<0)phi1+=TMath::Pi()*2.;
918       Float_t dPhi = phi1 - phi0;
919       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
920       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
921       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
922         ij1 = ij;
923         pT1 = ptJet;
924       }
925     }
926     // fill jet histos for kmax jets
927
928       Float_t phiJet = jet->Phi();
929       Float_t etaJet = jet->Eta();
930       if(phiJet<0)phiJet+=TMath::Pi()*2.;    
931       fh1TmpRho->Reset();
932       if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
933
934       fh1PtIn[iType][kMaxJets]->Fill(ptJet);
935       // fill leading jets...
936       AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
937       //      AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
938       Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
939
940       var1[1] = ptJet;
941       var1[4] = phiBin;
942       var1[5] = jet->EffectiveAreaCharged();
943
944       var2[1] = ptJet;
945       var2[3] = etaJet;
946       var2[4] = phiJet;
947       var2[5] = jet->EffectiveAreaCharged();
948       if(ij<kMaxJets){
949         if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
950         var1[0] = ij;
951         var2[0] = ij;
952         fhnJetPt[iType]->Fill(var1);
953         fhnJetPtQA[iType]->Fill(var2);
954       }
955       var1[0] = kMaxJets;// fill for all jets
956       var2[0] = kMaxJets;// fill for all jets
957       fhnJetPt[iType]->Fill(var1);
958       fhnJetPtQA[iType]->Fill(var2);
959       if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
960
961       if(particlesList.GetSize()&&ij<kMaxJets){
962         // Particles... correlated with jets...
963         for(int it = 0;it<particlesList.GetEntries();++it){
964           AliVParticle *part = (AliVParticle*)particlesList.At(it);
965           Float_t deltaR = jet->DeltaR(part);
966           if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
967         }
968         // fill the jet shapes
969       }// if we have particles
970   }// Jet Loop
971
972
973   // do something with dijets...
974   if(ij0>=0&&ij1>0){
975     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
976     Double_t ptJet0 = jet0->Pt();
977     Double_t phiJet0 = jet0->Phi();
978     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;       
979
980     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
981     Double_t ptJet1 = jet1->Pt();
982     Double_t phiJet1 = jet1->Phi();
983     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;       
984
985     Float_t deltaPhi = phiJet0 - phiJet1;
986     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
987     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();      
988     deltaPhi = TMath::Abs(deltaPhi);
989     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);      
990
991     Float_t asym = 9999;
992     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
993       fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
994       fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);        
995       fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);        
996       Float_t minv = 2.*(jet0->P()*jet1->P()-
997                          jet0->Px()*jet1->Px()- 
998                          jet0->Py()*jet1->Py()- 
999                          jet0->Pz()*jet1->Pz());    // assume mass == 0;
1000       if(minv<0)minv=0; // prevent numerical instabilities
1001       minv = TMath::Sqrt(minv);
1002       fh1DijetMinv[iType]->Fill(minv);            
1003   }
1004   
1005
1006
1007   // count down the jets above thrueshold
1008   Int_t nOver = nJets;
1009   if(nOver>0){
1010     TIterator *jetIter = jetsList.MakeIterator();
1011     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());  
1012     if(tmpJet){
1013       Float_t pt = tmpJet->Pt();
1014       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1015         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1016         while(pt<ptCut&&tmpJet){
1017           nOver--;
1018           tmpJet = (AliAODJet*)(jetIter->Next()); 
1019           if(tmpJet){
1020             pt = tmpJet->Pt();
1021           }
1022         }
1023         if(nOver<=0)break;
1024         fh2NJetsPt[iType]->Fill(ptCut,nOver);
1025       }
1026     }
1027     delete jetIter;
1028   }
1029 }
1030
1031 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1032
1033   if(fFlagJetType[iType]<=0)return;
1034   Int_t refMult = fMultRec;
1035   if(iType==kJetGen||iType==kJetGenFull){
1036     refMult = fMultGen;
1037
1038   }
1039
1040   // 
1041   Double_t var3[5]; // track number;p_{T};cent;#tracks;RP
1042   var3[2] = fCentrality;
1043   var3[3] = refMult;
1044   Double_t var4[5]; // track number;p_{T};cent;#eta;#phi
1045   var4[2] = fCentrality;
1046   Int_t nTrackOver = particlesList.GetSize();
1047   // do the same for tracks and jets
1048   if(nTrackOver>0){
1049     TIterator *trackIter = particlesList.MakeIterator();
1050     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());  
1051     Float_t pt = tmpTrack->Pt();
1052     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1053       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1054       while(pt<ptCut&&tmpTrack){
1055         nTrackOver--;
1056         tmpTrack = (AliVParticle*)(trackIter->Next()); 
1057         if(tmpTrack){
1058           pt = tmpTrack->Pt();
1059         }
1060       }
1061       if(nTrackOver<=0)break;
1062       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1063     }
1064
1065     
1066     trackIter->Reset();
1067     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1068     Float_t sumPt = 0;
1069
1070     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1071       Float_t tmpPt = tmpTrack->Pt();
1072       fh1PtTracksIn[iType]->Fill(tmpPt);
1073       fh1PtTracksInLow[iType]->Fill(tmpPt);
1074
1075       sumPt += tmpPt;
1076       Float_t tmpPhi = tmpTrack->Phi();
1077       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1078
1079       
1080       Float_t phiRP = tmpPhi-fRPAngle;
1081       if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1082       if(phiRP<0)phiRP += TMath::Pi();
1083       if(phiRP<0)phiRP += TMath::Pi();
1084       const float allPhi = -1./180.*TMath::Pi();
1085
1086       fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1087       fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1088
1089       fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1090       fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1091
1092       Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1093       var3[0] = 1;
1094       var3[1] = tmpPt;
1095       var3[4] = phiBin;
1096       
1097       var4[0] = 1;
1098       var4[1] = tmpPt;
1099       var4[3] = tmpTrack->Eta();
1100       var4[4] = tmpPhi;
1101
1102
1103       fhnTrackPt[iType]->Fill(var3);
1104       fhnTrackPtQA[iType]->Fill(var4);
1105
1106       if(tmpTrack==leading){
1107         var3[0] = 0;
1108         var4[0] = 0;
1109         fhnTrackPt[iType]->Fill(var3);
1110         fhnTrackPtQA[iType]->Fill(var4);
1111         continue;
1112       }
1113     }  
1114     fh1SumPtTrack[iType]->Fill(sumPt);
1115     delete trackIter;
1116   }
1117
1118 }
1119
1120
1121 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1122
1123
1124   // Fill al the matching histos
1125   // It is important that the acceptances for the mathing are as large as possible
1126   // to avoid false matches on the edge of acceptance
1127   // therefore we add some extra matching jets as overhead
1128
1129   static TArrayI aGenIndex(recJetsList.GetEntries());
1130   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1131
1132   static TArrayI aRecIndex(genJetsList.GetEntries());
1133   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1134
1135   if(fDebug){
1136     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1137     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1138   }
1139   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1140                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1141                                             aGenIndex,aRecIndex,fDebug);
1142
1143   if(fDebug){
1144     for(int i = 0;i< aGenIndex.GetSize();++i){ 
1145       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]); 
1146     }
1147     for(int i = 0;i< aRecIndex.GetSize();++i){
1148       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]); 
1149     }
1150   }
1151
1152   Double_t container[6];
1153
1154   // loop over generated jets
1155   // consider the 
1156   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1157     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1158     Double_t ptGen = genJet->Pt();
1159     Double_t phiGen = genJet->Phi();
1160     if(phiGen<0)phiGen+=TMath::Pi()*2.;    
1161     Double_t etaGen = genJet->Eta();
1162     container[3] = ptGen;
1163     container[4] = etaGen;
1164     container[5] = phiGen;
1165     fhnJetContainer->Fill(&container[3],kStep0);
1166     if(JetSelected(genJet)){
1167       fhnJetContainer->Fill(&container[3],kStep1);
1168       Int_t ir = aRecIndex[ig];
1169       if(ir>=0&&ir<recJetsList.GetEntries()){   
1170         fhnJetContainer->Fill(&container[3],kStep2);
1171         AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir); 
1172         if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep4);
1173         if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep3);
1174       }
1175     }
1176   }// loop over generated jets used for matching...
1177
1178
1179
1180   // loop over reconstructed jets
1181   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1182     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1183     Double_t etaRec = recJet->Eta();
1184     Double_t ptRec = recJet->Pt();
1185     Double_t phiRec = recJet->Phi();
1186     if(phiRec<0)phiRec+=TMath::Pi()*2.;    
1187     // do something with dijets...
1188     
1189     container[0] = ptRec;
1190     container[1] = etaRec;
1191     container[2] = phiRec;
1192
1193     fhnJetContainer->Fill(container,kStep0+kMaxStep);
1194     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1195   
1196     if(JetSelected(recJet)){
1197       fhnJetContainer->Fill(container,kStep1+kMaxStep);
1198       // Fill Correlation
1199       Int_t ig = aGenIndex[ir];
1200       if(ig>=0 && ig<genJetsList.GetEntries()){
1201         fhnJetContainer->Fill(container,kStep2+kMaxStep);
1202         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1203         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1204         Double_t ptGen  = genJet->Pt();
1205         Double_t phiGen = genJet->Phi();
1206         if(phiGen<0)phiGen+=TMath::Pi()*2.; 
1207         Double_t etaGen = genJet->Eta();
1208       
1209         container[3] = ptGen;
1210         container[4] = etaGen;
1211         container[5] = phiGen;
1212         // 
1213         // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1214         // 
1215         if(JetSelected(genJet))fhnJetContainer->Fill(container,kStep4+kMaxStep);
1216         fhnJetContainer->Fill(container,kStep3+kMaxStep);
1217         fhnCorrelation->Fill(container,0);
1218         if(ptGen>0){
1219           Float_t delta = (ptRec-ptGen)/ptGen;
1220           fh2RelPtFGen->Fill(ptGen,delta);
1221           fh2PtFGen->Fill(ptGen,ptRec);
1222         }
1223       } 
1224     }// loop over reconstructed jets
1225   }
1226   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1227 }
1228
1229
1230 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1231   //
1232   // Create the particle container for the correction framework manager and 
1233   // link it
1234   //
1235   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
1236   const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
1237   const Double_t kEtamin = -3.0, kEtamax = 3.0;
1238   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1239
1240   // can we neglect migration in eta and phi?
1241   // phi should be no problem since we cover full phi and are phi symmetric
1242   // eta migration is more difficult  due to needed acceptance correction
1243   // in limited eta range
1244
1245   //arrays for the number of bins in each dimension
1246   Int_t iBin[kNvar];
1247   iBin[0] = 125; //bins in pt
1248   iBin[1] =  1; //bins in eta 
1249   iBin[2] = 1; // bins in phi
1250
1251   //arrays for lower bounds :
1252   Double_t* binEdges[kNvar];
1253   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1254     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1255
1256   //values for bin lower bounds
1257   //  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);  
1258   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1259   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1260   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1261
1262
1263   fhnJetContainer = new AliTHn(Form("fahnJetContainer"),Form("AliTHn jet info"),kMaxStep*2,kNvar,iBin);
1264   for (int k=0; k<kNvar; k++) {
1265     fhnJetContainer->SetBinLimits(k,binEdges[k]);
1266   }
1267
1268   //create correlation matrix for unfolding
1269   Int_t thnDim[2*kNvar];
1270   for (int k=0; k<kNvar; k++) {
1271     //first half  : reconstructed 
1272     //second half : MC
1273     thnDim[k]      = iBin[k];
1274     thnDim[k+kNvar] = iBin[k];
1275   }
1276
1277   fhnCorrelation = new AliTHn("fahnCorrelation","AliTHn with correlations",1,2*kNvar,thnDim);
1278   for (int k=0; k<kNvar; k++) {
1279     fhnCorrelation->SetBinLimits(k,binEdges[k]);
1280     fhnCorrelation->SetBinLimits(k+kNvar,binEdges[k]);
1281   }
1282
1283   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1284     delete [] binEdges[ivar];
1285
1286
1287 }
1288
1289 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1290 {
1291 // Terminate analysis
1292 //
1293     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1294 }
1295
1296
1297 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1298
1299   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1300
1301   Int_t iCount = 0;
1302   if(type==kTrackAOD){
1303     AliAODEvent *aod = 0;
1304     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1305     else aod = AODEvent();
1306     if(!aod){
1307       return iCount;
1308     }
1309     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1310       AliAODTrack *tr = aod->GetTrack(it);
1311       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1312       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1313       if(tr->Pt()<fMinTrackPt)continue;
1314       if(fDebug>0){
1315         if(tr->Pt()>20){
1316           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1317           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1318           tr->Print();
1319           //    tr->Dump();
1320           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1321           if(fESD){
1322             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1323             esdTr->Print("");
1324             // esdTr->Dump();
1325           }
1326         }
1327       }
1328       list->Add(tr);
1329       iCount++;
1330     }
1331   }
1332   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1333     AliMCEvent* mcEvent = MCEvent();
1334     if(!mcEvent)return iCount;
1335     // we want to have alivpartilces so use get track
1336     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1337       if(!mcEvent->IsPhysicalPrimary(it))continue;
1338       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1339       if(part->Pt()<fMinTrackPt)continue;
1340       if(type == kTrackKineAll){
1341         list->Add(part);
1342         iCount++;
1343       }
1344       else if(type == kTrackKineCharged){
1345         if(part->Particle()->GetPDG()->Charge()==0)continue;
1346         list->Add(part);
1347         iCount++;
1348       }
1349     }
1350   }
1351   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1352     AliAODEvent *aod = 0;
1353     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1354     else aod = AODEvent();
1355     if(!aod)return iCount;
1356     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1357     if(!tca)return iCount;
1358     for(int it = 0;it < tca->GetEntriesFast();++it){
1359       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1360       if(!part)continue;
1361       if(part->Pt()<fMinTrackPt)continue;
1362       if(!part->IsPhysicalPrimary())continue;
1363       if(type == kTrackAODMCAll){
1364         list->Add(part);
1365         iCount++;
1366       }
1367       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1368         if(part->Charge()==0)continue;
1369         if(kTrackAODMCCharged){
1370           list->Add(part);
1371         }
1372         else {
1373           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1374           list->Add(part);
1375         }
1376         iCount++;
1377       }
1378     }
1379   }// AODMCparticle
1380   list->Sort();
1381   return iCount;
1382
1383 }
1384
1385
1386 Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1387     AliAODEvent *aod = 0;
1388     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1389     else aod = AODEvent();
1390     if(!aod){
1391       return 101;
1392     }
1393     return aod->GetHeader()->GetCentrality();
1394 }
1395
1396
1397
1398 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1399   Bool_t selected = false;
1400   
1401   if(!jet)return selected;
1402
1403   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1404     selected = kTRUE;
1405   }
1406   return selected;
1407
1408 }
1409
1410 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1411
1412   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1413   Int_t iCount = 0;
1414
1415   if(!jarray){
1416     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1417     return iCount;
1418   }
1419
1420
1421   for(int ij=0;ij<jarray->GetEntries();ij++){
1422     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1423     if(!jet)continue;
1424     if(type==0){
1425       // No cut at all, main purpose here is sorting      
1426       list->Add(jet);
1427       iCount++;
1428     }
1429     else if(type == 1){
1430       // eta cut
1431       if(JetSelected(jet)){
1432         list->Add(jet);
1433         iCount++;
1434       }
1435     }
1436   }
1437
1438   list->Sort();
1439   return iCount;
1440
1441 }
1442
1443
1444 Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1445   if(!jets)return 0;
1446
1447   Int_t refMult = 0;
1448   for(int ij = 0;ij < jets->GetEntries();++ij){
1449     AliAODJet* jet = (AliAODJet*)jets->At(ij);
1450     if(!jet)continue;
1451     TRefArray *refs = jet->GetRefTracks();
1452     if(!refs)continue;
1453     refMult += refs->GetEntries();
1454   }
1455   return refMult;
1456
1457 }
1458
1459
1460 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1461   if(!jet)return 0;
1462   TRefArray *refs = jet->GetRefTracks();
1463   if(!refs) return 0;
1464   AliVParticle *leading = 0;
1465   Float_t fMaxPt = 0;
1466   for(int i = 0;i<refs->GetEntries();i++){
1467     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1468     if(!tmp)continue;
1469     if(tmp->Pt()>fMaxPt){
1470       leading = tmp;
1471       fMaxPt = tmp->Pt();
1472     }
1473   }
1474   return leading;
1475 }
1476
1477
1478 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1479   if(!jet)return 0;
1480   if(!list) return 0;
1481   AliVParticle *leading = 0;
1482   Float_t fMaxPt = 0;
1483   for(int i = 0;i<list->GetEntries();i++){
1484     AliVParticle *tmp = (AliVParticle*)(list->At(i));
1485     if(!tmp)continue;
1486     if(jet->DeltaR(tmp)>r)continue;
1487     if(tmp->Pt()>fMaxPt){
1488       leading = tmp;
1489       fMaxPt = tmp->Pt();
1490     }
1491   }
1492   return leading;
1493 }
1494
1495 Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1496 {
1497     Int_t phibin=-1;
1498     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1499     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1500     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1501     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1502     return phibin;
1503 }
1504