7531f12de64033c7fe17f9abbd5f662e98de5f8f
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaParticleJetFinderCorrelation.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //_________________________________________________________________________
17 // Class for the analysis of particle (direct gamma) -jet (jet found with finder) correlations
18 //*-- Author: Gustavo Conesa (LNF-INFN)
19 //*-- Modified: Adam Matyja (INP PAN Cracow) 
20 //////////////////////////////////////////////////////////////////////////////
21
22
23 // --- ROOT system ---
24 #include "TH2F.h"
25 #include "TClonesArray.h"
26 #include "TClass.h"
27 //#include "Riostream.h"
28
29 //---- AliRoot system ----
30 #include "AliCaloTrackReader.h"
31 #include "AliAODJet.h"
32 #include "AliAnaParticleJetFinderCorrelation.h" 
33 #include "AliAODPWG4ParticleCorrelation.h"
34 #include "AliVTrack.h"
35 #include "AliAODCaloCluster.h"
36 #include "AliAODEvent.h"
37
38 //jets
39 #include "AliAODJetEventBackground.h"
40 #include "TRandom2.h"
41 //MC access
42 #include "AliStack.h"
43 #include "AliMCAnalysisUtils.h"
44 #include "AliAODMCParticle.h"
45 // --- Detectors ---
46 #include "AliEMCALGeometry.h"
47
48 ClassImp(AliAnaParticleJetFinderCorrelation)
49
50 //________________________________________________________________________
51 AliAnaParticleJetFinderCorrelation::AliAnaParticleJetFinderCorrelation() : 
52 AliAnaCaloTrackCorrBaseClass(),  
53   fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fRatioMaxCut(0.),  fRatioMinCut(0.), 
54   fConeSize(0.), fPtThresholdInCone(0.),fUseJetRefTracks(kTRUE),
55   fMakeCorrelationInHistoMaker(kFALSE), fSelectIsolated(kTRUE),
56   fJetConeSize(0.4),fJetMinPt(5),fJetAreaFraction(0.6),
57 //fNonStandardJetFromReader(kTRUE), 
58   fJetBranchName("jets"),
59   fBackgroundJetFromReader(kTRUE),
60   fBkgJetBranchName("jets"),
61   fGammaConeSize(0.3),fUseBackgroundSubtractionGamma(kFALSE),fSaveGJTree(kTRUE),
62   fMostEnergetic(kFALSE),fMostOpposite(kTRUE), fUseHistogramJetBkg(kTRUE),
63   fUseHistogramTracks(kTRUE),fUseHistogramJetTracks(kTRUE),fMCStudies(kFALSE),fGenerator(0),
64   fMomentum(),
65   fhDeltaEta(0), /*fhDeltaPhi(0),*/fhDeltaPhiCorrect(0),fhDeltaPhi0PiCorrect(0), fhDeltaPt(0), fhPtRatio(0), fhPt(0),
66   fhFFz(0),fhFFxi(0),fhFFpt(0),fhNTracksInCone(0),
67   fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetFFzCor(0),fhJetFFxiCor(0),
68   fhBkgFFz(),fhBkgFFxi(),fhBkgFFpt(),fhBkgNTracksInCone(),fhBkgSumPtInCone(),fhBkgSumPtnTracksInCone(),//<<---new
69   fhNjetsNgammas(0),fhCuts(0),
70   fhDeltaEtaBefore(0),fhDeltaPhiBefore(0),fhDeltaPtBefore(0),fhPtRatioBefore(0),
71   fhPtBefore(0),fhDeltaPhi0PiCorrectBefore(0),
72   fhJetPtBefore(0),fhJetPt(0),fhJetPtMostEne(0),fhJetPhi(0),fhJetEta(0),fhJetEtaVsPt(0),
73   fhJetPhiVsEta(0),fhJetEtaVsNpartInJet(0),fhJetEtaVsNpartInJetBkg(0),fhJetChBkgEnergyVsPt(0),fhJetChAreaVsPt(0),/*fhJetNjet(0),*/
74   fhTrackPhiVsEta(0),fhTrackAveTrackPt(0),fhJetNjetOverPtCut(),
75 /*fhJetChBkgEnergyVsPtEtaGt05(0),fhJetChBkgEnergyVsPtEtaLe05(0),fhJetChAreaVsPtEtaGt05(0),fhJetChAreaVsPtEtaLe05(0),*/
76   fhJetChBkgEnergyVsArea(0),fhJetRhoVsPt(0),fhJetRhoVsCentrality(0),//fhJetBkgRho(0),
77   fhJetNparticlesInJet(0),fhJetDeltaEtaDeltaPhi(0),fhJetDeltaEtaDeltaPhiAllTracks(0),
78   fhJetAveTrackPt(0),fhJetNtracksInJetAboveThr(),fhJetRatioNTrkAboveToNTrk(),fhJetNtrackRatioMostEne(),
79   fhJetNtrackRatioJet5GeV(),fhJetNtrackRatioLead5GeV(),
80   fhBkgJetBackground(),fhBkgJetSigma(),fhBkgJetArea(),fhPhotonPtMostEne(0),
81   fhPhotonAverageEnergy(0),fhPhotonRatioAveEneToMostEne(0),fhPhotonAverageEnergyMinus1(0),fhPhotonRatioAveEneMinus1ToMostEne(0),
82   fhPhotonNgammaMoreAverageToNgamma(0),fhPhotonNgammaMoreAverageMinus1ToNgamma(0),fhPhotonNgammaOverPtCut(),
83   fhPhotonBkgRhoVsNtracks(0),fhPhotonBkgRhoVsNclusters(0),fhPhotonBkgRhoVsCentrality(0),
84   fhPhotonBkgRhoVsNcells(0),fhPhotonPt(0),fhPhotonPtCorrected(0),fhPhotonPtCorrectedZoom(0),fhPhotonPtDiff(0),
85   fhPhotonPtDiffVsCentrality(0),fhPhotonPtDiffVsNcells(0),fhPhotonPtDiffVsNtracks(0),fhPhotonPtDiffVsNclusters(0),
86   fhPhotonSumPtInCone(0),fhPhotonSumPtCorrectInCone(0),fhPhotonSumPtChargedInCone(0),
87   fhSelectedJetPhiVsEta(0),fhSelectedJetChBkgEnergyVsPtJet(0),fhSelectedJetChAreaVsPtJet(0),fhSelectedJetNjet(0),fhSelectedNtracks(0),
88   fhSelectedTrackPhiVsEta(0),fhCuts2(0),
89   fhSelectedPhotonNLMVsPt(0),fhSelectedPhotonLambda0VsPt(0), fhRandomPhiEta(),
90   fhMCPhotonCuts(0),fhMCPhotonPt(0),fhMCPhotonEtaPhi(0),fhMCJetOrigin(0),
91   fhMCJetNPartVsPt(0),fhMCJetChNPartVsPt(0),fhMCJetNPart150VsPt(0),fhMCJetChNPart150VsPt(0),fhMCJetChNPart150ConeVsPt(0),
92   fhMCJetRatioChFull(0),fhMCJetRatioCh150Ch(0),
93   fhMCJetEtaPhi(0),fhMCJetChEtaPhi(0),fhMCJet150EtaPhi(0),fhMCJetCh150EtaPhi(0),fhMCJetCh150ConeEtaPhi(0),
94 fTreeGJ     (0),
95 fGamPt      (0),
96 fGamLambda0 (0),
97 fGamNLM     (0),
98 fGamSumPtCh (0),
99 fGamTime    (0),
100 fGamNcells  (0),
101 fGamEta     (0),
102 fGamPhi     (0),
103 fGamSumPtNeu(0),
104 fGamNtracks (0),
105 fGamNclusters(0),
106 fGamAvEne   (0),
107 fJetPhi     (0),
108 fJetEta     (0),
109 fJetPt      (0),
110 fJetBkgChEne(0),
111 fJetArea    (0),
112 fJetNtracks (0),
113 fJetNtracks1(0),
114 fJetNtracks2(0),
115 fJetRho(0),
116 fEventNumber(0),
117 fNtracks    (0),
118 fZvertex    (0),
119 fCentrality (0),
120 fIso(0),
121 fGamRho(0),
122 fMCGamPt        (0),
123 fMCGamEta       (0),
124 fMCGamPhi       (0),
125 fMCPartonType   (0),
126 fMCJetPt        (0),
127 fMCJetChPt      (0),
128 fMCJet150Pt     (0),
129 fMCJetCh150Pt   (0),
130 fMCJetNPart     (0),
131 fMCJetChNPart   (0),
132 fMCJet150NPart  (0),
133 fMCJetCh150NPart(0),
134 fMCJetEta       (0),
135 fMCJetPhi       (0),
136 fMCJetChEta     (0),
137 fMCJetChPhi     (0),
138 fMCJet150Eta    (0),
139 fMCJet150Phi    (0),
140 fMCJetCh150Eta  (0),
141 fMCJetCh150Phi  (0),
142 fMCJetCh150ConePt(0),
143 fMCJetCh150ConeNPart(0),
144 fMCJetCh150ConeEta(0),
145 fMCJetCh150ConePhi(0)
146
147 {
148   //Default Ctor
149   //printf("constructor\n");
150   
151   //Initialize parameters
152   InitParameters();
153   for(Int_t i=0;i<10;i++){
154     fhJetNjetOverPtCut[i] = 0;
155     fhPhotonNgammaOverPtCut[i] = 0;
156   }
157   fGenerator = new TRandom2();
158   fGenerator->SetSeed(0);
159 }
160
161 //___________________________________________________________________
162 AliAnaParticleJetFinderCorrelation::~AliAnaParticleJetFinderCorrelation(){
163   delete fGenerator;
164 }
165
166
167 //___________________________________________________________________
168 TList *  AliAnaParticleJetFinderCorrelation::GetCreateOutputObjects()
169 {  
170   // Create histograms to be saved in output file and 
171   // store them in fOutputContainer
172   //printf("GetCreateOutputObjects\n");    
173
174   TList * outputContainer = new TList() ; 
175   outputContainer->SetName("ParticleJetFinderHistos") ; 
176   
177   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();
178   //    Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
179   //    Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
180   Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();
181   //    Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
182   //    Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
183   Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
184   //    Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
185 //      Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();        
186   
187 //  fhDeltaPhi  = new TH2F("DeltaPhi","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-6,4); 
188 //  fhDeltaPhi->SetYTitle("#Delta #phi");
189 //  fhDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
190 //  outputContainer->Add(fhDeltaPhi);
191
192   fhDeltaPhiCorrect  = new TH2F("DeltaPhiCorrect","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5); 
193   fhDeltaPhiCorrect->SetYTitle("#Delta #phi");
194   fhDeltaPhiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
195   outputContainer->Add(fhDeltaPhiCorrect);
196
197   fhDeltaPhi0PiCorrect  = new TH2F("DeltaPhi0PiCorrect","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5); 
198   fhDeltaPhi0PiCorrect->SetYTitle("#Delta #phi");
199   fhDeltaPhi0PiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
200   outputContainer->Add(fhDeltaPhi0PiCorrect);
201
202
203   fhDeltaEta  = new TH2F("DeltaEta","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2); 
204   fhDeltaEta->SetYTitle("#Delta #eta");
205   fhDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
206   outputContainer->Add(fhDeltaEta);
207   
208   fhDeltaPt  = new TH2F("DeltaPt","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,150,-50,100); 
209   fhDeltaPt->SetYTitle("#Delta p_{T}");
210   fhDeltaPt->SetXTitle("p_{T trigger} (GeV/c)"); 
211   outputContainer->Add(fhDeltaPt);
212   
213   fhPtRatio  = new TH2F("PtRatio","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.); 
214   fhPtRatio->SetYTitle("ratio");
215   fhPtRatio->SetXTitle("p_{T trigger} (GeV/c)");
216   outputContainer->Add(fhPtRatio);
217   
218   fhPt  = new TH2F("Pt","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
219   fhPt->SetYTitle("p_{T jet}(GeV/c)");
220   fhPt->SetXTitle("p_{T trigger} (GeV/c)");
221   outputContainer->Add(fhPt);
222   
223   fhFFz  = new TH2F("FFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,2);  
224   fhFFz->SetYTitle("z");
225   fhFFz->SetXTitle("p_{T trigger}");
226   outputContainer->Add(fhFFz) ;
227         
228   fhFFxi  = new TH2F("FFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,10.); 
229   fhFFxi->SetYTitle("#xi");
230   fhFFxi->SetXTitle("p_{T trigger}");
231   outputContainer->Add(fhFFxi) ;
232   
233   fhFFpt  = new TH2F("FFpt","p_{T i charged} vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,50.); 
234   fhFFpt->SetYTitle("p_{T charged hadron}");
235   fhFFpt->SetXTitle("p_{T trigger}");
236   outputContainer->Add(fhFFpt) ;
237   
238   fhNTracksInCone  = new TH2F("NTracksInCone","Number of tracks in cone vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,150.); 
239   fhNTracksInCone->SetYTitle("p_{T charged hadron}");
240   fhNTracksInCone->SetXTitle("p_{T trigger}");
241   outputContainer->Add(fhNTracksInCone) ;
242   
243   //FF according to jet axis
244   fhJetFFz  = new TH2F("JetFFz","z = p_{T i charged}/p_{T jet} vs p_{T jet}",nptbins,ptmin,ptmax,200,0.,2);  
245   fhJetFFz->SetYTitle("z");
246   fhJetFFz->SetXTitle("p_{T jet}");
247   outputContainer->Add(fhJetFFz) ;
248         
249   fhJetFFxi  = new TH2F("JetFFxi","#xi = ln(p_{T jet}/p_{T i charged}) vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,10.); 
250   fhJetFFxi->SetYTitle("#xi");
251   fhJetFFxi->SetXTitle("p_{T jet}");
252   outputContainer->Add(fhJetFFxi) ;
253   
254   fhJetFFpt  = new TH2F("JetFFpt","p_{T i charged} vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,50.); 
255   fhJetFFpt->SetYTitle("p_{T charged hadron}");
256   fhJetFFpt->SetXTitle("p_{T jet}");
257   outputContainer->Add(fhJetFFpt) ;
258
259   fhJetFFzCor  = new TH2F("JetFFzCor","z = -cos(#alpha(jet,trig))*p_{T i charged}/p_{T jet} vs p_{T jet}",nptbins,ptmin,ptmax,200,0.,2);  
260   fhJetFFzCor->SetYTitle("z");
261   fhJetFFzCor->SetXTitle("p_{T jet}");
262   outputContainer->Add(fhJetFFzCor) ;
263         
264   fhJetFFxiCor  = new TH2F("JetFFxiCor","#xi = ln(p_{T jet}/(-cos(#alpha(jet,trig))*p_{T i charged})) vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,10.); 
265   fhJetFFxiCor->SetYTitle("#xi");
266   fhJetFFxiCor->SetXTitle("p_{T jet}");
267   outputContainer->Add(fhJetFFxiCor) ;
268
269
270   //background FF
271   fhBkgFFz[0]  = new TH2F("BkgFFzRC",  "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg RC"  ,nptbins,ptmin,ptmax,200,0.,2);  
272   fhBkgFFz[1]  = new TH2F("BkgFFzPCG", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg PCG" ,nptbins,ptmin,ptmax,200,0.,2);  
273   fhBkgFFz[2]  = new TH2F("BkgFFzPCJ", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg PCJ" ,nptbins,ptmin,ptmax,200,0.,2);  
274   fhBkgFFz[3]  = new TH2F("BkgFFzMP",  "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg MP"  ,nptbins,ptmin,ptmax,200,0.,2);  
275   fhBkgFFz[4]  = new TH2F("BkgFFzTest","z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg Test",nptbins,ptmin,ptmax,200,0.,2);  
276   for(Int_t i=0;i<5;i++){
277     fhBkgFFz[i]->SetYTitle("z");
278     fhBkgFFz[i]->SetXTitle("p_{T trigger}");
279     outputContainer->Add(fhBkgFFz[i]) ;
280   }
281
282   fhBkgFFxi[0]  = new TH2F("BkgFFxiRC",  "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg RC",  nptbins,ptmin,ptmax,100,0.,10.); 
283   fhBkgFFxi[1]  = new TH2F("BkgFFxiPCG", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,10.); 
284   fhBkgFFxi[2]  = new TH2F("BkgFFxiPCJ", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,10.); 
285   fhBkgFFxi[3]  = new TH2F("BkgFFxiMP",  "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg MP",  nptbins,ptmin,ptmax,100,0.,10.); 
286   fhBkgFFxi[4]  = new TH2F("BkgFFxiTest","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg Test",nptbins,ptmin,ptmax,100,0.,10.); 
287   for(Int_t i=0;i<5;i++){
288     fhBkgFFxi[i]->SetYTitle("#xi");
289     fhBkgFFxi[i]->SetXTitle("p_{T trigger}");
290     outputContainer->Add(fhBkgFFxi[i]) ;
291   }
292
293   fhBkgFFpt[0]  = new TH2F("BkgFFptRC",  "p_{T i charged} vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,50.); 
294   fhBkgFFpt[1]  = new TH2F("BkgFFptPCG", "p_{T i charged} vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,50.); 
295   fhBkgFFpt[2]  = new TH2F("BkgFFptPCJ", "p_{T i charged} vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,50.); 
296   fhBkgFFpt[3]  = new TH2F("BkgFFptMP",  "p_{T i charged} vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,50.); 
297   fhBkgFFpt[4]  = new TH2F("BkgFFptTest","p_{T i charged} vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,50.); 
298   for(Int_t i=0;i<5;i++){
299     fhBkgFFpt[i]->SetYTitle("p_{T charged hadron}");
300     fhBkgFFpt[i]->SetXTitle("p_{T trigger}");
301     outputContainer->Add(fhBkgFFpt[i]) ;
302   }
303
304   fhBkgNTracksInCone[0]  = new TH2F("BkgNTracksInConeRC",  "Number of tracks in cone vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,150.); 
305   fhBkgNTracksInCone[1]  = new TH2F("BkgNTracksInConePCG", "Number of tracks in cone vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,150.); 
306   fhBkgNTracksInCone[2]  = new TH2F("BkgNTracksInConePCJ", "Number of tracks in cone vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,150.); 
307   fhBkgNTracksInCone[3]  = new TH2F("BkgNTracksInConeMP",  "Number of tracks in cone vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,150.); 
308   fhBkgNTracksInCone[4]  = new TH2F("BkgNTracksInConeTest","Number of tracks in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,150.); 
309   for(Int_t i=0;i<5;i++){
310     fhBkgNTracksInCone[i]->SetYTitle("Number of tracks");
311     fhBkgNTracksInCone[i]->SetXTitle("p_{T trigger}");
312     outputContainer->Add(fhBkgNTracksInCone[i]) ;
313   }
314
315   fhBkgSumPtInCone[0]  = new TH2F("BkgSumPtInConeRC",  "Sum P_{T} in cone vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,100.);
316   fhBkgSumPtInCone[1]  = new TH2F("BkgSumPtInConePCG", "Sum P_{T} in cone vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,100.);
317   fhBkgSumPtInCone[2]  = new TH2F("BkgSumPtInConePCJ", "Sum P_{T} in cone vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,100.);
318   fhBkgSumPtInCone[3]  = new TH2F("BkgSumPtInConeMP",  "Sum P_{T} in cone vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,100.);
319   fhBkgSumPtInCone[4]  = new TH2F("BkgSumPtInConeTest","Sum P_{T} in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,100.);
320   for(Int_t i=0;i<5;i++){
321     fhBkgSumPtInCone[i]->SetYTitle("Sum P_{T}");
322     fhBkgSumPtInCone[i]->SetXTitle("p_{T trigger}");
323     outputContainer->Add(fhBkgSumPtInCone[i]) ;
324   }
325
326   fhBkgSumPtnTracksInCone[0]  = new TH2F("BkgSumPtnTracksInConeRC",  "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,20.);
327   fhBkgSumPtnTracksInCone[1]  = new TH2F("BkgSumPtnTracksInConePCG", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,20.);
328   fhBkgSumPtnTracksInCone[2]  = new TH2F("BkgSumPtnTracksInConePCJ", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,20.);
329   fhBkgSumPtnTracksInCone[3]  = new TH2F("BkgSumPtnTracksInConeMP",  "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,20.);
330   fhBkgSumPtnTracksInCone[4]  = new TH2F("BkgSumPtnTracksInConeTest","Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,20.);
331   for(Int_t i=0;i<5;i++){
332     fhBkgSumPtnTracksInCone[i]->SetYTitle("Sum p_{T}/Number of tracks");
333     fhBkgSumPtnTracksInCone[i]->SetXTitle("p_{T trigger}");
334     outputContainer->Add(fhBkgSumPtnTracksInCone[i]) ;
335   }
336
337
338   //temporary histograms
339   fhNjetsNgammas  = new TH2F("NjetsNgammas"," Number of jets vs number of gammas in event",20,0.,100.,10,0.,80.);  
340   fhNjetsNgammas->SetYTitle("Number of gammas");
341   fhNjetsNgammas->SetXTitle("Number of jets");
342   outputContainer->Add(fhNjetsNgammas) ;
343
344   fhCuts  = new TH1F("Cuts"," Cuts",10,0.,10.);  
345   fhCuts->SetYTitle("Counts");
346   fhCuts->SetXTitle("Cut number");
347   outputContainer->Add(fhCuts) ;
348
349   fhDeltaPhiBefore  = new TH2F("DeltaPhiBefore","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5); 
350   fhDeltaPhiBefore->SetYTitle("#Delta #phi");
351   fhDeltaPhiBefore->SetXTitle("p_{T trigger} (GeV/c)");
352   outputContainer->Add(fhDeltaPhiBefore);
353   
354   fhDeltaEtaBefore  = new TH2F("DeltaEtaBefore","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2); 
355   fhDeltaEtaBefore->SetYTitle("#Delta #eta");
356   fhDeltaEtaBefore->SetXTitle("p_{T trigger} (GeV/c)");
357   outputContainer->Add(fhDeltaEtaBefore);
358   
359   fhDeltaPtBefore  = new TH2F("DeltaPtBefore","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-50,50); 
360   fhDeltaPtBefore->SetYTitle("#Delta p_{T}");
361   fhDeltaPtBefore->SetXTitle("p_{T trigger} (GeV/c)"); 
362   outputContainer->Add(fhDeltaPtBefore);
363   
364   fhPtRatioBefore  = new TH2F("PtRatioBefore","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.); 
365   fhPtRatioBefore->SetYTitle("ratio");
366   fhPtRatioBefore->SetXTitle("p_{T trigger} (GeV/c)");
367   outputContainer->Add(fhPtRatioBefore);
368   
369   fhPtBefore  = new TH2F("PtBefore","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
370   fhPtBefore->SetYTitle("p_{T jet}(GeV/c)");
371   fhPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
372   outputContainer->Add(fhPtBefore);
373
374   fhDeltaPhi0PiCorrectBefore  = new TH2F("DeltaPhi0PiCorrectBefore","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5); 
375   fhDeltaPhi0PiCorrectBefore->SetYTitle("#Delta #phi");
376   fhDeltaPhi0PiCorrectBefore->SetXTitle("p_{T trigger} (GeV/c)");
377   outputContainer->Add(fhDeltaPhi0PiCorrectBefore);
378
379   //temporary jet histograms
380   fhJetPtBefore            = new TH1F("JetPtBefore","JetPtBefore",150,-50,100); 
381   fhJetPtBefore->SetYTitle("Counts");
382   fhJetPtBefore->SetXTitle("p_{T jet}(GeV/c)");
383   outputContainer->Add(fhJetPtBefore) ;
384
385   fhJetPt            = new TH1F("JetPt","JetPt",150,-50,100); 
386   fhJetPt->SetYTitle("Counts");
387   fhJetPt->SetXTitle("p_{T jet}(GeV/c)");
388   outputContainer->Add(fhJetPt) ;
389
390   fhJetPtMostEne            = new TH1F("JetPtMostEne","JetPtMostEne",150,0,150); 
391   fhJetPtMostEne->SetYTitle("Counts");
392   fhJetPtMostEne->SetXTitle("p_{T jet}(GeV/c)");
393   outputContainer->Add(fhJetPtMostEne) ;
394
395   fhJetPhi           = new TH1F("JetPhi","JetPhi",130,0,6.5); 
396   fhJetPhi->SetYTitle("Counts");
397   fhJetPhi->SetXTitle("#phi_{jet}");
398   outputContainer->Add(fhJetPhi) ;
399
400   fhJetEta           = new TH1F("JetEta","JetEta",100,-1,1); 
401   fhJetEta->SetYTitle("Counts");
402   fhJetEta->SetXTitle("#eta_{jet}");
403   outputContainer->Add(fhJetEta) ;
404
405   fhJetEtaVsPt      = new TH2F("JetEtaVsPt","JetEtaVsPt",100,0,100,50,-1,1);
406   fhJetEtaVsPt->SetYTitle("#eta_{jet}");
407   fhJetEtaVsPt->SetXTitle("p_{T,jet}(GeV/c)");
408   outputContainer->Add(fhJetEtaVsPt) ;
409
410   fhJetPhiVsEta      = new TH2F("JetPhiVsEta","JetPhiVsEta",65,0,6.5,50,-1,1); 
411   fhJetPhiVsEta->SetYTitle("#eta_{jet}");
412   fhJetPhiVsEta->SetXTitle("#phi_{jet}");
413   outputContainer->Add(fhJetPhiVsEta) ;
414
415   fhJetEtaVsNpartInJet= new TH2F("JetEtaVsNpartInJet","JetEtaVsNpartInJet",50,-1,1,100,0.,200.); 
416   fhJetEtaVsNpartInJet->SetYTitle("N_{tracks-in-jet}");
417   fhJetEtaVsNpartInJet->SetXTitle("#eta_{jet}");
418   outputContainer->Add(fhJetEtaVsNpartInJet) ;
419
420   fhJetEtaVsNpartInJetBkg= new TH2F("JetEtaVsNpartInJetBkg","JetEtaVsNpartInJetBkg",50,-1,1,100,0.,200.); 
421   fhJetEtaVsNpartInJetBkg->SetYTitle("N_{tracks-in-jet}");
422   fhJetEtaVsNpartInJetBkg->SetXTitle("#eta_{jet}");
423   outputContainer->Add(fhJetEtaVsNpartInJetBkg) ;
424
425   fhJetChBkgEnergyVsPt = new TH2F("JetBkgChEnergyVsPt","JetBkgChEnergyVsPt",100,0,100,90,0,90); 
426   fhJetChBkgEnergyVsPt->SetYTitle("Jet Bkg Energy (GeV)");
427   fhJetChBkgEnergyVsPt->SetXTitle("p_{T,jet} (GeV/c)");
428   outputContainer->Add(fhJetChBkgEnergyVsPt);
429   
430   fhJetChAreaVsPt      = new TH2F("JetChAreaVsPt","JetChAreaVsPt",100,0,100,50,0,1); 
431   fhJetChAreaVsPt->SetYTitle("Jet Area");
432   fhJetChAreaVsPt->SetXTitle("p_{T,jet} (GeV/c)");
433   outputContainer->Add(fhJetChAreaVsPt);
434   
435   if(IsHistogramTracks()){
436     fhTrackPhiVsEta      = new TH2F("TrackPhiVsEta","TrackPhiVsEta",65,0,6.5,50,-1,1); 
437     fhTrackPhiVsEta->SetYTitle("#eta_{track}");
438     fhTrackPhiVsEta->SetXTitle("#phi_{track}");
439     outputContainer->Add(fhTrackPhiVsEta) ;
440
441     fhTrackAveTrackPt      = new TH1F("TrackAveTrackPt","TrackAveTrackPt",45,0,1.5);
442     fhTrackAveTrackPt->SetYTitle("Counts");
443     fhTrackAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
444     outputContainer->Add(fhTrackAveTrackPt);
445   
446   }//end of IsHistogramTracks()
447
448   for(Int_t i=0;i<10;i++){
449     fhJetNjetOverPtCut[i]      = new TH1F(Form("JetNjetOverPtCut%d", i),Form("JetNjetOverPtCut%d", i),100,0,100);
450     fhJetNjetOverPtCut[i]->SetYTitle("Counts");
451     fhJetNjetOverPtCut[i]->SetXTitle("N_{jets} over threshold");
452     outputContainer->Add(fhJetNjetOverPtCut[i]);
453   }
454
455   fhJetChBkgEnergyVsArea = new TH2F("JetBkgChEnergyVsArea","JetBkgChEnergyVsArea",100,0,100,70,0,0.7); 
456   fhJetChBkgEnergyVsArea->SetXTitle("Jet Bkg Energy (GeV)");
457   fhJetChBkgEnergyVsArea->SetYTitle("Area");
458   outputContainer->Add(fhJetChBkgEnergyVsArea);
459
460   fhJetRhoVsPt           = new TH2F("JetRhoVsPt","JetRhoVsPt",100,0,100,100,0,150); 
461   fhJetRhoVsPt->SetYTitle("Rho");
462   fhJetRhoVsPt->SetXTitle("p_{T,jet} (GeV/c)");
463   outputContainer->Add(fhJetRhoVsPt);
464
465   if(IsHistogramJetBkg()){
466     fhJetRhoVsCentrality           = new TH2F("JetRhoVsCentrality","JetRhoVsCentrality",100,0,100,100,0,200);
467     fhJetRhoVsCentrality->SetYTitle("Rho");
468     fhJetRhoVsCentrality->SetXTitle("Centrality");
469     outputContainer->Add(fhJetRhoVsCentrality);
470   }
471
472   fhJetNparticlesInJet           = new TH1F("JetNparticlesInJet","JetNparticlesInJet",100,0,200);
473   fhJetNparticlesInJet->SetXTitle("N^{particles}");
474   fhJetNparticlesInJet->SetYTitle("N^{jets}");
475   outputContainer->Add(fhJetNparticlesInJet);
476
477   fhJetDeltaEtaDeltaPhi      = new TH2F("JetDeltaEtaDeltaPhi","#Delta #eta^{jet-track} vs. #Delta #phi^{jet-track} for jet tracks",100,-0.8,0.8,100,-0.8,0.8);
478   fhJetDeltaEtaDeltaPhi->SetXTitle("#Delta #eta^{jet-track}");
479   fhJetDeltaEtaDeltaPhi->SetYTitle("#Delta #phi^{jet-track}");
480   outputContainer->Add(fhJetDeltaEtaDeltaPhi );
481
482
483   fhJetDeltaEtaDeltaPhiAllTracks      = new TH2F("JetDeltaEtaDeltaPhiAllTracks","#Delta #eta^{jet-track} vs. #Delta #phi^{jet-track} for all tracks",100,-3.2,3.2,100,-3.2,3.2);
484   fhJetDeltaEtaDeltaPhiAllTracks->SetXTitle("#Delta #eta^{jet-track}");
485   fhJetDeltaEtaDeltaPhiAllTracks->SetYTitle("#Delta #phi^{jet-track}");
486   outputContainer->Add(fhJetDeltaEtaDeltaPhiAllTracks);
487
488
489   if(IsHistogramJetTracks()){
490     fhJetAveTrackPt           = new TH1F("JetAveTrackPt","JetAveTrackPt",45,0.,1.5);
491     fhJetAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
492     fhJetAveTrackPt->SetYTitle("Counts");
493     outputContainer->Add(fhJetAveTrackPt);
494     
495     for(Int_t i=0;i<6;i++){
496       if(i==0) fhJetNtracksInJetAboveThr[i]      = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,200);
497       else fhJetNtracksInJetAboveThr[i]      = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,100);
498       fhJetNtracksInJetAboveThr[i]->SetXTitle("p_{T,jet} (GeV/c)");
499       fhJetNtracksInJetAboveThr[i]->SetYTitle("N_{tracks} over threshold");
500       outputContainer->Add(fhJetNtracksInJetAboveThr[i]);
501     }
502     
503     for(Int_t i=0;i<5;i++){
504       fhJetRatioNTrkAboveToNTrk[i]      = new TH2F(Form("JetRatioNTrkAboveToNTrk%d", i),Form("JetRatioNTrkAboveToNTrk%d", i),100,0,100,40,0,1);
505       fhJetRatioNTrkAboveToNTrk[i]->SetXTitle("p_{T,jet} (GeV/c)");
506       fhJetRatioNTrkAboveToNTrk[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
507       outputContainer->Add(fhJetRatioNTrkAboveToNTrk[i]);
508       
509       fhJetNtrackRatioMostEne[i]      = new TH2F(Form("JetNtrackRatioMostEne%d", i),Form("JetNtrackRatioMostEne%d", i),100,0,100,40,0,1);
510       fhJetNtrackRatioMostEne[i]->SetXTitle("p_{T,jet} (GeV/c)");
511       fhJetNtrackRatioMostEne[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
512       outputContainer->Add(fhJetNtrackRatioMostEne[i]);
513       
514       fhJetNtrackRatioJet5GeV[i]      = new TH2F(Form("JetNtrackRatioJet5GeV%d", i),Form("JetNtrackRatioJet5GeV%d", i),100,0,100,40,0,1);
515       fhJetNtrackRatioJet5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
516       fhJetNtrackRatioJet5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
517       outputContainer->Add(fhJetNtrackRatioJet5GeV[i]);
518       
519       fhJetNtrackRatioLead5GeV[i]      = new TH2F(Form("JetNtrackRatioLead5GeV%d", i),Form("JetNtrackRatioLead5GeV%d", i),100,0,100,40,0,1);
520       fhJetNtrackRatioLead5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
521       fhJetNtrackRatioLead5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
522       outputContainer->Add(fhJetNtrackRatioLead5GeV[i]);
523     }
524   }//end of if IsHistogramJetTracks
525
526   //temporary background jets histograms
527   if(IsHistogramJetBkg()){
528     for(Int_t i=0;i<4;i++){
529       fhBkgJetBackground[i]      = new TH1F(Form("BkgJetBackground%d", i),Form("BkgJetBackground%d", i),100,0,200);
530       fhBkgJetBackground[i]->SetXTitle("<#rho> (GeV/c)");
531       fhBkgJetBackground[i]->SetYTitle("Counts");
532       outputContainer->Add(fhBkgJetBackground[i]);
533       
534       fhBkgJetSigma[i]      = new TH1F(Form("BkgJetSigma%d", i),Form("BkgJetSigma%d", i),100,0,50);
535       fhBkgJetSigma[i]->SetXTitle("#sigma (GeV/c)");
536       fhBkgJetSigma[i]->SetYTitle("Counts");
537       outputContainer->Add(fhBkgJetSigma[i]);
538       
539       fhBkgJetArea[i]      = new TH1F(Form("BkgJetArea%d", i),Form("BkgJetArea%d", i),100,0,1);
540       fhBkgJetArea[i]->SetXTitle("<A>");
541       fhBkgJetArea[i]->SetYTitle("Counts");
542       outputContainer->Add(fhBkgJetArea[i]);
543     }
544   }
545
546   //temporary photon histograms
547   fhPhotonPtMostEne = new TH1F("PhotonPtMostEne","PhotonPtMostEne",100,0,100);
548   fhPhotonPtMostEne->SetYTitle("Counts");
549   fhPhotonPtMostEne->SetXTitle("p_{T,#gamma} (GeV/c)");
550   outputContainer->Add(fhPhotonPtMostEne);
551
552 //  fhPhotonIndexMostEne = new TH1F("PhotonIndexMostEne","PhotonIndexMostEne",100,0,100);
553 //  fhPhotonIndexMostEne->SetYTitle("Counts");
554 //  fhPhotonIndexMostEne->SetXTitle("Index");
555 //  outputContainer->Add(fhPhotonIndexMostEne);
556
557   fhPhotonAverageEnergy = new TH1F("PhotonAverageEnergy","PhotonAverageEnergy",100,0,10);
558   fhPhotonAverageEnergy->SetYTitle("Counts");
559   fhPhotonAverageEnergy->SetXTitle("p_{T,#gamma} (GeV/c)");
560   outputContainer->Add(fhPhotonAverageEnergy);
561
562   fhPhotonRatioAveEneToMostEne = new TH1F("PhotonRatioAveEneToMostEne","PhotonRatioAveEneToMostEne",100,0,1);
563   fhPhotonRatioAveEneToMostEne->SetYTitle("Counts");
564   fhPhotonRatioAveEneToMostEne->SetXTitle("Ratio");
565   outputContainer->Add(fhPhotonRatioAveEneToMostEne);
566
567   fhPhotonAverageEnergyMinus1 = new TH1F("PhotonAverageEnergyMinus1","PhotonAverageEnergyMinus1",100,0,10);
568   fhPhotonAverageEnergyMinus1->SetYTitle("Counts");
569   fhPhotonAverageEnergyMinus1->SetXTitle("p_{T,#gamma} (GeV/c)");
570   outputContainer->Add(fhPhotonAverageEnergyMinus1);
571
572   fhPhotonRatioAveEneMinus1ToMostEne = new TH1F("PhotonRatioAveEneMinus1ToMostEne","PhotonRatioAveEneMinus1ToMostEne",100,0,1);
573   fhPhotonRatioAveEneMinus1ToMostEne->SetYTitle("Counts");
574   fhPhotonRatioAveEneMinus1ToMostEne->SetXTitle("Ratio");
575   outputContainer->Add(fhPhotonRatioAveEneMinus1ToMostEne);
576
577   fhPhotonNgammaMoreAverageToNgamma = new TH1F("PhotonNgammaMoreAverageToNgamma","PhotonNgammaMoreAverageToNgamma",100,0,1);
578   fhPhotonNgammaMoreAverageToNgamma->SetYTitle("Counts");
579   fhPhotonNgammaMoreAverageToNgamma->SetXTitle("Ratio");
580   outputContainer->Add(fhPhotonNgammaMoreAverageToNgamma);
581
582   fhPhotonNgammaMoreAverageMinus1ToNgamma = new TH1F("PhotonNgammaMoreAverageMinus1ToNgamma","PhotonNgammaMoreAverageMinus1ToNgamma",100,0,1);
583   fhPhotonNgammaMoreAverageMinus1ToNgamma->SetYTitle("Counts");
584   fhPhotonNgammaMoreAverageMinus1ToNgamma->SetXTitle("Ratio");
585   outputContainer->Add(fhPhotonNgammaMoreAverageMinus1ToNgamma);
586
587   for(Int_t i=0;i<10;i++){
588     fhPhotonNgammaOverPtCut[i]      = new TH1F(Form("PhotonNgammaOverPtCut%d",i),Form("PhotonNgammaOverPtCut%d",i),100,0,100);
589     fhPhotonNgammaOverPtCut[i]->SetYTitle("Counts");
590     fhPhotonNgammaOverPtCut[i]->SetXTitle("N_{#gamma} over threshold");
591     outputContainer->Add(fhPhotonNgammaOverPtCut[i]);
592   }
593
594   fhPhotonBkgRhoVsNtracks = new TH2F("PhotonBkgRhoVsNtracks","PhotonBkgRhoVsNtracks",200,0,2500,75,0,1.5);
595   //fhPhotonBkgRhoVsNtracks->SetXTitle("Counts");
596   fhPhotonBkgRhoVsNtracks->SetXTitle("Ntracks");
597   fhPhotonBkgRhoVsNtracks->SetYTitle("Rho");
598   outputContainer->Add(fhPhotonBkgRhoVsNtracks);
599
600   fhPhotonBkgRhoVsNclusters = new TH2F("PhotonBkgRhoVsNclusters","PhotonBkgRhoVsNclusters",50,0,100,75,0,1.5);
601   fhPhotonBkgRhoVsNclusters->SetXTitle("Nclusters");
602   fhPhotonBkgRhoVsNclusters->SetYTitle("Rho");
603   outputContainer->Add(fhPhotonBkgRhoVsNclusters);
604
605   fhPhotonBkgRhoVsCentrality = new TH2F("PhotonBkgRhoVsCentrality","PhotonBkgRhoVsCentrality",100,0,100,75,0,1.5);
606   fhPhotonBkgRhoVsCentrality->SetXTitle("Centrality");
607   fhPhotonBkgRhoVsCentrality->SetYTitle("Rho");
608   outputContainer->Add(fhPhotonBkgRhoVsCentrality);
609
610   fhPhotonBkgRhoVsNcells = new TH2F("PhotonBkgRhoVsNcells","PhotonBkgRhoVsNcells",100,0,200,75,0,1.5);
611   fhPhotonBkgRhoVsNcells->SetXTitle("N_{cells}");
612   fhPhotonBkgRhoVsNcells->SetYTitle("Rho");
613   outputContainer->Add(fhPhotonBkgRhoVsNcells);
614
615   fhPhotonPt = new TH1F("PhotonPt","PhotonPt",220,-10,100);
616   fhPhotonPt->SetXTitle("p_{T,#gamma} (GeV/c)");
617   fhPhotonPt->SetYTitle("Counts");
618   outputContainer->Add(fhPhotonPt);
619
620   fhPhotonPtCorrected = new TH1F("PhotonPtCorrected","PhotonPtCorrected",220,-10,100);
621   fhPhotonPtCorrected->SetXTitle("p_{T,#gamma} (GeV/c)");
622   fhPhotonPtCorrected->SetYTitle("Counts");
623   outputContainer->Add(fhPhotonPtCorrected);
624
625   fhPhotonPtDiff = new TH1F("PhotonPtDiff","PhotonPtDiff",50,0,10);
626   fhPhotonPtDiff->SetXTitle("p_{T,#gamma} (GeV/c)");
627   fhPhotonPtDiff->SetYTitle("Counts");
628   outputContainer->Add(fhPhotonPtDiff);
629
630   fhPhotonPtDiffVsNtracks = new TH2F("PhotonPtDiffVsNtracks","PhotonPtDiffVsNtracks",200,0,2500,50,0,5);
631   fhPhotonPtDiffVsNtracks->SetXTitle("N_{tracks}");
632   fhPhotonPtDiffVsNtracks->SetYTitle("<#rho^{#gamma}>*N_{cells}");
633   outputContainer->Add(fhPhotonPtDiffVsNtracks);
634
635   fhPhotonPtDiffVsNclusters = new TH2F("PhotonPtDiffVsNclusters","PhotonPtDiffVsNclusters",50,0,100,50,0,5);
636   fhPhotonPtDiffVsNclusters->SetXTitle("N_{clusters}");
637   fhPhotonPtDiffVsNclusters->SetYTitle("<#rho^{#gamma}>*N_{cells}");
638   outputContainer->Add(fhPhotonPtDiffVsNclusters);
639
640   fhPhotonPtDiffVsCentrality = new TH2F("PhotonPtDiffVsCentrality","PhotonPtDiffVsCentrality",100,0,100,50,0,5);
641   fhPhotonPtDiffVsCentrality->SetXTitle("Centrality");
642   fhPhotonPtDiffVsCentrality->SetYTitle("<#rho^{#gamma}>*N_{cells}");
643   outputContainer->Add(fhPhotonPtDiffVsCentrality);
644
645   fhPhotonPtDiffVsNcells = new TH2F("PhotonPtDiffVsNcells","PhotonPtDiffVsNcells",100,0,200,50,0,5);
646   fhPhotonPtDiffVsNcells->SetXTitle("N_{cells}");
647   fhPhotonPtDiffVsNcells->SetYTitle("<#rho^{#gamma}>*N_{cells}");
648   outputContainer->Add(fhPhotonPtDiffVsNcells);
649
650
651   fhPhotonPtCorrectedZoom = new TH1F("PhotonPtCorrectedZoom","PhotonPtCorrectedZoom",200,-5,5);
652   fhPhotonPtCorrectedZoom->SetXTitle("p_{T,#gamma} (GeV/c)");
653   fhPhotonPtCorrectedZoom->SetYTitle("Counts");
654   outputContainer->Add(fhPhotonPtCorrectedZoom);
655
656   fhPhotonSumPtInCone = new TH1F("PhotonSumPtInCone","PhotonSumPtInCone",100,0,100);
657   fhPhotonSumPtInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
658   fhPhotonSumPtInCone->SetYTitle("Counts");
659   outputContainer->Add(fhPhotonSumPtInCone);
660
661   fhPhotonSumPtCorrectInCone = new TH1F("PhotonSumPtCorrectInCone","PhotonSumPtCorrectInCone",100,-20,80);
662   fhPhotonSumPtCorrectInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
663   fhPhotonSumPtCorrectInCone->SetYTitle("Counts");
664   outputContainer->Add(fhPhotonSumPtCorrectInCone);
665
666   fhPhotonSumPtChargedInCone = new TH1F("PhotonSumPtChargedInCone","PhotonSumPtChargedInCone",100,0,100);
667   fhPhotonSumPtChargedInCone->SetXTitle("#Sigma p_{T,#gamma}^{ch} (GeV/c)");
668   fhPhotonSumPtChargedInCone->SetYTitle("Counts");
669   outputContainer->Add(fhPhotonSumPtChargedInCone);
670
671
672   //temporary jet histograms after selection
673   fhSelectedJetPhiVsEta      = new TH2F("SelectedJetSelectedPhiVsEta","SelectedJetPhiVsEta",65,0,6.5,50,-1,1); 
674   fhSelectedJetPhiVsEta->SetYTitle("#eta_{jet}");
675   fhSelectedJetPhiVsEta->SetXTitle("#phi_{jet}");
676   outputContainer->Add(fhSelectedJetPhiVsEta) ;
677
678   fhSelectedJetChBkgEnergyVsPtJet = new TH2F("SelectedJetBkgChEnergyVsPtJet","SelectedJetBkgChEnergyVsPtJet",100,0,100,90,0,90); 
679   fhSelectedJetChBkgEnergyVsPtJet->SetYTitle("Jet Bkg Energy (GeV)");
680   fhSelectedJetChBkgEnergyVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
681   outputContainer->Add(fhSelectedJetChBkgEnergyVsPtJet);
682   
683   fhSelectedJetChAreaVsPtJet      = new TH2F("SelectedJetChAreaVsPtJet","SelectedJetChAreaVsPtJet",100,0,100,50,0,1); 
684   fhSelectedJetChAreaVsPtJet->SetYTitle("Jet Area");
685   fhSelectedJetChAreaVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
686   outputContainer->Add(fhSelectedJetChAreaVsPtJet);
687
688   fhSelectedJetNjet      = new TH1F("SelectedJetNjet","SelectedJetNjet",100,0,100); 
689   fhSelectedJetNjet->SetYTitle("Counts");
690   fhSelectedJetNjet->SetXTitle("N_{jets} per event");
691   outputContainer->Add(fhSelectedJetNjet);
692
693   fhSelectedNtracks      = new TH1F("SelectedNtracks","SelectedNtracks",100,0,2000); 
694   fhSelectedNtracks->SetYTitle("Counts");
695   fhSelectedNtracks->SetXTitle("N_{tracks} per event");
696   outputContainer->Add(fhSelectedNtracks);
697
698   fhSelectedTrackPhiVsEta      = new TH2F("SelectedTrackPhiVsEta","SelectedTrackPhiVsEta",65,0,6.5,50,-1,1); 
699   fhSelectedTrackPhiVsEta->SetYTitle("#eta_{track}");
700   fhSelectedTrackPhiVsEta->SetXTitle("#phi_{track}");
701   outputContainer->Add(fhSelectedTrackPhiVsEta) ;
702
703   fhCuts2      = new TH1F("Cuts2","Cuts2",10,0.,10.); 
704   fhCuts2->SetYTitle("Counts");
705   fhCuts2->SetXTitle("Cut number");
706   outputContainer->Add(fhCuts2);
707
708   fhSelectedPhotonNLMVsPt      = new TH2F("SelectedPhotonNLMVsPt","SelectedPhotonNLMVsPt",100,0,100,10,0,10);
709   fhSelectedPhotonNLMVsPt->SetYTitle("NLM");
710   fhSelectedPhotonNLMVsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
711   outputContainer->Add(fhSelectedPhotonNLMVsPt);
712
713   fhSelectedPhotonLambda0VsPt      = new TH2F("SelectedPhotonLambda0VsPt","SelectedPhotonLambda0VsPt",100,0,100,50,0,5);
714   fhSelectedPhotonLambda0VsPt->SetYTitle("#lambda_{0}");
715   fhSelectedPhotonLambda0VsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
716   outputContainer->Add(fhSelectedPhotonLambda0VsPt);
717
718   //random
719   fhRandomPhiEta[0]  = new TH2F("RandomPhiEtaRC",  "RandomPhiEtaRC",            100,0,6.5,100,-1.,1.);
720   fhRandomPhiEta[1]  = new TH2F("RandomPhiEtaPCG", "RandomPhiEtaPerpConePhoton",100,0,6.5,100,-1.,1.);
721   fhRandomPhiEta[2]  = new TH2F("RandomPhiEtaPCJ", "RandomPhiEtaPerpConeJet",   100,0,6.5,100,-1.,1.);
722   fhRandomPhiEta[3]  = new TH2F("RandomPhiEtaMP",  "RandomPhiEtaMidPoint",      100,0,6.5,100,-1.,1.);
723   fhRandomPhiEta[4]  = new TH2F("RandomPhiEtaTest","RandomPhiEtaTest",          100,0,6.5,100,-1.,1.);
724   for(Int_t i=0;i<5;i++){
725     fhRandomPhiEta[i]->SetXTitle("#phi");
726     fhRandomPhiEta[i]->SetYTitle("#eta");
727     outputContainer->Add(fhRandomPhiEta[i]);
728   }
729
730   //MC generated photons and jets information  
731   if(fMCStudies) {
732     fhMCPhotonCuts      = new TH1F("MCPhotonCuts","MCPhotonCuts",10,0.,10.); 
733     fhMCPhotonCuts->SetYTitle("Counts");
734     fhMCPhotonCuts->SetXTitle("Cut number");
735     outputContainer->Add(fhMCPhotonCuts);
736     
737     fhMCPhotonPt      = new TH1F("MCPhotonPt","MCPhotonPt",100,0.,100.); 
738     fhMCPhotonPt->SetYTitle("Counts");
739     fhMCPhotonPt->SetXTitle("p_{T,#gamma}^{gen} (GeV/c)");
740     outputContainer->Add(fhMCPhotonPt);
741     
742     fhMCPhotonEtaPhi      = new TH2F("MCPhotonEtaPhi","MCPhotonEtaPhi",65,0,6.5,50,-1,1); 
743     fhMCPhotonEtaPhi->SetYTitle("#eta_{#gamma}");
744     fhMCPhotonEtaPhi->SetXTitle("#phi_{#gamma}");
745     outputContainer->Add(fhMCPhotonEtaPhi) ;
746     
747     fhMCJetOrigin      = new TH1F("MCJetOrigin","MCJetOrigin",35,-10.,25.); 
748     fhMCJetOrigin->SetYTitle("Counts");
749     fhMCJetOrigin->SetXTitle("PDG number");
750     outputContainer->Add(fhMCJetOrigin);
751     
752     fhMCJetNPartVsPt      = new TH2F("MCJetNPartVsPt","MCJetNPartVsPt",100,0.,100.,100,0.,100.); 
753     fhMCJetNPartVsPt->SetYTitle("N_{particles}");
754     fhMCJetNPartVsPt->SetXTitle("p_{T,full-jet}^{gen} (GeV/c)");
755     outputContainer->Add(fhMCJetNPartVsPt) ;
756     
757     fhMCJetChNPartVsPt      = new TH2F("MCJetChNPartVsPt","MCJetChNPartVsPt",100,0.,100.,100,0.,100.); 
758     fhMCJetChNPartVsPt->SetYTitle("N_{particles}");
759     fhMCJetChNPartVsPt->SetXTitle("p_{T,charged-jet}^{gen} (GeV/c)");
760     outputContainer->Add(fhMCJetChNPartVsPt) ;
761     
762     fhMCJetNPart150VsPt      = new TH2F("MCJetNPart150VsPt","MCJetNPart150VsPt",100,0.,100.,100,0.,100.); 
763     fhMCJetNPart150VsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
764     fhMCJetNPart150VsPt->SetXTitle("p_{T,full-jet}^{gen} (GeV/c)");
765     outputContainer->Add(fhMCJetNPart150VsPt) ;
766     
767     fhMCJetChNPart150VsPt      = new TH2F("MCJetChNPart150VsPt","MCJetChNPart150VsPt",100,0.,100.,100,0.,100.); 
768     fhMCJetChNPart150VsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
769     fhMCJetChNPart150VsPt->SetXTitle("p_{T,charged-jet}^{gen} (GeV/c)");
770     outputContainer->Add(fhMCJetChNPart150VsPt) ;
771     
772     fhMCJetChNPart150ConeVsPt      = new TH2F("MCJetChNPart150ConeVsPt","MCJetChNPart150ConeVsPt",100,0.,100.,50,0.,50.); 
773     fhMCJetChNPart150ConeVsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
774     fhMCJetChNPart150ConeVsPt->SetXTitle("p_{T,charged-jet,R=0.4}^{gen} (GeV/c)");
775     outputContainer->Add(fhMCJetChNPart150ConeVsPt) ;
776
777     fhMCJetRatioChFull      = new TH1F("MCJetRatioChFull","MCJetRatioChFull",100,0.,1.); 
778     fhMCJetRatioChFull->SetYTitle("Counts");
779     fhMCJetRatioChFull->SetXTitle("p_{T,charged-jet}^{gen}/p_{T,full-jet}^{gen}");
780     outputContainer->Add(fhMCJetRatioChFull);
781     
782     fhMCJetRatioCh150Ch      = new TH1F("MCJetRatioCh150Ch","MCJetRatioCh150Ch",100,0.,1.); 
783     fhMCJetRatioCh150Ch->SetYTitle("Counts");
784     fhMCJetRatioCh150Ch->SetXTitle("p_{T,charged-jet,pT>150MeV/c}^{gen}/p_{T,charged-jet}^{gen}");
785     outputContainer->Add(fhMCJetRatioCh150Ch);
786     
787     fhMCJetEtaPhi      = new TH2F("MCJetEtaPhi","MCJetEtaPhi",65,0,6.5,50,-1,1); 
788     fhMCJetEtaPhi->SetYTitle("#eta_{jet}");
789     fhMCJetEtaPhi->SetXTitle("#phi_{jet}");
790     outputContainer->Add(fhMCJetEtaPhi) ;
791     
792     fhMCJetChEtaPhi      = new TH2F("MCJetChEtaPhi","MCJetChEtaPhi",65,0,6.5,50,-1,1); 
793     fhMCJetChEtaPhi->SetYTitle("#eta_{jet}");
794     fhMCJetChEtaPhi->SetXTitle("#phi_{jet}");
795     outputContainer->Add(fhMCJetChEtaPhi) ;
796     
797     fhMCJet150EtaPhi      = new TH2F("MCJet150EtaPhi","MCJet150EtaPhi",65,0,6.5,50,-1,1); 
798     fhMCJet150EtaPhi->SetYTitle("#eta_{jet}");
799     fhMCJet150EtaPhi->SetXTitle("#phi_{jet}");
800     outputContainer->Add(fhMCJet150EtaPhi) ;
801     
802     fhMCJetCh150EtaPhi      = new TH2F("MCJetCh150EtaPhi","MCJetCh150EtaPhi",65,0,6.5,50,-1,1); 
803     fhMCJetCh150EtaPhi->SetYTitle("#eta_{jet}");
804     fhMCJetCh150EtaPhi->SetXTitle("#phi_{jet}");
805     outputContainer->Add(fhMCJetCh150EtaPhi) ;
806
807     fhMCJetCh150ConeEtaPhi      = new TH2F("MCJetCh150ConeEtaPhi","MCJetCh150ConeEtaPhi",65,0,6.5,50,-1,1); 
808     fhMCJetCh150ConeEtaPhi->SetYTitle("#eta_{jet}");
809     fhMCJetCh150ConeEtaPhi->SetXTitle("#phi_{jet}");
810     outputContainer->Add(fhMCJetCh150ConeEtaPhi) ;
811   }
812
813   //tree with data
814   if(fSaveGJTree){
815     fTreeGJ=new TTree("fTreeGJ","fTreeGJ");
816     fTreeGJ->Branch("fGamPt"      ,&fGamPt    ,"fGamPt/D");//! fGamPt
817     fTreeGJ->Branch("fGamLambda0" ,&fGamLambda0  ,"fGamLambda0/D");
818     fTreeGJ->Branch("fGamNLM"     ,&fGamNLM      ,"fGamNLM/I");
819     fTreeGJ->Branch("fGamSumPtCh" ,&fGamSumPtCh  ,"fGamSumPtCh/D");
820     fTreeGJ->Branch("fGamNtracks" ,&fGamNtracks  ,"fGamNtracks/I");
821     fTreeGJ->Branch("fGamTime"    ,&fGamTime     ,"fGamTime/D");
822     fTreeGJ->Branch("fGamNcells"  ,&fGamNcells   ,"fGamNcells/I");
823     fTreeGJ->Branch("fGamEta"     ,&fGamEta      ,"fGamEta/D");
824     fTreeGJ->Branch("fGamPhi"     ,&fGamPhi      ,"fGamPhi/D");
825     fTreeGJ->Branch("fGamSumPtNeu",&fGamSumPtNeu ,"fGamSumPtNeu/D");
826     fTreeGJ->Branch("fGamNclusters" ,&fGamNclusters  ,"fGamNclusters/I");
827     fTreeGJ->Branch("fGamAvEne"   ,&fGamAvEne    ,"fGamAvEne/D");
828     fTreeGJ->Branch("fJetPhi"     ,&fJetPhi      ,"fJetPhi/D");
829     fTreeGJ->Branch("fJetEta"     ,&fJetEta      ,"fJetEta/D");
830     fTreeGJ->Branch("fJetPt"      ,&fJetPt       ,"fJetPt/D");
831     fTreeGJ->Branch("fJetBkgChEne",&fJetBkgChEne ,"fJetBkgChEne/D");
832     fTreeGJ->Branch("fJetArea"    ,&fJetArea     ,"fJetArea/D");
833     fTreeGJ->Branch("fJetNtracks" ,&fJetNtracks  ,"fJetNtracks/I");
834     fTreeGJ->Branch("fJetNtracks1" ,&fJetNtracks1  ,"fJetNtracks1/I");
835     fTreeGJ->Branch("fJetNtracks2" ,&fJetNtracks2  ,"fJetNtracks2/I");
836     fTreeGJ->Branch("fJetRho" ,&fJetRho  ,"fJetRho/D");
837     fTreeGJ->Branch("fEventNumber",&fEventNumber ,"fEventNumber/I");
838     fTreeGJ->Branch("fNtracks"    ,&fNtracks     ,"fNtracks/I");
839     fTreeGJ->Branch("fZvertex"    ,&fZvertex     ,"fZvertex/D");
840     fTreeGJ->Branch("fCentrality" ,&fCentrality  ,"fCentrality/D");
841     fTreeGJ->Branch("fIso" ,&fIso  ,"fIso/O");
842     fTreeGJ->Branch("fGamRho" ,&fGamRho  ,"fGamRho/D");
843
844     fTreeGJ->Branch("fMCGamPt"         ,&fMCGamPt         ,"fMCGamPt/D"        );
845     fTreeGJ->Branch("fMCGamEta"        ,&fMCGamEta        ,"fMCGamEta/D"       );
846     fTreeGJ->Branch("fMCGamPhi"        ,&fMCGamPhi        ,"fMCGamPhi/D"       );
847     fTreeGJ->Branch("fMCPartonType"    ,&fMCPartonType    ,"fMCPartonType/I"   );
848     fTreeGJ->Branch("fMCJetPt"         ,&fMCJetPt         ,"fMCJetPt/D"        );
849     fTreeGJ->Branch("fMCJetChPt"       ,&fMCJetChPt       ,"fMCJetChPt/D"      );
850     fTreeGJ->Branch("fMCJet150Pt"      ,&fMCJet150Pt      ,"fMCJet150Pt/D"     );
851     fTreeGJ->Branch("fMCJetCh150Pt"    ,&fMCJetCh150Pt    ,"fMCJetCh150Pt/D"   );
852     fTreeGJ->Branch("fMCJetNPart"      ,&fMCJetNPart      ,"fMCJetNPart/I"     );
853     fTreeGJ->Branch("fMCJetChNPart"    ,&fMCJetChNPart    ,"fMCJetChNPart/I"   );
854     fTreeGJ->Branch("fMCJet150NPart"   ,&fMCJet150NPart   ,"fMCJet150NPart/I"  );
855     fTreeGJ->Branch("fMCJetCh150NPart" ,&fMCJetCh150NPart ,"fMCJetCh150NPart/I");
856     fTreeGJ->Branch("fMCJetEta"        ,&fMCJetEta        ,"fMCJetEta/D"       );
857     fTreeGJ->Branch("fMCJetPhi"        ,&fMCJetPhi        ,"fMCJetPhi/D"       );
858     fTreeGJ->Branch("fMCJetChEta"      ,&fMCJetChEta      ,"fMCJetChEta/D"     );
859     fTreeGJ->Branch("fMCJetChPhi"      ,&fMCJetChPhi      ,"fMCJetChPhi/D"     );
860     fTreeGJ->Branch("fMCJet150Eta"     ,&fMCJet150Eta     ,"fMCJet150Eta/D"    );
861     fTreeGJ->Branch("fMCJet150Phi"     ,&fMCJet150Phi     ,"fMCJet150Phi/D"    );
862     fTreeGJ->Branch("fMCJetCh150Eta"   ,&fMCJetCh150Eta   ,"fMCJetCh150Eta/D"  );
863     fTreeGJ->Branch("fMCJetCh150Phi"   ,&fMCJetCh150Phi   ,"fMCJetCh150Phi/D"  );  
864
865     fTreeGJ->Branch("fMCJetCh150ConePt"    ,&fMCJetCh150ConePt    ,"fMCJetCh150ConePt/D"  );    
866     fTreeGJ->Branch("fMCJetCh150ConeNPart" ,&fMCJetCh150ConeNPart ,"fMCJetCh150ConeNPart/I");
867     fTreeGJ->Branch("fMCJetCh150ConeEta"   ,&fMCJetCh150ConeEta   ,"fMCJetCh150ConeEta/D"  );   
868     fTreeGJ->Branch("fMCJetCh150ConePhi"   ,&fMCJetCh150ConePhi   ,"fMCJetCh150ConePhi/D"  );   
869     
870     outputContainer->Add(fTreeGJ);
871   }
872
873   return outputContainer;
874
875 }
876
877 //_______________________________________________________
878 void AliAnaParticleJetFinderCorrelation::InitParameters()
879 {
880   //printf("InitParameters\n");
881   //Initialize the parameters of the analysis.
882   SetInputAODName("PWG4Particle");
883   AddToHistogramsName("AnaJetFinderCorr_");
884
885   fDeltaPhiMinCut    = 2.6 ;
886   fDeltaPhiMaxCut    = 3.7 ; 
887   fRatioMaxCut       = 1.2 ; 
888   fRatioMinCut       = 0.3 ; 
889   fConeSize          = 0.4 ;
890   fPtThresholdInCone = 0.5 ;
891   fUseJetRefTracks   = kFALSE ;
892   fMakeCorrelationInHistoMaker = kFALSE ;
893   fSelectIsolated = kTRUE;
894   fJetConeSize = 0.4 ;
895   fJetMinPt = 15. ; //GeV/c
896   fJetAreaFraction = 0.6 ;
897   fJetBranchName = "jets";
898   fBkgJetBranchName = "jets";
899   fGammaConeSize = 0.4;
900   fUseBackgroundSubtractionGamma = kFALSE;
901   fSaveGJTree = kTRUE;
902   fMostEnergetic = kFALSE;
903   fMostOpposite = kTRUE;
904   fUseHistogramJetBkg = kTRUE;
905   fUseHistogramTracks = kTRUE;
906   fUseHistogramJetTracks = kTRUE;
907   fMCStudies = kFALSE;
908 }
909
910 //__________________________________________________________________
911 Int_t  AliAnaParticleJetFinderCorrelation::SelectJet(AliAODPWG4Particle * particle, TClonesArray *aodRecJets)
912 {
913   //Input for jets is TClonesArray *aodRecJets
914   //Returns the index of the jet that is opposite to the particle
915   //printf(" SelectJet ");
916   
917   Double_t particlePt=particle->Pt();
918   if(fUseBackgroundSubtractionGamma) {
919     Int_t clusterIDtmp = particle->GetCaloLabel(0) ;
920     Int_t nCells=0;
921     AliVCluster *cluster=0;
922     if(!(clusterIDtmp<0) ){
923       Int_t iclustmp = -1;
924       TObjArray* clusters = GetEMCALClusters();
925       cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
926       nCells = cluster->GetNCells();
927     }
928     particlePt-=(fGamRho*nCells);
929   }
930   if(particlePt<=0) {
931     //printf("Particle with negative  or 0 pt\n");
932     return -2;
933   }
934   
935   Int_t njets = aodRecJets->GetEntriesFast();
936   AliAODJet * jet = 0 ;
937   Int_t index = -1;
938   Double_t dphiprev= 10000;
939   Double_t particlePhi=particle->Phi();
940   Double_t deltaPhi=-10000.;// in the range (0; 2*pi)
941   Double_t jetPt=0.;
942   
943   for(Int_t ijet = 0; ijet < njets ; ijet++){
944     jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
945     if(!jet)
946     {
947       AliInfo("Jet not in container");
948       continue;
949     }
950     fhCuts2->Fill(2.,1.);
951     jetPt=jet->Pt();
952     if(jetPt<fJetMinPt) continue;
953     fhCuts2->Fill(3.,1.);
954     if(fBackgroundJetFromReader ){
955       jetPt-= (fJetRho * jet->EffectiveAreaCharged() );
956     }
957     if(jetPt<0.) continue;
958     //put jet eta requirement here |eta_jet|<0.9-jet_cone_size
959     fhCuts2->Fill(4.,1.);
960     if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
961     fhCuts2->Fill(5.,1.);
962     if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
963     fhCuts2->Fill(6.,1.);
964     //printf("jet found\n");
965     Double_t deltaPhi0pi  = TMath::Abs(particle->Phi()-jet->Phi());
966     Double_t ratio = jetPt/particlePt;
967     
968     deltaPhi = particlePhi - jet->Phi() ;
969     if ( deltaPhi0pi > TMath::Pi() ) deltaPhi0pi = 2. * TMath::Pi() - deltaPhi0pi ;
970     if(deltaPhi<0) deltaPhi +=(TMath::Pi()*2.);
971     
972     fhDeltaPtBefore ->Fill(particlePt, particlePt - jetPt);
973     fhDeltaPhiBefore->Fill(particlePt, deltaPhi);
974     fhDeltaEtaBefore->Fill(particlePt, particle->Eta() - jet->Eta());
975     fhPtRatioBefore ->Fill(particlePt, jetPt/particlePt);
976     fhPtBefore      ->Fill(particlePt, jetPt);
977     
978     fhDeltaPhi0PiCorrectBefore->Fill(particlePt, deltaPhi0pi);//good
979     
980     if(GetDebug() > 5)
981       printf("AliAnaParticleJetFinderCorrelation::SelectJet() - Jet %d, Ratio pT %2.3f, Delta phi %2.3f\n",ijet,ratio,deltaPhi);
982     
983     if((deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) &&
984        (ratio > fRatioMinCut) && (ratio < fRatioMaxCut)  &&
985        (TMath::Abs(deltaPhi-TMath::Pi()) < TMath::Abs(dphiprev-TMath::Pi())  )){//In case there is more than one jet select the most opposite.
986       dphiprev = deltaPhi;
987       index = ijet ;
988     }//Selection cuts
989   }//AOD Jet loop
990   
991   return index ;
992   
993 }
994
995 //__________________________________________________________________
996 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD()
997 {
998   //Particle-Jet Correlation Analysis, fill AODs
999   //  printf("I use MakeAnalysisFillAOD\n");
1000   //Get the event, check if there are AODs available, if not, skip this analysis
1001   AliAODEvent * event = NULL;
1002   
1003   //  cout<<"GetReader()->GetOutputEvent() "<<GetReader()->GetOutputEvent()<<endl;
1004   //  cout<<"GetReader()->GetDataType() "<<GetReader()->GetDataType() <<endl;
1005   //  cout<<"AliCaloTrackReader::kAOD "<<AliCaloTrackReader::kAOD<<endl;
1006   //  cout<<"GetReader()->GetInputEvent() "<<GetReader()->GetInputEvent()<<endl;
1007   
1008   if(GetReader()->GetOutputEvent())
1009   {
1010     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1011   }
1012   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1013   {
1014     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1015   }
1016   else
1017   {
1018     if(GetDebug() > 3) AliInfo("There are no jets available for this analysis");
1019     return;
1020   }
1021   
1022   if(!GetInputAODBranch() || !event)
1023   {
1024     AliFatal(Form("No input particles in AOD with name branch < %s > \n",
1025                   GetInputAODName().Data()));
1026     return; // Trick coverity
1027   }
1028   
1029   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
1030   {
1031     AliFatal(Form("Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s>",
1032                   GetInputAODBranch()->GetClass()->GetName()));
1033     return; // Trick coverity
1034   }
1035   
1036   //
1037   // non-standard jets
1038   //
1039   Int_t nJets=-1;
1040   TClonesArray *aodRecJets = 0;
1041   //if(IsNonStandardJetFromReader()){//jet branch from reader
1042   if(GetDebug() > 3) AliInfo(Form("GetNonStandardJets function (from reader) is called"));
1043   aodRecJets = GetNonStandardJets();
1044   if(GetDebug() > 3) AliInfo(Form("aodRecJets %p",aodRecJets));
1045   if(aodRecJets==0x0)
1046     {
1047       if(GetDebug() > 3) event->Print();
1048       AliFatal("List of jets is null");
1049       return;
1050     }
1051   nJets=aodRecJets->GetEntries();
1052   if(GetDebug() > 3) printf("nJets %d\n",nJets);
1053   //}
1054   //end of new part
1055   
1056   if(nJets==0) {
1057     //printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1058     return;
1059   }
1060   
1061   //
1062   //Background jets
1063   //
1064   AliAODJetEventBackground* aodBkgJets = 0;
1065   if(IsBackgroundJetFromReader()){//jet branch from reader
1066     if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
1067     aodBkgJets = GetBackgroundJets();
1068     if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
1069     if(aodBkgJets==0x0){
1070       if(GetDebug() > 3) event->Print();
1071       AliFatal("No jet background found\n");
1072       return; // Trick coverity        
1073     }
1074     if(GetDebug() > 3) aodBkgJets->Print("c");
1075   }
1076   
1077   Double_t rhoEvent=0.;
1078   if(aodBkgJets && IsBackgroundJetFromReader() ) {
1079     rhoEvent = aodBkgJets->GetBackground(2);//hardcoded
1080   }
1081   fJetRho = rhoEvent;
1082   
1083   
1084   Int_t ntrig =  GetInputAODBranch()->GetEntriesFast() ;
1085   if(GetDebug() > 3){
1086     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Begin jet finder  correlation analysis, fill AODs \n");
1087     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", ntrig);
1088     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In standard jet branch aod entries %d\n", event->GetNJets());
1089     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In non standard jet branch aod entries %d\n", nJets);
1090   }
1091   
1092   //if(nJets==0)   return;//to speed up
1093   //  cout<<"ntrig po return "<<ntrig<<endl;
1094   
1095   //
1096   //calculate average cell energy without most energetic photon
1097   //
1098   Double_t medianPhotonRho=0.;
1099   TObjArray* clusters = GetEMCALClusters();
1100   Int_t clusterIDtmp;
1101   Int_t iclustmp = -1;
1102   AliVCluster *cluster=0;
1103   
1104   if(IsBackgroundSubtractionGamma()){
1105     //
1106     // Find most energetic photon without background subtraction
1107     //
1108     Double_t maxPt=0.;
1109     Int_t maxIndex=-1;
1110     AliAODPWG4ParticleCorrelation* particlecorr =0;
1111     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1112       particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1113       if(particlecorr->Pt() > maxPt) {
1114         maxPt = particlecorr->Pt();
1115         maxIndex = iaod;
1116       }
1117     }
1118     
1119     //
1120     // calculate background energy per cell
1121     //
1122     Int_t numberOfcells=0;
1123     if(ntrig>1){
1124       Double_t *photonRhoArr=new Double_t[ntrig-1];
1125       Int_t photonRhoArrayIndex=0;
1126       for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1127         particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1128         if(iaod==maxIndex) continue;
1129         clusterIDtmp = particlecorr->GetCaloLabel(0) ;
1130         if(clusterIDtmp < 0) continue;
1131         cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1132         photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1133         numberOfcells+=cluster->GetNCells();
1134         photonRhoArrayIndex++;
1135       }
1136       if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1137       delete [] photonRhoArr;
1138     }
1139   }//end of if background calculation for gamma
1140   fGamRho = medianPhotonRho;
1141   
1142   
1143   //
1144   //take most energetic photon and most energetic jet and combine
1145   //
1146   if(fMostEnergetic){
1147     //
1148     // most energetic photon with background subtraction
1149     //
1150     Double_t mostEnePhotonPt=0.;
1151     Int_t indexMostEnePhoton=-1;
1152     AliAODPWG4ParticleCorrelation* particle =0;
1153     Double_t ptCorrect=0.;
1154     Int_t nCells=0;
1155     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1156       particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1157       clusterIDtmp = particle->GetCaloLabel(0) ;
1158       if(!(clusterIDtmp<0)){
1159         cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1160         nCells = cluster->GetNCells();
1161       }
1162       ptCorrect = particle->Pt() - medianPhotonRho * nCells;
1163       if( ptCorrect > mostEnePhotonPt ){
1164         mostEnePhotonPt = ptCorrect;
1165         indexMostEnePhoton = iaod ;
1166       }
1167     }
1168     //    printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Photon with index %d selected \n",indexMostEnePhoton);
1169     //
1170     // most energetic jet with background subtraction
1171     //
1172     Double_t mostEneJetPt=0.;
1173     Int_t indexMostEneJet=-1;
1174     AliAODJet * jet = 0 ;
1175     //Double_t ptCorrect=0.;
1176     for(Int_t ijet = 0; ijet < nJets ; ijet++){
1177       jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1178       if(!jet) continue;
1179       if(jet->Pt()<fJetMinPt) continue;
1180       if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
1181       if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
1182       ptCorrect = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();
1183       if(ptCorrect > mostEneJetPt){
1184         mostEneJetPt = ptCorrect;
1185         indexMostEneJet = ijet;
1186       }
1187     }
1188     //    printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",indexMostEneJet);
1189     
1190     //
1191     // assign jet to photon
1192     //
1193     if(indexMostEneJet>=0 && indexMostEnePhoton>=0){
1194       particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(indexMostEnePhoton));
1195       jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(indexMostEneJet));
1196       if(jet)particle->SetRefJet(jet);
1197     }
1198   }//end of take most energetic photon and most ene. jet after bkg subtraction
1199   
1200   if(fMostOpposite){
1201     //Bool_t isJetFound=kFALSE;//new
1202     //Loop on stored AOD particles, trigger
1203     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1204       AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1205       
1206       //Correlate with jets
1207       Int_t ijet = SelectJet(particle,aodRecJets);//input for jets is TClonesArray
1208       if(ijet > -1){
1209         //isJetFound=kTRUE;
1210         if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",ijet);
1211         AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1212         if(jet)particle->SetRefJet(jet);
1213         //printf("Most opposite found\n");
1214       }
1215     } // input aod loop
1216     //  if(GetReader()->WriteDeltaAODToFile() && isJetFound) WriteJetsToOutputBranch(aodRecJets);
1217   }//end of take most opposite photon and jet after bkg subtraction
1218   
1219   
1220   if(GetDebug() > 1 ) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
1221
1222
1223 //__________________________________________________________________
1224 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms()
1225 {
1226   //Particle-Jet Correlation Analysis, fill histograms
1227   if(GetDebug() > 3 ) {
1228     printf("I use MakeAnalysisFillHistograms\n");
1229     printf("ntracks before iso %d\n",GetCTSTracks()->GetEntriesFast() );
1230   }
1231
1232   //Get the event, check if there are AODs available, if not, skip this analysis
1233   AliAODEvent * event = NULL;
1234   
1235   //printf("GetReader()->GetOutputEvent() %d\n",GetReader()->GetOutputEvent() );
1236   //printf("GetReader()->GetDataType() %d\n",GetReader()->GetDataType() );
1237   //printf("AliCaloTrackReader::kAOD %d\n",AliCaloTrackReader::kAOD );
1238   //printf("GetReader()->GetInputEvent() %d\n",GetReader()->GetInputEvent() );
1239   
1240   if(GetReader()->GetOutputEvent())
1241   {
1242     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1243   }
1244   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1245   {
1246     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1247   }
1248   else {
1249     if(GetDebug() > 3) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - There are no jets available for this analysis\n");
1250     return;
1251   }
1252   
1253   if(!GetInputAODBranch() || !event){
1254
1255     AliFatal(Form("No input particles in AOD with name branch < %s > \n",
1256                   GetInputAODName().Data()));
1257     return; // Trick coverity        
1258   }
1259   
1260   Int_t nJets=-1;
1261   TClonesArray *aodRecJets = 0;
1262   //if(IsNonStandardJetFromReader()){//branch read in reader from reader
1263   if (GetDebug() > 3) printf("GetNonStandardJets function (from reader) is called\n");
1264   aodRecJets = GetNonStandardJets();
1265   if(aodRecJets==0x0){
1266     if(GetDebug() > 3) event->Print();
1267     AliFatal("Jets container not found\n");
1268     return; // trick coverity
1269   }
1270   nJets=aodRecJets->GetEntries();
1271   //}
1272   if(nJets==0) {
1273     //    printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1274     GetReader()->FillInputNonStandardJets();
1275     aodRecJets = GetNonStandardJets();
1276     if(aodRecJets) nJets=aodRecJets->GetEntries();
1277     //    printf("nJets = %d\n",nJets);
1278     return;
1279   }
1280   
1281   //
1282   //Background jets
1283   //
1284   AliAODJetEventBackground* aodBkgJets = 0;
1285   if(IsBackgroundJetFromReader()){//jet branch from reader
1286     if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
1287     aodBkgJets = GetBackgroundJets();
1288     if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
1289     if(aodBkgJets==0x0){
1290       if(GetDebug() > 3) event->Print();
1291       AliFatal("No jet background container found");
1292       return; // trick coverity  
1293     }
1294     if(GetDebug() > 3) aodBkgJets->Print("c");
1295   }
1296   
1297   
1298   //
1299   // only background jets informations
1300   //
1301   //  Float_t pTback = 0;
1302   Double_t rhoEvent=0.;
1303   if(aodBkgJets) {
1304     if(IsBackgroundJetFromReader() ) rhoEvent = aodBkgJets->GetBackground(2);
1305     if(IsHistogramJetBkg()) {
1306       fhJetRhoVsCentrality->Fill(GetEventCentrality(),rhoEvent);
1307       for(Int_t i=0;i<4;i++){
1308         fhBkgJetBackground[i]->Fill(aodBkgJets->GetBackground(i));
1309         fhBkgJetSigma[i]->Fill(aodBkgJets->GetSigma(i));
1310         fhBkgJetArea[i]->Fill(aodBkgJets->GetMeanarea(i));
1311       }
1312     }//end of if fill HistogramJetBkg
1313   }//end if aodBkgJets exists
1314   
1315   //
1316   //only track information
1317   //
1318   Int_t nCTSTracks = GetCTSTracks()->GetEntriesFast();
1319   AliAODTrack *aodtrack;
1320   Int_t itrack = 0;
1321   if(IsHistogramTracks()) {
1322     Double_t sumTrackPt=0;
1323     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1324       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1325       if(!aodtrack) continue;
1326       fhTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1327       sumTrackPt+=aodtrack->Pt();
1328     }
1329     if(nCTSTracks)
1330       fhTrackAveTrackPt->Fill(sumTrackPt/nCTSTracks);
1331   }//end if IsHistogramTracks
1332   
1333   //
1334   //only jet informations
1335   //
1336   AliAODJet * jettmp = 0 ;
1337   Double_t jetPttmp = 0.;
1338   Double_t var1tmp = 0.;
1339   Double_t var2tmp = 0.;
1340   //  fhJetNjet->Fill(nJets);
1341   Double_t ptMostEne=0;
1342   //  Int_t indexMostEne=-1;
1343   Int_t nJetsOverThreshold[10]={nJets,0,0,0,0,0,0,0,0,0};
1344   Int_t iCounter=0;
1345   Double_t sumJetTrackPt=0.;
1346   Int_t sumNJetTrack=0;
1347   Int_t nTracksInJet=0;
1348   Int_t itrk=0;
1349   for(Int_t ijet = 0; ijet < nJets ; ijet++){
1350     jettmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1351     if(!jettmp) continue;
1352     fhJetPtBefore->Fill(jettmp->Pt());
1353     jetPttmp  = jettmp->Pt() - rhoEvent * jettmp->EffectiveAreaCharged();//<<---changed here
1354     
1355     //calculate number of tracks above 1,2,3,4,5 GeV in jet
1356     AliVTrack* jettrack = 0x0 ;
1357     
1358     Int_t nTrackThrGeV[5]={0,0,0,0,0};
1359     nTracksInJet=(jettmp->GetRefTracks())->GetEntriesFast();
1360     fhJetNparticlesInJet->Fill(nTracksInJet);
1361     if(nTracksInJet==0) continue;
1362     sumNJetTrack+=nTracksInJet;
1363     for(itrack=0;itrack<nTracksInJet;itrack++){
1364       jettrack=(AliVTrack *) ((jettmp->GetRefTracks())->At(itrack));
1365       if(!jettrack) continue;
1366       
1367       fhJetDeltaEtaDeltaPhi->Fill(jettmp->Eta()-jettrack->Eta(),jettmp->Phi()-jettrack->Phi());
1368       sumJetTrackPt+=jettrack->Pt();
1369       if(IsHistogramJetTracks()){
1370         if(jettrack->Pt()>1.) nTrackThrGeV[0]++;
1371         if(jettrack->Pt()>2.) nTrackThrGeV[1]++;
1372         if(jettrack->Pt()>3.) nTrackThrGeV[2]++;
1373         if(jettrack->Pt()>4.) nTrackThrGeV[3]++;
1374         if(jettrack->Pt()>5.) nTrackThrGeV[4]++;
1375       }//end of if IsHistogramJetTracks
1376     }//end of jet track loop
1377     //printf("jetPt %f ntrks %d ntrks>1,2,3,4,5GeV in jet %d, %d, %d, %d, %d\n",jetPttmp,nTracksInJet,nTrackThrGeV[0],nTrackThrGeV[1],nTrackThrGeV[2],nTrackThrGeV[3],nTrackThrGeV[4]);
1378     
1379     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1380       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1381       if(aodtrack) fhJetDeltaEtaDeltaPhiAllTracks->Fill(jettmp->Eta()-aodtrack->Eta(),jettmp->Phi()-aodtrack->Phi());
1382     }
1383     
1384     
1385     if(IsHistogramJetTracks()){
1386       fhJetNtracksInJetAboveThr[0]->Fill(jetPttmp,nTracksInJet);//all jets
1387       
1388       for(itrk=0;itrk<5;itrk++) {
1389         fhJetNtracksInJetAboveThr[itrk+1]->Fill(jetPttmp,nTrackThrGeV[itrk]);//all jets
1390         fhJetRatioNTrkAboveToNTrk[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);//all jets
1391       }
1392       if(ijet==0){//most ene jet
1393         for(itrk=0;itrk<5;itrk++)
1394           fhJetNtrackRatioMostEne[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1395       }
1396       if(jetPttmp>5){//jet with pt>5GeV/c
1397         for(itrk=0;itrk<5;itrk++)
1398           fhJetNtrackRatioJet5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1399       }
1400       if(nTrackThrGeV[4]>0){//jet with leading particle pt>5GeV/c
1401         for(itrk=0;itrk<5;itrk++)
1402           fhJetNtrackRatioLead5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1403       }
1404     }//end of if IsHistogramJetTracks
1405     
1406     
1407     Double_t effectiveChargedBgEnergy=(IsBackgroundJetFromReader()?rhoEvent * jettmp->EffectiveAreaCharged():jettmp->ChargedBgEnergy());
1408     
1409     
1410     fhJetChBkgEnergyVsArea->Fill(effectiveChargedBgEnergy,jettmp->EffectiveAreaCharged());
1411     //if(jettmp->EffectiveAreaCharged()>0)
1412     fhJetRhoVsPt->Fill(jetPttmp,jettmp->ChargedBgEnergy()*jettmp->EffectiveAreaCharged());
1413     
1414     if(jetPttmp>ptMostEne) {
1415       ptMostEne = jetPttmp;
1416       //indexMostEne=ijet;
1417     }
1418     fhJetPt->Fill(jetPttmp);
1419     fhJetChBkgEnergyVsPt->Fill(jetPttmp,effectiveChargedBgEnergy);
1420     fhJetChAreaVsPt->Fill(jetPttmp,jettmp->EffectiveAreaCharged());
1421     if(GetDebug()>5) printf("ChargedBgEnergy %f EffectiveAreaCharged %f\n", jettmp->ChargedBgEnergy(),jettmp->EffectiveAreaCharged());
1422     for(iCounter=1;iCounter<10;iCounter++){
1423       if(jetPttmp>iCounter) nJetsOverThreshold[iCounter]++;
1424     }
1425     
1426     var1tmp  = jettmp->Phi();
1427     var2tmp  = jettmp->Eta();
1428     fhJetPhi->Fill(var1tmp);
1429     fhJetEta->Fill(var2tmp);
1430     fhJetPhiVsEta->Fill(var1tmp,var2tmp);
1431     fhJetEtaVsPt->Fill(jetPttmp,var2tmp);
1432     fhJetEtaVsNpartInJet->Fill(var2tmp,nTracksInJet);
1433     if(jetPttmp>0)
1434       fhJetEtaVsNpartInJetBkg->Fill(var2tmp,nTracksInJet);
1435     
1436   }//end of jet loop
1437   if(IsHistogramJetTracks()){
1438     if(sumNJetTrack>0){
1439       //printf("average track pt %f\n",sumJetTrackPt/sumNJetTrack);
1440       fhJetAveTrackPt->Fill(sumJetTrackPt/sumNJetTrack);
1441     }
1442   }//end of if IsHistogramJetTracks
1443   
1444   
1445   fhJetPtMostEne->Fill(ptMostEne);
1446   for(iCounter=0;iCounter<10;iCounter++){
1447     fhJetNjetOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1448     nJetsOverThreshold[iCounter]=0;
1449   }
1450   
1451   //end of only jet part
1452   
1453   //
1454   // Photons
1455   //
1456   Int_t ntrig   =  GetInputAODBranch()->GetEntriesFast() ;
1457   if(GetDebug() > 1){
1458     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Begin jet finder  correlation analysis, fill histograms \n");
1459     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", ntrig);
1460     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In jet output branch aod entries %d\n", event->GetNJets());
1461   }
1462   fhNjetsNgammas->Fill(nJets,ntrig);
1463   
1464   //if(nJets==0)   return;//to speed up
1465   //printf("ntrig %d, njets %d\n",ntrig,nJets);
1466   
1467   //
1468   //Fill only temporary photon histograms
1469   //
1470   Double_t maxPt=0.;
1471   Int_t maxIndex=-1;
1472   Double_t tmpPt=0.;
1473   Double_t sumPt=0.;
1474   nJetsOverThreshold[0]=ntrig;
1475   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1476     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1477     tmpPt = particlecorr->Pt();
1478     if(tmpPt>maxPt) {
1479       maxPt = tmpPt;
1480       maxIndex = iaod;
1481     }
1482     for(iCounter=1;iCounter<10;iCounter++){
1483       if(tmpPt>iCounter) nJetsOverThreshold[iCounter]++;
1484     }
1485     sumPt+=tmpPt;
1486   }
1487   for(iCounter=0;iCounter<10;iCounter++){
1488     fhPhotonNgammaOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1489     //    nJetsOverThreshold[iCounter]=0;
1490   }
1491   
1492   fGamAvEne=0;
1493   TObjArray* clusters = GetEMCALClusters();
1494   //printf("calculate median bkg energy for photons ");
1495   Double_t medianPhotonRho=0.;
1496   Int_t clusterID;
1497   Int_t iclustmp = -1;
1498   Int_t numberOfcells=0;
1499   AliVCluster *cluster = 0;
1500   if(ntrig>1){
1501     Double_t *photonRhoArr=new Double_t[ntrig-1];
1502     fhPhotonPtMostEne->Fill(maxPt);
1503     //    fhPhotonIndexMostEne->Fill(indexMaxPt);
1504     fhPhotonAverageEnergy->Fill(sumPt/ntrig);
1505     fhPhotonRatioAveEneToMostEne->Fill(sumPt/(ntrig*maxPt));
1506     fhPhotonAverageEnergyMinus1->Fill((sumPt-maxPt)/(ntrig-1));
1507     fGamAvEne=(sumPt-maxPt)/(ntrig-1);
1508     fhPhotonRatioAveEneMinus1ToMostEne->Fill((sumPt-maxPt)/((ntrig-1)*maxPt));
1509     
1510     Int_t counterGamma=0;
1511     Int_t counterGammaMinus1=0;
1512     
1513     Int_t photonRhoArrayIndex=0;
1514     //TObjArray* clusterstmp = GetEMCALClusters();
1515     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1516       AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1517       if( particlecorr->Pt() > sumPt/ntrig ) counterGamma++;
1518       if( particlecorr->Pt() > (sumPt-maxPt)/(ntrig-1) ) counterGammaMinus1++;
1519       
1520       if(iaod==maxIndex) continue;
1521       clusterID = particlecorr->GetCaloLabel(0) ;
1522       if(clusterID < 0) continue;
1523       cluster = FindCluster(clusters,clusterID,iclustmp);
1524       photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1525       numberOfcells+=cluster->GetNCells();
1526       photonRhoArrayIndex++;
1527     }
1528     if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1529     delete [] photonRhoArr;
1530     fhPhotonNgammaMoreAverageToNgamma->Fill((Double_t)counterGamma / (Double_t)ntrig);
1531     fhPhotonNgammaMoreAverageMinus1ToNgamma->Fill((Double_t)counterGammaMinus1 / (Double_t)ntrig);
1532   }
1533   //printf("median = %f\n",medianPhotonRho);
1534   fhPhotonBkgRhoVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),medianPhotonRho);
1535   fhPhotonBkgRhoVsNclusters->Fill(ntrig,medianPhotonRho);
1536   fhPhotonBkgRhoVsCentrality->Fill(GetEventCentrality(),medianPhotonRho);
1537   fhPhotonBkgRhoVsNcells->Fill(numberOfcells,medianPhotonRho);
1538   
1539   
1540   AliVCluster *cluster2 = 0;
1541   Double_t photon2Corrected=0;
1542   Double_t sumPtTmp=0.;
1543   Double_t sumPtCorrectTmp=0.;
1544   AliVTrack* trackTmp = 0x0 ;
1545   TVector3 p3Tmp;
1546   
1547   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1548     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1549     clusterID = particlecorr->GetCaloLabel(0) ;
1550     if(clusterID < 0) continue;
1551     cluster = FindCluster(clusters,clusterID,iclustmp);
1552     fhPhotonPt->Fill(particlecorr->Pt());
1553     fhPhotonPtCorrected->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1554     fhPhotonPtDiff->Fill(cluster->GetNCells() * medianPhotonRho);
1555     fhPhotonPtDiffVsCentrality->Fill(GetEventCentrality(),cluster->GetNCells() * medianPhotonRho);
1556     fhPhotonPtDiffVsNcells->Fill(numberOfcells,cluster->GetNCells() * medianPhotonRho);
1557     fhPhotonPtDiffVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),cluster->GetNCells() * medianPhotonRho);
1558     fhPhotonPtDiffVsNclusters->Fill(ntrig,cluster->GetNCells() * medianPhotonRho);
1559     
1560     fhPhotonPtCorrectedZoom->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1561     
1562     //test: sum_pt in the cone 0.3 for each photon
1563     //should be: random fake gamma from MB
1564     //is: each gamma for EMCEGA
1565     sumPtTmp=0.;
1566     sumPtCorrectTmp=0.;
1567     
1568     for(Int_t iaod2 = 0; iaod2 < ntrig ; iaod2++){
1569       if(iaod==iaod2) continue;
1570       AliAODPWG4ParticleCorrelation* particlecorr2 =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod2));
1571       clusterID = particlecorr2->GetCaloLabel(0) ;
1572       if(clusterID < 0) continue;
1573       cluster2 = FindCluster(clusters,clusterID,iclustmp);
1574       photon2Corrected = particlecorr2->Pt() - cluster2->GetNCells() * medianPhotonRho;
1575       
1576       //if(Pt()<0.5) continue; //<<hardcoded here //FIXME
1577       if( TMath::Sqrt((particlecorr->Eta()-particlecorr2->Eta())*(particlecorr->Eta()-particlecorr2->Eta()) +
1578                       (particlecorr->Phi()-particlecorr2->Phi())*(particlecorr->Phi()-particlecorr2->Phi()) )<fGammaConeSize ){//if(/*cone is correct*/){
1579         sumPtTmp+= particlecorr2->Pt();
1580         sumPtCorrectTmp+=photon2Corrected;
1581       }
1582     }
1583     fhPhotonSumPtInCone->Fill(sumPtTmp);
1584     fhPhotonSumPtCorrectInCone->Fill(sumPtCorrectTmp);
1585     
1586     //test: sum_pt in the cone 0.3 for each track
1587     //should be: random fake gamma from MB
1588     //is: each gamma for EMCEGA
1589     sumPtTmp=0.;
1590     for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++){
1591       trackTmp = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1592       p3Tmp.SetXYZ(trackTmp->Px(),trackTmp->Py(),trackTmp->Pz());
1593       if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1594                       (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1595         sumPtTmp+=p3Tmp.Pt();
1596       }
1597     }//end of loop over tracks
1598     fhPhotonSumPtChargedInCone->Fill(sumPtTmp);
1599   }
1600   
1601   //End of Fill temporary photon histograms
1602   
1603   //
1604   // Apply background subtraction for photons
1605   //
1606   fGamRho = medianPhotonRho;
1607   if(!IsBackgroundSubtractionGamma()) medianPhotonRho=0;
1608   
1609   
1610   //
1611   //Get vertex for cluster momentum calculation <<----new here
1612   //
1613   Double_t vertex[] = {0,0,0} ; //vertex ;
1614   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1615     GetReader()->GetVertex(vertex);
1616   fZvertex = vertex[2];
1617   
1618   //
1619   //Loop on stored AOD particles, trigger
1620   //
1621   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1622     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1623     fhCuts->Fill(0);
1624     fhCuts2->Fill(0.,(Double_t)nJets);
1625     if(GetDebug() > 5) printf("OnlyIsolated %d  !particlecorr->IsIsolated() %d \n",OnlyIsolated(), !particlecorr->IsIsolated());
1626     
1627     if(OnlyIsolated() && !particlecorr->IsIsolated()) continue;
1628     fhCuts->Fill(1);
1629     fhCuts2->Fill(1.,nJets);
1630     
1631     if(nJets>0) {
1632       fhCuts->Fill(2);
1633     }
1634     
1635     //Recover the jet correlated, found previously.
1636     AliAODJet   * jet = particlecorr->GetJet();
1637     //If correlation not made before, do it now.
1638     if(fMakeCorrelationInHistoMaker){
1639       //Correlate with jets
1640       Int_t ijet = SelectJet(particlecorr,aodRecJets);//input for jets is TClonesArray
1641       if(ijet > -1){
1642         if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Jet with index %d selected \n",ijet);
1643         //jet = event->GetJet(ijet);
1644         jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1645         
1646        if(jet) particlecorr->SetRefJet(jet);
1647         
1648       }
1649     }
1650     
1651     if (!jet) continue ;
1652     fhCuts->Fill(3);
1653     fhCuts2->Fill(7.,1.);
1654
1655     //check MC genereted information
1656     if(fMCStudies) FindMCgenInfo();
1657
1658     //
1659     //Fill Correlation Histograms
1660     //
1661     clusterID = particlecorr->GetCaloLabel(0) ;
1662     if(!(clusterID<0)){
1663       cluster = FindCluster(clusters,clusterID,iclustmp);
1664       //fill tree variables
1665       fGamNcells = cluster->GetNCells();
1666     }
1667     Double_t ptTrig = particlecorr->Pt() - medianPhotonRho * fGamNcells;//<<---changed here
1668     Double_t ptJet = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();//<<---changed here
1669     Double_t phiTrig = particlecorr->Phi();
1670     Double_t phiJet = jet->Phi();
1671     Double_t etaTrig = particlecorr->Eta();
1672     Double_t etaJet = jet->Eta();
1673     Double_t deltaPhi=phiTrig-phiJet;
1674     if(deltaPhi<0)deltaPhi+=(TMath::Pi()*2.);
1675     //printf("pT trigger %2.3f, pT jet %2.3f, Delta phi %2.3f, Delta eta %2.3f, Delta pT %2.3f, ratio %2.3f \n",
1676     //  ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
1677     fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet);
1678     //    fhDeltaPhi->Fill(ptTrig, phiTrig-phiJet);//need to be shifted by 2pi
1679     
1680     fhDeltaPhiCorrect->Fill(ptTrig, deltaPhi);//correct
1681     
1682     Double_t deltaPhiCorrect = TMath::Abs( particlecorr->Phi() - jet->Phi() );
1683     if ( deltaPhiCorrect > TMath::Pi() ) deltaPhiCorrect = 2. * TMath::Pi() - deltaPhiCorrect ;
1684     fhDeltaPhi0PiCorrect->Fill(ptTrig, deltaPhiCorrect);
1685     
1686     fhDeltaEta->Fill(ptTrig, etaTrig-etaJet);
1687     fhPtRatio ->Fill(ptTrig, ptJet/ptTrig);
1688     fhPt      ->Fill(ptTrig, ptJet);
1689     
1690     fhSelectedJetPhiVsEta->Fill(phiJet,etaJet);
1691     fhSelectedJetChBkgEnergyVsPtJet->Fill(ptJet,(IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy()) );
1692     fhSelectedJetChAreaVsPtJet->Fill(ptJet,jet->EffectiveAreaCharged());
1693     fhSelectedJetNjet->Fill(nJets);
1694     fhSelectedNtracks->Fill(GetCTSTracks()->GetEntriesFast());//to be checked
1695     fhSelectedPhotonNLMVsPt->Fill(ptTrig,particlecorr->GetNLM());
1696     
1697     
1698     if(clusterID < 0 ){
1699       fhSelectedPhotonLambda0VsPt->Fill(ptTrig,-1);
1700       //fill tree variables
1701       fGamLambda0  = -1;
1702       fGamTime = -1;
1703       fGamNcells = 0;
1704       fGamSumPtNeu=0;
1705     }
1706     else{
1707       //Int_t iclus = -1;
1708       //      TObjArray* clusters = GetEMCALClusters();
1709       //cluster = FindCluster(clusters,clusterID,iclustmp);
1710       Double_t lambda0=cluster->GetM02();
1711       fhSelectedPhotonLambda0VsPt->Fill(ptTrig,lambda0);
1712       //fill tree variables
1713       fGamLambda0  = lambda0;
1714       fGamTime = cluster->GetTOF();
1715       //fGamNcells = cluster->GetNCells();
1716       
1717       fGamSumPtNeu=0;
1718       fGamNclusters=0;
1719       //TVector3 p3Tmp;
1720       //Double_t scalarProduct=0;
1721       //Double_t vectorLength=particlecorr->P();
1722       for(Int_t icalo=0; icalo <clusters->GetEntriesFast(); icalo++){
1723         AliVCluster* calo = (AliVCluster *) clusters->At(icalo);
1724         if(clusterID==calo->GetID()) continue;//the same cluster as trigger
1725         calo->GetMomentum(fMomentum,vertex) ;//Assume that come from vertex in straight line
1726         //printf("min pt %f\n",GetMinPt());
1727         if(fMomentum.Pt()<GetMinPt()) continue; //<<hardcoded here //FIXME 0.5 check if correct
1728         p3Tmp.SetXYZ(fMomentum.Px(),fMomentum.Py(),fMomentum.Pz());
1729         //calculate sum pt in the cone
1730         if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1731                         (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1732           //scalarProduct = particlecorr->Px()*fMomentum.Px() + particlecorr->Py()*fMomentum.Py() + particlecorr->Pz()*fMomentum.Pz();
1733           //scalarProduct/=fMomentum.P();
1734           //scalarProduct/=vectorLength;
1735           //if(scalarProduct>TMath::Cos(0.3)) {//FIXME photon radius
1736           fGamSumPtNeu+=fMomentum.Pt();
1737           fGamNclusters++;
1738         }
1739       }
1740     }
1741     
1742     //sum pt of charged tracks in the gamma isolation cone
1743     //starts here
1744     fGamSumPtCh=0;
1745     fGamNtracks=0;
1746     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1747       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1748       if(!aodtrack) continue;
1749       fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());//fill histogram here
1750       //      if(aodtrack->Pt()<0.15) continue;//hardcoded
1751       if(aodtrack->Pt()<fPtThresholdInCone) continue;
1752       if(!aodtrack->IsHybridGlobalConstrainedGlobal()) continue;
1753       if(TMath::Sqrt((particlecorr->Phi() - aodtrack->Phi())*(particlecorr->Phi() - aodtrack->Phi()) +
1754                      (particlecorr->Eta() - aodtrack->Eta())*(particlecorr->Eta() - aodtrack->Eta()) ) <fGammaConeSize ) {
1755         fGamSumPtCh+=aodtrack->Pt();
1756         fGamNtracks++;
1757       }
1758     }
1759     //ends here
1760     
1761     //    for(Int_t itrack = 0; itrack < nCTSTracks ; itrack++){
1762     //      aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1763     //      fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1764     //    }
1765     
1766     //
1767     // Background Fragmentation function
1768     //
1769     TVector3 gammaVector,jetVector;
1770     gammaVector.SetXYZ(particlecorr->Px(),particlecorr->Py(),particlecorr->Pz());
1771     jetVector.SetXYZ(jet->Px(),jet->Py(),jet->Pz());
1772     CalculateBkg(gammaVector,jetVector,vertex,1);//jet perp
1773     CalculateBkg(gammaVector,jetVector,vertex,2);//RC
1774     CalculateBkg(gammaVector,jetVector,vertex,3);//mid point
1775     CalculateBkg(gammaVector,jetVector,vertex,4);//gamma perp
1776     //CalculateBkg(gammaVector,jetVector,vertex,5);/test
1777     Double_t angleJetGam = gammaVector.Angle(jetVector);
1778     //printf("angleJetGam %f\n",angleJetGam*180/TMath::Pi());
1779     
1780     //
1781     // Fragmentation function
1782     //
1783     Float_t      rad = 0, pt = 0, eta = 0, phi = 0;
1784     Int_t        npartcone = 0;
1785     TVector3 p3;
1786     
1787     Int_t ntracks =  0;
1788     if(GetDebug()>3){
1789       printf ("fUseJetRefTracks %d\n",fUseJetRefTracks );
1790       printf ("jet->GetRefTracks() %p\n",jet->GetRefTracks());
1791       printf ("GetCTSTracks() %p\n",GetCTSTracks() );
1792     }
1793     
1794     if(!fUseJetRefTracks)
1795       ntracks =GetCTSTracks()->GetEntriesFast();
1796     else //If you want to use jet tracks from JETAN
1797       ntracks =  (jet->GetRefTracks())->GetEntriesFast();
1798     
1799     if(GetDebug()>3)    printf ("ntracks %d\n",ntracks);
1800     AliVTrack* track = 0x0 ;
1801     for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
1802       if(!fUseJetRefTracks)
1803         track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1804       else //If you want to use jet tracks from JETAN
1805         track = (AliVTrack *) ((jet->GetRefTracks())->At(ipr));
1806       
1807       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1808       pt    = p3.Pt();
1809       eta  = p3.Eta();
1810       phi  = p3.Phi() ;
1811       if(phi < 0) phi+=TMath::TwoPi();
1812       
1813       //Check if there is any particle inside cone with pt larger than  fPtThreshold
1814       rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
1815       if(rad < fConeSize  && pt > fPtThresholdInCone){
1816         //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
1817         npartcone++;
1818         fhFFz ->Fill(ptTrig, pt/ptTrig);
1819         fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt));
1820         fhFFpt->Fill(ptTrig, pt);
1821         
1822         //according to jet axis
1823         fhJetFFz ->Fill(ptJet, pt/ptJet);
1824         fhJetFFxi->Fill(ptJet, TMath::Log(ptJet/pt));
1825         fhJetFFpt->Fill(ptJet, pt);
1826         
1827         
1828         if(TMath::Cos(angleJetGam)<0 && ptJet!=0 && pt!=0 ){
1829           fhJetFFzCor ->Fill(ptJet, -pt*TMath::Cos(angleJetGam)/ptJet);
1830           fhJetFFxiCor->Fill(ptJet, TMath::Log(ptJet/(-pt*TMath::Cos(angleJetGam))));
1831         }
1832       }
1833     }//Tracks
1834     fhNTracksInCone->Fill(ptTrig, npartcone);
1835     //fill tree here for each photon-jet (isolation required usually)
1836     
1837     fGamPt      = ptTrig;
1838     //fGamLambda0  = ;//filled earlier
1839     fGamNLM      = particlecorr->GetNLM();
1840     //fGamSumPtCh  = ;//filled earlier
1841     //fGamTime     = particlecorr->GetTOF();//filled earlier
1842     //fGamNcells   = particlecorr->GetNCells();//filled earlier
1843     fGamEta      = etaTrig;
1844     fGamPhi      = phiTrig;
1845     //fGamSumPtNeu = ;//filled earlier
1846     //fGamNtracks  = ;//filled earlier
1847     //fGamNclusters= ;//filled earlier
1848     //fGamAvEne    = ;//filled earlier
1849     fJetPhi      = phiJet;
1850     fJetEta      = etaJet;
1851     fJetPt       = ptJet;
1852     fJetBkgChEne = (IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy());
1853     fJetArea     = jet->EffectiveAreaCharged();
1854     fJetNtracks  = (jet->GetRefTracks())->GetEntriesFast();
1855     fEventNumber = 0;
1856     fNtracks     = GetCTSTracks()->GetEntriesFast();
1857     fCentrality  = GetEventCentrality();
1858     fIso         = particlecorr->IsIsolated();
1859     
1860     Int_t nTrk1GeV=0;
1861     Int_t nTrk2GeV=0;
1862     for(itrack=0;itrack < fJetNtracks;itrack++){
1863       track = (AliVTrack *) ((jet->GetRefTracks())->At(itrack));
1864       if(track->Pt()>1.) nTrk1GeV++;
1865       if(track->Pt()>2.) nTrk2GeV++;
1866     }
1867     
1868     fJetNtracks1 = nTrk1GeV;
1869     fJetNtracks2 = nTrk2GeV;
1870     
1871     if(fSaveGJTree) fTreeGJ->Fill();
1872   }//AOD trigger particle loop
1873   if(GetDebug() > 1) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
1874 }
1875
1876
1877 //__________________________________________________________________
1878 void AliAnaParticleJetFinderCorrelation::Print(const Option_t * opt) const
1879 {
1880   
1881   //Print some relevant parameters set for the analysis
1882   if(! opt)
1883     return;
1884   
1885   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1886   AliAnaCaloTrackCorrBaseClass::Print(" ");
1887
1888   printf("Phi trigger-jet        <     %3.2f\n", fDeltaPhiMaxCut) ; 
1889   printf("Phi trigger-jet        >     %3.2f\n", fDeltaPhiMinCut) ;
1890   printf("pT Ratio trigger/jet   <     %3.2f\n", fRatioMaxCut) ; 
1891   printf("pT Ratio trigger/jet   >     %3.2f\n", fRatioMinCut) ;
1892   printf("fConeSize              =     %3.2f\n", fConeSize) ; 
1893   printf("fPtThresholdInCone     =     %3.2f\n", fPtThresholdInCone) ;
1894   printf("fUseJetRefTracks         =     %d\n",    fUseJetRefTracks) ;
1895   printf("fMakeCorrelationInHistoMaker     =     %d\n",    fMakeCorrelationInHistoMaker) ;      
1896   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
1897   printf("Reconstructed jet cone size = %3.2f\n", fJetConeSize) ;
1898   printf("Reconstructed jet minimum pt before background subtraction = %3.2f\n", fJetMinPt) ;
1899   printf("Reconstructed jet minimum area fraction = %3.2f\n", fJetAreaFraction) ;
1900
1901   //if(!fNonStandardJetFromReader){
1902   printf("fJetBranchName =   %s\n", fJetBranchName.Data()) ;
1903   //}
1904   if(!fBackgroundJetFromReader){
1905     printf("fBkgJetBranchName =   %s\n", fBkgJetBranchName.Data()) ;
1906   }
1907
1908   printf("Isolation cone size = %3.2f\n", fGammaConeSize) ;
1909   printf("fUseBackgroundSubtractionGamma = %d\n",fUseBackgroundSubtractionGamma);
1910   printf("fSaveGJTree = %d\n",fSaveGJTree);
1911   printf("fMostEnergetic = %d\n",fMostEnergetic);
1912   printf("fMostOpposite = %d\n",fMostOpposite);
1913
1914   printf("fUseHistogramJetBkg = %d\n",fUseHistogramJetBkg);
1915   printf("fUseHistogramTracks = %d\n",fUseHistogramTracks);
1916   printf("fUseHistogramJetTracks = %d\n",fUseHistogramJetTracks);
1917   printf("fMCStudies = %d\n",fMCStudies);
1918
1919
1920
1921 //__________________________________________________________________
1922 void AliAnaParticleJetFinderCorrelation::CalculateBkg(TVector3 gamma, TVector3 jet,Double_t vertex[3],Int_t type=2){
1923   //
1924   // calculate background for fragmentation function and fill histograms
1925   // 1. 90 degrees from jet axis in random place = perpendicular cone
1926   // 2. Random cone not belonging to jet cone nor photon cone
1927   // 3. In the middle point from jet and photon momentum vectors
1928   // 4. 90 degrees from photon direction in random place = perpendicular cone 2
1929
1930   //
1931   // implementation of 2 works, 1 and 4 works
1932   //
1933   Double_t gammaPt  = gamma.Pt();
1934   Double_t gammaEta = gamma.Eta();
1935   Double_t gammaPhi = gamma.Phi();
1936   Double_t jetEta   = jet.Eta();
1937   Double_t jetPhi   = jet.Phi();
1938
1939   //refference direction of background
1940   Double_t refEta=0.,refPhi=0.;
1941   Double_t rad = 0,rad2 = 0.;
1942   if(type==1){//perpendicular to jet axis
1943     //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
1944
1945     Double_t xVar;
1946     Double_t yVar;
1947     Double_t newX=0.;
1948     Double_t newY=0.;
1949     Double_t newZ=0.;
1950     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1951     Double_t jx=jet.Px();
1952     Double_t jy=jet.Py();
1953     Double_t jz=jet.Pz();
1954     //if(jz==0)printf("problem\n");
1955     //X axis 
1956     Double_t Xx=1.0-vertex[0];
1957     Double_t Xy=-1.0*vertex[1];
1958     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1959     //Y axis
1960     Double_t Yx=jy*Xz-jz*Xy;
1961     Double_t Yy=jz*Xx-jx*Xz;
1962     Double_t Yz=jx*Xy-jy*Xx;
1963     //Determinant
1964     Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
1965     if(det==0)printf("problem det==0\n");
1966     Double_t detX = 0.;
1967     Double_t detY = 0.;
1968     Double_t detZ = 0.;
1969
1970 //    Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
1971 //    printf("scalar jet.P o X: %f\n",tmpScalar);
1972 //    tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
1973 //    printf("scalar jet.P o Y: %f\n",tmpScalar);
1974 //    tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
1975 //    printf("scalar X o Y: %f\n",tmpScalar);
1976
1977     TVector3 perp;
1978     //randomise
1979     do{
1980       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1981       //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
1982       xVar=TMath::Cos(refPhi);
1983       yVar=TMath::Sin(refPhi);
1984       //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
1985       //zVar=0 in new surface frame
1986       detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
1987       detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
1988       detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
1989
1990       newX=detX/det;
1991       newY=detY/det;
1992       newZ=detZ/det;
1993
1994       perp.SetXYZ(newX,newY,newZ);
1995       refEta = perp.Eta();
1996       refPhi = perp.Phi();//output in <-pi, pi> range
1997       if(refPhi<0)refPhi+=2*TMath::Pi();
1998       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1999       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2000       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2001     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2002     fhRandomPhiEta[2]->Fill(refPhi,refEta);
2003
2004   }
2005   else if(type==2){//random cone
2006     //randomise
2007     do{
2008       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2009       refEta=fGenerator->Uniform(-(0.9-fJetConeSize),0.9-fJetConeSize);
2010       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2011       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2012       //check if reference is not within jet cone or gamma cone +0.4
2013       //example: FF fConeSize=0.4, fJetConeSize=0.4, fIsoGammaConeSize=0.3
2014     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize);
2015     //photon:0.7=0.4+0.3; jets:0.8=0.4 +0.4 //rad<fConeSize+fJetConeSize rad2<fConeSize+0.3
2016     fhRandomPhiEta[0]->Fill(refPhi,refEta);
2017   }
2018   else if(type==4){//perpendicular to photon axis
2019     Double_t xVar;
2020     Double_t yVar;
2021     Double_t newX=0.;
2022     Double_t newY=0.;
2023     Double_t newZ=0.;
2024     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
2025     Double_t jx=gamma.Px();
2026     Double_t jy=gamma.Py();
2027     Double_t jz=gamma.Pz();
2028     //if(jz==0)printf("problem\n");
2029     //X axis 
2030     Double_t Xx=1.0-vertex[0];
2031     Double_t Xy=-1.0*vertex[1];
2032     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2033     //Y axis
2034     Double_t Yx=jy*Xz-jz*Xy;
2035     Double_t Yy=jz*Xx-jx*Xz;
2036     Double_t Yz=jx*Xy-jy*Xx;
2037     //Determinant
2038     Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
2039     if(det==0)printf("problem det==0\n");
2040     Double_t detX = 0.;
2041     Double_t detY = 0.;
2042     Double_t detZ = 0.;
2043
2044 //    Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
2045 //    printf("scalar jet.P o X: %f\n",tmpScalar);
2046 //    tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
2047 //    printf("scalar jet.P o Y: %f\n",tmpScalar);
2048 //    tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
2049 //    printf("scalar X o Y: %f\n",tmpScalar);
2050
2051     TVector3 perp;
2052     //randomise
2053     do{
2054       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2055       //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
2056       xVar=TMath::Cos(refPhi);
2057       yVar=TMath::Sin(refPhi);
2058       //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
2059       //zVar=0 in new surface frame
2060       detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
2061       detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
2062       detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
2063
2064       newX=detX/det;
2065       newY=detY/det;
2066       newZ=detZ/det;
2067
2068       perp.SetXYZ(newX,newY,newZ);
2069       refEta = perp.Eta();
2070       refPhi = perp.Phi();//output in <-pi, pi> range
2071       if(refPhi<0)refPhi+=2*TMath::Pi();
2072       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2073       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2074       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2075     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2076     fhRandomPhiEta[1]->Fill(refPhi,refEta);
2077
2078   }
2079   else if(type==3){//mid point
2080
2081     Double_t jx=jet.Px();
2082     Double_t jy=jet.Py();
2083     Double_t jz=jet.Pz();
2084     //    if(jz==0)printf("problem\n");
2085     Double_t gx=gamma.Px();
2086     Double_t gy=gamma.Py();
2087     Double_t gz=gamma.Pz();
2088
2089     Double_t cosAlpha=(jx*gx+jy*gy+jz*gz)/(jet.Mag()*gamma.Mag());
2090     Double_t cosinus=TMath::Sqrt((cosAlpha+1.)/2.);
2091     //perpendicular axis
2092     Double_t Zx=gy*jz-gz*jy;
2093     Double_t Zy=gz*jx-gx*jz;
2094     Double_t Zz=gx*jy-gy*jx;
2095
2096     //Determinant
2097     Double_t det = Zx*gy*jz + Zy*gz*jx + Zz*gx*jy - Zz*gy*jx - Zy*gx*jz - Zx*gz*jy;
2098
2099     Double_t newX=0.;
2100     Double_t newY=0.;
2101     Double_t newZ=0.;
2102     if(det!=0) {
2103       Double_t detX =            -Zy*gz*cosinus +Zz*cosinus*jy + Zz*gy*cosinus - Zy*cosinus*jz;
2104       Double_t detY = Zx*cosinus*jz          - Zz*gx*cosinus - Zz*cosinus*jx            + Zx*gz*cosinus;
2105       Double_t detZ = -Zx*gy*cosinus + Zy*cosinus*jx + Zy*gx*cosinus - Zx*cosinus*jy;
2106
2107       newX=detX/det;
2108       newY=detY/det;
2109       newZ=detZ/det;
2110     }
2111
2112     TVector3 perp;
2113     perp.SetXYZ(newX,newY,newZ);
2114     refEta = perp.Eta();
2115     refPhi = perp.Phi();//output in <-pi, pi> range
2116     if(refPhi<0)refPhi+=2*TMath::Pi();
2117     rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2118     rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2119       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2120
2121     if (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize) fhRandomPhiEta[3]->Fill(refPhi,refEta);
2122   }
2123   else if(type==5){//tmp                                                                                                                                                   
2124     //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);                                                                                                         
2125
2126     Double_t xVar;
2127     Double_t newX=0.;
2128     Double_t newY=0.;
2129     Double_t newZ=0.;
2130     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]                                                                                                                          
2131     Double_t jx=jet.Px();
2132     Double_t jy=jet.Py();
2133     Double_t jz=jet.Pz();
2134     //    if(jz==0)printf("problem\n");
2135     //X axis                                                                                                                                                               
2136     Double_t Xx=1.0-vertex[0];
2137     Double_t Xy=-1.0*vertex[1];
2138     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2139     //Y axis                                                                                                                                                               
2140     Double_t Yx=jy*Xz-jz*Xy;
2141     Double_t Yy=jz*Xx-jx*Xz;
2142     Double_t Yz=jx*Xy-jy*Xx;
2143
2144     // X and Y length                                                                                                                                                      
2145     Double_t Xlength=TMath::Sqrt(Xx*Xx+Xy*Xy+Xz*Xz);
2146     Double_t Ylength=TMath::Sqrt(Yx*Yx+Yy*Yy+Yz*Yz);
2147     Double_t ratio=Ylength/Xlength;
2148
2149     TVector3 perp;
2150     //randomise                                                                                                                                                            
2151     do{
2152       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2153       xVar=TMath::Tan(refPhi)/ratio;
2154       newX=xVar*Yx+Xx;
2155       newY=xVar*Yy+Xy;
2156       newZ=xVar*Yz+Xz;
2157
2158       perp.SetXYZ(newX,newY,newZ);
2159       refEta = perp.Eta();
2160       refPhi = perp.Phi();//output in <-pi, pi> range                                                                                                                      
2161       if(refPhi<0)refPhi+=2*TMath::Pi();
2162       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2163       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2164       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2165     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2166     fhRandomPhiEta[4]->Fill(refPhi,refEta);
2167   }
2168
2169
2170
2171   //calculate FF in background
2172   Int_t ntracks =  0;
2173   ntracks =GetCTSTracks()->GetEntriesFast();
2174   AliVTrack* track = 0x0 ;
2175   TVector3 p3;
2176
2177   Double_t pt = 0, eta = 0, phi = 0;
2178   Int_t  npartcone = 0;
2179   Double_t sumPt=0.;
2180   for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
2181     track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
2182     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
2183     pt   = p3.Pt();
2184     if(pt<fPtThresholdInCone) {//0.150
2185       //printf("problem: track pt < %f MeV/c \n",fPtThresholdInCone);
2186       continue;
2187     }
2188     eta  = p3.Eta() ;
2189     phi  = p3.Phi() ;
2190     if(phi < 0) phi+=TMath::TwoPi();
2191     //Check if there is any particle inside cone with pt larger than  fPtThreshold
2192     rad = TMath::Sqrt((eta-refEta)*(eta-refEta) + (phi-refPhi)*(phi-refPhi));
2193     if(rad < fConeSize  && pt > fPtThresholdInCone){    
2194       //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
2195         npartcone++;
2196         sumPt+=pt;
2197         if(type==1){//perp jet
2198           fhBkgFFz[1] ->Fill(gammaPt, pt/gammaPt);
2199           fhBkgFFxi[1]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2200           fhBkgFFpt[1]->Fill(gammaPt, pt);
2201         }
2202         else if(type==2){//RC
2203           fhBkgFFz[0] ->Fill(gammaPt, pt/gammaPt);
2204           fhBkgFFxi[0]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2205           fhBkgFFpt[0]->Fill(gammaPt, pt);
2206         }
2207         else if(type==3){//mid point
2208           fhBkgFFz[3] ->Fill(gammaPt, pt/gammaPt);
2209           fhBkgFFxi[3]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2210           fhBkgFFpt[3]->Fill(gammaPt, pt);
2211         }
2212         else if(type==4){//perp jet
2213           fhBkgFFz[2] ->Fill(gammaPt, pt/gammaPt);
2214           fhBkgFFxi[2]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2215           fhBkgFFpt[2]->Fill(gammaPt, pt);
2216         }
2217         else if(type==5){//test
2218           fhBkgFFz[4] ->Fill(gammaPt, pt/gammaPt);
2219           fhBkgFFxi[4]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2220           fhBkgFFpt[4]->Fill(gammaPt, pt);
2221         }
2222
2223
2224     }
2225   }//end of loop over tracks
2226   Double_t sumOverTracks=0.;
2227   if(npartcone!=0) sumOverTracks = sumPt/npartcone;
2228   if(type==1) {
2229     fhBkgNTracksInCone[1]->Fill(gammaPt, npartcone);
2230     fhBkgSumPtInCone[1]->Fill(gammaPt,sumPt);
2231     fhBkgSumPtnTracksInCone[1]->Fill(gammaPt,sumOverTracks);
2232   }
2233   else if(type==2) {
2234     fhBkgNTracksInCone[0]->Fill(gammaPt, npartcone);
2235     fhBkgSumPtInCone[0]->Fill(gammaPt,sumPt);
2236     fhBkgSumPtnTracksInCone[0]->Fill(gammaPt,sumOverTracks);
2237   }
2238   else if(type==3) {
2239     fhBkgNTracksInCone[3]->Fill(gammaPt, npartcone);
2240     fhBkgSumPtInCone[3]->Fill(gammaPt,sumPt);
2241     fhBkgSumPtnTracksInCone[3]->Fill(gammaPt,sumOverTracks);
2242   }
2243   else if(type==4) {
2244     fhBkgNTracksInCone[2]->Fill(gammaPt, npartcone);
2245     fhBkgSumPtInCone[2]->Fill(gammaPt,sumPt);
2246     fhBkgSumPtnTracksInCone[2]->Fill(gammaPt,sumOverTracks);
2247   }
2248   else if(type==5) {
2249     fhBkgNTracksInCone[4]->Fill(gammaPt, npartcone);
2250     fhBkgSumPtInCone[4]->Fill(gammaPt,sumPt);
2251     fhBkgSumPtnTracksInCone[4]->Fill(gammaPt,sumOverTracks);
2252   }
2253 }
2254
2255
2256
2257
2258
2259 //__________________________________________________________________
2260 void AliAnaParticleJetFinderCorrelation::FindMCgenInfo(){
2261   //
2262   // Find information about photon and (quark or gluon) on generated level 
2263   //
2264
2265   //frequently used variables
2266   Int_t pdg    = 0 ;
2267   Int_t mother = -1 ; 
2268   Int_t absID  = 0 ;
2269
2270   //Double_t photonY   = -100 ;
2271   //Double_t photonE   = -1 ;
2272   Double_t photonPt  = -1 ;
2273   Double_t photonPhi =  100 ;
2274   Double_t photonEta = -1 ;
2275   Bool_t   inacceptance = kFALSE;
2276   AliAODMCParticle * primTmp = NULL;
2277
2278   //jet counters
2279   Int_t nParticlesInJet=0;
2280   Int_t nChargedParticlesInJet=0;
2281   Int_t nParticlesInJet150=0;
2282   Int_t nChargedParticlesInJet150=0;
2283   Int_t nChargedParticlesInJet150Cone=0;
2284
2285   Double_t eneParticlesInJet=0.;
2286   Double_t eneChargedParticlesInJet=0.;
2287   Double_t eneParticlesInJet150=0.;
2288   Double_t eneChargedParticlesInJet150=0.;
2289   Double_t eneChargedParticlesInJet150Cone=0.;
2290
2291   Double_t pxParticlesInJet=0.;
2292   Double_t pxChargedParticlesInJet=0.;
2293   Double_t pxParticlesInJet150=0.;
2294   Double_t pxChargedParticlesInJet150=0.;
2295   Double_t pxChargedParticlesInJet150Cone=0.;
2296
2297   Double_t pyParticlesInJet=0.;
2298   Double_t pyChargedParticlesInJet=0.;
2299   Double_t pyParticlesInJet150=0.;
2300   Double_t pyChargedParticlesInJet150=0.;
2301   Double_t pyChargedParticlesInJet150Cone=0.;
2302
2303   Double_t etaParticlesInJet=0.;
2304   Double_t etaChargedParticlesInJet=0.;
2305   Double_t etaParticlesInJet150=0.;
2306   Double_t etaChargedParticlesInJet150=0.;
2307   Double_t etaChargedParticlesInJet150Cone=0.;
2308
2309   Double_t phiParticlesInJet=0.;
2310   Double_t phiChargedParticlesInJet=0.;
2311   Double_t phiParticlesInJet150=0.;
2312   Double_t phiChargedParticlesInJet150=0.;
2313   Double_t phiChargedParticlesInJet150Cone=0.;
2314
2315   Double_t ptParticlesInJet=0.;
2316   Double_t ptChargedParticlesInJet=0.;
2317   Double_t ptParticlesInJet150=0.;
2318   Double_t ptChargedParticlesInJet150=0.;
2319   Double_t ptChargedParticlesInJet150Cone=0.;
2320
2321   Double_t coneJet=0.;
2322   Double_t coneChargedJet=0.;
2323   Double_t coneJet150=0.;
2324   Double_t coneChargedJet150=0.;
2325
2326   std::vector<Int_t> jetParticleIndex;
2327
2328   if(GetReader()->ReadStack()) {//ESD
2329      if(GetDebug()>3) printf("I use stack\n");
2330   }//end Stack 
2331   else if(GetReader()->ReadAODMCParticles()) {//AOD
2332     TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
2333     if(mcparticles){
2334       //jet origin
2335       //index =6 and 7 is hard scattering (jet-quark or photon)
2336       primTmp = (AliAODMCParticle *) mcparticles->At(6);
2337       pdg=primTmp->GetPdgCode();
2338        if(GetDebug()>3) printf("id 6 pdg %d, pt %f ",pdg,primTmp->Pt() );
2339       if(TMath::Abs(pdg)<=6 ||pdg==21) {
2340         fhMCJetOrigin->Fill(pdg);
2341         fMCPartonType=pdg;
2342       }
2343       primTmp = (AliAODMCParticle *) mcparticles->At(7);
2344       pdg=primTmp->GetPdgCode();
2345        if(GetDebug()>3) printf("id 7 pdg %d, pt %f\n",pdg,primTmp->Pt() );
2346       if(TMath::Abs(pdg)<=6 ||pdg==21) {
2347         fhMCJetOrigin->Fill(pdg);
2348         fMCPartonType=pdg;
2349       }
2350       //end of jet origin
2351
2352       Int_t nprim = mcparticles->GetEntriesFast();
2353       for(Int_t i=0; i < nprim; i++) {
2354         if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
2355         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
2356         pdg = prim->GetPdgCode();
2357         mother=prim->GetMother();
2358         //photon=22, gluon=21, quarks=(1,...,6), antiquarks=(-1,...,-6)
2359         if(pdg == 22){//photon
2360           fhMCPhotonCuts->Fill(0);
2361           if(prim->GetStatus()!=1) continue;
2362           fhMCPhotonCuts->Fill(1);
2363            if(GetDebug()>5) printf("id %d, prim %d, physPrim %d, status %d\n",i,prim->IsPrimary(),prim->IsPhysicalPrimary(),prim->GetStatus());
2364           while(mother>7){
2365             primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2366             mother=primTmp->GetMother();
2367           }
2368           if(mother<6)continue;
2369           fhMCPhotonCuts->Fill(2);
2370           primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2371           if(primTmp->GetPdgCode()!=22)continue;
2372           fhMCPhotonCuts->Fill(3);
2373
2374           //Get photon kinematics
2375           photonPt  = prim->Pt() ;
2376           photonPhi = prim->Phi() ;
2377           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
2378           photonEta = prim->Eta() ;
2379           fhMCPhotonPt->Fill(photonPt);
2380           fhMCPhotonEtaPhi->Fill(photonPhi,photonEta);
2381           
2382           //Check if photons hit the Calorimeter
2383           fMomentum.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
2384           inacceptance = kFALSE;
2385           if(GetCaloUtils()->IsEMCALGeoMatrixSet()){
2386             fhMCPhotonCuts->Fill(4);
2387             //check if in EMCAL
2388             if(GetEMCALGeometry()) {
2389               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
2390               if(absID >= 0) inacceptance = kTRUE;
2391               if(GetDebug() > 3) printf("In EMCAL Real acceptance? %d\n",inacceptance);
2392             }
2393             else{
2394               if(GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) inacceptance = kTRUE ;
2395               if(GetDebug() > 3) printf("In EMCAL fiducial cut acceptance? %d\n",inacceptance);
2396             }
2397           }else{//no EMCAL nor EMCALGeoMatrixSet
2398             printf("not EMCALGeoMatrix set\n");
2399           }//end of check if EMCAL
2400           if(inacceptance)fhMCPhotonCuts->Fill(5);
2401           if(GetDebug() > 5)
2402             printf("Photon Energy %f, Pt %f\n",prim->E(),prim->Pt());
2403           fMCGamPt=photonPt;
2404           fMCGamEta=photonEta;
2405           fMCGamPhi=photonPhi;
2406         }//end of check if photon
2407         else{//not photon
2408           if(prim->GetStatus()!=1) continue;
2409           if(GetDebug() > 5)
2410             printf("id %d, prim %d, physPrim %d, status %d, pdg %d, E %f ",i,prim->IsPrimary(),prim->IsPhysicalPrimary(),prim->GetStatus(),prim->GetPdgCode(),prim->E());
2411           while(mother>7){
2412             primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2413             mother=primTmp->GetMother();
2414             if(GetDebug() > 5) printf("next mother %d ",mother);
2415           }
2416           if(GetDebug() > 5) printf("\n");
2417           if(mother<6)continue;//soft part
2418           primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2419           pdg=primTmp->GetPdgCode();
2420           if( !(TMath::Abs(pdg)<=6 || pdg==21) ) continue;//origin not hard q or g
2421
2422           //jetParticleIndex.Add(&i);
2423           jetParticleIndex.push_back(i);
2424
2425         nParticlesInJet++;
2426         eneParticlesInJet+=prim->E();
2427         pxParticlesInJet+=prim->Px();
2428         pyParticlesInJet+=prim->Py();
2429         etaParticlesInJet+=(prim->E()*prim->Eta());
2430         photonPhi = prim->Phi() ;
2431         if(photonPhi < 0) photonPhi+=TMath::TwoPi();
2432         phiParticlesInJet+=(prim->E()*photonPhi);
2433         
2434         if(prim->Charge()!=0) {
2435           nChargedParticlesInJet++;
2436           eneChargedParticlesInJet+=prim->E();
2437           pxChargedParticlesInJet+=prim->Px();
2438           pyChargedParticlesInJet+=prim->Py();
2439           etaChargedParticlesInJet+=(prim->E()*prim->Eta());
2440           phiChargedParticlesInJet+=(prim->E()*photonPhi);
2441         }
2442         if(prim->Pt()>0.150 && TMath::Abs(prim->Eta())<0.9 ) {//in TPC acceptance :)
2443           nParticlesInJet150++;
2444           eneParticlesInJet150+=prim->E();
2445           pxParticlesInJet150+=prim->Px();
2446           pyParticlesInJet150+=prim->Py();
2447           etaParticlesInJet150+=(prim->E()*prim->Eta());
2448           phiParticlesInJet150+=(prim->E()*photonPhi);
2449         }
2450         if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(prim->Eta())<0.9 ) {//in TPC acceptance :)
2451           nChargedParticlesInJet150++;
2452           eneChargedParticlesInJet150+=prim->E();
2453           pxChargedParticlesInJet150+=prim->Px();
2454           pyChargedParticlesInJet150+=prim->Py();
2455           etaChargedParticlesInJet150+=(prim->E()*prim->Eta());
2456           phiChargedParticlesInJet150+=(prim->E()*photonPhi);
2457         }
2458
2459         }//end of check pdg
2460       }//end of loop over primaries
2461
2462       if(eneParticlesInJet != 0.) {
2463         etaParticlesInJet/=eneParticlesInJet ;
2464         phiParticlesInJet/=eneParticlesInJet ;
2465       }
2466       if(eneChargedParticlesInJet != 0) {
2467         etaChargedParticlesInJet/=eneChargedParticlesInJet;
2468         phiChargedParticlesInJet/=eneChargedParticlesInJet;
2469       }
2470       if(eneParticlesInJet150 != 0) {
2471         etaParticlesInJet150/=eneParticlesInJet150;
2472         phiParticlesInJet150/=eneParticlesInJet150;
2473       }
2474       if(eneChargedParticlesInJet150 != 0) {
2475         etaChargedParticlesInJet150/=eneChargedParticlesInJet150;
2476         phiChargedParticlesInJet150/=eneChargedParticlesInJet150;
2477       }
2478
2479       ptParticlesInJet=TMath::Sqrt(pxParticlesInJet*pxParticlesInJet+pyParticlesInJet*pyParticlesInJet);
2480       ptChargedParticlesInJet=TMath::Sqrt(pxChargedParticlesInJet*pxChargedParticlesInJet+pyChargedParticlesInJet*pyChargedParticlesInJet);
2481       ptParticlesInJet150=TMath::Sqrt(pxParticlesInJet150*pxParticlesInJet150+pyParticlesInJet150*pyParticlesInJet150);
2482       ptChargedParticlesInJet150=TMath::Sqrt(pxChargedParticlesInJet150*pxChargedParticlesInJet150+pyChargedParticlesInJet150*pyChargedParticlesInJet150);
2483
2484       Double_t distance=0.;
2485       Double_t eta=0.;
2486       Double_t phi=0.;
2487       Double_t mostPtCharged=0.;
2488       Int_t mostmostPtChargedId=-1;
2489       std::vector<Int_t>::iterator it;
2490       for( it=jetParticleIndex.begin(); it!=jetParticleIndex.end(); ++it ){
2491         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(*it);
2492         eta = prim->Eta();
2493         phi = prim->Phi();
2494         if(phi < 0) phi+=TMath::TwoPi();
2495         //full jet
2496         distance=TMath::Sqrt((eta-etaParticlesInJet)*(eta-etaParticlesInJet)+(phi-phiParticlesInJet)*(phi-phiParticlesInJet));
2497         if(distance>coneJet) coneJet=distance;
2498         //charged jet
2499         distance=TMath::Sqrt((eta-etaChargedParticlesInJet)*(eta-etaChargedParticlesInJet)+(phi-phiChargedParticlesInJet)*(phi-phiChargedParticlesInJet));
2500         if(distance>coneChargedJet) coneChargedJet=distance;
2501         //
2502         distance=TMath::Sqrt((eta-etaParticlesInJet150)*(eta-etaParticlesInJet150)+(phi-phiParticlesInJet150)*(phi-phiParticlesInJet150));
2503         if(distance>coneJet150 && TMath::Abs(eta)<0.9 ) coneJet150=distance;
2504         //
2505         distance=TMath::Sqrt((eta-etaChargedParticlesInJet150)*(eta-etaChargedParticlesInJet150)+(phi-phiChargedParticlesInJet150)*(phi-phiChargedParticlesInJet150));
2506         if(distance>coneChargedJet150 && TMath::Abs(eta)<0.9) coneChargedJet150=distance;
2507
2508         if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(eta)<0.9) {
2509           if(prim->Pt()>mostPtCharged) {
2510             mostPtCharged=prim->Pt();
2511             mostmostPtChargedId=(*it);
2512           }
2513         }
2514
2515         if(distance<=0.4){
2516           if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(eta)<0.9) {
2517             nChargedParticlesInJet150Cone++;
2518             eneChargedParticlesInJet150Cone+=prim->E();
2519             pxChargedParticlesInJet150Cone+=prim->Px();
2520             pyChargedParticlesInJet150Cone+=prim->Py();
2521             etaChargedParticlesInJet150Cone+=(prim->E()*eta);
2522             phiChargedParticlesInJet150Cone+=(prim->E()*phi);
2523           }
2524
2525         }
2526
2527       }//end of loop over jet particle indexes
2528       if(eneChargedParticlesInJet150Cone != 0) {
2529         etaChargedParticlesInJet150Cone/=eneChargedParticlesInJet150Cone;
2530         phiChargedParticlesInJet150Cone/=eneChargedParticlesInJet150Cone;
2531       }
2532       ptChargedParticlesInJet150Cone=TMath::Sqrt(pxChargedParticlesInJet150Cone*pxChargedParticlesInJet150Cone+pyChargedParticlesInJet150Cone*pyChargedParticlesInJet150Cone);
2533       if(nChargedParticlesInJet150>0 && nChargedParticlesInJet150Cone<1){//no particles in cone? take the most energetic one
2534         nChargedParticlesInJet150Cone=1;
2535         etaChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Eta();
2536         phiChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Phi();
2537         ptChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Pt();
2538       }
2539
2540
2541     }//mc array exists and data is MC
2542   }// read AOD MC
2543
2544   jetParticleIndex.clear();
2545
2546
2547   //printouts
2548   if(GetDebug() > 3) {
2549     printf("cone full %f, charged %f, full150 %f, charged150 %f\n",coneJet,coneChargedJet,coneJet150,coneChargedJet150);
2550     printf("Npart %d, NchPart %d, Npart(pt>150M) %d, NchPart(pt>150M) %d, NchPart(pt>150M)Cone %d\n",nParticlesInJet,nChargedParticlesInJet,nParticlesInJet150,nChargedParticlesInJet150,nChargedParticlesInJet150Cone);
2551     printf("Etot %f, Ech %f, E(pt>150M) %f, Ech(pt>150M) %f\n",eneParticlesInJet,eneChargedParticlesInJet,eneParticlesInJet150,eneChargedParticlesInJet150);
2552     printf("pt %f, ptch %f, pt(pt>150M) %f,ptch(pt>150M) %f,ptch(pt>150M)Cone %f\n",ptParticlesInJet,ptChargedParticlesInJet,ptParticlesInJet150,ptChargedParticlesInJet150,ptChargedParticlesInJet150Cone);
2553     printf("eta/phi tot %f/%f, ch %f/%f, tot150 %f/%f,  ch150 %f/%f, ch150cone %f/%f\n",etaParticlesInJet,phiParticlesInJet,etaChargedParticlesInJet,phiChargedParticlesInJet,etaParticlesInJet150,phiParticlesInJet150,etaChargedParticlesInJet150,phiChargedParticlesInJet150,etaChargedParticlesInJet150Cone,phiChargedParticlesInJet150Cone);
2554   }
2555   //fill histograms
2556   if(ptParticlesInJet) fhMCJetRatioChFull->Fill(ptChargedParticlesInJet/ptParticlesInJet);
2557   if(ptChargedParticlesInJet) fhMCJetRatioCh150Ch->Fill(ptChargedParticlesInJet150/ptChargedParticlesInJet);
2558
2559   fhMCJetNPartVsPt     ->Fill(ptParticlesInJet,nParticlesInJet);
2560   fhMCJetChNPartVsPt   ->Fill(ptChargedParticlesInJet,nChargedParticlesInJet);
2561   fhMCJetNPart150VsPt  ->Fill(ptParticlesInJet150,nParticlesInJet150);
2562   fhMCJetChNPart150VsPt->Fill(ptChargedParticlesInJet150,nChargedParticlesInJet150);
2563   fhMCJetChNPart150ConeVsPt->Fill(ptChargedParticlesInJet150Cone,nChargedParticlesInJet150Cone);
2564
2565   fhMCJetEtaPhi->Fill(phiParticlesInJet,etaParticlesInJet);
2566   fhMCJetChEtaPhi->Fill(phiChargedParticlesInJet,etaChargedParticlesInJet);
2567   fhMCJet150EtaPhi->Fill(phiParticlesInJet150,etaParticlesInJet150);
2568   fhMCJetCh150EtaPhi->Fill(phiChargedParticlesInJet150,etaChargedParticlesInJet150);
2569   fhMCJetCh150ConeEtaPhi->Fill(phiChargedParticlesInJet150Cone,etaChargedParticlesInJet150Cone);
2570
2571   //fill tree
2572   fMCJetPt      = ptParticlesInJet;
2573   fMCJetChPt    = ptChargedParticlesInJet;      
2574   fMCJet150Pt   = ptParticlesInJet150;     
2575   fMCJetCh150Pt = ptChargedParticlesInJet150;   
2576   fMCJetNPart      = nParticlesInJet;     
2577   fMCJetChNPart    = nChargedParticlesInJet;   
2578   fMCJet150NPart   = nParticlesInJet150;  
2579   fMCJetCh150NPart = nChargedParticlesInJet150;
2580   fMCJetEta      = etaParticlesInJet          ;
2581   fMCJetPhi      = phiParticlesInJet          ;
2582   fMCJetChEta    = etaChargedParticlesInJet   ;
2583   fMCJetChPhi    = phiChargedParticlesInJet   ;
2584   fMCJet150Eta   = etaParticlesInJet150       ;
2585   fMCJet150Phi   = phiParticlesInJet150       ;
2586   fMCJetCh150Eta = etaChargedParticlesInJet150;
2587   fMCJetCh150Phi = phiChargedParticlesInJet150;
2588
2589   fMCJetCh150ConePt    = ptChargedParticlesInJet150Cone;
2590   fMCJetCh150ConeNPart = nChargedParticlesInJet150Cone;
2591   fMCJetCh150ConeEta   = etaChargedParticlesInJet150Cone;
2592   fMCJetCh150ConePhi   = phiChargedParticlesInJet150Cone;
2593
2594 }//end of method FindMCgenInfo