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