hjet mass ana update + bug fix
[u/mrichter/AliRoot.git] / PWGJE / 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 #include <TROOT.h>
23 #include <TRandom.h>
24 #include <TRandom3.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 <TProfile2D.h>
35 #include <TList.h>
36 #include <TLorentzVector.h>
37 #include <TClonesArray.h>
38 #include <TRefArray.h>
39 #include  "TDatabasePDG.h"
40
41 #include "AliAnalysisTaskJetSpectrum2.h"
42 #include "AliAnalysisManager.h"
43 #include "AliJetFinder.h"
44 #include "AliTHn.h"
45 #include "AliJetHeader.h"
46 #include "AliJetReader.h"
47 #include "AliJetReaderHeader.h"
48 #include "AliUA1JetHeaderV1.h"
49 #include "AliESDEvent.h"
50 #include "AliAODEvent.h"
51 #include "AliAODHandler.h"
52 #include "AliAODTrack.h"
53 #include "AliAODJet.h"
54 #include "AliAODJetEventBackground.h"
55 #include "AliAODMCParticle.h"
56 #include "AliMCEventHandler.h"
57 #include "AliMCEvent.h"
58 #include "AliStack.h"
59 #include "AliGenPythiaEventHeader.h"
60 #include "AliJetKineReaderHeader.h"
61 #include "AliGenCocktailEventHeader.h"
62 #include "AliInputEventHandler.h"
63 #include "AliCFContainer.h"
64
65 #include "AliAnalysisHelperJetTasks.h"
66
67 ClassImp(AliAnalysisTaskJetSpectrum2)
68
69 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): 
70 AliAnalysisTaskSE(),
71   fJetHeaderRec(0x0),
72   fJetHeaderGen(0x0),
73   fAODIn(0x0),
74   fAODOut(0x0),
75   fAODExtension(0x0),
76   fhnJetContainer(0x0),
77   fhnCorrelation(0x0),
78   fhnEvent(0x0),
79   f1PtScale(0x0),
80   fBranchRec("jets"),
81   fBranchGen(""),
82   fBranchBkgRec(""), 
83   fBranchBkgGen(""), 
84   fNonStdFile(""),
85   fRandomizer(0x0),
86   fUseAODJetInput(kFALSE),
87   fUseAODTrackInput(kFALSE),
88   fUseAODMCInput(kFALSE),
89   fUseGlobalSelection(kFALSE),
90   fUseExternalWeightOnly(kFALSE),
91   fLimitGenJetEta(kFALSE),
92   fDoMatching(kFALSE),
93   fNMatchJets(5),
94   fNRPBins(3),
95   fTRP(0),
96   fDebug(0),
97   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
98   fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
99   fFilterMask(0),
100   fEventSelectionMask(0),
101   fNTrigger(0),
102   fTriggerBit(0x0),
103   fNAcceptance(0),
104   fNBinsLeadingTrackPt(10),
105   fNBinsMult(20),
106   fAnalysisType(0),
107   fTrackTypeRec(kTrackUndef),
108   fTrackTypeGen(kTrackUndef),
109   fEventClass(0),
110   fRPMethod(0),
111   fAvgTrials(1),
112   fExternalWeight(1),    
113   fJetRecEtaWindow(0.5),
114   fTrackRecEtaWindow(0.5),
115   fMinJetPt(0),
116   fMinTrackPt(0.15),
117   fDeltaPhiWindow(90./180.*TMath::Pi()),
118   fAcceptancePhiMin(0x0),
119   fAcceptancePhiMax(0x0),
120   fAcceptanceEtaMin(0x0),
121   fAcceptanceEtaMax(0x0),
122   fCentrality(100),
123   fRPAngle(0),
124   fMultRec(0),
125   fMultGen(0),
126   fTriggerName(0x0),
127   fh1Xsec(0x0),
128   fh1Trials(0x0),
129   fh1AvgTrials(0x0),
130   fh1PtHard(0x0),
131   fh1PtHardNoW(0x0),  
132   fh1PtHardTrials(0x0),
133   fh1ZVtx(0x0),
134   fh1RP(0x0),
135   fh1Centrality(0x0),
136   fh1TmpRho(0x0),
137   fh2MultRec(0x0),
138   fh2MultGen(0x0),
139   fh2RPCentrality(0x0),
140   fh2PtFGen(0x0),
141   fh2deltaPt1Pt2(0x0),
142   fh2RelPtFGen(0x0),
143   fh3RelPtFGenLeadTrkPt(0x0),
144   fh1EvtSelection(0),
145   fMaxVertexZ(100.),
146   fMinNcontributors(0),
147   fRejectPileup(0),
148   fHistList(0x0)  
149 {
150
151   for(int ij = 0;ij <kJetTypes;++ij){    
152     fFlagJetType[ij] = 1; // default = on
153     fh1NJets[ij] = 0;
154     fh1SumPtTrack[ij] = 0;
155     fh1PtJetsIn[ij] = 0;
156     fh1PtJetsInRej[ij] = 0;
157     fh1PtJetsInBest[ij] = 0;
158     fh1PtTracksIn[ij] = 0;
159     fh1PtTracksInLow[ij] = 0;
160     fh2NJetsPt[ij]  = 0;
161     fh2NTracksPt[ij]  = 0;
162     fp2MultRPPhiTrackPt[ij] = 0;
163     fp2CentRPPhiTrackPt[ij] = 0;
164     fhnJetPt[ij] = 0;
165     fhnJetPtBest[ij] = 0;
166     fhnJetPtRej[ij] = 0;
167     fhnJetPtQA[ij] = 0;
168     fhnTrackPt[ij] = 0;
169     fhnTrackPtQA[ij] = 0;
170     for(int i = 0;i <= kMaxJets;++i){
171       fh2LTrackPtJetPt[ij][i] = 0;
172       fh1PtIn[ij][i] = 0;
173     }
174
175     fh1DijetMinv[ij] = 0;      
176     fh2DijetDeltaPhiPt[ij] = 0;  
177     fh2DijetAsymPt[ij] = 0; 
178     fh2DijetPt2vsPt1[ij] = 0;
179     fh2DijetDifvsSum[ij] = 0;
180
181   }
182 }
183
184 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
185   AliAnalysisTaskSE(name),
186   fJetHeaderRec(0x0),
187   fJetHeaderGen(0x0),
188   fAODIn(0x0),
189   fAODOut(0x0),
190   fAODExtension(0x0),
191   fhnJetContainer(0x0),
192   fhnCorrelation(0x0),
193   fhnEvent(0x0),
194   f1PtScale(0x0),
195   fBranchRec("jets"),
196   fBranchGen(""),
197   fBranchBkgRec(""),
198   fBranchBkgGen(""),
199   fNonStdFile(""),
200   fRandomizer(0x0),
201   fUseAODJetInput(kFALSE),
202   fUseAODTrackInput(kFALSE),
203   fUseAODMCInput(kFALSE),
204   fUseGlobalSelection(kFALSE),
205   fUseExternalWeightOnly(kFALSE),
206   fLimitGenJetEta(kFALSE),
207   fDoMatching(kFALSE),
208   fNMatchJets(5),
209   fNRPBins(3),
210   fTRP(0),
211   fDebug(0),
212   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
213   fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
214   fFilterMask(0),
215   fEventSelectionMask(0),
216   fNTrigger(0),
217   fTriggerBit(0x0),
218   fNAcceptance(0),
219   fNBinsLeadingTrackPt(10),
220   fNBinsMult(20),
221   fAnalysisType(0),
222   fTrackTypeRec(kTrackUndef),
223   fTrackTypeGen(kTrackUndef),
224   fEventClass(0),
225   fRPMethod(0),
226   fAvgTrials(1),
227   fExternalWeight(1),    
228   fJetRecEtaWindow(0.5),
229   fTrackRecEtaWindow(0.5),
230   fMinJetPt(0),
231   fMinTrackPt(0.15),
232   fDeltaPhiWindow(90./180.*TMath::Pi()),
233   fAcceptancePhiMin(0x0),
234   fAcceptancePhiMax(0x0),
235   fAcceptanceEtaMin(0x0),
236   fAcceptanceEtaMax(0x0),
237   fCentrality(100),
238   fRPAngle(0),
239   fMultRec(0),
240   fMultGen(0),
241   fTriggerName(0x0),
242   fh1Xsec(0x0),
243   fh1Trials(0x0),
244   fh1AvgTrials(0x0),
245   fh1PtHard(0x0),
246   fh1PtHardNoW(0x0),  
247   fh1PtHardTrials(0x0),
248   fh1ZVtx(0x0),
249   fh1RP(0x0),
250   fh1Centrality(0x0),
251   fh1TmpRho(0x0),
252   fh2MultRec(0x0),
253   fh2MultGen(0x0),
254   fh2RPCentrality(0x0),
255   fh2PtFGen(0x0),
256   fh2deltaPt1Pt2(0x0),
257   fh2RelPtFGen(0x0),
258   fh3RelPtFGenLeadTrkPt(0x0),
259   fh1EvtSelection(0),
260   fMaxVertexZ(100.),
261   fMinNcontributors(0),
262   fRejectPileup(0),
263   fHistList(0x0)
264 {
265
266   for(int ij = 0;ij <kJetTypes;++ij){    
267     fFlagJetType[ij] = 1; // default = on
268     fh1NJets[ij] = 0;
269     fh1SumPtTrack[ij] = 0;
270     fh1PtJetsIn[ij] = 0;
271     fh1PtJetsInRej[ij] = 0;
272     fh1PtJetsInBest[ij] = 0;
273     fh1PtTracksIn[ij] = 0;
274     fh1PtTracksInLow[ij] = 0;
275     fh2NJetsPt[ij]  = 0;
276     fh2NTracksPt[ij]  = 0;
277     fp2MultRPPhiTrackPt[ij] = 0;
278     fp2CentRPPhiTrackPt[ij] = 0;
279     fhnJetPt[ij] = 0;
280     fhnJetPtBest[ij] = 0;
281     fhnJetPtRej[ij] = 0;
282     fhnJetPtQA[ij] = 0;
283     fhnTrackPt[ij] = 0;
284     fhnTrackPtQA[ij] = 0;
285     for(int i = 0;i <= kMaxJets;++i){
286       fh2LTrackPtJetPt[ij][i] = 0;
287       fh1PtIn[ij][i] = 0;
288     }
289
290     fh1DijetMinv[ij] = 0;      
291     fh2DijetDeltaPhiPt[ij] = 0;  
292     fh2DijetAsymPt[ij] = 0; 
293     fh2DijetPt2vsPt1[ij] = 0;
294     fh2DijetDifvsSum[ij] = 0;
295   } 
296
297   DefineOutput(1, TList::Class());  
298 }
299
300
301
302 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
303 {
304
305
306
307   //
308   // Implemented Notify() to read the cross sections
309   // and number of trials from pyxsec.root
310   // 
311   
312   // Fetch the aod also from the input in,
313   // have todo it in notify
314   
315   
316   fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
317   //  assume that the AOD is in the general output...
318   fAODOut  = AODEvent();
319   
320   if(fNonStdFile.Length()!=0){
321     // case that we have an AOD extension we need can fetch the jets from the extended output
322     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
323     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);    
324     if(!fAODExtension){
325       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
326     }
327   }
328
329
330   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
331   Float_t xsection = 0;
332   Float_t ftrials  = 1;
333
334   fAvgTrials = 1;
335   if(tree){
336     TFile *curfile = tree->GetCurrentFile();
337     if (!curfile) {
338       Error("Notify","No current file");
339       return kFALSE;
340     }
341     if(!fh1Xsec||!fh1Trials){
342       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
343       return kFALSE;
344     }
345     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
346     fh1Xsec->Fill("<#sigma>",xsection);
347     // construct a poor man average trials 
348     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
349     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
350     fh1Trials->Fill("#sum{ntrials}",ftrials); 
351   }  
352
353   if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
354
355   return kTRUE;
356 }
357
358 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
359 {
360
361   
362   // Connect the AOD
363
364   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
365   OpenFile(1);
366   if(!fHistList)fHistList = new TList(); 
367   PostData(1, fHistList); // post data in any case once
368
369   if(!fRandomizer)fRandomizer = new TRandom3(0);
370
371   fHistList->SetOwner(kTRUE);
372   Bool_t oldStatus = TH1::AddDirectoryStatus(); 
373   TH1::AddDirectory(kFALSE);
374
375
376   
377   // event npsparse cent, mult
378   const Int_t nBinsSparse0 = 4;
379   const Int_t nBins0[nBinsSparse0] = {     100, 500,fNTrigger,125};
380   const Double_t xmin0[nBinsSparse0]  = {    0,   0, -0.5,-2};
381   const Double_t xmax0[nBinsSparse0]  = {  100,5000,fNTrigger-0.5,248};
382       
383
384   fhnEvent = new THnSparseF("fhnEvent",";cent;mult:trigger;#rho",nBinsSparse0,
385                             nBins0,xmin0,xmax0);
386   fHistList->Add(fhnEvent);
387
388   if(fDoMatching){
389     MakeJetContainer();
390     fHistList->Add(fhnCorrelation);
391     fHistList->Add(fhnJetContainer);
392   }
393   //
394   //  Histogram
395     
396
397   fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 7, -0.5, 6.5);
398   fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
399   fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
400   fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
401   fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
402   fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
403   fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
404   fh1EvtSelection->GetXaxis()->SetBinLabel(7,"pileup: rejected");
405
406   fHistList->Add(fh1EvtSelection);
407  
408   const Int_t nBinPt = 150;
409   Double_t binLimitsPt[nBinPt+1];
410   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
411     if(iPt == 0){
412       binLimitsPt[iPt] = -50.0;
413     }
414     else {// 1.0
415       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
416     }
417   }
418   const Int_t nBinPhi = 90;
419   Double_t binLimitsPhi[nBinPhi+1];
420   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
421     if(iPhi==0){
422       binLimitsPhi[iPhi] = 0;
423     }
424     else{
425       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
426     }
427   }
428
429   const Int_t nBinEta = 40;
430   Double_t binLimitsEta[nBinEta+1];
431   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
432     if(iEta==0){
433       binLimitsEta[iEta] = -2.0;
434     }
435     else{
436       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
437     }
438   }
439
440
441   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
442   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
443   fHistList->Add(fh1Xsec);
444   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
445   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
446   fHistList->Add(fh1Trials);
447   fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
448   fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
449   fHistList->Add(fh1AvgTrials);
450   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
451   fHistList->Add(fh1PtHard);
452   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
453   fHistList->Add(fh1PtHardNoW);
454   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
455   fHistList->Add(fh1PtHardTrials);
456   
457   fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
458   fHistList->Add(fh1ZVtx);
459
460
461   fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
462   fHistList->Add(fh1RP);
463
464   fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",103,-1,102);
465   fHistList->Add(fh1Centrality);
466
467   fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",500,0,5000,500,0.,5000);
468   fHistList->Add(fh2MultRec);
469   fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,0,5000,500,0.,5000);
470   fHistList->Add(fh2MultGen);
471
472
473   fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
474   fHistList->Add(fh2RPCentrality);
475
476   fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
477   fHistList->Add(fh2PtFGen);
478
479   //  fh2deltaPt1Pt2(0x0),
480   fh2deltaPt1Pt2 = new TH3F("fh2deltaPt1Pt2",Form("%s vs. %s;delta;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
481   fHistList->Add(fh2deltaPt1Pt2);
482
483   fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
484   fHistList->Add(fh2RelPtFGen);
485
486   const Int_t nBinsLeadingTrackPt = 10;
487   const Double_t binArrayLeadingTrackPt[nBinsLeadingTrackPt+1] = {0.,1.,2.,3.,4.,5.,6.,8.,10.,12.,200.}; //store pT of leading track in jet
488
489   const Int_t nBinDeltaPt = 241;
490   Double_t binLimitsDeltaPt[nBinDeltaPt+1];
491   for(Int_t iDeltaPt = 0;iDeltaPt <= nBinDeltaPt;iDeltaPt++){
492     if(iDeltaPt == 0){
493       binLimitsDeltaPt[iDeltaPt] = -2.41;
494     }
495     else {
496       binLimitsDeltaPt[iDeltaPt] =  binLimitsDeltaPt[iDeltaPt-1] + 0.02;
497     }
498   }
499
500   fh3RelPtFGenLeadTrkPt = new TH3F("fh3RelPtFGenLeadTrkPt",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen};p_{T}^{leading track}",nBinPt,binLimitsPt,nBinDeltaPt,binLimitsDeltaPt,nBinsLeadingTrackPt,binArrayLeadingTrackPt);
501   fHistList->Add(fh3RelPtFGenLeadTrkPt);  
502
503   for(int ij = 0;ij <kJetTypes;++ij){    
504     TString cAdd = "";
505     TString cJetBranch = "";
506     if(ij==kJetRec){
507       cAdd = "Rec";
508       cJetBranch = fBranchRec.Data();
509     }
510     else if (ij==kJetGen){
511       cAdd = "Gen";
512       cJetBranch = fBranchGen.Data();
513     }
514     else if (ij==kJetRecFull){
515       cAdd = "RecFull";
516       cJetBranch = fBranchRec.Data();
517     }
518     else if (ij==kJetGenFull){
519       cAdd = "GenFull";
520       cJetBranch = fBranchGen.Data();
521     }
522
523     if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
524     if(!fFlagJetType[ij])continue;
525
526     fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
527     fHistList->Add(fh1NJets[ij]);
528     
529     fh1PtJetsIn[ij]  = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
530     fHistList->Add(fh1PtJetsIn[ij]);
531
532     fh1PtJetsInRej[ij]  = new TH1F(Form("fh1PtJets%sInRej",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
533     fHistList->Add(fh1PtJetsInRej[ij]);
534
535     fh1PtJetsInBest[ij]  = new TH1F(Form("fh1PtJets%sInBest",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
536     fHistList->Add(fh1PtJetsInBest[ij]);
537     
538     fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
539     fHistList->Add(fh1PtTracksIn[ij]);
540     
541     fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),100,0.,5.);
542     fHistList->Add(fh1PtTracksInLow[ij]);
543     
544     fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),1000,0.,3000.);
545     fHistList->Add(fh1SumPtTrack[ij]);
546     
547     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);
548     fHistList->Add(fh2NJetsPt[ij]);
549     
550     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);
551     fHistList->Add(fh2NTracksPt[ij]);
552
553     
554     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");
555     fHistList->Add(fp2MultRPPhiTrackPt[ij]);
556     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");
557     fHistList->Add(fp2CentRPPhiTrackPt[ij]);    
558
559     // Bins:  Jet number: pTJet, cent, mult, RP, Area, trigger, leading track pT total bins = 4.5M
560     const Int_t nBinsSparse1 = 9;
561     const Int_t nBinsArea = 10;
562     Int_t nbinr5=120;
563     Int_t nbinr5min=-50;
564   
565     if(cJetBranch.Contains("KT05")){
566      nbinr5=132;
567      nbinr5min=-80;}
568          
569     Int_t nBins1[nBinsSparse1] = {     kMaxJets+1,nbinr5, 10,  fNBinsMult,    fNRPBins, nBinsArea,fNTrigger,fNBinsLeadingTrackPt,fNAcceptance+1};
570     if(cJetBranch.Contains("RandomCone")){
571       nBins1[1] = 600;
572       nBins1[5] = 1;
573     }
574     const Double_t xmin1[nBinsSparse1]  = {        -0.5, static_cast<Double_t>(nbinr5min),  0,   0,        -0.5, 0.,         -0.5,  0.,           -0.5,};
575     const Double_t xmax1[nBinsSparse1]  = {kMaxJets+0.5,250,100,4000,fNRPBins-0.5,1.0,fNTrigger-0.5,200.,fNAcceptance+0.5};
576      
577     const Double_t binArrayArea[nBinsArea+1] = {xmin1[5],0.07,0.1,0.15,0.2,0.25,0.3,0.4,0.5,0.6,xmax1[5]};
578
579
580     fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leading track p_{T};acceptance bin",nBinsSparse1,nBins1,xmin1,xmax1);
581     fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
582     if(fNBinsLeadingTrackPt>1) fhnJetPt[ij]->SetBinEdges(7,binArrayLeadingTrackPt);
583     fHistList->Add(fhnJetPt[ij]);
584
585
586     // Bins:  pTJet, cent, trigger, 
587     const Int_t nBinsSparse1b = 3;
588     Int_t nBins1b[nBinsSparse1b] = {120, 10,fNTrigger};
589     const Double_t xmin1b[nBinsSparse1b]  = {-50,  0,-0.5};
590     const Double_t xmax1b[nBinsSparse1b]  = {250,100,fNTrigger-0.5};
591
592     fhnJetPtBest[ij] = new THnSparseF(Form("fhnJetPtBest%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b);
593     fHistList->Add(fhnJetPtBest[ij]);
594
595     fhnJetPtRej[ij] = new THnSparseF(Form("fhnJetPtRej%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b);
596     fHistList->Add(fhnJetPtRej[ij]);
597     
598     // Bins:  Jet number: pTJet, cent, eta, phi, Area, trigger, acceptance, signed pT leading
599     const Int_t nBinsSparse2 = 9;
600     Int_t nBins2[nBinsSparse2] = {     kMaxJets+1, 60,   8,  6, 90, 100,fNTrigger,fNAcceptance+1,20};
601     if(cJetBranch.Contains("RandomCone")){
602       nBins2[5] = 1;
603     }
604     const Double_t xmin2[nBinsSparse2]  = {        -0.5, -50,  0,-0.9,              0,  0.,        -0.5,            -0.5, -100};
605     const Double_t xmax2[nBinsSparse2]  = {kMaxJets+0.5, 250, 80, 0.9, 2.*TMath::Pi(),1.0,fNTrigger-0.5,fNAcceptance+0.5,  100};
606     fhnJetPtQA[ij] = new THnSparseF(Form("fhnJetPtQA%s",cAdd.Data()),";jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading",nBinsSparse2,nBins2,xmin2,xmax2);
607     //    fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
608     fHistList->Add(fhnJetPtQA[ij]);
609     
610     // Bins:track number  pTtrack, cent, mult, RP.  charge total bins = 224 k
611     const Int_t nBinsSparse3 = 7;
612     const Int_t nRPBinsSparse3 = 3;
613     const Int_t nBins3[nBinsSparse3] = {       2,    100,     10,   1,    nRPBinsSparse3,fNTrigger,2};
614     const Double_t xmin3[nBinsSparse3]  = { -0.5,      0,   0,      0,        -0.5,-0.5,-1.5};
615     const Double_t xmax3[nBinsSparse3]  = { 1.5,     200, 100,   4000,nRPBinsSparse3-0.5,fNTrigger-0.5,1.5};  
616     
617     // change the binning of the pT axis:
618     Double_t *xPt3 = new Double_t[nBins3[1]+1];
619     xPt3[0] = 0.;
620     for(int i = 1; i<=nBins3[1];i++){
621       if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.05; // 1 - 40
622       else if(xPt3[i-1]<4)xPt3[i] = xPt3[i-1] + 0.2; // 41 - 50
623       else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 0.5; // 50 - 62
624       else if(xPt3[i-1]<15)xPt3[i] = xPt3[i-1] +  1.; // 62 - 67
625       else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] + 2.; // 67 - 72
626       else if(xPt3[i-1]<60)xPt3[i] = xPt3[i-1] + 5; // 72 - 76
627       else xPt3[i] = xPt3[i-1] + 10; // 76 - 100 = 140 
628     }
629     
630     fhnTrackPt[ij] = new THnSparseF(Form("fhnTrackPt%s",cAdd.Data()),";track number;p_{T};cent;#tracks;RP;trigger",nBinsSparse3,nBins3,xmin3,xmax3);
631     fhnTrackPt[ij]->SetBinEdges(1,xPt3);
632     fHistList->Add(fhnTrackPt[ij]);
633     delete [] xPt3;
634
635     // Track QA bins track nr, pTrack, cent, eta, phi bins 5.4 M
636     const Int_t nBinsSparse4 = 6;
637     const Int_t nBins4[nBinsSparse4] =    {    2, 50,  10,   20, 180,2};
638     const Double_t xmin4[nBinsSparse4]  = { -0.5,  0,   0, -1.0,  0.,-1.5};
639     const Double_t xmax4[nBinsSparse4]  = {  1.5,150, 100,  1.0,2.*TMath::Pi(),1.5};  
640     
641     // change the binning ot the pT axis:
642     Double_t *xPt4 = new Double_t[nBins4[1]+1];
643     xPt4[0] = 0.;
644     for(int i = 1; i<=nBins4[1];i++){
645       if(xPt4[i-1]<2)xPt4[i] = xPt4[i-1] + 0.1;
646       else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 0.5;
647       else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] +  1.;
648       else if(xPt4[i-1]<30)xPt4[i] = xPt4[i-1] +  2.5;
649       else xPt4[i] = xPt4[i-1] + 5.;
650     }
651     fhnTrackPtQA[ij] = new THnSparseF(Form("fhnTrackPtQA%s",cAdd.Data()),";track number;p_{T};cent;#eta;#phi;sign",nBinsSparse4,nBins4,xmin4,xmax4);
652     fhnTrackPtQA[ij]->SetBinEdges(1,xPt4);
653     fHistList->Add(fhnTrackPtQA[ij]);
654     delete [] xPt4;
655     
656     for(int i = 0;i <= kMaxJets;++i){
657       fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
658       fHistList->Add(fh1PtIn[ij][i]);
659
660
661       if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
662       fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
663                                          Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
664                                               cAdd.Data()),
665                                          200,0.,200.,nBinPt,binLimitsPt);
666       fHistList->Add(fh2LTrackPtJetPt[ij][i]);
667     }
668
669
670     fh1DijetMinv[ij]                = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
671     fHistList->Add(fh1DijetMinv[ij]);
672
673     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);
674     fHistList->Add(fh2DijetDeltaPhiPt[ij]);
675
676     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);
677     fHistList->Add(fh2DijetAsymPt[ij]);
678
679     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.);
680     fHistList->Add(fh2DijetPt2vsPt1[ij]);
681     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.);
682     fHistList->Add( fh2DijetDifvsSum[ij]);
683   }   
684   // =========== Switch on Sumw2 for all histos ===========
685   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
686     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
687     if (h1){
688       h1->Sumw2();
689       continue;
690     }
691     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
692     if(hn)hn->Sumw2();
693   }
694   TH1::AddDirectory(oldStatus);
695 }
696
697 void AliAnalysisTaskJetSpectrum2::Init()
698 {
699   //
700   // Initialization
701   //
702
703   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
704
705 }
706
707 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
708
709
710   Bool_t selected = kTRUE;
711
712   if(fUseGlobalSelection&&fEventSelectionMask==0){
713     selected = AliAnalysisHelperJetTasks::Selected();
714   }
715   else if(fUseGlobalSelection&&fEventSelectionMask>0){
716     selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
717   }
718
719   if(!selected){
720     // no selection by the service task, we continue
721     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
722     fh1EvtSelection->Fill(1.); 
723     PostData(1, fHistList);
724     return;
725   }
726
727   if(fEventClass>0){
728     selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
729   }
730
731   if(!selected){
732     // no selection by the service task, we continue
733     if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d  Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
734     fh1EvtSelection->Fill(2.);
735     PostData(1, fHistList);
736     return;
737   }
738
739
740   static AliAODEvent* aod = 0;
741
742   // take all other information from the aod we take the tracks from
743   if(!aod){
744     if(fUseAODTrackInput)aod = fAODIn;
745     else aod = fAODOut;
746   }
747
748
749   //
750   // Execute analysis for current event
751   //
752   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);  
753   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
754   if(!aodH){
755     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
756     return;
757   }
758
759   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
760   TClonesArray *aodRecJets = 0;
761
762   if(fAODOut&&!aodRecJets){
763     aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
764   }
765   if(fAODExtension&&!aodRecJets){
766     aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
767   }
768   if(fAODIn&&!aodRecJets){
769     aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
770   }
771
772
773
774   if(!aodRecJets){
775     if(fDebug){
776
777       Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
778       if(fAODIn){
779         Printf("Input AOD >>>>");
780         fAODIn->Print();    
781       }
782       if(fAODExtension){
783         Printf("AOD Extension >>>>");
784         fAODExtension->Print();
785       }
786       if(fAODOut){
787         Printf("Output AOD >>>>");
788         fAODOut->Print();    
789       }
790     }
791     return;
792   }
793
794   TClonesArray *aodGenJets = 0;
795   if(fBranchGen.Length()>0){
796     if(fAODOut&&!aodGenJets){
797       aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
798     }
799     if(fAODExtension&&!aodGenJets){
800       aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
801     }
802     if(fAODIn&&!aodGenJets){
803       aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
804     }
805
806     if(!aodGenJets){
807       Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
808       return;
809     }
810   }
811
812   TClonesArray *aodBackRecJets = 0;
813   if(fBranchBkgRec.Length()>0){
814     if(fAODOut&&!aodBackRecJets){
815       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
816     }
817     if(fAODExtension&&!aodBackRecJets){
818       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
819     }
820     if(fAODIn&&!aodBackRecJets){
821       aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
822     }
823
824     if(!aodBackRecJets){
825       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
826       return;
827     }
828   }
829
830
831   TClonesArray *aodBackGenJets = 0;
832
833   if(fBranchBkgGen.Length()>0){
834     if(fAODOut&&!aodBackGenJets){
835       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
836     }
837     if(fAODExtension&&!aodBackGenJets){
838       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
839     }
840     if(fAODIn&&!aodBackGenJets){
841       aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
842     }
843
844     if(!aodBackGenJets){
845       Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
846       return;
847     }
848   }
849
850   AliAODVertex* primVtx = aod->GetPrimaryVertex();
851   Int_t nTracksPrim = primVtx->GetNContributors();
852   
853   
854   if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
855   if(nTracksPrim < fMinNcontributors){
856     if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__); 
857     fh1EvtSelection->Fill(3.);
858     PostData(1, fHistList);
859     return;
860   }
861   
862   if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
863     if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ()); 
864     fh1EvtSelection->Fill(4.);
865     PostData(1, fHistList);
866     return; 
867   }
868   
869   TString primVtxName(primVtx->GetName());
870
871   if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
872     if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
873     fh1EvtSelection->Fill(5.);
874     PostData(1, fHistList);
875     return;
876   }
877
878   if(fRejectPileup && AliAnalysisHelperJetTasks::IsPileUp()){
879     if (fDebug > 1) Printf("%s:%d SPD pileup: event REJECTED...",(char*)__FILE__,__LINE__);
880     fh1EvtSelection->Fill(6.);
881     PostData(1, fHistList);
882     return;
883   }
884
885   if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__); 
886   fh1EvtSelection->Fill(0.);
887
888  
889   // new Scheme
890   // first fill all the pure  histograms, i.e. full jets 
891   // in case of the correaltion limit to the n-leading jets
892
893   // reconstructed
894
895   
896   // generated
897
898
899   // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
900
901
902
903   TList genJetsList;         // full acceptance
904   TList genJetsListCut;      // acceptance cut
905   TList recJetsList;         // full acceptance
906   TList recJetsListCut;      // acceptance cut
907
908
909   GetListOfJets(&recJetsList,aodRecJets,0);
910   GetListOfJets(&recJetsListCut,aodRecJets,1);
911
912   if(fBranchGen.Length()>0){
913     GetListOfJets(&genJetsList,aodGenJets,0);
914     GetListOfJets(&genJetsListCut,aodGenJets,1);
915   }
916
917   Double_t eventW = 1;
918   Double_t ptHard = 0; 
919   Double_t nTrials = 1; // Trials for MC trigger 
920   fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials); 
921
922   // Getting some global properties
923   fCentrality = GetCentrality();
924   if(fCentrality<=0)fCentrality = 0;
925   fh1Centrality->Fill(fCentrality);
926
927
928   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
929     // this is the part we only use when we have MC information
930     AliMCEvent* mcEvent = MCEvent();
931     //    AliStack *pStack = 0; 
932     if(!mcEvent){
933       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
934       return;
935     }
936     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
937     if(pythiaGenHeader){
938       nTrials = pythiaGenHeader->Trials();
939       ptHard  = pythiaGenHeader->GetPtHard();
940       int iProcessType = pythiaGenHeader->ProcessType();
941       // 11 f+f -> f+f
942       // 12 f+barf -> f+barf
943       // 13 f+barf -> g+g
944       // 28 f+g -> f+g
945       // 53 g+g -> f+barf
946       // 68 g+g -> g+g
947       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
948       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
949     } 
950   }// (fAnalysisType&kMCESD)==kMCESD)
951
952
953   // we simply fetch the tracks/mc particles as a list of AliVParticles
954
955   TList recParticles;
956   TList genParticles;
957
958   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
959
960   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
961   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
962   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
963
964   //  CalculateReactionPlaneAngle(&recParticles);
965   fRPAngle = 0;
966   
967   if(fRPMethod==0)fRPAngle = ((AliVAODHeader*)aod->GetHeader())->GetEventplane();
968   else if(fRPMethod==1||fRPMethod==2){
969     AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(aod->GetHeader());
970     if(!aodheader) AliFatal("Not a standard AOD");
971     fRPAngle = aodheader->GetQTheta(fRPMethod);
972   }
973   fh1RP->Fill(fRPAngle);
974   fh2RPCentrality->Fill(fCentrality,fRPAngle);
975   // Event control and counting ...  
976   // MC
977   fh1PtHard->Fill(ptHard,eventW);
978   fh1PtHardNoW->Fill(ptHard,1);
979   fh1PtHardTrials->Fill(ptHard,nTrials);
980
981   // Real
982   if(aod->GetPrimaryVertex()){// No vtx for pure MC
983     fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
984   }
985
986
987   Int_t recMult1 = recParticles.GetEntries();
988   Int_t genMult1 = genParticles.GetEntries();
989
990   Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
991   Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
992
993   fh2MultRec->Fill(recMult1,recMult2);
994   fh2MultGen->Fill(genMult1,genMult2);
995   fMultRec = recMult1;
996   if(fMultRec<=0)fMultRec = recMult2;
997   fMultGen = genMult1;
998   if(fMultGen<=0)fMultGen = genMult2;
999
1000   Double_t var0[4] = {0,};
1001   var0[0] = fCentrality;
1002   var0[1] = fMultRec;
1003   for(int it=0;it<fNTrigger;it++){
1004     if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1005       var0[2] = it;
1006       var0[3] = GetRho(recJetsList);
1007       fhnEvent->Fill(var0);
1008     }
1009   }
1010   // the loops for rec and gen should be indentical... pass it to a separate
1011   // function ...
1012   // Jet Loop
1013   // Track Loop
1014   // Jet Jet Loop
1015   // Jet Track loop
1016
1017   FillJetHistos(recJetsListCut,recParticles,kJetRec);
1018   FillJetHistos(recJetsList,recParticles,kJetRecFull);
1019   FillTrackHistos(recParticles,kJetRec);
1020
1021   FillJetHistos(genJetsListCut,genParticles,kJetGen);
1022   FillJetHistos(genJetsList,genParticles,kJetGenFull);
1023   FillTrackHistos(genParticles,kJetGen);
1024
1025   // Here follows the jet matching part
1026   // delegate to a separated method?
1027
1028   if(fDoMatching){
1029     FillMatchHistos(recJetsList,genJetsList);
1030   }
1031
1032   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1033   PostData(1, fHistList);
1034 }
1035
1036 void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
1037
1038   if(iType>=kJetTypes){
1039     return;
1040   }
1041   if(!fFlagJetType[iType])return;
1042
1043   Int_t refMult = fMultRec;
1044   if(iType==kJetGen||iType==kJetGenFull){
1045     refMult = fMultGen;
1046   }
1047
1048   Int_t nJets = jetsList.GetEntries(); 
1049   fh1NJets[iType]->Fill(nJets);
1050
1051   if(nJets<=0)return;
1052   
1053   Float_t pT1 = 0;
1054   Float_t phi0 = 0;
1055   Float_t phi1 = 0;
1056   Int_t ij0 = -1;
1057   Int_t ij1 = -1;
1058
1059   Double_t var1[9] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt
1060   Double_t var1b[3] = {0,}; // p_{T,jet};cent;trigger;
1061
1062   var1[2] = fCentrality; 
1063   var1b[1] = fCentrality; 
1064   var1[3] = refMult;
1065
1066   
1067
1068   Double_t var2[9] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area;trigger
1069   var2[2] = fCentrality;
1070
1071   for(int ij = 0;ij < nJets;ij++){
1072     AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
1073     Float_t ptJet = jet->Pt();
1074     
1075     
1076     if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
1077
1078     var1b[0] = ptJet;
1079     if(jet->Trigger()&fJetTriggerBestMask){
1080       fh1PtJetsInBest[iType]->Fill(ptJet);
1081       for(int it = 0;it <fNTrigger;it++){
1082         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1083           var1b[2] = it;
1084           fhnJetPtBest[iType]->Fill(var1b);
1085         }
1086       }
1087     }
1088     if(jet->Trigger()&fJetTriggerExcludeMask){
1089       fh1PtJetsInRej[iType]->Fill(ptJet);
1090       for(int it = 0;it <fNTrigger;it++){
1091         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1092           var1b[2] = it;
1093           fhnJetPtRej[iType]->Fill(var1b);
1094         }
1095       }
1096       continue;
1097     }
1098     fh1PtJetsIn[iType]->Fill(ptJet);
1099     
1100     // find the dijets assume sorting and acceptance cut...
1101     if(ij==0){
1102       ij0 = ij;
1103       phi0 = jet->Phi();
1104       if(phi0<0)phi0+=TMath::Pi()*2.;
1105     }
1106     else if(ptJet>pT1){
1107       // take only the backward hemisphere??                                                        
1108       phi1 = jet->Phi();
1109       if(phi1<0)phi1+=TMath::Pi()*2.;
1110       Float_t dPhi = phi1 - phi0;
1111       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1112       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1113       if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
1114         ij1 = ij;
1115         pT1 = ptJet;
1116       }
1117     }
1118     // fill jet histos for kmax jets
1119     Float_t phiJet = jet->Phi();
1120     Float_t etaJet = jet->Eta();
1121     if(phiJet<0)phiJet+=TMath::Pi()*2.;    
1122     fh1TmpRho->Reset();
1123     if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
1124
1125     fh1PtIn[iType][kMaxJets]->Fill(ptJet);
1126     // fill leading jets...
1127     AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
1128     //AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
1129     Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
1130     Double_t ptLead = jet->GetPtLeading(); //pT of leading jet
1131     AliVParticle *tt = (AliVParticle*)particlesList.At(0);
1132     Float_t ttphi=tt->Phi();
1133     if(ttphi<0)ttphi+=TMath::Pi()*2.;  
1134     Float_t ttpt=tt->Pt();
1135     Int_t phiBintt = GetPhiBin(ttphi-fRPAngle);
1136     Double_t dphitrigjet=RelativePhi(ttphi,phiJet);
1137     if(fTRP==1){
1138       if(TMath::Abs(dphitrigjet)<TMath::Pi()-0.6) continue; 
1139       var1[1] = ptJet;
1140       var1[4] = phiBintt;
1141       var1[5] = jet->EffectiveAreaCharged();
1142       var1[7] = ttpt;
1143       var1[8] = CheckAcceptance(phiJet,etaJet);
1144     }
1145     
1146     if(fTRP==0){
1147       var1[1] = ptJet;
1148       var1[4] = phiBin;
1149       var1[5] = jet->EffectiveAreaCharged();
1150       var1[7] = ptLead;
1151       var1[8] = CheckAcceptance(phiJet,etaJet);
1152     }
1153     
1154     //jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading
1155     var2[1] = ptJet;
1156     var2[3] = etaJet;
1157     var2[4] = phiJet;
1158     var2[5] = jet->EffectiveAreaCharged();
1159     var2[7] = var1[8];
1160     var2[8] = (leadTrack?leadTrack->Charge()*ptLead:0);//pT of leading jet x charge
1161
1162     if(ij<kMaxJets){
1163       fh2LTrackPtJetPt[iType][ij]->Fill(jet->GetPtLeading(),ptJet);
1164       var1[0] = ij;
1165       var2[0] = ij;
1166       for(int it = 0;it <fNTrigger;it++){
1167         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1168           var1[6] = it;
1169           var2[6] = it;
1170           fhnJetPt[iType]->Fill(var1);
1171           fhnJetPtQA[iType]->Fill(var2);
1172         }
1173       }
1174     }
1175     var1[0] = kMaxJets;// fill for all jets
1176     var2[0] = kMaxJets;// fill for all jets
1177     for(int it = 0;it <fNTrigger;it++){
1178       if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1179         var1[6] = it;
1180         var2[6] = it;
1181         fhnJetPt[iType]->Fill(var1);
1182         fhnJetPtQA[iType]->Fill(var2);
1183       }
1184     }
1185
1186     fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet);
1187
1188     if(particlesList.GetSize()&&ij<kMaxJets){
1189       // Particles... correlated with jets...
1190       for(int it = 0;it<particlesList.GetEntries();++it){
1191         AliVParticle *part = (AliVParticle*)particlesList.At(it);
1192         Float_t deltaR = jet->DeltaR(part);
1193         if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1194       }
1195       // fill the jet shapes
1196     }// if we have particles
1197   }// Jet Loop
1198
1199
1200   // do something with dijets...
1201   if(ij0>=0&&ij1>0){
1202     AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1203     Double_t ptJet0 = jet0->Pt();
1204     Double_t phiJet0 = jet0->Phi();
1205     if(phiJet0<0)phiJet0+=TMath::Pi()*2.;       
1206
1207     AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1208     Double_t ptJet1 = jet1->Pt();
1209     Double_t phiJet1 = jet1->Phi();
1210     if(phiJet1<0)phiJet1+=TMath::Pi()*2.;       
1211
1212     Float_t deltaPhi = phiJet0 - phiJet1;
1213     if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1214     if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();      
1215     deltaPhi = TMath::Abs(deltaPhi);
1216     fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);      
1217
1218     Float_t asym = 9999;
1219     if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1220     fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1221     fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);        
1222     fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);        
1223     Float_t minv = 2.*(jet0->P()*jet1->P()-
1224                        jet0->Px()*jet1->Px()- 
1225                        jet0->Py()*jet1->Py()- 
1226                        jet0->Pz()*jet1->Pz());    // assume mass == 0;
1227     if(minv<0)minv=0; // prevent numerical instabilities
1228     minv = TMath::Sqrt(minv);
1229     fh1DijetMinv[iType]->Fill(minv);            
1230   }
1231   
1232
1233
1234   // count down the jets above thrueshold
1235   Int_t nOver = nJets;
1236   if(nOver>0){
1237     TIterator *jetIter = jetsList.MakeIterator();
1238     AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());  
1239     if(tmpJet){
1240       Float_t pt = tmpJet->Pt();
1241       for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1242         Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1243         while(pt<ptCut&&tmpJet){
1244           nOver--;
1245           tmpJet = (AliAODJet*)(jetIter->Next()); 
1246           if(tmpJet){
1247             pt = tmpJet->Pt();
1248           }
1249         }
1250         if(nOver<=0)break;
1251         fh2NJetsPt[iType]->Fill(ptCut,nOver);
1252       }
1253     }
1254     delete jetIter;
1255   }
1256 }
1257
1258 void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1259
1260   if(fFlagJetType[iType]<=0)return;
1261   Int_t refMult = fMultRec;
1262   if(iType==kJetGen||iType==kJetGenFull){
1263     refMult = fMultGen;
1264
1265   }
1266
1267   // 
1268   Double_t var3[7]; // track number;p_{T};cent;#tracks;RP;cahrge
1269   var3[2] = fCentrality;
1270   var3[3] = refMult;
1271   Double_t var4[6]; // track number;p_{T};cent;#eta;#phi;sign
1272   var4[2] = fCentrality;
1273   Int_t nTrackOver = particlesList.GetSize();
1274   // do the same for tracks and jets
1275   if(nTrackOver>0){
1276     TIterator *trackIter = particlesList.MakeIterator();
1277     AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());  
1278     Float_t pt = tmpTrack->Pt();
1279     for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1280       Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1281       while(pt<ptCut&&tmpTrack){
1282         nTrackOver--;
1283         tmpTrack = (AliVParticle*)(trackIter->Next()); 
1284         if(tmpTrack){
1285           pt = tmpTrack->Pt();
1286         }
1287       }
1288       if(nTrackOver<=0)break;
1289       fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1290     }
1291
1292     
1293     trackIter->Reset();
1294     AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1295     Float_t sumPt = 0;
1296
1297     while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1298       Float_t tmpPt = tmpTrack->Pt();
1299       fh1PtTracksIn[iType]->Fill(tmpPt);
1300       fh1PtTracksInLow[iType]->Fill(tmpPt);
1301
1302       sumPt += tmpPt;
1303       Float_t tmpPhi = tmpTrack->Phi();
1304       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1305
1306       
1307       Float_t phiRP = tmpPhi-fRPAngle;
1308       if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1309       if(phiRP<0)phiRP += TMath::Pi();
1310       if(phiRP<0)phiRP += TMath::Pi();
1311       const float allPhi = -1./180.*TMath::Pi();
1312
1313       if(tmpPt<100){
1314         fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1315         fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1316         
1317         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1318         fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1319       }
1320       Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1321       var3[0] = 1;
1322       var3[1] = tmpPt;
1323       var3[4] = phiBin;
1324       var3[6] = tmpTrack->Charge();
1325       
1326       var4[0] = 1;
1327       var4[1] = tmpPt;
1328       var4[3] = tmpTrack->Eta();
1329       var4[4] = tmpPhi;
1330       var4[5] = tmpTrack->Charge();
1331
1332       
1333       for(int it = 0;it <fNTrigger;it++){
1334         if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1335           var3[0] = 1;
1336           var3[5] = it;
1337           fhnTrackPt[iType]->Fill(var3);
1338           if(tmpTrack==leading){
1339             var3[0] = 0;         
1340             fhnTrackPt[iType]->Fill(var3);   
1341           }
1342         }
1343       }
1344
1345
1346
1347
1348       fhnTrackPtQA[iType]->Fill(var4);
1349       if(tmpTrack==leading){
1350         var4[0] = 0;
1351         fhnTrackPtQA[iType]->Fill(var4);
1352         continue;
1353       }
1354
1355
1356
1357
1358       
1359     }  
1360     fh1SumPtTrack[iType]->Fill(sumPt);
1361     delete trackIter;
1362   }
1363
1364 }
1365
1366
1367 void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1368
1369
1370   // Fill al the matching histos
1371   // It is important that the acceptances for the mathing are as large as possible
1372   // to avoid false matches on the edge of acceptance
1373   // therefore we add some extra matching jets as overhead
1374
1375   static TArrayI aGenIndex(recJetsList.GetEntries());
1376   if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1377
1378   static TArrayI aRecIndex(genJetsList.GetEntries());
1379   if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1380
1381   if(fDebug){
1382     Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1383     Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1384   }
1385   AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1386                                             &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1387                                             aGenIndex,aRecIndex,fDebug);
1388
1389   if(fDebug){
1390     for(int i = 0;i< aGenIndex.GetSize();++i){ 
1391       if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]); 
1392     }
1393     for(int i = 0;i< aRecIndex.GetSize();++i){
1394       if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]); 
1395     }
1396   }
1397
1398   Double_t container[8];
1399   Double_t containerRec[4];
1400   Double_t containerGen[4];
1401
1402   // loop over generated jets
1403   // consider the 
1404   for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1405     AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1406     Double_t ptGen = genJet->Pt();
1407     Double_t phiGen = genJet->Phi();
1408     if(phiGen<0)phiGen+=TMath::Pi()*2.;    
1409     Double_t etaGen = genJet->Eta();
1410     containerGen[0] = ptGen;
1411     containerGen[1] = etaGen;
1412     containerGen[2] = phiGen;
1413     containerGen[3] = genJet->GetPtLeading();
1414
1415     fhnJetContainer->Fill(containerGen,kStep0); //all generated jets
1416
1417     if(JetSelected(genJet))
1418       fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window
1419
1420     Int_t ir = aRecIndex[ig];
1421     if(ir>=0&&ir<recJetsList.GetEntries()){   
1422       AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir); 
1423
1424       fhnJetContainer->Fill(containerGen,kStep2); // all generated jets with reconstructed partner
1425
1426       if(JetSelected(genJet)){  
1427         fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner
1428         
1429         if(JetSelected(recJet)) {
1430           fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window
1431           
1432         }
1433         containerGen[3] = recJet->GetPtLeading();
1434         fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner with leading track on reconstructed level
1435       }
1436     }
1437   }// loop over generated jets used for matching...
1438
1439
1440
1441   // loop over reconstructed jets
1442   for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1443     AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1444     Double_t etaRec = recJet->Eta();
1445     Double_t ptRec = recJet->Pt();
1446     Double_t phiRec = recJet->Phi();
1447     if(phiRec<0)phiRec+=TMath::Pi()*2.;    
1448     
1449     containerRec[0] = ptRec;
1450     containerRec[1] = etaRec;
1451     containerRec[2] = phiRec;
1452     containerRec[3] = recJet->GetPtLeading();
1453
1454     container[0] = ptRec;
1455     container[1] = etaRec;
1456     container[2] = phiRec;
1457     container[3] = recJet->GetPtLeading();
1458
1459     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1460   
1461     if(JetSelected(recJet)){
1462       fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window
1463       // Fill Correlation
1464       Int_t ig = aGenIndex[ir];
1465       if(ig>=0 && ig<genJetsList.GetEntries()){
1466         if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1467         AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1468         Double_t ptGen  = genJet->Pt();
1469         Double_t phiGen = genJet->Phi();
1470         if(phiGen<0)phiGen+=TMath::Pi()*2.; 
1471         Double_t etaGen = genJet->Eta();
1472
1473         containerGen[0] = ptGen;
1474         containerGen[1] = etaGen;
1475         containerGen[2] = phiGen;
1476         containerGen[3] = genJet->GetPtLeading();
1477
1478         container[4] = ptGen;
1479         container[5] = etaGen;
1480         container[6] = phiGen;
1481         container[7] = genJet->GetPtLeading();
1482
1483         fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner
1484
1485         fhnCorrelation->Fill(container);
1486         
1487         Float_t delta = (ptRec-ptGen)/ptGen;
1488         fh2RelPtFGen->Fill(ptGen,delta);
1489         fh3RelPtFGenLeadTrkPt->Fill(ptGen,delta,recJet->GetPtLeading());
1490         fh2PtFGen->Fill(ptGen,ptRec);
1491
1492         fh2deltaPt1Pt2->Fill(ptRec-ptGen,ptRec,ptGen);
1493         
1494       } 
1495     }// loop over reconstructed jets
1496   }
1497   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1498 }
1499
1500
1501 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1502   //
1503   // Create the particle container for the correction framework manager and 
1504   // link it
1505   //
1506   const Int_t kNvar   = 4 ; //number of variables on the grid:pt,eta, phi
1507   const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
1508   const Double_t kEtamin = -1.0, kEtamax = 1.0;
1509   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1510   const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.;
1511
1512   // can we neglect migration in eta and phi?
1513   // phi should be no problem since we cover full phi and are phi symmetric
1514   // eta migration is more difficult  due to needed acceptance correction
1515   // in limited eta range
1516
1517   //arrays for the number of bins in each dimension
1518   Int_t iBin[kNvar];
1519   iBin[0] = 125; //bins in pt
1520   iBin[1] =  4; //bins in eta 
1521   iBin[2] =  1; // bins in phi
1522   iBin[3] =  10; // bins in leading track Pt
1523
1524   //arrays for lower bounds :
1525   Double_t* binEdges[kNvar];
1526   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1527     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1528
1529   //values for bin lower bounds
1530   //  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);  
1531   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1532   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1533   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1534   binEdges[3][0]= kPtLeadingTrackPtMin;
1535   binEdges[3][1]= 1.;
1536   binEdges[3][2]= 2.;
1537   binEdges[3][3]= 3.;
1538   binEdges[3][4]= 4.;
1539   binEdges[3][5]= 5.;
1540   binEdges[3][6]= 6.;
1541   binEdges[3][7]= 8.;
1542   binEdges[3][8]= 10.;
1543   binEdges[3][9]= 12.;
1544   binEdges[3][10]= kPtLeadingTrackPtMax;
1545
1546   fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin);
1547   for (int k=0; k<kNvar; k++) {
1548     fhnJetContainer->SetBinLimits(k,binEdges[k]);
1549   }
1550
1551   //create correlation matrix for unfolding
1552   Int_t thnDim[2*kNvar];
1553   for (int k=0; k<kNvar; k++) {
1554     //first half  : reconstructed 
1555     //second half : MC
1556     thnDim[k]      = iBin[k];
1557     thnDim[k+kNvar] = iBin[k];
1558   }
1559
1560   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim);
1561   for (int k=0; k<kNvar; k++) {
1562     fhnCorrelation->SetBinEdges(k,binEdges[k]);
1563     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1564   }
1565
1566   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1567     delete [] binEdges[ivar];
1568
1569
1570 }
1571
1572 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1573 {
1574   // Terminate analysis
1575   //
1576   if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1577 }
1578
1579
1580 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1581
1582   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1583
1584   Int_t iCount = 0;
1585   if(type==kTrackAOD){
1586     AliAODEvent *aod = 0;
1587     if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1588     else aod = AODEvent();
1589     if(!aod){
1590       return iCount;
1591     }
1592     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1593       AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
1594       if(!tr) AliFatal("Not a standard AOD");
1595       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1596       if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1597       if(tr->Pt()<fMinTrackPt)continue;
1598       if(fDebug>0){
1599         if(tr->Pt()>20){
1600           Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1601           Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1602           tr->Print();
1603           //    tr->Dump();
1604           AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1605           if(fESD){
1606             AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1607             esdTr->Print("");
1608             // esdTr->Dump();
1609           }
1610         }
1611       }
1612       list->Add(tr);
1613       iCount++;
1614     }
1615   }
1616   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1617     AliMCEvent* mcEvent = MCEvent();
1618     if(!mcEvent)return iCount;
1619     // we want to have alivpartilces so use get track
1620     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1621       if(!mcEvent->IsPhysicalPrimary(it))continue;
1622       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1623       if(part->Pt()<fMinTrackPt)continue;
1624       if(type == kTrackKineAll){
1625         list->Add(part);
1626         iCount++;
1627       }
1628       else if(type == kTrackKineCharged){
1629         if(part->Particle()->GetPDG()->Charge()==0)continue;
1630         list->Add(part);
1631         iCount++;
1632       }
1633     }
1634   }
1635   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1636     AliAODEvent *aod = 0;
1637     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1638     else aod = AODEvent();
1639     if(!aod)return iCount;
1640     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1641     if(!tca)return iCount;
1642     for(int it = 0;it < tca->GetEntriesFast();++it){
1643       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1644       if(!part)continue;
1645       if(part->Pt()<fMinTrackPt)continue;
1646       if(!part->IsPhysicalPrimary())continue;
1647       if(type == kTrackAODMCAll){
1648         list->Add(part);
1649         iCount++;
1650       }
1651       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1652         if(part->Charge()==0)continue;
1653         if(kTrackAODMCCharged){
1654           list->Add(part);
1655         }
1656         else {
1657           if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1658           list->Add(part);
1659         }
1660         iCount++;
1661       }
1662     }
1663   }// AODMCparticle
1664   list->Sort();
1665   return iCount;
1666
1667 }
1668
1669
1670 Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1671   AliAODEvent *aod = 0;
1672   if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1673   else aod = AODEvent();
1674   if(!aod){
1675     return 101;
1676   }
1677   return ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
1678 }
1679
1680
1681
1682 Bool_t  AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1683   Bool_t selected = false;
1684   
1685   if(!jet)return selected;
1686
1687   if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1688     selected = kTRUE;
1689   }
1690   return selected;
1691
1692 }
1693
1694 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1695
1696   if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1697   Int_t iCount = 0;
1698
1699   if(!jarray){
1700     Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1701     return iCount;
1702   }
1703
1704
1705   for(int ij=0;ij<jarray->GetEntries();ij++){
1706     AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1707     if(!jet)continue;
1708     if(type==0){
1709       // No cut at all, main purpose here is sorting      
1710       list->Add(jet);
1711       iCount++;
1712     }
1713     else if(type == 1){
1714       // eta cut
1715       if(JetSelected(jet)){
1716         list->Add(jet);
1717         iCount++;
1718       }
1719     }
1720   }
1721
1722   list->Sort();
1723   return iCount;
1724
1725 }
1726
1727
1728 Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1729   if(!jets)return 0;
1730
1731   Int_t refMult = 0;
1732   for(int ij = 0;ij < jets->GetEntries();++ij){
1733     AliAODJet* jet = (AliAODJet*)jets->At(ij);
1734     if(!jet)continue;
1735     TRefArray *refs = jet->GetRefTracks();
1736     if(!refs)continue;
1737     refMult += refs->GetEntries();
1738   }
1739   return refMult;
1740
1741 }
1742
1743
1744 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1745   if(!jet)return 0;
1746   TRefArray *refs = jet->GetRefTracks();
1747   if(!refs) return 0;
1748   AliVParticle *leading = 0;
1749   Float_t fMaxPt = 0;
1750   for(int i = 0;i<refs->GetEntries();i++){
1751     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1752     if(!tmp)continue;
1753     if(tmp->Pt()>fMaxPt){
1754       leading = tmp;
1755       fMaxPt = tmp->Pt();
1756     }
1757   }
1758   return leading;
1759 }
1760
1761
1762 AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1763   if(!jet)return 0;
1764   if(!list) return 0;
1765   AliVParticle *leading = 0;
1766   Float_t fMaxPt = 0;
1767   for(int i = 0;i<list->GetEntries();i++){
1768     AliVParticle *tmp = (AliVParticle*)(list->At(i));
1769     if(!tmp)continue;
1770     if(jet->DeltaR(tmp)>r)continue;
1771     if(tmp->Pt()>fMaxPt){
1772       leading = tmp;
1773       fMaxPt = tmp->Pt();
1774     }
1775   }
1776   return leading;
1777 }
1778
1779
1780 Double_t AliAnalysisTaskJetSpectrum2::RelativePhi(Double_t mphi,Double_t vphi){
1781
1782   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1783   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1784   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1785   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1786   double dphi = mphi-vphi;
1787   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1788   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1789   return dphi;//dphi in [-Pi, Pi]
1790 }
1791
1792
1793
1794
1795
1796
1797 Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1798 {
1799   Int_t phibin=-1;
1800   if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1801   Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1802   phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1803   if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1804   return phibin;
1805 }
1806
1807 void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
1808   //
1809   // set number of trigger bins
1810   //
1811   if(n>0){
1812     fNTrigger = n;
1813     delete [] fTriggerName;
1814     fTriggerName = new TString [fNTrigger];
1815     delete [] fTriggerBit;fTriggerBit = 0;
1816     fTriggerBit = new UInt_t [fNTrigger];
1817   }
1818   else{
1819     fNTrigger = 0;
1820   }
1821 }
1822
1823
1824 void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
1825   //
1826   // set trigger bin
1827   //
1828   if(i<fNTrigger){
1829     fTriggerBit[i] = it;
1830     fTriggerName[i] =  c;     
1831   } 
1832 }
1833
1834
1835 void AliAnalysisTaskJetSpectrum2::SetNAcceptance(Int_t n){
1836   //
1837   // Set number of acceptance bins
1838   //
1839
1840   if(n>0){
1841     fNAcceptance = n;
1842     delete [] fAcceptancePhiMin; 
1843     delete [] fAcceptancePhiMax;
1844     delete [] fAcceptanceEtaMin;
1845     delete [] fAcceptanceEtaMax;
1846
1847     fAcceptancePhiMin = new Float_t[fNAcceptance];
1848     fAcceptancePhiMax = new Float_t[fNAcceptance];
1849     fAcceptanceEtaMin = new Float_t[fNAcceptance];
1850     fAcceptanceEtaMax = new Float_t[fNAcceptance];
1851   }
1852   else{
1853     fNTrigger = 0;
1854   }
1855 }
1856
1857 void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){
1858   //
1859   // Set acceptance bins
1860   //
1861
1862   if(i<fNAcceptance){
1863     Printf("%s:%d Setting acceptance %d phi %3.3f-%3.3f eta %3.3f-%3.3f",(char*)__FILE__,__LINE__,i,phiMin,phiMax,etaMin,etaMax); 
1864
1865     fAcceptancePhiMin[i] = phiMin; 
1866     fAcceptancePhiMax[i] = phiMax;
1867     fAcceptanceEtaMin[i] = etaMin;
1868     fAcceptanceEtaMax[i] = etaMax;
1869   }
1870   else{
1871     fNTrigger = 0;
1872   }
1873 }
1874
1875
1876 Int_t AliAnalysisTaskJetSpectrum2::CheckAcceptance(Float_t phi,Float_t eta){
1877
1878   //
1879   // return aceptnace bin do no use ovelapping bins
1880   //
1881
1882   for(int i = 0;i<fNAcceptance;i++){    
1883     if(phi<fAcceptancePhiMin[i])continue;
1884     if(phi>fAcceptancePhiMax[i])continue;
1885     if(eta<fAcceptanceEtaMin[i])continue;
1886     if(eta>fAcceptanceEtaMax[i])continue;
1887     return i;
1888   }
1889   // catch the rest
1890   return fNAcceptance;
1891 }
1892
1893 Float_t AliAnalysisTaskJetSpectrum2::GetRho(TList &list){
1894
1895   // invert the correction
1896   AliAODJet *jet = (AliAODJet*)list.At(0);  // highest pt jet
1897   if(!jet)return -1;
1898   if(jet->EffectiveAreaCharged()<=0)return -1;
1899   Float_t rho = jet->ChargedBgEnergy()/jet->EffectiveAreaCharged();
1900   return rho;
1901  
1902 }
1903
1904
1905
1906 AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1907   // 
1908   delete [] fTriggerBit;
1909   delete [] fTriggerName;
1910
1911   delete [] fAcceptancePhiMin; 
1912   delete [] fAcceptancePhiMax;
1913   delete [] fAcceptanceEtaMin;
1914   delete [] fAcceptanceEtaMax;
1915
1916
1917
1918 }