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