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