adding requirement on tzero or vzero if needed
[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      clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
1511      fAODJetBackgroundOut->SetBackground(0,bkg1,sigma1,meanarea1);
1512
1513      //     fh1BiARandomCones[0]->Fill(omCone-(bkg1*areaRandomCone));    
1514      //  fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));    
1515      
1516      clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
1517      fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1518      //  fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));    
1519      //   fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));    
1520
1521    }
1522   }
1523
1524   if(fStoreRhoLeadingTrackCorr) {
1525     vector<fastjet::PseudoJet> jets3=sortedJets;
1526     if(jets3.size()>2) jets3.erase(jets3.begin(),jets3.begin()+2); 
1527
1528     Double_t bkg2=0;
1529     Double_t sigma2=0.;
1530     Double_t meanarea2=0.;
1531
1532     clustSeq.get_median_rho_and_sigma(jets3, range, false, bkg2, sigma2, meanarea2, true);
1533     fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
1534
1535     //Get leading track in event
1536     AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1537
1538     fh3CentvsRhoLeadingTrackPt->Fill(cent,bkg2,leading->Pt());
1539     fh3CentvsSigmaLeadingTrackPt->Fill(cent,sigma2,leading->Pt());
1540     fh3MultvsRhoLeadingTrackPt->Fill(nCh,bkg2,leading->Pt());
1541     fh3MultvsSigmaLeadingTrackPt->Fill(nCh,sigma2,leading->Pt());
1542      
1543
1544     //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
1545     //exclude 2 leading jets from event
1546     //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R  (Near side to leading track)
1547     //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
1548     //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R)  (Away side to leading track)
1549     //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R)   (Orthogonal to leading track)
1550
1551     //Assuming R=0.2 for background clusters
1552
1553     Double_t bkg2Q[4]      = {0.};
1554     Double_t sigma2Q[4]    = {0.};
1555     Double_t meanarea2Q[4] = {0.};
1556
1557     //Get phi of leading track in event
1558     Float_t phiLeadingTrack = leading->Phi();
1559     Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
1560
1561     //Quadrant1 - near side
1562     phiMin = phiLeadingTrack - phiStep;
1563     phiMax = phiLeadingTrack + phiStep;
1564     fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
1565     clustSeq.get_median_rho_and_sigma(jets3, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
1566
1567     //Quadrant2 - orthogonal
1568     Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
1569     if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
1570     phiMin = phiQ2 - phiStep;
1571     phiMax = phiQ2 + phiStep;
1572     fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
1573     clustSeq.get_median_rho_and_sigma(jets3, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
1574
1575     //Quadrant3 - away side
1576     Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
1577     if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
1578     phiMin = phiQ3 - phiStep;
1579     phiMax = phiQ3 + phiStep;
1580     fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
1581     clustSeq.get_median_rho_and_sigma(jets3, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
1582
1583     //Quadrant4 - othogonal
1584     Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
1585     if(phiQ4>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
1586     phiMin = phiQ4 - phiStep;
1587     phiMax = phiQ4 + phiStep;
1588     fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
1589     clustSeq.get_median_rho_and_sigma(jets3, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
1590
1591     fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
1592     fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
1593     fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
1594     fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
1595
1596     fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
1597     fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
1598     fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
1599     fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
1600
1601     fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
1602     fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
1603     fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
1604     fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
1605
1606     fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
1607     fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
1608     fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
1609     fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
1610
1611     fh3CentvsDeltaRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0]-bkg2,leading->Pt());
1612     fh3CentvsDeltaRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1]-bkg2,leading->Pt());
1613     fh3CentvsDeltaRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2]-bkg2,leading->Pt());
1614     fh3CentvsDeltaRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3]-bkg2,leading->Pt());
1615
1616   }
1617    
1618
1619   
1620   
1621  
1622   // fill track information
1623   Int_t nTrackOver = recParticles.GetSize();
1624   // do the same for tracks and jets
1625
1626   if(nTrackOver>0){
1627    TIterator *recIter = recParticles.MakeIterator();
1628    AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());  
1629    Float_t pt = tmpRec->Pt();
1630
1631    //    Printf("Leading track p_t %3.3E",pt);
1632    for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1633      Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1634      while(pt<ptCut&&tmpRec){
1635        nTrackOver--;
1636        tmpRec = (AliVParticle*)(recIter->Next()); 
1637        if(tmpRec){
1638          pt = tmpRec->Pt();
1639        }
1640      }
1641      if(nTrackOver<=0)break;
1642      fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1643    }
1644    
1645    recIter->Reset();
1646    AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1647    Float_t phi = leading->Phi();
1648    if(phi<0)phi+=TMath::Pi()*2.;    
1649    Float_t eta = leading->Eta();
1650    pt  = leading->Pt();
1651    while((tmpRec = (AliVParticle*)(recIter->Next()))){
1652      Float_t tmpPt = tmpRec->Pt();
1653      Float_t tmpEta = tmpRec->Eta();
1654      fh1PtTracksRecIn->Fill(tmpPt);
1655      fh2TrackEtaPt->Fill(tmpEta,tmpPt);
1656      if(tmpRec==leading){
1657        fh1PtTracksLeadingRecIn->Fill(tmpPt);
1658        fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt);
1659        continue;
1660      }
1661       // correlation
1662      Float_t tmpPhi =  tmpRec->Phi();
1663      
1664      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1665      Float_t dPhi = phi - tmpPhi;
1666      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1667      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1668      Float_t dEta = eta - tmpRec->Eta();
1669      fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1670      fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1671      fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
1672    }  
1673    delete recIter;
1674  }
1675
1676  // find the random jets
1677
1678  fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
1679   
1680  // fill the jet information from random track
1681  const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
1682  const vector <fastjet::PseudoJet> &sortedJetsRan    = sorted_by_pt(inclusiveJetsRan);
1683
1684   fh1NJetsRecRan->Fill(sortedJetsRan.size());
1685
1686  // loop over all jets an fill information, first on is the leading jet
1687
1688  Int_t nRecOverRan = inclusiveJetsRan.size();
1689  Int_t nRecRan     = inclusiveJetsRan.size();
1690
1691  if(inclusiveJetsRan.size()>0){
1692    AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
1693    Float_t pt = leadingJet.Pt();
1694    
1695    Int_t iCount = 0;
1696    TLorentzVector vecarearanb;
1697
1698    for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
1699      Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
1700       while(pt<ptCut&&iCount<nRecRan){
1701         nRecOverRan--;
1702         iCount++;
1703         if(iCount<nRecRan){
1704           pt = sortedJetsRan[iCount].perp();
1705         }
1706       }
1707       if(nRecOverRan<=0)break;
1708       fh2NRecJetsPtRan->Fill(ptCut,nRecOverRan);
1709     }
1710     Float_t phi = leadingJet.Phi();
1711     if(phi<0)phi+=TMath::Pi()*2.;    
1712     pt = leadingJet.Pt();
1713
1714     // correlation of leading jet with random tracks
1715
1716     for(unsigned int ip = 0; ip < inputParticlesRecRan.size();ip++)
1717       { 
1718         Float_t tmpPt = inputParticlesRecRan[ip].perp();
1719         // correlation
1720         Float_t tmpPhi =  inputParticlesRecRan[ip].phi();
1721         if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1722         Float_t dPhi = phi - tmpPhi;
1723         if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1724         if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
1725         fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
1726         fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt);
1727       }  
1728
1729     Int_t nAodOutJetsRan = 0;
1730      AliAODJet *aodOutJetRan = 0;
1731     for(int j = 0; j < nRecRan;j++){
1732       AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E());
1733       Float_t tmpPt = tmpRec.Pt();
1734       fh1PtJetsRecInRan->Fill(tmpPt);
1735       // Fill Spectra with constituents
1736       const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
1737       fh1NConstRecRan->Fill(constituents.size());
1738       fh2NConstPtRan->Fill(tmpPt,constituents.size());
1739       fh2PtNchRan->Fill(nCh,tmpPt);
1740       fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size());
1741
1742
1743      if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
1744        aodOutJetRan =  new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
1745        Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
1746        aodOutJetRan->GetRefTracks()->Clear();
1747        aodOutJetRan->SetEffArea(arearan,0);
1748        fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);  
1749        vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
1750        aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
1751
1752      }
1753
1754       // correlation
1755       Float_t tmpPhi =  tmpRec.Phi();
1756       if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
1757
1758       if(j==0){
1759         fh1PtJetsLeadingRecInRan->Fill(tmpPt);
1760         fh1NConstLeadingRecRan->Fill(constituents.size());
1761         fh2NConstLeadingPtRan->Fill(tmpPt,constituents.size());
1762         continue;
1763       }
1764     }  
1765
1766      
1767     if(fAODJetBackgroundOut){
1768      Double_t bkg3=0.;
1769      Double_t sigma3=0.;
1770      Double_t meanarea3=0.;
1771      clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true);
1772      fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
1773      //     float areaRandomCone = rRandomCone2 *TMath::Pi();         
1774      /*
1775      fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));    
1776      fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));    
1777      */
1778     }
1779
1780
1781
1782  }
1783
1784
1785  // do the event selection if activated
1786  if(fJetTriggerPtCut>0){
1787    bool select = false;
1788    Float_t minPt = fJetTriggerPtCut;
1789    /*
1790    // hard coded for now ...
1791    // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
1792    if(cent<10)minPt = 50;
1793    else if(cent<30)minPt = 42;
1794    else if(cent<50)minPt = 28;
1795    else if(cent<80)minPt = 18;
1796    */
1797    float rho = 0;
1798    if(externalBackground)rho = externalBackground->GetBackground(2);
1799    if(fTCAJetsOut){
1800      for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
1801        AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
1802        Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
1803        if(ptSub>=minPt){
1804          select = true;
1805          break;
1806        }
1807      }
1808    }   
1809  
1810    if(select){
1811      static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1812      fh1CentralitySelect->Fill(cent);
1813      fh1ZSelect->Fill(zVtx);
1814      aodH->SetFillAOD(kTRUE);
1815    }
1816  }
1817  if (fDebug > 2){
1818    if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
1819    if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
1820    if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
1821    if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
1822  }
1823  PostData(1, fHistList);
1824 }
1825
1826 void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
1827 {
1828   //
1829   // Terminate analysis
1830   //
1831     if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
1832
1833     if(fMomResH1Fit) delete fMomResH1Fit;
1834     if(fMomResH2Fit) delete fMomResH2Fit;
1835     if(fMomResH3Fit) delete fMomResH3Fit;
1836     
1837 }
1838
1839
1840 Int_t  AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
1841
1842   //
1843   // get list of tracks/particles for different types
1844   //
1845
1846   if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1847
1848   Int_t iCount = 0;
1849   if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
1850     if(type!=kTrackAODextraonly) {
1851       AliAODEvent *aod = 0;
1852       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1853       else aod = AODEvent();
1854       if(!aod){
1855         if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1856         return iCount;
1857       }
1858
1859       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1860         AliAODTrack *tr = aod->GetTrack(it);
1861         Bool_t bGood = false;
1862         if(fFilterType == 0)bGood = true;
1863         else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1864         else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1865         if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
1866           if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());     
1867           continue;
1868         }
1869         if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
1870           if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());     
1871           continue;
1872         }
1873         if(tr->Pt()<fTrackPtCut){
1874           if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());      
1875           continue;
1876         }
1877         if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());        
1878         list->Add(tr);
1879         iCount++;
1880       }
1881     }
1882     if(type==kTrackAODextra || type==kTrackAODextraonly) {
1883       AliAODEvent *aod = 0;
1884       if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1885       else aod = AODEvent();
1886       
1887       if(!aod){
1888         return iCount;
1889       }
1890       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
1891       if(!aodExtraTracks)return iCount;
1892       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
1893         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
1894         if (!track) continue;
1895
1896         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
1897         if(!trackAOD)continue;
1898         Bool_t bGood = false;
1899         if(fFilterType == 0)bGood = true;
1900         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
1901         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
1902         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
1903         if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
1904         if(trackAOD->Pt()<fTrackPtCut) continue;
1905         list->Add(trackAOD);
1906         iCount++;
1907       }
1908     }
1909   }
1910   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
1911     AliMCEvent* mcEvent = MCEvent();
1912     if(!mcEvent)return iCount;
1913     // we want to have alivpartilces so use get track
1914     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1915       if(!mcEvent->IsPhysicalPrimary(it))continue;
1916       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1917       if(type == kTrackKineAll){
1918         if(part->Pt()<fTrackPtCut)continue;
1919         list->Add(part);
1920         iCount++;
1921       }
1922       else if(type == kTrackKineCharged){
1923         if(part->Particle()->GetPDG()->Charge()==0)continue;
1924         if(part->Pt()<fTrackPtCut)continue;
1925         list->Add(part);
1926         iCount++;
1927       }
1928     }
1929   }
1930   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1931     AliAODEvent *aod = 0;
1932     if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1933     else aod = AODEvent();
1934     if(!aod)return iCount;
1935     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1936     if(!tca)return iCount;
1937     for(int it = 0;it < tca->GetEntriesFast();++it){
1938       AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
1939       if(!part->IsPhysicalPrimary())continue;
1940       if(type == kTrackAODMCAll){
1941         if(part->Pt()<fTrackPtCut)continue;
1942         list->Add(part);
1943         iCount++;
1944       }
1945       else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1946         if(part->Charge()==0)continue;
1947         if(part->Pt()<fTrackPtCut)continue;
1948         if(kTrackAODMCCharged){
1949           list->Add(part);
1950         }
1951         else {
1952           if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
1953           list->Add(part);
1954         }
1955         iCount++;
1956       }
1957     }
1958   }// AODMCparticle
1959   list->Sort();
1960   return iCount;
1961 }
1962
1963 void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
1964
1965   TFile *f = TFile::Open(fPathTrPtResolution.Data());
1966
1967   if(!f)return;
1968
1969   TProfile *fProfPtPtSigma1PtGlobSt     = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
1970   TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
1971   TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
1972
1973   SetSmearResolution(kTRUE);
1974   SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
1975   
1976
1977 }
1978
1979 void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
1980
1981   TFile *f = TFile::Open(fPathTrEfficiency.Data());
1982   if(!f)return;
1983
1984   TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
1985   TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
1986   TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
1987
1988   SetDiceEfficiency(kTRUE);
1989
1990   if(fChangeEfficiencyFraction>0.) {
1991
1992     TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
1993     
1994     for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
1995       Double_t content = hEffPosGlobSt->GetBinContent(bin);
1996       hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
1997     }
1998
1999     SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2000
2001   }
2002   else
2003     SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
2004
2005 }
2006
2007 void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
2008
2009   //
2010   // set mom res profiles
2011   //
2012
2013   if(fMomResH1) delete fMomResH1;
2014   if(fMomResH2) delete fMomResH2;
2015   if(fMomResH3) delete fMomResH3;
2016
2017   fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
2018   fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
2019   fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
2020
2021 }
2022
2023 void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
2024   //
2025   // set tracking efficiency histos
2026   //
2027
2028   fhEffH1 = (TH1*)h1->Clone("fhEffH1");
2029   fhEffH2 = (TH1*)h2->Clone("fhEffH2");
2030   fhEffH3 = (TH1*)h3->Clone("fhEffH3");
2031 }
2032
2033 Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
2034
2035   //
2036   // Get smearing on generated momentum
2037   //
2038
2039   //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
2040
2041   TProfile *fMomRes = 0x0;
2042   if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
2043   if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
2044   if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
2045
2046   if(!fMomRes) {
2047     return 0.;
2048   }
2049
2050
2051   Double_t smear = 0.;
2052
2053   if(pt>20.) {
2054     if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
2055     if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
2056     if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
2057   }
2058   else {
2059
2060     Int_t bin = fMomRes->FindBin(pt);
2061
2062     smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
2063
2064   }
2065
2066   if(fMomRes) delete fMomRes;
2067   
2068   return smear;
2069 }
2070
2071 void AliAnalysisTaskJetCluster::FitMomentumResolution() {
2072   //
2073   // Fit linear function on momentum resolution at high pT
2074   //
2075
2076   if(!fMomResH1Fit && fMomResH1) {
2077     fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
2078     fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
2079     fMomResH1Fit ->SetRange(5.,100.);
2080   }
2081
2082   if(!fMomResH2Fit && fMomResH2) {
2083     fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
2084     fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
2085     fMomResH2Fit ->SetRange(5.,100.);
2086   }
2087
2088   if(!fMomResH3Fit && fMomResH3) {
2089     fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
2090     fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
2091     fMomResH3Fit ->SetRange(5.,100.);
2092   }
2093
2094 }
2095
2096 /*
2097 Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
2098   for(int i = 0; i < particles.GetEntries(); i++){
2099     AliVParticle *vp = (AliVParticle*)particles.At(i);
2100     // Carefull energy is not well determined in real data, should not matter for p_T scheme?
2101     fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E());
2102     jInp.set_user_index(i);
2103     inputParticles.push_back(jInp);
2104   }
2105
2106   return 0;
2107
2108 }
2109 */