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