Add a protection to avoid crash in QA (Julian)
[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetCluster.cxx
1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets   
4 // *******************************************
5
6
7 /**************************************************************************
8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9  *                                                                        *
10  * Author: The ALICE Off-line Project.                                    *
11  * Contributors are mentioned in the code where appropriate.              *
12  *                                                                        *
13  * Permission to use, copy, modify and distribute this software and its   *
14  * documentation strictly for non-commercial purposes is hereby granted   *
15  * without fee, provided that the above copyright notice appears in all   *
16  * copies and that both the copyright notice and this permission notice   *
17  * appear in the supporting documentation. The authors make no claims     *
18  * about the suitability of this software for any purpose. It is          *
19  * provided "as is" without express or implied warranty.                  *
20  **************************************************************************/
21
22  
23 #include <TROOT.h>
24 #include <TRandom3.h>
25 #include <TSystem.h>
26 #include <TInterpreter.h>
27 #include <TChain.h>
28 #include <TRefArray.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 <TF1.h>
36 #include <TList.h>
37 #include <TLorentzVector.h>
38 #include <TClonesArray.h>
39 #include  "TDatabasePDG.h"
40
41 #include "AliAnalysisTaskJetCluster.h"
42 #include "AliAnalysisManager.h"
43 #include "AliJetFinder.h"
44 #include "AliJetHeader.h"
45 #include "AliJetReader.h"
46 #include "AliESDEvent.h"
47 #include "AliAODEvent.h"
48 #include "AliAODHandler.h"
49 #include "AliAODExtension.h"
50 #include "AliAODTrack.h"
51 #include "AliAODJet.h"
52 #include "AliAODMCParticle.h"
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
55 #include "AliStack.h"
56 #include "AliGenPythiaEventHeader.h"
57 #include "AliJetKineReaderHeader.h"
58 #include "AliGenCocktailEventHeader.h"
59 #include "AliInputEventHandler.h"
60 #include "AliAODJetEventBackground.h"
61
62 #include "fastjet/PseudoJet.hh"
63 #include "fastjet/ClusterSequenceArea.hh"
64 #include "fastjet/AreaDefinition.hh"
65 #include "fastjet/JetDefinition.hh"
66 // get info on how fastjet was configured
67 #include "fastjet/config.h"
68
69 using std::vector;
70
71 ClassImp(AliAnalysisTaskJetCluster)
72
73 AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
74   //
75   // Destructor
76   //
77
78   delete fRef;
79   delete fRandom;
80
81   if(fTCAJetsOut)fTCAJetsOut->Delete();
82   delete fTCAJetsOut;
83   
84   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
85   delete fTCAJetsOutRan;
86   
87   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
88   delete fTCARandomConesOut;
89   
90   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
91   delete fTCARandomConesOutRan;
92   
93   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
94   delete fAODJetBackgroundOut;
95 }
96
97 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): 
98   AliAnalysisTaskSE(),
99   fAOD(0x0),
100   fAODExtension(0x0),
101   fRef(new TRefArray),
102   fUseAODTrackInput(kFALSE),
103   fUseAODMCInput(kFALSE),
104   fUseBackgroundCalc(kFALSE),
105   fEventSelection(kFALSE),     
106   fRequireVZEROAC(kFALSE),     
107   fRequireTZEROvtx(kFALSE),
108   fFilterMask(0),
109   fFilterMaskBestPt(0),
110   fFilterType(0),
111   fJetTypes(kJet),
112   fTrackTypeRec(kTrackUndef),
113   fTrackTypeGen(kTrackUndef),  
114   fNSkipLeadingRan(0),
115   fNSkipLeadingCone(0),
116   fNRandomCones(0),
117   fAvgTrials(1),
118   fExternalWeight(1),
119   fTrackEtaWindow(0.9),    
120   fRecEtaWindow(0.5),
121   fTrackPtCut(0.),                                                      
122   fJetOutputMinPt(0.150),
123   fMaxTrackPtInJet(100.),
124   fJetTriggerPtCut(0),
125   fVtxZCut(8),
126   fVtxR2Cut(1),
127   fCentCutUp(0),
128   fCentCutLo(0),
129   fStoreRhoLeadingTrackCorr(kFALSE),
130   fNonStdBranch(""),
131   fBackgroundBranch(""),
132   fNonStdFile(""),
133   fMomResH1(0x0),
134   fMomResH2(0x0),
135   fMomResH3(0x0),
136   fMomResH1Fit(0x0),
137   fMomResH2Fit(0x0),
138   fMomResH3Fit(0x0),
139   fhEffH1(0x0),
140   fhEffH2(0x0),
141   fhEffH3(0x0),
142   fUseTrPtResolutionSmearing(kFALSE),
143   fUseDiceEfficiency(kFALSE),
144   fDiceEfficiencyMinPt(-1.),
145   fUseTrPtResolutionFromOADB(kFALSE),
146   fUseTrEfficiencyFromOADB(kFALSE),
147   fPathTrPtResolution(""),
148   fPathTrEfficiency(""),
149   fChangeEfficiencyFraction(0.),
150   fEfficiencyFixed(1.),
151   fRparam(1.0), 
152   fAlgorithm(fastjet::kt_algorithm),
153   fStrategy(fastjet::Best),
154   fRecombScheme(fastjet::BIpt_scheme),
155   fAreaType(fastjet::active_area), 
156   fGhostArea(0.01),
157   fActiveAreaRepeats(1),
158   fGhostEtamax(1.5),
159   fTCAJetsOut(0x0),
160   fTCAJetsOutRan(0x0),
161   fTCARandomConesOut(0x0),
162   fTCARandomConesOutRan(0x0),
163   fAODJetBackgroundOut(0x0),
164   fRandom(0),
165   fh1Xsec(0x0),
166   fh1Trials(0x0),
167   fh1PtHard(0x0),
168   fh1PtHardNoW(0x0),  
169   fh1PtHardTrials(0x0),
170   fh1NJetsRec(0x0),
171   fh1NConstRec(0x0),
172   fh1NConstLeadingRec(0x0),
173   fh1PtJetsRecIn(0x0),
174   fh1PtJetsLeadingRecIn(0x0),
175   fh1PtJetConstRec(0x0),
176   fh1PtJetConstLeadingRec(0x0),
177   fh1PtTracksRecIn(0x0),
178   fh1PtTracksLeadingRecIn(0x0),
179   fh1NJetsRecRan(0x0),
180   fh1NConstRecRan(0x0),
181   fh1PtJetsLeadingRecInRan(0x0),
182   fh1NConstLeadingRecRan(0x0),
183   fh1PtJetsRecInRan(0x0),
184   fh1PtTracksGenIn(0x0),
185   fh1Nch(0x0),
186   fh1CentralityPhySel(0x0), 
187   fh1Centrality(0x0), 
188   fh1CentralitySelect(0x0),
189   fh1ZPhySel(0x0), 
190   fh1Z(0x0), 
191   fh1ZSelect(0x0),
192   fh2NRecJetsPt(0x0),
193   fh2NRecTracksPt(0x0),
194   fh2NConstPt(0x0),
195   fh2NConstLeadingPt(0x0),
196   fh2JetPhiEta(0x0),
197   fh2LeadingJetPhiEta(0x0),
198   fh2JetEtaPt(0x0),
199   fh2LeadingJetEtaPt(0x0),
200   fh2TrackEtaPt(0x0),
201   fh2LeadingTrackEtaPt(0x0),
202   fh2JetsLeadingPhiEta(0x0),
203   fh2JetsLeadingPhiPt(0x0),
204   fh2TracksLeadingPhiEta(0x0),
205   fh2TracksLeadingPhiPt(0x0),
206   fh2TracksLeadingJetPhiPt(0x0),
207   fh2JetsLeadingPhiPtW(0x0),
208   fh2TracksLeadingPhiPtW(0x0),
209   fh2TracksLeadingJetPhiPtW(0x0),
210   fh2NRecJetsPtRan(0x0),
211   fh2NConstPtRan(0x0),
212   fh2NConstLeadingPtRan(0x0),
213   fh2PtNch(0x0),
214   fh2PtNchRan(0x0),
215   fh2PtNchN(0x0),
216   fh2PtNchNRan(0x0),
217   fh2TracksLeadingJetPhiPtRan(0x0),
218   fh2TracksLeadingJetPhiPtWRan(0x0),
219   fh3CentvsRhoLeadingTrackPt(0x0),
220   fh3CentvsSigmaLeadingTrackPt(0x0),
221   fh3MultvsRhoLeadingTrackPt(0x0),
222   fh3MultvsSigmaLeadingTrackPt(0x0),
223   fh3CentvsRhoLeadingTrackPtQ1(0x0),
224   fh3CentvsRhoLeadingTrackPtQ2(0x0),
225   fh3CentvsRhoLeadingTrackPtQ3(0x0),
226   fh3CentvsRhoLeadingTrackPtQ4(0x0),
227   fh3CentvsSigmaLeadingTrackPtQ1(0x0),
228   fh3CentvsSigmaLeadingTrackPtQ2(0x0),
229   fh3CentvsSigmaLeadingTrackPtQ3(0x0),
230   fh3CentvsSigmaLeadingTrackPtQ4(0x0),
231   fh3MultvsRhoLeadingTrackPtQ1(0x0),
232   fh3MultvsRhoLeadingTrackPtQ2(0x0),
233   fh3MultvsRhoLeadingTrackPtQ3(0x0),
234   fh3MultvsRhoLeadingTrackPtQ4(0x0),
235   fh3MultvsSigmaLeadingTrackPtQ1(0x0),
236   fh3MultvsSigmaLeadingTrackPtQ2(0x0),
237   fh3MultvsSigmaLeadingTrackPtQ3(0x0),
238   fh3MultvsSigmaLeadingTrackPtQ4(0x0),
239   fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
240   fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
241   fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
242   fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
243   fh2PtGenPtSmeared(0x0),
244   fp1Efficiency(0x0),
245   fp1PtResolution(0x0),
246   fHistList(0x0)  
247 {
248   //
249   // Constructor
250   //
251
252   for(int i = 0;i<3;i++){
253     fh1BiARandomCones[i] = 0;
254     fh1BiARandomConesRan[i] = 0;
255   }
256   for(int i = 0;i<kMaxCent;i++){
257     fh2JetsLeadingPhiPtC[i] = 0;     
258     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
259     fh2TracksLeadingJetPhiPtC[i] = 0;
260     fh2TracksLeadingJetPhiPtWC[i] = 0;
261   }
262 }
263
264 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
265   AliAnalysisTaskSE(name),
266   fAOD(0x0),
267   fAODExtension(0x0),
268   fRef(new TRefArray),
269   fUseAODTrackInput(kFALSE),
270   fUseAODMCInput(kFALSE),
271   fUseBackgroundCalc(kFALSE),
272   fEventSelection(kFALSE),
273   fRequireVZEROAC(kFALSE),     
274   fRequireTZEROvtx(kFALSE), 
275   fFilterMask(0),
276   fFilterMaskBestPt(0),
277   fFilterType(0),
278   fJetTypes(kJet),
279   fTrackTypeRec(kTrackUndef),
280   fTrackTypeGen(kTrackUndef),
281   fNSkipLeadingRan(0),
282   fNSkipLeadingCone(0),
283   fNRandomCones(0),
284   fAvgTrials(1),
285   fExternalWeight(1),    
286   fTrackEtaWindow(0.9),    
287   fRecEtaWindow(0.5),
288   fTrackPtCut(0.),                                                      
289   fJetOutputMinPt(0.150),
290   fMaxTrackPtInJet(100.),
291   fJetTriggerPtCut(0),
292   fVtxZCut(8),
293   fVtxR2Cut(1),
294   fCentCutUp(0),
295   fCentCutLo(0),
296   fStoreRhoLeadingTrackCorr(kFALSE),
297   fNonStdBranch(""),
298   fBackgroundBranch(""),
299   fNonStdFile(""),
300   fMomResH1(0x0),
301   fMomResH2(0x0),
302   fMomResH3(0x0),
303   fMomResH1Fit(0x0),
304   fMomResH2Fit(0x0),
305   fMomResH3Fit(0x0),
306   fhEffH1(0x0),
307   fhEffH2(0x0),
308   fhEffH3(0x0),
309   fUseTrPtResolutionSmearing(kFALSE),
310   fUseDiceEfficiency(kFALSE),
311   fDiceEfficiencyMinPt(-1.),
312   fUseTrPtResolutionFromOADB(kFALSE),
313   fUseTrEfficiencyFromOADB(kFALSE),
314   fPathTrPtResolution(""),
315   fPathTrEfficiency(""),
316   fChangeEfficiencyFraction(0.),
317   fEfficiencyFixed(1.),
318   fRparam(1.0), 
319   fAlgorithm(fastjet::kt_algorithm),
320   fStrategy(fastjet::Best),
321   fRecombScheme(fastjet::BIpt_scheme),
322   fAreaType(fastjet::active_area), 
323   fGhostArea(0.01),
324   fActiveAreaRepeats(1),
325   fGhostEtamax(1.5),
326   fTCAJetsOut(0x0),
327   fTCAJetsOutRan(0x0),
328   fTCARandomConesOut(0x0),
329   fTCARandomConesOutRan(0x0),
330   fAODJetBackgroundOut(0x0),
331   fRandom(0),
332   fh1Xsec(0x0),
333   fh1Trials(0x0),
334   fh1PtHard(0x0),
335   fh1PtHardNoW(0x0),  
336   fh1PtHardTrials(0x0),
337   fh1NJetsRec(0x0),
338   fh1NConstRec(0x0),
339   fh1NConstLeadingRec(0x0),
340   fh1PtJetsRecIn(0x0),
341   fh1PtJetsLeadingRecIn(0x0),
342   fh1PtJetConstRec(0x0),
343   fh1PtJetConstLeadingRec(0x0),
344   fh1PtTracksRecIn(0x0),
345   fh1PtTracksLeadingRecIn(0x0),
346   fh1NJetsRecRan(0x0),
347   fh1NConstRecRan(0x0),
348   fh1PtJetsLeadingRecInRan(0x0),
349   fh1NConstLeadingRecRan(0x0),
350   fh1PtJetsRecInRan(0x0),
351   fh1PtTracksGenIn(0x0),
352   fh1Nch(0x0),
353   fh1CentralityPhySel(0x0), 
354   fh1Centrality(0x0), 
355   fh1CentralitySelect(0x0),
356   fh1ZPhySel(0x0), 
357   fh1Z(0x0), 
358   fh1ZSelect(0x0),
359   fh2NRecJetsPt(0x0),
360   fh2NRecTracksPt(0x0),
361   fh2NConstPt(0x0),
362   fh2NConstLeadingPt(0x0),
363   fh2JetPhiEta(0x0),
364   fh2LeadingJetPhiEta(0x0),
365   fh2JetEtaPt(0x0),
366   fh2LeadingJetEtaPt(0x0),
367   fh2TrackEtaPt(0x0),
368   fh2LeadingTrackEtaPt(0x0),
369   fh2JetsLeadingPhiEta(0x0),
370   fh2JetsLeadingPhiPt(0x0),
371   fh2TracksLeadingPhiEta(0x0),
372   fh2TracksLeadingPhiPt(0x0),
373   fh2TracksLeadingJetPhiPt(0x0),
374   fh2JetsLeadingPhiPtW(0x0),
375   fh2TracksLeadingPhiPtW(0x0),
376   fh2TracksLeadingJetPhiPtW(0x0),
377   fh2NRecJetsPtRan(0x0),
378   fh2NConstPtRan(0x0),
379   fh2NConstLeadingPtRan(0x0),
380   fh2PtNch(0x0),
381   fh2PtNchRan(0x0),
382   fh2PtNchN(0x0),
383   fh2PtNchNRan(0x0),
384   fh2TracksLeadingJetPhiPtRan(0x0),
385   fh2TracksLeadingJetPhiPtWRan(0x0),
386   fh3CentvsRhoLeadingTrackPt(0x0),
387   fh3CentvsSigmaLeadingTrackPt(0x0),
388   fh3MultvsRhoLeadingTrackPt(0x0),
389   fh3MultvsSigmaLeadingTrackPt(0x0),
390   fh3CentvsRhoLeadingTrackPtQ1(0x0),
391   fh3CentvsRhoLeadingTrackPtQ2(0x0),
392   fh3CentvsRhoLeadingTrackPtQ3(0x0),
393   fh3CentvsRhoLeadingTrackPtQ4(0x0),
394   fh3CentvsSigmaLeadingTrackPtQ1(0x0),
395   fh3CentvsSigmaLeadingTrackPtQ2(0x0),
396   fh3CentvsSigmaLeadingTrackPtQ3(0x0),
397   fh3CentvsSigmaLeadingTrackPtQ4(0x0),
398   fh3MultvsRhoLeadingTrackPtQ1(0x0),
399   fh3MultvsRhoLeadingTrackPtQ2(0x0),
400   fh3MultvsRhoLeadingTrackPtQ3(0x0),
401   fh3MultvsRhoLeadingTrackPtQ4(0x0),
402   fh3MultvsSigmaLeadingTrackPtQ1(0x0),
403   fh3MultvsSigmaLeadingTrackPtQ2(0x0),
404   fh3MultvsSigmaLeadingTrackPtQ3(0x0),
405   fh3MultvsSigmaLeadingTrackPtQ4(0x0),
406   fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0),
407   fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0),
408   fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0),
409   fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0),
410   fh2PtGenPtSmeared(0x0),
411   fp1Efficiency(0x0),
412   fp1PtResolution(0x0),
413   fHistList(0x0)
414 {
415   //
416   // named ctor
417   //
418
419   for(int i = 0;i<3;i++){
420     fh1BiARandomCones[i] = 0;
421     fh1BiARandomConesRan[i] = 0;
422   }
423   for(int i = 0;i<kMaxCent;i++){
424     fh2JetsLeadingPhiPtC[i] = 0;     
425     fh2JetsLeadingPhiPtWC[i] = 0;      //! jet correlation with leading jet
426     fh2TracksLeadingJetPhiPtC[i] = 0;
427     fh2TracksLeadingJetPhiPtWC[i] = 0;
428   }
429  DefineOutput(1, TList::Class());  
430 }
431
432
433
434 Bool_t AliAnalysisTaskJetCluster::Notify()
435 {
436   //
437   // Implemented Notify() to read the cross sections
438   // and number of trials from pyxsec.root
439   // 
440   return kTRUE;
441 }
442
443 void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
444 {
445
446   //
447   // Create the output container
448   //
449
450   fRandom = new TRandom3(0);
451
452
453   // Connect the AOD
454
455
456   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
457
458   
459
460   if(fNonStdBranch.Length()!=0)
461     {
462       // only create the output branch if we have a name
463       // Create a new branch for jets...
464       //  -> cleared in the UserExec....
465       // here we can also have the case that the brnaches are written to a separate file
466       
467       if(fJetTypes&kJet){
468         fTCAJetsOut = new TClonesArray("AliAODJet", 0);
469         fTCAJetsOut->SetName(fNonStdBranch.Data());
470         AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
471       }
472
473       if(fJetTypes&kJetRan){
474         fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
475         fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
476         if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
477           if(  fEfficiencyFixed < 1.)
478             fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dEffFixed%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
479           else
480             fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
481         }
482         AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
483       }
484
485       if(fUseBackgroundCalc){
486         if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
487           fAODJetBackgroundOut = new AliAODJetEventBackground();
488           fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
489           if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
490             if(  fEfficiencyFixed < 1.)
491               fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dEffFixed%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.)));
492             else
493               fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
494           }
495           AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());  
496         }
497       }
498       // create the branch for the random cones with the same R 
499       TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
500       if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) {
501         if(  fEfficiencyFixed < 1.)
502           cName = Form("%sDetector%d%dEffFixed%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.),fNSkipLeadingCone);
503         else
504           cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
505       }
506       if(fNRandomCones>0){
507         if(fJetTypes&kRC){
508           if(!AODEvent()->FindListObject(cName.Data())){
509             fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
510             fTCARandomConesOut->SetName(cName.Data());
511             AddAODBranch("TClonesArray",&fTCARandomConesOut,fNonStdFile.Data());
512           }
513         }
514         // create the branch with the random for the random cones on the random event
515         if(fJetTypes&kRCRan){
516           cName = Form("%sRandomCone_random",fNonStdBranch.Data());
517           if(!AODEvent()->FindListObject(cName.Data())){
518             fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
519             fTCARandomConesOutRan->SetName(cName.Data());
520             AddAODBranch("TClonesArray",&fTCARandomConesOutRan,fNonStdFile.Data());
521           }
522         }
523       }
524     
525       if(fNonStdFile.Length()!=0){
526         // 
527         // case that we have an AOD extension we need to fetch the jets from the extended output
528         // we identify the extension aod event by looking for the branchname
529         AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
530         // case that we have an AOD extension we need can fetch the background maybe from the extended output                                                                  
531         fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
532       }
533     }
534
535
536   if(!fHistList)fHistList = new TList();
537   fHistList->SetOwner();
538   PostData(1, fHistList); // post data in any case once
539
540   Bool_t oldStatus = TH1::AddDirectoryStatus();
541   TH1::AddDirectory(kFALSE);
542
543   //
544   //  Histogram
545     
546   const Int_t nBinPt = 100;
547   Double_t binLimitsPt[nBinPt+1];
548   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
549     if(iPt == 0){
550       binLimitsPt[iPt] = 0.0;
551     }
552     else {// 1.0
553       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.0;
554     }
555   }
556   
557   const Int_t nBinPhi = 90;
558   Double_t binLimitsPhi[nBinPhi+1];
559   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
560     if(iPhi==0){
561       binLimitsPhi[iPhi] = -1.*TMath::Pi();
562     }
563     else{
564       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
565     }
566   }
567
568
569
570   const Int_t nBinEta = 40;
571   Double_t binLimitsEta[nBinEta+1];
572   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
573     if(iEta==0){
574       binLimitsEta[iEta] = -2.0;
575     }
576     else{
577       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
578     }
579   }
580
581   const int nChMax = 5000;
582
583   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
584   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
585
586   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
587   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
588
589
590   fh1NJetsRec = new TH1F("fh1NJetsRec","N reconstructed jets",120,-0.5,119.5);
591   fh1NJetsRecRan = new TH1F("fh1NJetsRecRan","N reconstructed jets",120,-0.5,119.5);
592
593   fh1NConstRec = new TH1F("fh1NConstRec","# jet constituents",120,-0.5,119.5);
594   fh1NConstRecRan = new TH1F("fh1NConstRecRan","# jet constituents",120,-0.5,119.5);
595   fh1NConstLeadingRec = new TH1F("fh1NConstLeadingRec","jet constituents",120,-0.5,119.5);
596   fh1NConstLeadingRecRan = new TH1F("fh1NConstLeadingRecRan","jet constituents",120,-0.5,119.5);
597
598
599   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
600   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
601   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
602
603   fh1PtJetsRecIn  = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
604   fh1PtJetsRecInRan  = new TH1F("fh1PtJetsRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
605   fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
606   fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
607   fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
608   fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
609   fh1PtTracksRecIn  = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
610   fh1PtTracksLeadingRecIn  = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
611   fh1PtTracksGenIn  = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
612   fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
613
614   fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
615   fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
616   fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
617
618   fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
619   fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
620   fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
621
622   fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
623   fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
624   fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
625   // 
626
627
628   fh2NConstPt = new TH2F("fh2NConstPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
629   fh2NConstLeadingPt = new TH2F("fh2NConstLeadingPt","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
630   fh2NConstPtRan = new TH2F("fh2NConstPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
631   fh2NConstLeadingPtRan = new TH2F("fh2NConstLeadingPtRan","Number of constituents ;p_{T} (GeV/c);N",nBinPt,binLimitsPt,50,-0.5,49.5);
632
633   fh2PtNch = new TH2F("fh2PtNch","p_T of cluster vs. multiplicity; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
634   fh2PtNchRan = new TH2F("fh2PtNchRan","p_T of cluster vs. multiplicity ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
635   fh2PtNchN = new TH2F("fh2PtNchN","p_T of cluster vs. multiplicity N weighted; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
636   fh2PtNchNRan = new TH2F("fh2PtNchNRan","p_T of cluster vs. multiplicity N weighted ran; N_{ch};p_{T} (GeV/c);",nChMax,-0.5,nChMax-0.5,nBinPt,binLimitsPt);
637
638
639
640   fh2JetPhiEta  = new TH2F("fh2JetPhiEta","eta vs phi all jets;#phi;#eta",
641                            nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
642   fh2LeadingJetPhiEta  = new TH2F("fh2LeadingJetPhiEta","eta vs phi leading jets;#phi;#eta",
643                                   nBinPhi,0.,2.*TMath::Pi(),nBinEta,binLimitsEta);
644
645   fh2JetEtaPt  = new TH2F("fh2JetEtaPt","pt vs eta all jets;#eta;p_{T}",
646                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
647   fh2LeadingJetEtaPt  = new TH2F("fh2LeadingJetEtaPt","pT vs eta leading jets;#eta;p_{T}",
648                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
649
650   fh2TrackEtaPt  = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}",
651                           nBinEta,binLimitsEta,nBinPt,binLimitsPt);
652   fh2LeadingTrackEtaPt  = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}",
653                                  nBinEta,binLimitsEta,nBinPt,binLimitsPt);
654
655
656
657   fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
658                                 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
659   fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
660                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
661   fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
662                                     nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
663   fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
664                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
665   fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
666                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
667   fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
668                                  nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
669
670   fh2JetsLeadingPhiPtW      = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)",
671                                 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
672   fh2TracksLeadingPhiPtW    = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)",
673                                     nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
674
675   fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
676                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
677   fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
678                                        nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
679
680   if(fStoreRhoLeadingTrackCorr) {
681     fh3CentvsRhoLeadingTrackPt = new TH3F("fh3CentvsRhoLeadingTrackPt","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.,100,0.,100.);
682     fh3CentvsSigmaLeadingTrackPt = new TH3F("fh3CentvsSigmaLeadingTrackPt","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.,100,0.,100.);
683     fh3MultvsRhoLeadingTrackPt = new TH3F("fh3MultvsRhoLeadingTrackPt","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.,100,0.,100.);
684     fh3MultvsSigmaLeadingTrackPt = new TH3F("fh3MultvsSigmaLeadingTrackPt","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.,100,0.,100.);
685
686
687     fh3CentvsRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
688     fh3CentvsRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
689     fh3CentvsRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
690     fh3CentvsRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
691
692     fh3CentvsSigmaLeadingTrackPtQ1 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ1","centrality vs background density Q1; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
693     fh3CentvsSigmaLeadingTrackPtQ2 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ2","centrality vs background density Q2; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
694     fh3CentvsSigmaLeadingTrackPtQ3 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ3","centrality vs background density Q3; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
695     fh3CentvsSigmaLeadingTrackPtQ4 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ4","centrality vs background density Q4; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
696
697     fh3MultvsRhoLeadingTrackPtQ1 = new TH3F("fh3MultvsRhoLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
698     fh3MultvsRhoLeadingTrackPtQ2 = new TH3F("fh3MultvsRhoLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
699     fh3MultvsRhoLeadingTrackPtQ3 = new TH3F("fh3MultvsRhoLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
700     fh3MultvsRhoLeadingTrackPtQ4 = new TH3F("fh3MultvsRhoLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
701
702     fh3MultvsSigmaLeadingTrackPtQ1 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
703     fh3MultvsSigmaLeadingTrackPtQ2 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
704     fh3MultvsSigmaLeadingTrackPtQ3 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
705     fh3MultvsSigmaLeadingTrackPtQ4 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
706
707
708     fh3CentvsDeltaRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
709     fh3CentvsDeltaRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
710     fh3CentvsDeltaRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
711     fh3CentvsDeltaRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.);
712
713     fHistList->Add(fh3CentvsRhoLeadingTrackPt);
714     fHistList->Add(fh3CentvsSigmaLeadingTrackPt);
715     fHistList->Add(fh3MultvsRhoLeadingTrackPt);
716     fHistList->Add(fh3MultvsSigmaLeadingTrackPt);
717
718     fHistList->Add(fh3CentvsRhoLeadingTrackPtQ1);
719     fHistList->Add(fh3CentvsRhoLeadingTrackPtQ2);
720     fHistList->Add(fh3CentvsRhoLeadingTrackPtQ3);
721     fHistList->Add(fh3CentvsRhoLeadingTrackPtQ4);
722
723     fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ1);
724     fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ2);
725     fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ3);
726     fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ4);
727
728     fHistList->Add(fh3MultvsRhoLeadingTrackPtQ1);
729     fHistList->Add(fh3MultvsRhoLeadingTrackPtQ2);
730     fHistList->Add(fh3MultvsRhoLeadingTrackPtQ3);
731     fHistList->Add(fh3MultvsRhoLeadingTrackPtQ4);
732
733     fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ1);
734     fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ2);
735     fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ3);
736     fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ4);
737
738     fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ1);
739     fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ2);
740     fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ3);
741     fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ4);
742
743   }
744   
745   //Detector level effects histos
746   fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
747
748   fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
749   fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
750
751   fHistList->Add(fh2PtGenPtSmeared);
752   fHistList->Add(fp1Efficiency);
753   fHistList->Add(fp1PtResolution);
754
755   if(fNRandomCones>0&&fUseBackgroundCalc){
756     for(int i = 0;i<3;i++){
757       fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
758       fh1BiARandomConesRan[i] =  new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
759     }
760   }
761
762   for(int i = 0;i < kMaxCent;i++){
763     fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
764     fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
765     fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
766     fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
767   }
768
769   const Int_t saveLevel = 3; // large save level more histos
770   if(saveLevel>0){
771     fHistList->Add(fh1Xsec);
772     fHistList->Add(fh1Trials);
773
774     fHistList->Add(fh1NJetsRec);
775     fHistList->Add(fh1NConstRec);
776     fHistList->Add(fh1NConstLeadingRec);
777     fHistList->Add(fh1PtJetsRecIn);
778     fHistList->Add(fh1PtJetsLeadingRecIn);
779     fHistList->Add(fh1PtTracksRecIn);
780     fHistList->Add(fh1PtTracksLeadingRecIn);
781     fHistList->Add(fh1PtJetConstRec);
782     fHistList->Add(fh1PtJetConstLeadingRec);
783     fHistList->Add(fh1NJetsRecRan);
784     fHistList->Add(fh1NConstRecRan);
785     fHistList->Add(fh1PtJetsLeadingRecInRan);
786     fHistList->Add(fh1NConstLeadingRecRan);
787     fHistList->Add(fh1PtJetsRecInRan);
788     fHistList->Add(fh1Nch);
789     fHistList->Add(fh1Centrality);
790     fHistList->Add(fh1CentralitySelect);
791     fHistList->Add(fh1CentralityPhySel);
792     fHistList->Add(fh1Z);
793     fHistList->Add(fh1ZSelect);
794     fHistList->Add(fh1ZPhySel);
795     if(fNRandomCones>0&&fUseBackgroundCalc){
796       for(int i = 0;i<3;i++){
797         fHistList->Add(fh1BiARandomCones[i]);
798         fHistList->Add(fh1BiARandomConesRan[i]);
799       }
800     }
801   for(int i = 0;i < kMaxCent;i++){
802     fHistList->Add(fh2JetsLeadingPhiPtC[i]);
803     fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
804     fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
805     fHistList->Add(fh2TracksLeadingJetPhiPtWC[i]);
806   }
807
808     fHistList->Add(fh2NRecJetsPt);
809     fHistList->Add(fh2NRecTracksPt);
810     fHistList->Add(fh2NConstPt);
811     fHistList->Add(fh2NConstLeadingPt);
812     fHistList->Add(fh2PtNch);
813     fHistList->Add(fh2PtNchRan);
814     fHistList->Add(fh2PtNchN);
815     fHistList->Add(fh2PtNchNRan);
816     fHistList->Add(fh2JetPhiEta);
817     fHistList->Add(fh2LeadingJetPhiEta);
818     fHistList->Add(fh2JetEtaPt);
819     fHistList->Add(fh2LeadingJetEtaPt);
820     fHistList->Add(fh2TrackEtaPt);
821     fHistList->Add(fh2LeadingTrackEtaPt);
822     fHistList->Add(fh2JetsLeadingPhiEta );
823     fHistList->Add(fh2JetsLeadingPhiPt);
824     fHistList->Add(fh2TracksLeadingPhiEta);
825     fHistList->Add(fh2TracksLeadingPhiPt);
826     fHistList->Add(fh2TracksLeadingJetPhiPt);
827     fHistList->Add(fh2JetsLeadingPhiPtW);
828     fHistList->Add(fh2TracksLeadingPhiPtW);
829     fHistList->Add(fh2TracksLeadingJetPhiPtW);
830     fHistList->Add(fh2NRecJetsPtRan);
831     fHistList->Add(fh2NConstPtRan);
832     fHistList->Add(fh2NConstLeadingPtRan);
833     fHistList->Add(fh2TracksLeadingJetPhiPtRan);
834     fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
835   }
836
837   // =========== Switch on Sumw2 for all histos ===========
838   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
839     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
840     if (h1){
841       h1->Sumw2();
842       continue;
843     }
844     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
845     if(hn)hn->Sumw2();
846   }
847   TH1::AddDirectory(oldStatus);
848 }
849
850 void AliAnalysisTaskJetCluster::LocalInit()
851 {
852   //
853   // Initialization
854   //
855
856   if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
857
858   if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
859   if(fUseTrEfficiencyFromOADB)   LoadTrEfficiencyRootFileFromOADB();
860
861
862   FitMomentumResolution();
863
864 }
865
866 void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
867 {
868
869   // handle and reset the output jet branch 
870
871   if(fTCAJetsOut)fTCAJetsOut->Delete();
872   if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
873   if(fTCARandomConesOut)fTCARandomConesOut->Delete();
874   if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
875   if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
876
877   AliAODJetEventBackground* externalBackground = 0;
878   if(!externalBackground&&fBackgroundBranch.Length()){
879     externalBackground =  (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
880     if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
881     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
882   }
883   //
884   // Execute analysis for current event
885   //
886   AliESDEvent *fESD = 0;
887   if(fUseAODTrackInput){    
888     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
889     if(!fAOD){
890       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
891       return;
892     }
893     // fetch the header
894   }
895   else{
896     //  assume that the AOD is in the general output...
897     fAOD  = AODEvent();
898     if(!fAOD){
899       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
900       return;
901     }
902     if(fDebug>0){
903       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
904     }
905   }
906
907   //Check if information is provided detector level effects
908   if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
909   if(!fhEffH1 || !fhEffH2 || !fhEffH3 )      fUseDiceEfficiency = kFALSE;
910   if(  fEfficiencyFixed < 1. )               fUseDiceEfficiency = kTRUE;  
911
912   Bool_t selectEvent =  false;
913   Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
914
915   Float_t cent = 0;
916   Float_t zVtx  = 0;
917   Int_t cenClass = -1;
918   if(fAOD){
919     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
920     TString vtxTitle(vtxAOD->GetTitle());
921     zVtx = vtxAOD->GetZ();
922
923     cent = fAOD->GetHeader()->GetCentrality();
924     if(cent<10)cenClass = 0;
925     else if(cent<30)cenClass = 1;
926     else if(cent<50)cenClass = 2;
927     else if(cent<80)cenClass = 3;
928     if(physicsSelection){
929       fh1CentralityPhySel->Fill(cent);
930       fh1ZPhySel->Fill(zVtx);
931     }
932
933     if(fEventSelection){
934       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
935         Float_t yvtx = vtxAOD->GetY();
936         Float_t xvtx = vtxAOD->GetX();
937         Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
938         if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
939           if(physicsSelection){
940             selectEvent = true;
941           }
942         }
943       }
944       if(fCentCutUp>0){
945         if(cent<fCentCutLo||cent>fCentCutUp){
946           selectEvent = false;
947         }
948       }
949     }else{
950       selectEvent = true;
951     }
952   }
953
954
955   Bool_t T0 = false;
956   Bool_t V0 = false;
957   const AliAODVZERO  *vzero = fAOD->GetVZEROData();
958   if(vzero){
959     if((vzero->GetTriggerChargeA()>0)&&(vzero->GetTriggerChargeC()>0)){
960       V0 = true;
961     }
962   }
963   
964   const AliAODTZERO  *tzero = fAOD->GetTZEROData();
965   if(tzero){
966     if(TMath::Abs(tzero->GetT0VertexRaw())<100){
967       T0 = true;
968     }
969   }
970   
971   if(fRequireVZEROAC&&fRequireTZEROvtx)selectEvent = selectEvent&&V0&&T0;
972   else if(fRequireTZEROvtx)selectEvent = selectEvent&&T0;
973   else if(fRequireVZEROAC)selectEvent = selectEvent&&V0;
974
975
976   if(!selectEvent){
977     PostData(1, fHistList);
978     return;
979   }
980   fh1Centrality->Fill(cent);  
981   fh1Z->Fill(zVtx);
982   fh1Trials->Fill("#sum{ntrials}",1);
983   
984
985   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
986
987   // ==== General variables needed
988
989
990
991   // we simply fetch the tracks/mc particles as a list of AliVParticles
992
993   TList recParticles;
994   TList genParticles;
995
996   Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
997   Float_t nCh = recParticles.GetEntries(); 
998   fh1Nch->Fill(nCh);
999   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
1000   nT = GetListOfTracks(&genParticles,fTrackTypeGen);
1001   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
1002
1003   // find the jets....
1004
1005   vector<fastjet::PseudoJet> inputParticlesRec;
1006   vector<fastjet::PseudoJet> inputParticlesRecRan;
1007   
1008   // Generate the random cones
1009   
1010   AliAODJet vTmpRan(1,0,0,1);
1011   for(int i = 0; i < recParticles.GetEntries(); i++){
1012     AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1013
1014     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
1015     // we take total momentum here
1016
1017     if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1018       //Add particles to fastjet in case we are not running toy model
1019       fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1020       jInp.set_user_index(i);
1021       inputParticlesRec.push_back(jInp);
1022     }
1023     else if(fUseDiceEfficiency) {
1024
1025       // Dice to decide if particle is kept or not - toy  model for efficiency
1026       //
1027       Double_t sumEff = 0.;
1028       Double_t pT = 0.;
1029       Double_t eff[3] = {0.};
1030       Int_t cat[3] = {0};
1031
1032       Double_t rnd = fRandom->Uniform(1.);
1033       if(  fEfficiencyFixed<1. ) {
1034         sumEff = fEfficiencyFixed;
1035       } else {
1036
1037         pT = vp->Pt();
1038         Double_t pTtmp = pT;
1039         if(pT>10.) pTtmp = 10.;
1040         Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
1041         Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
1042         Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
1043
1044         //Sort efficiencies from large to small
1045         if(eff1>eff2 && eff1>eff3) { 
1046           eff[0] = eff1; 
1047           cat[0] = 1;
1048           if(eff2>eff3) {
1049             eff[1] = eff2;
1050             eff[2] = eff3; 
1051             cat[1] = 2;
1052             cat[2] = 3;
1053           } else {
1054             eff[1] = eff3;
1055             eff[2] = eff2; 
1056             cat[1] = 3;
1057             cat[2] = 2;
1058           }
1059         }
1060         else if(eff2>eff1 && eff2>eff3) {
1061           eff[0] = eff2;
1062           cat[0] = 2;
1063           if(eff1>eff3) {
1064             eff[1] = eff1;
1065             eff[2] = eff3; 
1066             cat[1] = 1;
1067             cat[2] = 3;
1068           } else {
1069             eff[1] = eff3;
1070             eff[2] = eff1; 
1071             cat[1] = 3;
1072             cat[2] = 1;
1073           }
1074         }
1075         else if(eff3>eff1 && eff3>eff2) {
1076           eff[0] = eff3;
1077           cat[0] = 3;
1078           if(eff1>eff2) {
1079             eff[1] = eff1;
1080             eff[2] = eff2; 
1081             cat[1] = 1;
1082             cat[2] = 2;
1083           } else {
1084             eff[1] = eff2;
1085             eff[2] = eff1; 
1086             cat[1] = 2;
1087             cat[2] = 1;
1088           }
1089         }
1090         
1091         sumEff = eff[0]+eff[1]+eff[2];
1092       }
1093       fp1Efficiency->Fill(vp->Pt(),sumEff);
1094       if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
1095
1096       if(fUseTrPtResolutionSmearing) {
1097         //Smear momentum of generated particle
1098         Double_t smear = 1.;
1099         //Select hybrid track category
1100         if(rnd<=eff[2]) 
1101           smear = GetMomentumSmearing(cat[2],pT);
1102         else if(rnd<=(eff[2]+eff[1])) 
1103           smear = GetMomentumSmearing(cat[1],pT);
1104         else 
1105           smear = GetMomentumSmearing(cat[0],pT);
1106
1107         fp1PtResolution->Fill(vp->Pt(),smear);
1108
1109         Double_t sigma = vp->Pt()*smear;
1110         Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
1111         
1112         Double_t phi   = vp->Phi();
1113         Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
1114         Double_t pX    = pTrec * TMath::Cos(phi);
1115         Double_t pY    = pTrec * TMath::Sin(phi);
1116         Double_t pZ    = pTrec/TMath::Tan(theta);
1117         Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
1118
1119         fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
1120
1121         fastjet::PseudoJet jInp(pX,pY,pZ,p);
1122         jInp.set_user_index(i);
1123         inputParticlesRec.push_back(jInp);
1124
1125       }
1126       else {
1127         fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
1128         jInp.set_user_index(i);
1129         inputParticlesRec.push_back(jInp);
1130
1131       }
1132
1133     }
1134
1135     // the randomized input changes eta and phi, but keeps the p_T
1136     if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1137       Double_t pT = vp->Pt();
1138       Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
1139       Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
1140       
1141       Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));  
1142       Double_t pZ = pT/TMath::Tan(theta);
1143
1144       Double_t pX = pT * TMath::Cos(phi);
1145       Double_t pY = pT * TMath::Sin(phi);
1146       Double_t p  = TMath::Sqrt(pT*pT+pZ*pZ); 
1147       fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
1148
1149       jInpRan.set_user_index(i);
1150       inputParticlesRecRan.push_back(jInpRan);
1151       vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
1152     }
1153
1154     // fill the tref array, only needed when we write out jets
1155     if(fTCAJetsOut){
1156       if(i == 0){
1157         fRef->Delete(); // make sure to delete before placement new...
1158         if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
1159           new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
1160         } 
1161       }
1162       if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp);  //TRefArray does not work with toy model ...
1163     }
1164   }// recparticles
1165
1166   if(inputParticlesRec.size()==0){
1167     if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
1168     PostData(1, fHistList);
1169     return;
1170   }
1171   
1172   // run fast jet
1173   // employ setters for these...
1174
1175  
1176   // now create the object that holds info about ghosts                        
1177   /*
1178   if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
1179     // reduce CPU time...
1180     fGhostArea = 0.5; 
1181     fActiveAreaRepeats = 0; 
1182   }
1183   */
1184
1185   fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
1186   fastjet::AreaType areaType =   fastjet::active_area;
1187   fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
1188   fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
1189   fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
1190   
1191   //range where to compute background
1192   Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0;
1193   phiMin = 0;
1194   phiMax = 2*TMath::Pi();
1195   rapMax = fGhostEtamax - fRparam;
1196   rapMin = - fGhostEtamax + fRparam;
1197   fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
1198  
1199
1200   const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
1201   const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
1202
1203  
1204   fh1NJetsRec->Fill(sortedJets.size());
1205
1206  // loop over all jets an fill information, first on is the leading jet
1207
1208   Int_t nRecOver = inclusiveJets.size();
1209   Int_t nRec     = inclusiveJets.size();
1210   if(inclusiveJets.size()>0){
1211     AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
1212     Double_t area = clustSeq.area(sortedJets[0]);
1213     leadingJet.SetEffArea(area,0);
1214     Float_t pt = leadingJet.Pt();
1215     Int_t nAodOutJets = 0;
1216     Int_t nAodOutTracks = 0;
1217     AliAODJet *aodOutJet = 0;
1218
1219     Int_t iCount = 0;
1220     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1221       Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1222       while(pt<ptCut&&iCount<nRec){
1223         nRecOver--;
1224         iCount++;
1225         if(iCount<nRec){
1226           pt = sortedJets[iCount].perp();
1227         }
1228       }
1229       if(nRecOver<=0)break;
1230       fh2NRecJetsPt->Fill(ptCut,nRecOver);
1231     }
1232     Float_t phi = leadingJet.Phi();
1233     if(phi<0)phi+=TMath::Pi()*2.;    
1234     Float_t eta = leadingJet.Eta();
1235     Float_t pTback = 0;
1236     if(externalBackground){
1237       // carefull has to be filled in a task before
1238       // todo, ReArrange to the botom 
1239      pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
1240     }
1241     pt = leadingJet.Pt() - pTback;
1242     // correlation of leading jet with tracks
1243     TIterator *recIter = recParticles.MakeIterator();
1244     recIter->Reset();
1245     AliVParticle *tmpRecTrack = 0;
1246     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
1247       Float_t tmpPt = tmpRecTrack->Pt();
1248       // correlation
1249       Float_t tmpPhi =  tmpRecTrack->Phi();     
1250       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1251       Float_t dPhi = phi - tmpPhi;
1252       if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1253       if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1254       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
1255       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
1256       if(cenClass>=0){
1257         fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
1258         fh2TracksLeadingJetPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1259       }
1260
1261     }  
1262     
1263    
1264     TLorentzVector vecareab;
1265     for(int j = 0; j < nRec;j++){
1266       AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
1267       aodOutJet = 0;
1268       nAodOutTracks = 0;
1269       Float_t tmpPt = tmpRec.Pt();  
1270       
1271       if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
1272         aodOutJet =  new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
1273         aodOutJet->GetRefTracks()->Clear();
1274         Double_t area1 = clustSeq.area(sortedJets[j]);
1275         aodOutJet->SetEffArea(area1,0);
1276         fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);  
1277         vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());     
1278         aodOutJet->SetVectorAreaCharged(&vecareab);
1279       }
1280
1281
1282       Float_t tmpPtBack = 0;
1283       if(externalBackground){
1284         // carefull has to be filled in a task before
1285        // todo, ReArrange to the botom
1286         tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
1287       }
1288       tmpPt = tmpPt - tmpPtBack;
1289       if(tmpPt<0)tmpPt = 0; // avoid negative weights...
1290       
1291       fh1PtJetsRecIn->Fill(tmpPt);
1292       // Fill Spectra with constituentsemacs
1293       const vector<fastjet::PseudoJet> &constituents = clustSeq.constituents(sortedJets[j]);
1294
1295       fh1NConstRec->Fill(constituents.size());
1296       fh2PtNch->Fill(nCh,tmpPt);
1297       fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
1298       fh2NConstPt->Fill(tmpPt,constituents.size());
1299       // loop over constiutents and fill spectrum
1300    
1301       AliVParticle *partLead = 0;
1302       Float_t ptLead = -1;
1303
1304       for(unsigned int ic = 0; ic < constituents.size();ic++){
1305         AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
1306         if(!part) continue;
1307         fh1PtJetConstRec->Fill(part->Pt());
1308         if(aodOutJet){
1309           if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
1310           if(part->Pt()>fMaxTrackPtInJet){
1311             aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1312           }
1313         }
1314         if(part->Pt()>ptLead){
1315           ptLead = part->Pt();
1316           partLead = part;
1317         }
1318         if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
1319       }
1320
1321       AliAODTrack *aodT = 0;
1322       if(partLead){
1323         if(aodOutJet){
1324           //set pT of leading constituent of jet
1325           aodOutJet->SetPtLeading(partLead->Pt());
1326           aodT = dynamic_cast<AliAODTrack*>(partLead);
1327           if(aodT){
1328             if(aodT->TestFilterBit(fFilterMaskBestPt)){
1329               aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
1330             }
1331           }
1332         }
1333       }
1334     
1335      // correlation
1336      Float_t tmpPhi =  tmpRec.Phi();
1337      Float_t tmpEta =  tmpRec.Eta();
1338      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;        
1339      if(j==0){
1340        fh1PtJetsLeadingRecIn->Fill(tmpPt);
1341        fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
1342        fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
1343        fh1NConstLeadingRec->Fill(constituents.size());
1344        fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
1345        continue;
1346      }
1347      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
1348      fh2JetEtaPt->Fill(tmpEta,tmpPt);
1349      Float_t dPhi = phi - tmpPhi;
1350      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1351      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1352      Float_t dEta = eta - tmpRec.Eta();
1353      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1354      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
1355      if(cenClass>=0){
1356        fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
1357        fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
1358      }
1359      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1360     }// loop over reconstructed jets
1361    delete recIter;
1362
1363
1364
1365    // Add the random cones...
1366    if(fNRandomCones>0&&fTCARandomConesOut){       
1367      // create a random jet within the acceptance
1368      Double_t etaMax = fTrackEtaWindow - fRparam;
1369      Int_t nCone = 0;
1370      Int_t nConeRan = 0;
1371      Double_t pTC = 1; // small number
1372      for(int ir = 0;ir < fNRandomCones;ir++){
1373        Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
1374        Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
1375        // massless jet
1376        Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1377        Double_t pZC = pTC/TMath::Tan(thetaC);
1378        Double_t pXC = pTC * TMath::Cos(phiC);
1379        Double_t pYC = pTC * TMath::Sin(phiC);
1380        Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1381        AliAODJet tmpRecC (pXC,pYC,pZC, pC); 
1382        bool skip = false;
1383        for(int jj = 0; jj < TMath::Min(nRec,fNSkipLeadingCone);jj++){// test for overlap with leading jets
1384          AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
1385          if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
1386            skip = true;
1387            break;
1388          }
1389        }
1390        // test for overlap with previous cones to avoid double counting
1391        for(int iic = 0;iic<ir;iic++){
1392          AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
1393          if(iicone){
1394            if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
1395              skip = true;
1396              break;
1397            }
1398          }
1399        }
1400        if(skip)continue;
1401        tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
1402        tmpRecC.SetPtLeading(-1.);
1403        if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
1404        if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
1405      }// loop over random cones creation
1406
1407   
1408      // loop over the reconstructed particles and add up the pT in the random cones
1409      // maybe better to loop over randomized particles not in the real jets...
1410      // but this by definition brings dow average energy in the whole  event
1411      AliAODJet vTmpRanR(1,0,0,1);
1412      for(int i = 0; i < recParticles.GetEntries(); i++){
1413        AliVParticle *vp = (AliVParticle*)recParticles.At(i);
1414        if(fTCARandomConesOut){
1415          for(int ir = 0;ir < fNRandomCones;ir++){
1416            AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);  
1417            if(jC&&jC->DeltaR(vp)<fRparam){
1418              if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1419              jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
1420              if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
1421            }
1422          }  
1423        }// add up energy in cone
1424
1425        // the randomized input changes eta and phi, but keeps the p_T
1426        if(i>=fNSkipLeadingRan){// eventually skip the leading particles
1427          Double_t pTR = vp->Pt();
1428          Double_t etaR = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
1429          Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
1430          
1431          Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));  
1432          Double_t pZR = pTR/TMath::Tan(thetaR);
1433          
1434          Double_t pXR = pTR * TMath::Cos(phiR);
1435          Double_t pYR = pTR * TMath::Sin(phiR);
1436          Double_t pR  = TMath::Sqrt(pTR*pTR+pZR*pZR); 
1437          vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
1438          if(fTCARandomConesOutRan){
1439            for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1440              AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);  
1441              if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
1442                if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
1443                jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
1444                if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
1445              }
1446            }  
1447          }
1448        }
1449      }// loop over recparticles
1450     
1451      Float_t jetArea = fRparam*fRparam*TMath::Pi();
1452      if(fTCARandomConesOut){
1453        for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
1454          // rescale the momentum vectors for the random cones
1455          
1456          AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
1457          if(rC){
1458            Double_t etaC = rC->Eta();
1459            Double_t phiC = rC->Phi();
1460            // massless jet, unit vector
1461            pTC = rC->ChargedBgEnergy();
1462            if(pTC<=0)pTC = 0.001; // for almost empty events
1463            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1464            Double_t pZC = pTC/TMath::Tan(thetaC);
1465            Double_t pXC = pTC * TMath::Cos(phiC);
1466            Double_t pYC = pTC * TMath::Sin(phiC);
1467            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1468            rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
1469            rC->SetBgEnergy(0,0);
1470            rC->SetEffArea(jetArea,0);
1471          }
1472        }
1473      }
1474      if(fTCARandomConesOutRan){
1475        for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
1476          AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
1477          // same wit random
1478          if(rC){
1479            Double_t etaC = rC->Eta();
1480            Double_t phiC = rC->Phi();
1481            // massless jet, unit vector
1482            pTC = rC->ChargedBgEnergy();
1483            if(pTC<=0)pTC = 0.001;// for almost empty events
1484            Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));  
1485            Double_t pZC = pTC/TMath::Tan(thetaC);
1486            Double_t pXC = pTC * TMath::Cos(phiC);
1487            Double_t pYC = pTC * TMath::Sin(phiC);
1488            Double_t pC  = TMath::Sqrt(pTC*pTC+pZC*pZC); 
1489            rC->SetPxPyPzE(pXC,pYC,pZC, pC); 
1490            rC->SetBgEnergy(0,0);
1491            rC->SetEffArea(jetArea,0);
1492          }
1493        }
1494      }
1495    }// if(fNRandomCones
1496   
1497    //background estimates:all bckg jets(0) & wo the 2 hardest(1)
1498
1499    if(fAODJetBackgroundOut){
1500      vector<fastjet::PseudoJet> jets2=sortedJets;
1501      if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); 
1502
1503      Double_t bkg1=0;
1504      Double_t sigma1=0.;
1505      Double_t meanarea1=0.;
1506      Double_t bkg2=0;
1507      Double_t sigma2=0.;
1508      Double_t meanarea2=0.;
1509
1510      Double_t bkg4=0;
1511      Double_t sigma4=0.;
1512      Double_t meanarea4=0.;
1513
1514      clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1515      fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1516
1517      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));    
1518      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));      
1519        clustSeq.get_median_rho_and_sigma(sortedJets, range, true, bkg4, sigma4, meanarea4, true);
1520      fAODJetBackgroundOut->SetBackground(3,bkg4,sigma4,meanarea4);
1521
1522      //////////////////////
1523      
1524      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1525      fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1526      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));    
1527      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));    
1528
1529    }
1530   }
1531
1532   if(fStoreRhoLeadingTrackCorr) {
1533     vector<fastjet::PseudoJet> jets3=sortedJets;
1534     if(jets3.size()>2) jets3.erase(jets3.begin(),jets3.begin()+2); 
1535
1536     Double_t bkg2=0;
1537     Double_t sigma2=0.;
1538     Double_t meanarea2=0.;
1539
1540     clustSeq.get_median_rho_and_sigma(jets3, range, false, bkg2, sigma2, meanarea2, true);
1541     fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1542
1543     //Get leading track in event
1544     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1545
1546     fh3CentvsRhoLeadingTrackPt->Fill(cent,bkg2,leading->Pt());
1547     fh3CentvsSigmaLeadingTrackPt->Fill(cent,sigma2,leading->Pt());
1548     fh3MultvsRhoLeadingTrackPt->Fill(nCh,bkg2,leading->Pt());
1549     fh3MultvsSigmaLeadingTrackPt->Fill(nCh,sigma2,leading->Pt());
1550      
1551
1552     //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1553     //exclude 2 leading jets from event
1554     //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R  (Near side to leading track)
1555     //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1556     //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R)  (Away side to leading track)
1557     //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R)   (Orthogonal to leading track)
1558
1559     //Assuming R=0.2 for background clusters
1560
1561     Double_t bkg2Q[4]      = {0.};
1562     Double_t sigma2Q[4]    = {0.};
1563     Double_t meanarea2Q[4] = {0.};
1564
1565     //Get phi of leading track in event
1566     Float_t phiLeadingTrack = leading->Phi();
1567     Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
1568
1569     //Quadrant1 - near side
1570     phiMin = phiLeadingTrack - phiStep;
1571     phiMax = phiLeadingTrack + phiStep;
1572     fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1573     clustSeq.get_median_rho_and_sigma(jets3, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1574
1575     //Quadrant2 - orthogonal
1576     Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
1577     if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
1578     phiMin = phiQ2 - phiStep;
1579     phiMax = phiQ2 + phiStep;
1580     fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1581     clustSeq.get_median_rho_and_sigma(jets3, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1582
1583     //Quadrant3 - away side
1584     Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
1585     if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
1586     phiMin = phiQ3 - phiStep;
1587     phiMax = phiQ3 + phiStep;
1588     fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1589     clustSeq.get_median_rho_and_sigma(jets3, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1590
1591     //Quadrant4 - othogonal
1592     Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
1593     if(phiQ4>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
1594     phiMin = phiQ4 - phiStep;
1595     phiMax = phiQ4 + phiStep;
1596     fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1597     clustSeq.get_median_rho_and_sigma(jets3, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1598
1599     fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
1600     fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
1601     fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
1602     fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
1603
1604     fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
1605     fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
1606     fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
1607     fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
1608
1609     fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
1610     fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
1611     fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
1612     fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
1613
1614     fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
1615     fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
1616     fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
1617     fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
1618
1619     fh3CentvsDeltaRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0]-bkg2,leading->Pt());
1620     fh3CentvsDeltaRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1]-bkg2,leading->Pt());
1621     fh3CentvsDeltaRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2]-bkg2,leading->Pt());
1622     fh3CentvsDeltaRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3]-bkg2,leading->Pt());
1623
1624   }
1625    
1626
1627   
1628   
1629  
1630   // fill track information
1631   Int_t nTrackOver = recParticles.GetSize();
1632   // do the same for tracks and jets
1633
1634   if(nTrackOver>0){
1635    TIterator *recIter = recParticles.MakeIterator();
1636    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());  
1637    Float_t pt = tmpRec->Pt();
1638
1639    //    Printf("Leading track p_t %3.3E",pt);
1640    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1641      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1642      while(pt<ptCut&&tmpRec){
1643        nTrackOver--;
1644        tmpRec = (AliVParticle*)(recIter->Next()); 
1645        if(tmpRec){
1646          pt = tmpRec->Pt();
1647        }
1648      }
1649      if(nTrackOver<=0)break;
1650      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1651    }
1652    
1653    recIter->Reset();
1654    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1655    Float_t phi = leading->Phi();
1656    if(phi<0)phi+=TMath::Pi()*2.;    
1657    Float_t eta = leading->Eta();
1658    pt  = leading->Pt();
1659    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1660      Float_t tmpPt = tmpRec->Pt();
1661      Float_t tmpEta = tmpRec->Eta();
1662      fh1PtTracksRecIn->Fill(tmpPt);
1663      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1664      if(tmpRec==leading){
1665        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1666        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1667        continue;
1668      }
1669       // correlation
1670      Float_t tmpPhi =  tmpRec->Phi();
1671      
1672      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1673      Float_t dPhi = phi - tmpPhi;
1674      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1675      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1676      Float_t dEta = eta - tmpRec->Eta();
1677      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1678      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1679      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1680    }  
1681    delete recIter;
1682  }
1683
1684  // find the random jets
1685
1686  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1687   
1688  // fill the jet information from random track
1689  const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1690  const vector <fastjet::PseudoJet> &sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1691
1692   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1693
1694  // loop over all jets an fill information, first on is the leading jet
1695
1696  Int_t nRecOverRan = inclusiveJetsRan.size();
1697  Int_t nRecRan     = inclusiveJetsRan.size();
1698
1699  if(inclusiveJetsRan.size()>0){
1700    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1701    Float_t pt = leadingJet.Pt();
1702    
1703    Int_t iCount = 0;
1704    TLorentzVector vecarearanb;
1705
1706    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1707      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1708       while(pt<ptCut&&iCount<nRecRan){
1709         nRecOverRan--;
1710         iCount++;
1711         if(iCount<nRecRan){
1712           pt = sortedJetsRan[iCount].perp();
1713         }
1714       }
1715       if(nRecOverRan<=0)break;
1716       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1717     }
1718     Float_t phi = leadingJet.Phi();
1719     if(phi<0)phi+=TMath::Pi()*2.;    
1720     pt = leadingJet.Pt();
1721
1722     // correlation of leading jet with random tracks
1723
1724     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1725       { 
1726         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1727         // correlation
1728         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1729         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1730         Float_t dPhi = phi - tmpPhi;
1731         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1732         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1733         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1734         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1735       }  
1736
1737     Int_t nAodOutJetsRan = 0;
1738      AliAODJet *aodOutJetRan = 0;
1739     for(int j = 0; j < nRecRan;j++){
1740       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1741       Float_t tmpPt = tmpRec.Pt();
1742       fh1PtJetsRecInRan->Fill(tmpPt);
1743       // Fill Spectra with constituents
1744       const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1745       fh1NConstRecRan->Fill(constituents.size());
1746       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1747       fh2PtNchRan->Fill(nCh,tmpPt);
1748       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1749
1750
1751      if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1752        aodOutJetRan =  new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1753        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1754        aodOutJetRan->GetRefTracks()->Clear();
1755        aodOutJetRan->SetEffArea(arearan,0);
1756        fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);  
1757        vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1758        aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1759
1760      }
1761
1762       // correlation
1763       Float_t tmpPhi =  tmpRec.Phi();
1764       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1765
1766       if(j==0){
1767         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1768         fh1NConstLeadingRecRan->Fill(constituents.size());
1769         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1770         continue;
1771       }
1772     }  
1773
1774      
1775     if(fAODJetBackgroundOut){
1776      Double_t bkg3=0.;
1777      Double_t sigma3=0.;
1778      Double_t meanarea3=0.;
1779      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1780      fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1781      //     float areaRandomCone = rRandomCone2 *TMath::Pi();         
1782      /*
1783      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));    
1784      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));    
1785      */
1786     }
1787
1788
1789
1790  }
1791
1792
1793  // do the event selection if activated
1794  if(fJetTriggerPtCut>0){
1795    bool select = false;
1796    Float_t minPt = fJetTriggerPtCut;
1797    /*
1798    // hard coded for now ...
1799    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1800    if(cent<10)minPt = 50;
1801    else if(cent<30)minPt = 42;
1802    else if(cent<50)minPt = 28;
1803    else if(cent<80)minPt = 18;
1804    */
1805    float rho = 0;
1806    if(externalBackground)rho = externalBackground->GetBackground(2);
1807    if(fTCAJetsOut){
1808      for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1809        AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1810        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1811        if(ptSub>=minPt){
1812          select = true;
1813          break;
1814        }
1815      }
1816    }   
1817  
1818    if(select){
1819      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1820      fh1CentralitySelect->Fill(cent);
1821      fh1ZSelect->Fill(zVtx);
1822      aodH->SetFillAOD(kTRUE);
1823    }
1824  }
1825  if (fDebug > 2){
1826    if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1827    if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1828    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1829    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1830  }
1831  PostData(1, fHistList);
1832 }
1833
1834 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1835 {
1836   //
1837   // Terminate analysis
1838   //
1839     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1840
1841     if(fMomResH1Fit) delete fMomResH1Fit;
1842     if(fMomResH2Fit) delete fMomResH2Fit;
1843     if(fMomResH3Fit) delete fMomResH3Fit;
1844     
1845 }
1846
1847
1848 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1849
1850   //
1851   // get list of tracks/particles for different types
1852   //
1853
1854   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1855
1856   Int_t iCount = 0;
1857   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1858     if(type!=kTrackAODextraonly) {
1859       AliAODEvent *aod = 0;
1860       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1861       else aod = AODEvent();
1862       if(!aod){
1863         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1864         return iCount;
1865       }
1866
1867       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1868         AliAODTrack *tr = aod->GetTrack(it);
1869         Bool_t bGood = false;
1870         if(fFilterType == 0)bGood = true;
1871         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1872         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1873         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1874           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());     
1875           continue;
1876         }
1877         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1878           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());     
1879           continue;
1880         }
1881         if(tr->Pt()<fTrackPtCut){
1882           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());      
1883           continue;
1884         }
1885         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());        
1886         list->Add(tr);
1887         iCount++;
1888       }
1889     }
1890     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1891       AliAODEvent *aod = 0;
1892       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1893       else aod = AODEvent();
1894       
1895       if(!aod){
1896         return iCount;
1897       }
1898       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1899       if(!aodExtraTracks)return iCount;
1900       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1901         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1902         if (!track) continue;
1903
1904         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1905         if(!trackAOD)continue;
1906         Bool_t bGood = false;
1907         if(fFilterType == 0)bGood = true;
1908         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1909         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1910         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1911         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1912         if(trackAOD->Pt()<fTrackPtCut) continue;
1913         list->Add(trackAOD);
1914         iCount++;
1915       }
1916     }
1917   }
1918   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1919     AliMCEvent* mcEvent = MCEvent();
1920     if(!mcEvent)return iCount;
1921     // we want to have alivpartilces so use get track
1922     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1923       if(!mcEvent->IsPhysicalPrimary(it))continue;
1924       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1925       if(type == kTrackKineAll){
1926         if(part->Pt()<fTrackPtCut)continue;
1927         list->Add(part);
1928         iCount++;
1929       }
1930       else if(type == kTrackKineCharged){
1931         if(part->Particle()->GetPDG()->Charge()==0)continue;
1932         if(part->Pt()<fTrackPtCut)continue;
1933         list->Add(part);
1934         iCount++;
1935       }
1936     }
1937   }
1938   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1939     AliAODEvent *aod = 0;
1940     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1941     else aod = AODEvent();
1942     if(!aod)return iCount;
1943     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1944     if(!tca)return iCount;
1945     for(int it = 0;it < tca->GetEntriesFast();++it){
1946       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1947       if(!part->IsPhysicalPrimary())continue;
1948       if(type == kTrackAODMCAll){
1949         if(part->Pt()<fTrackPtCut)continue;
1950         list->Add(part);
1951         iCount++;
1952       }
1953       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1954         if(part->Charge()==0)continue;
1955         if(part->Pt()<fTrackPtCut)continue;
1956         if(kTrackAODMCCharged){
1957           list->Add(part);
1958         }
1959         else {
1960           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1961           list->Add(part);
1962         }
1963         iCount++;
1964       }
1965     }
1966   }// AODMCparticle
1967   list->Sort();
1968   return iCount;
1969 }
1970
1971 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
1972
1973   TFile *f = TFile::Open(fPathTrPtResolution.Data());
1974
1975   if(!f)return;
1976
1977   TProfile *fProfPtPtSigma1PtGlobSt     = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
1978   TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
1979   TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
1980
1981   SetSmearResolution(kTRUE);
1982   SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
1983   
1984
1985 }
1986
1987 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
1988
1989   TFile *f = TFile::Open(fPathTrEfficiency.Data());
1990   if(!f)return;
1991
1992   TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
1993   TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
1994   TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
1995
1996   SetDiceEfficiency(kTRUE);
1997
1998   if(fChangeEfficiencyFraction>0.) {
1999
2000     TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
2001     
2002     for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
2003       Double_t content = hEffPosGlobSt->GetBinContent(bin);
2004       hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
2005     }
2006
2007     SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2008
2009   }
2010   else
2011     SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2012
2013 }
2014
2015 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
2016
2017   //
2018   // set mom res profiles
2019   //
2020
2021   if(fMomResH1) delete fMomResH1;
2022   if(fMomResH2) delete fMomResH2;
2023   if(fMomResH3) delete fMomResH3;
2024
2025   fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
2026   fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
2027   fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
2028
2029 }
2030
2031 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
2032   //
2033   // set tracking efficiency histos
2034   //
2035
2036   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
2037   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
2038   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
2039 }
2040
2041 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
2042
2043   //
2044   // Get smearing on generated momentum
2045   //
2046
2047   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
2048
2049   TProfile *fMomRes = 0x0;
2050   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
2051   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
2052   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
2053
2054   if(!fMomRes) {
2055     return 0.;
2056   }
2057
2058
2059   Double_t smear = 0.;
2060
2061   if(pt>20.) {
2062     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
2063     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
2064     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
2065   }
2066   else {
2067
2068     Int_t bin = fMomRes->FindBin(pt);
2069
2070     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
2071
2072   }
2073
2074   if(fMomRes) delete fMomRes;
2075   
2076   return smear;
2077 }
2078
2079 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
2080   //
2081   // Fit linear function on momentum resolution at high pT
2082   //
2083
2084   if(!fMomResH1Fit && fMomResH1) {
2085     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
2086     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
2087     fMomResH1Fit ->SetRange(5.,100.);
2088   }
2089
2090   if(!fMomResH2Fit && fMomResH2) {
2091     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
2092     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
2093     fMomResH2Fit ->SetRange(5.,100.);
2094   }
2095
2096   if(!fMomResH3Fit && fMomResH3) {
2097     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
2098     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
2099     fMomResH3Fit ->SetRange(5.,100.);
2100   }
2101
2102 }
2103
2104 /*
2105 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
2106   for(int i = 0; i < particles.GetEntries(); i++){
2107     AliVParticle *vp = (AliVParticle*)particles.At(i);
2108     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
2109     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
2110     jInp.set_user_index(i);
2111     inputParticles.push_back(jInp);
2112   }
2113
2114   return 0;
2115
2116 }
2117 */