6a902249dae889afb5a38ea9cfc3b0c7ef06466a
[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   
270 // =========================== All cases ===================================
271   
272 if (fOutput->GetNJets()>=1)
273 {
274         Float_t theta = 2.0*atan(exp(-fOutput->GetParton(0)->Eta()));
275         Float_t et    = fOutput->GetParton(0)->Energy() * TMath::Sin(theta);
276         fhInputOutput->Fill(et,fOutput->GetJet(0)->Energy());
277         if (fOutput->GetNJets()>1)
278         {
279                 for (Int_t counter = 0; counter<fOutput->GetNJets();counter++)  
280                 {         
281                         if (counter==0)
282                         {                 
283                                 jethighest = fOutput->GetJet(0);
284                                 jetsecond  = fOutput->GetJet(1);
285                         }
286                         if (counter>0)
287                         {
288                                 Float_t energyhighest = jethighest->Energy();
289                                 Float_t energysecond = jetsecond->Energy();
290                                 
291                                 if ((fOutput->GetJet(counter))->Energy()>energyhighest) 
292                                 {
293                                         jetsecond=jethighest;
294                                         jethighest=fOutput->GetJet(counter);
295                                 }else if ((fOutput->GetJet(counter))->Energy()>energysecond) 
296                                 {
297                                         jetsecond=fOutput->GetJet(counter);
298                                 }
299               
300                         }
301                 }
302         }else
303         {
304                 jethighest=fOutput->GetJet(0);
305                 jetsecond=0;
306         }
307         if ( 95.0 < jethighest->Energy() && jethighest->Energy() < 105.0  )
308         {
309                 doesJetMeetBinCriteria = 1;
310                 fhRecoBinJetEt->Fill(jethighest->Energy());
311                 fhRecoBinInputJetEt->Fill(et);
312         }
313                 
314 }
315
316 if (fOutput->GetNJets()>1)
317 {       
318 //========================= CASE 2 ===========================  
319   Int_t nPartons = fOutput->GetNPartons();
320   fhNJets2->Fill(fOutput->GetNJets());
321   AliEMCALParton* parton;
322
323   // End finding highest and second highest and continue
324   fhJetEt2->Fill(jethighest->Energy());
325   fhJetEta2->Fill(jethighest->Eta()  );
326   fhJetPhi2->Fill(jethighest->Phi() );
327   if (nPartons ==0) return;
328   parton = fOutput->GetParton(0);
329   
330   fhPartonEta2->Fill( parton->Eta() );
331   fhPartonPhi2->Fill( parton->Phi() );
332
333   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
334   fhEtaDiff2->Fill( jethighest->Eta() - parton->Eta() );
335   fhPhiDiff2->Fill( jethighest->Phi() - parton->Phi() );
336   fhEtaPhiSpread2->Fill(jethighest->Eta()-parton->Eta(),jethighest->Phi() - parton->Phi());
337   fhJetEtSecond2->Fill(jetsecond->Energy()); 
338   fhJetEtRatio2->Fill(jetsecond->Energy()/jethighest->Energy()); 
339   fhEtaPhiDist2->Fill( TMath::Sqrt((jethighest->Eta() - jetsecond->Eta())*(jethighest->Eta() - jetsecond->Eta())
340                       + (jethighest->Phi() - jetsecond->Phi())*(jethighest->Phi() - jetsecond->Phi())     ));  
341   /* 
342   Float_t *pt,*phi,*eta;
343   Int_t *pdg;
344   pt  = new Float_t[parton->GetNTracks()];
345   eta = new Float_t[parton->GetNTracks()];
346   phi = new Float_t[parton->GetNTracks()];
347   pdg = new Int_t[parton->GetNTracks()];*/
348
349
350
351  Float_t pt[2000];
352  Float_t eta[2000];
353  Float_t phi[2000];
354  Int_t pdg[2000];
355   
356   parton->GetTrackList(pt,eta,phi,pdg);
357   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ ) 
358   {
359       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
360            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue; 
361       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
362       Double_t rt = 2.0*atan(exp(-parton->Eta()));
363       Double_t ctt = cos(tt);
364       Double_t crt = cos(rt);
365       Double_t stt = sin(tt);
366       Double_t srt = sin(rt);
367       Double_t ctp = cos(phi[iT]);
368       Double_t crp = cos(parton->Phi());
369       Double_t stp = sin(phi[iT]);
370       Double_t srp = sin(parton->Phi());
371       //Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
372           Double_t alpha;   
373           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
374           {
375                   alpha = 0.0;
376           }else
377           {
378                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
379           }
380       Double_t correctp = pt[iT]/stt;   
381       fhPartonPL2->Fill( correctp*cos(alpha));
382       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +  
383       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 ) 
384               fhPartonJT2->Fill( correctp*sin(alpha));
385       if (fNominalEnergy == 0.0) {
386               fhPartonFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
387       }else 
388       {
389               fhPartonFragmFcn2->Fill(correctp*sin(tt)/fNominalEnergy);
390       }           
391       if (doesJetMeetBinCriteria)
392       {
393               fhRecoBinPartonPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
394       }
395   }// loop over tracks
396
397 /*
398   pt  = new Float_t[jet->NTracks()];
399   eta = new Float_t[jet->NTracks()];
400   phi = new Float_t[jet->NTracks()];
401   pdg = new Int_t[jet->NTracks()];*/
402   jethighest->TrackList(pt,eta,phi,pdg);
403   for(Int_t iT=0; iT< jethighest->NTracks() ; iT++ )
404   {
405           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
406           Double_t rt = 2.0*atan(exp(-jethighest->Eta()));
407           Double_t ctt = cos(tt);                   
408           Double_t crt = cos(rt);                         
409           Double_t stt = sin(tt);                               
410           Double_t srt = sin(rt);
411           Double_t ctp = cos(phi[iT]);
412           Double_t crp = cos(jethighest->Phi());
413           Double_t stp = sin(phi[iT]);
414           Double_t srp = sin(jethighest->Phi());
415          // Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
416           Double_t alpha;   
417           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
418           {
419                   alpha = 0.0;
420           }else
421           {
422                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
423           }
424           Double_t correctp = pt[iT]/stt;
425           fhJetPL2->Fill( correctp*cos(alpha));      
426           if ( (jethighest->Eta()-eta[iT])*(jethighest->Eta()-eta[iT]) +  
427                           (jethighest->Phi()-phi[iT])*(jethighest->Phi()-phi[iT]) < 0.2*0.2 )   
428                   fhJetJT2->Fill( correctp*sin(alpha));
429           if (fNominalEnergy==0.0){
430                   fhFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
431           } else
432           {
433                   fhFragmFcn2->Fill(  correctp*sin(tt)/fNominalEnergy );
434           }
435           if (doesJetMeetBinCriteria)
436           {
437                   fhRecoBinPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
438           }
439   }// loop over tracks
440   }
441
442   if (fOutput->GetNJets()==1)
443   {
444
445 //========================= CASE 1 ===========================  
446   Int_t nPartons = fOutput->GetNPartons();
447   if (fOutput->GetNJets()!=1) return;
448   AliEMCALParton* parton;
449   AliEMCALJet* jet; 
450   jet = fOutput->GetJet(0);
451   fhJetEt->Fill(jet->Energy());
452   fhJetEta->Fill(jet->Eta()  );
453   fhJetPhi->Fill(jet->Phi() );
454   if (nPartons ==0) return;
455   parton = fOutput->GetParton(0);
456   
457   fhPartonEta->Fill( parton->Eta() );
458   fhPartonPhi->Fill( parton->Phi() );
459
460   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
461   fhEtaDiff->Fill( jet->Eta() - parton->Eta() );
462   fhPhiDiff->Fill( jet->Phi() - parton->Phi() );
463   fhEtaPhiSpread->Fill(jet->Eta()-parton->Eta(),jet->Phi() - parton->Phi());
464  /* 
465   Float_t *pt,*phi,*eta;
466   Int_t *pdg;
467   pt  = new Float_t[parton->GetNTracks()];
468   eta = new Float_t[parton->GetNTracks()];
469   phi = new Float_t[parton->GetNTracks()];
470   pdg = new Int_t[parton->GetNTracks()];*/
471
472
473
474  Float_t pt[2000];
475  Float_t eta[2000];
476  Float_t phi[2000];
477  Int_t pdg[2000];
478   
479   parton->GetTrackList(pt,eta,phi,pdg);
480   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ ) 
481   {
482       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
483            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue; 
484       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
485       Double_t rt = 2.0*atan(exp(-parton->Eta()));
486       Double_t ctt = cos(tt);
487       Double_t crt = cos(rt);
488       Double_t stt = sin(tt);
489       Double_t srt = sin(rt);
490       Double_t ctp = cos(phi[iT]);
491       Double_t crp = cos(parton->Phi());
492       Double_t stp = sin(phi[iT]);
493       Double_t srp = sin(parton->Phi());
494      // Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
495           Double_t alpha;   
496           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
497           {
498                   alpha = 0.0;
499           }else
500           {
501                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);           }
502       Double_t correctp = pt[iT]/stt;   
503       fhPartonPL->Fill( correctp*cos(alpha));
504       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +  
505       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 ) 
506               fhPartonJT->Fill( correctp*sin(alpha));
507       if (fNominalEnergy == 0.0) {
508               fhPartonFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
509       }else 
510       {
511               fhPartonFragmFcn->Fill(correctp*sin(tt)/fNominalEnergy);
512       }
513       if (doesJetMeetBinCriteria)
514       {
515               fhRecoBinPartonPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
516       }
517   }// loop over tracks
518
519 /*
520   pt  = new Float_t[jet->NTracks()];
521   eta = new Float_t[jet->NTracks()];
522   phi = new Float_t[jet->NTracks()];
523   pdg = new Int_t[jet->NTracks()];*/
524   jet->TrackList(pt,eta,phi,pdg);
525   for(Int_t iT=0; iT< jet->NTracks() ; iT++ )
526   {
527           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
528           Double_t rt = 2.0*atan(exp(-jet->Eta()));
529           Double_t ctt = cos(tt);                   
530           Double_t crt = cos(rt);                         
531           Double_t stt = sin(tt);                               
532           Double_t srt = sin(rt);
533           Double_t ctp = cos(phi[iT]);
534           Double_t crp = cos(jet->Phi());
535           Double_t stp = sin(phi[iT]);
536           Double_t srp = sin(jet->Phi());
537           //Info("plots","acos(%1.16f)\nstt=%f\npt=%f",crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt,stt,pt[iT]);
538           //Info("plots","diff to 1 %f",1.0-crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
539           Double_t alpha;   
540           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
541           {
542                   alpha = 0.0;
543           }else
544           {
545                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
546           }
547           Double_t correctp = pt[iT]/stt;
548           fhJetPL->Fill( correctp*cos(alpha));      
549           if ( (jet->Eta()-eta[iT])*(jet->Eta()-eta[iT]) +  
550                           (jet->Phi()-phi[iT])*(jet->Phi()-phi[iT]) < 0.2*0.2 )   
551                   fhJetJT->Fill( correctp*sin(alpha));
552           if (fNominalEnergy==0.0){
553                   fhFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
554           } else
555           {
556                   fhFragmFcn->Fill(  correctp*sin(tt)/fNominalEnergy );
557           }
558           if (doesJetMeetBinCriteria)
559           {
560                   fhRecoBinPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
561           }
562   }// loop over tracks
563   }
564 }
565
566