Modified plots and made jet finder use SDigits
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALJetFinderPlots.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 /* $Id$ */
18
19
20 //_________________________________________________________________________
21 //  Class for Filling JetFinder Plots
22 // --
23 //*-- Author: Mark Horner (LBL/UCT)
24 // --
25 // --
26
27
28 #include "TMath.h"
29 #include "AliEMCALJetFinderPlots.h"
30
31 ClassImp(AliEMCALJetFinderPlots)
32         
33 AliEMCALJetFinderPlots::AliEMCALJetFinderPlots()
34 {
35         // Constructor to initialise variables
36   fInitialised = kFALSE;
37   fNominalEnergy = 0.0;
38   fConeRadius = 0.3;
39   fDebug = 0;
40   fOutput=0;
41     fhFragmFcn=0;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
42     fhPartonFragmFcn=0;// = new TH1F("hPartonFragmFcn","Fragmentation Function",100,0,1);
43    fhPartonJT=0;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
44     fhPartonPL=0;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
45     fhJetJT=0;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
46     fhJetPL=0;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
47     fhJetEt=0;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
48     fhJetEta=0;// = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
49     fhJetPhi=0;// = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
50     fhPartonEta=0;// = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
51     fhPartonPhi=0;// = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
52     fhEtaDiff=0;// = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
53     fhPhiDiff=0;//  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
54     fhNJets=0;// = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
55   fhEtaPhiSpread=0;         
56
57 fhFragmFcn2=0;  // ("hFragmFcn2","Fragmentation Function",100,0,1);
58 fhPartonFragmFcn2=0;// ("hFragmFcn2","Parton Fragmentation Function",100,0,1);
59 fhPartonJT2=0;  // ("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
60 fhPartonPL2=0;  // ("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
61 fhJetJT2=0;     // ("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
62 fhJetPL2=0;     // ("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
63 fhJetEt2=0;     // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
64 fhJetEta2=0;    // ("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
65 fhJetPhi2=0;    // ("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
66 fhPartonEta2=0; // ("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
67 fhPartonPhi2=0; // ("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
68 fhEtaDiff2=0;   // ("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
69 fhPhiDiff2=0;   // ("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
70 fhEtaPhiSpread2=0;      // ("hEtaPhiSpread2","#eta - #phi Distribution 
71                                                 //of Reconstructed Jets",192,-0.7,0.7,288,pi/3,pi);
72 fhNJets2=0;     // ("hNJets2","N Reconstructed jets",11,-0.5,10.5);
73 fhJetEtSecond2=0; //("hJetEtSecond2","E_{T}^{reco}",250,0.,250.); 
74 fhJetEtRatio2=0;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
75 fhEtaPhiDist2=0;  //("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
76 fhInputOutput=0;
77 //     TH2F                            *fhInputOutput;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
78      
79 fhRecoBinPt=0;         // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
80 fhRecoBinPtNoBg=0;             // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
81 fhRecoBinPartonPt=0;    // ("fhRecoBinPartonPt","Input Pt Distribution",100,0,1);
82 fhRecoBinJetEt=0;       // ("fhRecoJetEt","E_{T}^{reco}",250,0.,250.);
83 fhRecoBinInputJetEt=0;  // ("fhRecoInputJetEt","E_{T}^{reco}",250,0.,250.);
84 fhJetPT =0;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
85 fhPartonPT =0;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
86 fhJetPT2 =0;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
87 fhPartonPT2 =0;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
88 fhRecoBinFragmFcn =0;//new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
89 fhRecoBinFragmFcnNoBg =0;//new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
90 fhRecoBinPartonFragmFcn =0;// new TH1F("fhRecoBinPartonFragmFcn","Input Bin Fragm Fcn Distribution",100,0,1);
91 fhJetInvE=0;// = new TH1F("fhJetInvE","#frac{1}{E_{R}}",100,0,1);
92 fhJetInvE2=0;// = new TH1F("fhJetInvE","#frac{1}{E_{R}}",100,0,1);
93 fScaleFactor = 1.0/0.6731;
94 fhBackHisto=0;
95
96 }
97
98 void AliEMCALJetFinderPlots::InitPlots()
99 {
100 //========================= CASE 1 =======================================      
101     fhFragmFcn = new TH1F("hFragmFcn","Fragmentation Function",200,0,2);
102     fhFragmFcn->Sumw2();
103     fhJetPT = new TH1F("hJetPT","P_{T} Distribution",200,0,200);
104     fhJetPT->Sumw2();
105     fhPartonPT = new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,200);
106     fhPartonPT->Sumw2();
107     fhPartonFragmFcn = new TH1F("hPartonFragmFcn","Parton Fragmentation Function",200,0,2);
108     fhPartonFragmFcn->Sumw2();
109     fhPartonJT = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
110     fhPartonJT->Sumw2();
111     fhPartonPL = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
112     fhPartonPL->Sumw2();
113     fhJetJT = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
114     fhJetJT->Sumw2();
115     fhJetPL = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
116     fhJetPL->Sumw2();
117     fhJetEt = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
118     fhJetEt->Sumw2();
119     fhJetEtDiff = new TH1F("hJetEtDiff","E_{T}^{reco}-E_{T}^{Parton}",250,-124.,125.);
120     fhJetEtDiff->Sumw2();
121     fhJetEta = new      TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
122     fhJetEta->Sumw2();
123     fhJetPhi = new      TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
124     fhJetPhi->Sumw2();
125     fhPartonEta = new   TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
126     fhPartonEta->Sumw2();
127     fhPartonPhi = new   TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
128     fhPartonPhi->Sumw2();
129     fhEtaDiff = new     TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
130     fhEtaDiff->Sumw2();
131     fhPhiDiff  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
132     fhPhiDiff->Sumw2();
133     fhNJets = new       TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
134     fhNJets->Sumw2();
135     fhEtaPhiSpread = new TH2F("hEtaPhiSpread","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
136     fhEtaPhiSpread->Sumw2();    
137   fhNJets->SetXTitle("N_{jets}^{reco}/event");
138   fhNJets->SetYTitle("N_{events}");
139
140   //Jet properties
141   fhJetEt->SetFillColor(16);
142   fhJetEt->SetXTitle("E_{T}^{reco}");
143   
144   fhJetEta->SetFillColor(16);
145   fhJetEta->SetXTitle("#eta_{jet}^{reco}");
146   
147   fhJetPhi->SetFillColor(16);
148   fhJetPhi->SetXTitle("#phi_{jet}^{reco}");
149   
150   fhPartonEta->SetFillColor(16);
151   fhPartonEta->SetXTitle("#eta_{parton}");
152   
153   fhPartonPhi->SetFillColor(16);
154   fhPartonPhi->SetXTitle("#phi_{parton}");
155
156   fhPartonPL->SetXTitle("p (GeV/c)");
157   fhPartonJT->SetXTitle("p (GeV/c)");
158
159   fhPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{parton}");
160
161   //Jet component properties
162
163   fhJetPL->SetXTitle("p (GeV/c)");
164   fhJetJT->SetXTitle("p (GeV/c)");
165   fhFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
166   fhPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
167
168   fhEtaDiff->SetXTitle("#eta_{jet}^{reco}-#eta_{jet}^{input}");
169   fhPhiDiff->SetXTitle("#phi_{jet}^{reco}-#phi_{jet}^{input}");
170   fhEtaPhiSpread->SetXTitle("#eta");
171   fhEtaPhiSpread->SetYTitle("#phi");
172
173 //======================= CASE 2 ======================================
174   
175
176 fhFragmFcn2             = new TH1F("hFragmFcn2","Fragmentation Function",200,0,2);
177 fhFragmFcn2->Sumw2();
178     fhJetPT2 = new TH1F("hJetPT2","P_{T} Distribution",200,0,200);
179     fhJetPT2->Sumw2();
180     fhPartonPT2 = new TH1F("hPartonPT2","Parton P_{T} Distribution",200,0,1);
181     fhPartonPT2->Sumw2();
182 fhPartonFragmFcn2       = new TH1F("hPartonFragmFcn2","Parton Fragmentation Function",200,0,2);
183 fhPartonFragmFcn2->Sumw2();
184 fhPartonJT2             = new TH1F("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
185 fhPartonJT2->Sumw2();
186 fhPartonPL2             = new TH1F("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
187 fhPartonPL2->Sumw2();
188 fhJetJT2                        = new TH1F("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
189 fhJetJT2->Sumw2();
190 fhJetPL2                        = new TH1F("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
191 fhJetPL2->Sumw2();
192 fhJetEt2                        = new TH1F("hJetEt2","E_{T}^{reco}",250,0.,250.);
193 fhJetEt2->Sumw2();
194     fhJetEtDiff2 = new TH1F("hJetEtDiff2","E_{T}^{reco}-E_{T}^{Parton}",250,-124.,125.);
195     fhJetEtDiff2->Sumw2();
196 fhJetEta2               = new TH1F("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
197 fhJetEta2->Sumw2();
198 fhJetPhi2               = new TH1F("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
199 fhJetPhi2->Sumw2();
200 fhPartonEta2            = new TH1F("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
201 fhPartonEta2->Sumw2();
202 fhPartonPhi2            = new TH1F("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
203 fhPartonPhi2->Sumw2();
204 fhEtaDiff2              = new TH1F("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
205 fhEtaDiff2->Sumw2();
206 fhPhiDiff2              = new TH1F("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
207 fhPhiDiff2->Sumw2();
208 fhEtaPhiSpread2 = new TH2F("hEtaPhiSpread2","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
209 fhEtaPhiSpread2->Sumw2();
210 fhNJets2                        = new TH1F("hNJets2","N Reconstructed jets",11,-0.5,10.5);
211 fhNJets2->Sumw2();
212 fhJetEtSecond2          = new TH1F("hJetEtSecond2","E_{T}^{reco}",250,0.,250.); 
213 fhJetEtSecond2->Sumw2();
214 fhJetEtRatio2           = new TH1F("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
215 fhJetEtRatio2->Sumw2();
216 fhEtaPhiDist2           = new TH1F("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
217 fhEtaPhiDist2->Sumw2();
218
219 fhInputOutput=  new TH2F("hInputOutput","Input and Reconstruction Correlations;Input;Output",200,0,200,200,0,200);  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
220
221 //============================== Reconstruction Bin Comparison  ============================================
222
223 fhRecoBinPt =new TH1F("fhRecoBinPt","Reconstructed Pt Distribution",200,0,200);
224 fhRecoBinPt->Sumw2();
225 fhRecoBinPtNoBg =new TH1F("fhRecoBinPtNoBg","Reconstructed Pt Distribution Background Subtracted",200,0,200);
226 fhRecoBinPtNoBg->Sumw2();
227 fhRecoBinPartonPt = new TH1F("fhRecoBinPartonPt","Input Pt Distribution",200,0,200);
228 fhRecoBinPartonPt->Sumw2();
229 fhRecoBinFragmFcn =new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",200,0,2);
230 fhRecoBinFragmFcn->Sumw2();
231 fhRecoBinFragmFcnNoBg =new TH1F("fhRecoBinFragmFcnNoBg","Reconstructed Frag. Fcn With Background Removed",200,0,2);
232 fhRecoBinFragmFcnNoBg->Sumw2();
233 fhRecoBinPartonFragmFcn = new TH1F("fhRecoBinPartonFragmFcn","Input Bin Fragm Fcn Distribution",200,0,2);
234 fhRecoBinPartonFragmFcn->Sumw2();
235 fhRecoBinJetEt = new TH1F("fhRecoJetEt","E_{T}^{reco}",250,0.,250.);
236 fhRecoBinJetEt->Sumw2();
237 fhRecoBinInputJetEt = new TH1F("fhRecoInputJetEt","E_{T}^{reco}",250,0.,250.);
238 fhRecoBinInputJetEt->Sumw2();
239         
240
241 fhJetInvE = new TH1F("fhJetInvE","#frac{1}{E_{R}}",100,0,1);
242 fhJetInvE->Sumw2();
243 fhJetInvE2 = new TH1F("fhJetInvE2","#frac{1}{E_{R}}",100,0,1);
244 fhJetInvE2->Sumw2();
245                 
246
247
248   fInitialised = kTRUE; 
249   
250 }
251
252 AliEMCALJetFinderPlots::~AliEMCALJetFinderPlots()
253 {
254         // To ensure that all requested memory is returned
255 delete    fhFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
256 delete    fhPartonFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
257 delete   fhPartonJT;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
258 delete    fhPartonPL;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
259 delete    fhJetJT;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
260 delete    fhJetPL;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
261 delete    fhJetEt;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
262 delete  fhJetEtDiff;    // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
263 delete    fhJetEta;// = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
264 delete    fhJetPhi;// = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
265 delete    fhPartonEta;// = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
266 delete    fhPartonPhi;// = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
267 delete    fhEtaDiff;// = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
268 delete    fhPhiDiff;//  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
269 delete    fhNJets;// = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
270 delete    fhEtaPhiSpread;           
271
272         delete                          fhFragmFcn2;    // ("hFragmFcn2","Fragmentation Function",100,0,1);
273         delete                          fhPartonFragmFcn2;// ("hFragmFcn2","Parton Fragmentation Function",100,0,1);
274         delete                          fhPartonJT2;    // ("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
275         delete                          fhPartonPL2;    // ("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
276         delete                          fhJetJT2;       // ("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
277         delete                          fhJetPL2;       // ("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
278         delete                          fhJetEt2;       // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
279         delete                          fhJetEtDiff2;   // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
280         delete                          fhJetEta2;      // ("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
281         delete                          fhJetPhi2;      // ("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
282         delete                          fhPartonEta2;   // ("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
283         delete                          fhPartonPhi2;   // ("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
284         delete                          fhEtaDiff2;     // ("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
285         delete                          fhPhiDiff2;     // ("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
286         delete                          fhEtaPhiSpread2;        // ("hEtaPhiSpread2","#eta - #phi Distribution 
287                                                         //of Reconstructed Jets",192,-0.7,0.7,288,pi/3,pi);
288         delete                          fhNJets2;       // ("hNJets2","N Reconstructed jets",11,-0.5,10.5);
289         delete                          fhJetEtSecond2; //("hJetEtSecond2","E_{T}^{reco}",250,0.,250.); 
290         delete                          fhJetEtRatio2;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
291         delete                          fhEtaPhiDist2;  //("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
292
293         delete fhRecoBinPt;          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
294         delete fhRecoBinPtNoBg;          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
295         delete fhRecoBinPartonPt;    // ("fhRecoBinPartonPt","Input Pt Distribution",100,0,1);
296         delete fhRecoBinJetEt;       // ("fhRecoJetEt","E_{T}^{reco}",250,0.,250.);
297         delete fhRecoBinInputJetEt;  // ("fhRecoInputJetEt","E_{T}^{reco}",250,0.,250.);
298                                         
299         delete fhJetPT ;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
300         delete fhPartonPT ;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
301         delete fhJetPT2 ;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
302         delete fhPartonPT2 ;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
303         delete fhRecoBinFragmFcn;//new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
304         delete fhRecoBinFragmFcnNoBg;//new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
305         delete fhRecoBinPartonFragmFcn;// new TH1F("fhRecoBinPartonFragmFcn","Input Bin Fragm Fcn Distribution",100,0,1);
306     delete fhJetInvE;// = new TH1F("fhJetInvE","#frac{1}{E_{R}}",100,0,1);
307     delete fhJetInvE2;// = new TH1F("fhJetInvE","#frac{1}{E_{R}}",100,0,1);
308
309 }       
310
311 void AliEMCALJetFinderPlots::FillFromOutput(AliEMCALJetFinderOutput* output)
312 {
313         // Fill histograms from an output object
314 if (!fInitialised) InitPlots(); 
315   fOutput = output;
316   if (!fOutput) return;
317 // Make some temporary histograms to make sure we subtract 
318   // background properly
319 /*
320 tempFragmFcnNoBg =new TH1F("tempFragmFcnNoBg","Reconstructed Frag. Fcn With Background Removed",200,0,2);
321 tempPtNoBg =new TH1F("tempPtNoBg","Reconstructed Frag. Fcn With Background Removed",200,0,200);
322 tempFragmFcnNoBg->Fill(count/(jethighest->Energy()*fScaleFactor),-fhBackHisto->GetBinContent(count));
323 tempPtNoBg->AddBinContent(count,-fhBackHisto->GetBinContent(count));
324 */
325   
326   fhNJets->Fill(fOutput->GetNJets());
327   Bool_t doesJetMeetBinCriteria = 0;
328   AliEMCALJet* jethighest=0; 
329   AliEMCALJet* jetsecond=0; 
330   //  Find Highest and Second Highest Jet 
331   // (NB!!!!!!!)  Pointing into the EMCAL!!!!!!
332   
333 // =========================== All cases ===================================
334  
335
336         // I will make a little array of jet indices for which jets are in 
337         // the EMCAL then my counter can loop from 0 below - but it will 
338         // be the index of the array of applicable jets
339     
340   Int_t appjet[4];
341   Int_t numappjet=0;
342   
343   for (Int_t appc=0;appc<fOutput->GetNJets();appc++)
344   { // Check all jets for applicability
345           Float_t eta = fOutput->GetJet(appc)->Eta();
346           Float_t phi = fOutput->GetJet(appc)->Phi();
347           if (eta > -0.7 && eta < 0.7 && phi > 1./3.*TMath::Pi() && phi < TMath::Pi())
348           { // Then jet is applicable
349                   appjet[numappjet]=appc;                               
350                   numappjet++;
351           }
352   }
353         
354   
355 Float_t et=0;
356 if (numappjet >=1)
357 {
358         Float_t theta = 2.0*atan(exp(-fOutput->GetParton(0)->Eta()));
359         et    = fOutput->GetParton(0)->Energy() * TMath::Sin(theta);
360         if (fOutput->GetNJets()>1)
361         {
362                 for (Int_t counter = 0; counter<numappjet;counter++)  
363                 {         
364                         if (counter==0)
365                         {                 
366                                 jethighest = fOutput->GetJet(appjet[0]);
367                                 jetsecond  = fOutput->GetJet(appjet[1]);
368                         }
369                         if (counter>0)
370                         {
371                                 Float_t energyhighest = jethighest->Energy();
372                                 Float_t energysecond = jetsecond->Energy();
373                                 
374                                 if ((fOutput->GetJet(appjet[counter]))->Energy()>energyhighest) 
375                                 {
376                                         jetsecond=jethighest;
377                                         jethighest=fOutput->GetJet(appjet[counter]);
378                                 }else if ((fOutput->GetJet(appjet[counter]))->Energy()>energysecond) 
379                                 {
380                                         jetsecond=fOutput->GetJet(appjet[counter]);
381                                 }
382               
383                         }
384                 }
385         }else
386         {
387                 Float_t eta = fOutput->GetJet(0)->Eta();
388                 Float_t phi = fOutput->GetJet(0)->Phi();
389                 if (eta > -0.7 && eta < 0.7 && phi > 1./3.*TMath::Pi() && phi < TMath::Pi())
390                 { // Then jet is applicable
391                         jethighest=fOutput->GetJet(0);
392                         jetsecond=0;
393                 }else
394                 {
395                         Error("FillFromOutput","There is only one jet and it isn't in the area of applicability");
396                 }
397                                                 
398         }
399         if ( 95.0 < jethighest->Energy()*fScaleFactor && jethighest->Energy()*fScaleFactor < 105.0  )
400         {
401                 doesJetMeetBinCriteria = 1;
402                 fhRecoBinJetEt->Fill(jethighest->Energy()*fScaleFactor);
403                 fhRecoBinInputJetEt->Fill(et);
404         }
405         fhInputOutput->Fill(et,jethighest->Energy());
406                 
407 }
408
409 if (numappjet > 1)
410 {       
411 //========================= CASE 2 ===========================  
412   Int_t nPartons = fOutput->GetNPartons();
413   fhNJets2->Fill(fOutput->GetNJets());
414   AliEMCALParton* parton;
415
416   // End finding highest and second highest and continue
417   fhJetEt2->Fill(jethighest->Energy()*fScaleFactor);
418   fhJetEtDiff2->Fill(jethighest->Energy()*fScaleFactor-et);
419   fhJetInvE2->Fill(1.0/(jethighest->Energy()*fScaleFactor));
420   fhJetEta2->Fill(jethighest->Eta()  );
421   fhJetPhi2->Fill(jethighest->Phi() );
422   if (nPartons ==0) return;
423   parton = fOutput->GetParton(0);
424   
425   fhPartonEta2->Fill( parton->Eta() );
426   fhPartonPhi2->Fill( parton->Phi() );
427
428   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
429   fhEtaDiff2->Fill( jethighest->Eta() - parton->Eta() );
430   fhPhiDiff2->Fill( jethighest->Phi() - parton->Phi() );
431   fhEtaPhiSpread2->Fill(jethighest->Eta()-parton->Eta(),jethighest->Phi() - parton->Phi());
432   fhJetEtSecond2->Fill(jetsecond->Energy()*fScaleFactor); 
433   fhJetEtRatio2->Fill(jetsecond->Energy()/jethighest->Energy()); 
434   fhEtaPhiDist2->Fill( TMath::Sqrt((jethighest->Eta() - jetsecond->Eta())*(jethighest->Eta() - jetsecond->Eta())
435                       + (jethighest->Phi() - jetsecond->Phi())*(jethighest->Phi() - jetsecond->Phi())     ));  
436   /* 
437   Float_t *pt,*phi,*eta;
438   Int_t *pdg;
439   pt  = new Float_t[parton->GetNTracks()];
440   eta = new Float_t[parton->GetNTracks()];
441   phi = new Float_t[parton->GetNTracks()];
442   pdg = new Int_t[parton->GetNTracks()];*/
443
444
445
446  Float_t pt[2000];
447  Float_t eta[2000];
448  Float_t phi[2000];
449  Int_t pdg[2000];
450   
451   parton->GetTrackList(pt,eta,phi,pdg);
452   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ ) 
453   {
454       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
455            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue; 
456       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
457       Double_t rt = 2.0*atan(exp(-parton->Eta()));
458       Double_t ctt = cos(tt);
459       Double_t crt = cos(rt);
460       Double_t stt = sin(tt);
461       Double_t srt = sin(rt);
462       Double_t ctp = cos(phi[iT]);
463       Double_t crp = cos(parton->Phi());
464       Double_t stp = sin(phi[iT]);
465       Double_t srp = sin(parton->Phi());
466       //Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
467           Double_t alpha;   
468           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
469           {
470                   alpha = 0.0;
471           }else
472           {
473                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
474           }
475       Double_t correctp = pt[iT]/stt;   
476       fhPartonPL2->Fill( correctp*cos(alpha));
477       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +  
478       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 ) 
479               fhPartonJT2->Fill( correctp*sin(alpha));
480       fhPartonPT2->Fill(correctp*sin(tt));
481       if (fNominalEnergy == 0.0) {
482               fhPartonFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
483       }else 
484       {
485               fhPartonFragmFcn2->Fill(correctp*sin(tt)/fNominalEnergy);
486       }           
487       if (doesJetMeetBinCriteria)
488       {
489               fhRecoBinPartonPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
490       }
491   }// loop over tracks
492
493 /*
494   pt  = new Float_t[jet->NTracks()];
495   eta = new Float_t[jet->NTracks()];
496   phi = new Float_t[jet->NTracks()];
497   pdg = new Int_t[jet->NTracks()];*/
498   jethighest->TrackList(pt,eta,phi,pdg);
499   for(Int_t iT=0; iT< jethighest->NTracks() ; iT++ )
500   {
501           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
502           Double_t rt = 2.0*atan(exp(-jethighest->Eta()));
503           Double_t ctt = cos(tt);                   
504           Double_t crt = cos(rt);                         
505           Double_t stt = sin(tt);                               
506           Double_t srt = sin(rt);
507           Double_t ctp = cos(phi[iT]);
508           Double_t crp = cos(jethighest->Phi());
509           Double_t stp = sin(phi[iT]);
510           Double_t srp = sin(jethighest->Phi());
511          // Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
512           Double_t alpha;   
513           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
514           {
515                   alpha = 0.0;
516           }else
517           {
518                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
519           }
520           Double_t correctp = pt[iT]/stt;
521           fhJetPL2->Fill( correctp*cos(alpha));      
522           if ( (jethighest->Eta()-eta[iT])*(jethighest->Eta()-eta[iT]) +  
523                           (jethighest->Phi()-phi[iT])*(jethighest->Phi()-phi[iT]) < 0.2*0.2 )   
524                   fhJetJT2->Fill( correctp*sin(alpha));
525           fhJetPT2->Fill(correctp*sin(tt));
526           if (fNominalEnergy==0.0){
527                   fhFragmFcn2->Fill(  correctp*sin(tt)/(jethighest->Energy()*fScaleFactor)  );
528           } else
529           {
530                   fhFragmFcn2->Fill(  correctp*sin(tt)/fNominalEnergy );
531           }
532           if (doesJetMeetBinCriteria)
533           {
534                   fhRecoBinPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
535                   fhRecoBinPtNoBg->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
536                   fhRecoBinFragmFcn->Fill(  correctp*sin(tt)/(jethighest->Energy()*fScaleFactor)  ); // This is the jet fragmentation function
537                   fhRecoBinFragmFcnNoBg->Fill(  correctp*sin(tt)/(jethighest->Energy()*fScaleFactor)  ); // This is the jet fragmentation function
538           }
539   }// loop over tracks
540   }
541
542   if (numappjet == 1)
543   {
544
545 //========================= CASE 1 ===========================  
546   Int_t nPartons = fOutput->GetNPartons();
547   if (fOutput->GetNJets()!=1) return;
548   AliEMCALParton* parton;
549   AliEMCALJet* jet; 
550   jet = jethighest;//fOutput->GetJet(0);
551   fhJetEt->Fill(jet->Energy()*fScaleFactor);
552   fhJetEtDiff->Fill(jethighest->Energy()*fScaleFactor-et);
553   fhJetInvE->Fill(1.0/(jethighest->Energy()*fScaleFactor));
554   fhJetEta->Fill(jet->Eta()  );
555   fhJetPhi->Fill(jet->Phi() );
556   if (nPartons ==0) return;
557   parton = fOutput->GetParton(0);
558   
559   fhPartonEta->Fill( parton->Eta() );
560   fhPartonPhi->Fill( parton->Phi() );
561
562   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
563   fhEtaDiff->Fill( jet->Eta() - parton->Eta() );
564   fhPhiDiff->Fill( jet->Phi() - parton->Phi() );
565   fhEtaPhiSpread->Fill(jet->Eta()-parton->Eta(),jet->Phi() - parton->Phi());
566  /* 
567   Float_t *pt,*phi,*eta;
568   Int_t *pdg;
569   pt  = new Float_t[parton->GetNTracks()];
570   eta = new Float_t[parton->GetNTracks()];
571   phi = new Float_t[parton->GetNTracks()];
572   pdg = new Int_t[parton->GetNTracks()];*/
573
574
575
576  Float_t pt[2000];
577  Float_t eta[2000];
578  Float_t phi[2000];
579  Int_t pdg[2000];
580   
581   parton->GetTrackList(pt,eta,phi,pdg);
582   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ ) 
583   {
584       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
585            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue; 
586       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
587       Double_t rt = 2.0*atan(exp(-parton->Eta()));
588       Double_t ctt = cos(tt);
589       Double_t crt = cos(rt);
590       Double_t stt = sin(tt);
591       Double_t srt = sin(rt);
592       Double_t ctp = cos(phi[iT]);
593       Double_t crp = cos(parton->Phi());
594       Double_t stp = sin(phi[iT]);
595       Double_t srp = sin(parton->Phi());
596      // Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
597           Double_t alpha;   
598           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
599           {
600                   alpha = 0.0;
601           }else
602           {
603                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);           }
604       Double_t correctp = pt[iT]/stt;   
605       fhPartonPL->Fill( correctp*cos(alpha));
606       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +  
607       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 ) 
608               fhPartonJT->Fill( correctp*sin(alpha));
609       if (fNominalEnergy == 0.0) {
610               fhPartonFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
611           fhPartonPT->Fill(correctp*sin(tt));
612       }else 
613       {
614               fhPartonFragmFcn->Fill(correctp*sin(tt)/fNominalEnergy);
615       }
616       if (doesJetMeetBinCriteria)
617       {
618               fhRecoBinPartonPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
619       }
620   }// loop over tracks
621
622 /*
623   pt  = new Float_t[jet->NTracks()];
624   eta = new Float_t[jet->NTracks()];
625   phi = new Float_t[jet->NTracks()];
626   pdg = new Int_t[jet->NTracks()];*/
627   jet->TrackList(pt,eta,phi,pdg);
628   for(Int_t iT=0; iT< jet->NTracks() ; iT++ )
629   {
630           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
631           Double_t rt = 2.0*atan(exp(-jet->Eta()));
632           Double_t ctt = cos(tt);                   
633           Double_t crt = cos(rt);                         
634           Double_t stt = sin(tt);                               
635           Double_t srt = sin(rt);
636           Double_t ctp = cos(phi[iT]);
637           Double_t crp = cos(jet->Phi());
638           Double_t stp = sin(phi[iT]);
639           Double_t srp = sin(jet->Phi());
640           //Info("plots","acos(%1.16f)\nstt=%f\npt=%f",crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt,stt,pt[iT]);
641           //Info("plots","diff to 1 %f",1.0-crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
642           Double_t alpha;   
643           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
644           {
645                   alpha = 0.0;
646           }else
647           {
648                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
649           }
650           Double_t correctp = pt[iT]/stt;
651           fhJetPL->Fill( correctp*cos(alpha));      
652           if ( (jet->Eta()-eta[iT])*(jet->Eta()-eta[iT]) +  
653                (jet->Phi()-phi[iT])*(jet->Phi()-phi[iT]) < 0.2*0.2 )   
654                   fhJetJT->Fill( correctp*sin(alpha));
655           fhJetPT->Fill(correctp*sin(tt));
656           if (fNominalEnergy==0.0){
657                   fhFragmFcn->Fill(  correctp*sin(tt)/(jet->Energy()*fScaleFactor)  ); // This is the jet fragmentation function
658           } else
659           {
660                   fhFragmFcn->Fill(  correctp*sin(tt)/fNominalEnergy );
661           }
662           if (doesJetMeetBinCriteria)
663           {
664                   fhRecoBinPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
665                   fhRecoBinPtNoBg->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
666                   fhRecoBinFragmFcn->Fill(  correctp*sin(tt)/(jet->Energy()*fScaleFactor)  ); // This is the jet fragmentation function
667                   fhRecoBinFragmFcnNoBg->Fill(  correctp*sin(tt)/(jet->Energy()*fScaleFactor)  ); // This is the jet fragmentation function
668           }
669   }// loop over tracks
670   }
671
672 if (numappjet>=1 && fhBackHisto != 0 && doesJetMeetBinCriteria)
673   {
674     for (Int_t count=1;count<=100;count++)
675     {
676         fhRecoBinFragmFcnNoBg->Fill( ((Float_t)count)/(jethighest->Energy()*fScaleFactor),-fhBackHisto->GetBinContent(count));
677         fhRecoBinPtNoBg->AddBinContent(count,-fhBackHisto->GetBinContent(count));
678     } 
679   }
680
681
682 }
683
684