Made more robust
[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
77 }
78
79 void AliEMCALJetFinderPlots::InitPlots()
80 {
81 //========================= CASE 1 =======================================      
82     fhFragmFcn = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
83     fhFragmFcn->Sumw2();
84     fhPartonFragmFcn = new TH1F("hPartonFragmFcn","Parton Fragmentation Function",100,0,1);
85     fhPartonFragmFcn->Sumw2();
86     fhPartonJT = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
87     fhPartonJT->Sumw2();
88     fhPartonPL = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
89     fhPartonPL->Sumw2();
90     fhJetJT = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
91     fhJetJT->Sumw2();
92     fhJetPL = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
93     fhJetPL->Sumw2();
94     fhJetEt = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
95     fhJetEt->Sumw2();
96     fhJetEta = new      TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
97     fhJetEta->Sumw2();
98     fhJetPhi = new      TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
99     fhJetPhi->Sumw2();
100     fhPartonEta = new   TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
101     fhPartonEta->Sumw2();
102     fhPartonPhi = new   TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
103     fhPartonPhi->Sumw2();
104     fhEtaDiff = new     TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
105     fhEtaDiff->Sumw2();
106     fhPhiDiff  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
107     fhPhiDiff->Sumw2();
108     fhNJets = new       TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
109     fhNJets->Sumw2();
110     fhEtaPhiSpread = new TH2F("hEtaPhiSpread","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
111     fhEtaPhiSpread->Sumw2();    
112   fhNJets->SetXTitle("N_{jets}^{reco}/event");
113   fhNJets->SetYTitle("N_{events}");
114
115   //Jet properties
116   fhJetEt->SetFillColor(16);
117   fhJetEt->SetXTitle("E_{T}^{reco}");
118   
119   fhJetEta->SetFillColor(16);
120   fhJetEta->SetXTitle("#eta_{jet}^{reco}");
121   
122   fhJetPhi->SetFillColor(16);
123   fhJetPhi->SetXTitle("#phi_{jet}^{reco}");
124   
125   fhPartonEta->SetFillColor(16);
126   fhPartonEta->SetXTitle("#eta_{parton}");
127   
128   fhPartonPhi->SetFillColor(16);
129   fhPartonPhi->SetXTitle("#phi_{parton}");
130
131   fhPartonPL->SetXTitle("p (GeV/c)");
132   fhPartonJT->SetXTitle("p (GeV/c)");
133
134   fhPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{parton}");
135
136   //Jet component properties
137
138   fhJetPL->SetXTitle("p (GeV/c)");
139   fhJetJT->SetXTitle("p (GeV/c)");
140   fhFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
141   fhPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
142
143   fhEtaDiff->SetXTitle("#eta_{jet}^{reco}-#eta_{jet}^{input}");
144   fhPhiDiff->SetXTitle("#phi_{jet}^{reco}-#phi_{jet}^{input}");
145   fhEtaPhiSpread->SetXTitle("#eta");
146   fhEtaPhiSpread->SetYTitle("#phi");
147
148 //======================= CASE 2 ======================================
149   
150
151 fhFragmFcn2             = new TH1F("hFragmFcn2","Fragmentation Function",100,0,1);
152 fhFragmFcn2->Sumw2();
153 fhPartonFragmFcn2       = new TH1F("hPartonFragmFcn2","Parton Fragmentation Function",100,0,1);
154 fhPartonFragmFcn2->Sumw2();
155 fhPartonJT2             = new TH1F("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
156 fhPartonJT2->Sumw2();
157 fhPartonPL2             = new TH1F("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
158 fhPartonPL2->Sumw2();
159 fhJetJT2                        = new TH1F("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
160 fhJetJT2->Sumw2();
161 fhJetPL2                        = new TH1F("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
162 fhJetPL2->Sumw2();
163 fhJetEt2                        = new TH1F("hJetEt2","E_{T}^{reco}",250,0.,250.);
164 fhJetEt2->Sumw2();
165 fhJetEta2               = new TH1F("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
166 fhJetEta2->Sumw2();
167 fhJetPhi2               = new TH1F("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
168 fhJetPhi2->Sumw2();
169 fhPartonEta2            = new TH1F("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
170 fhPartonEta2->Sumw2();
171 fhPartonPhi2            = new TH1F("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
172 fhPartonPhi2->Sumw2();
173 fhEtaDiff2              = new TH1F("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
174 fhEtaDiff2->Sumw2();
175 fhPhiDiff2              = new TH1F("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
176 fhPhiDiff2->Sumw2();
177 fhEtaPhiSpread2 = new TH2F("hEtaPhiSpread2","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
178 fhEtaPhiSpread2->Sumw2();
179 fhNJets2                        = new TH1F("hNJets2","N Reconstructed jets",11,-0.5,10.5);
180 fhNJets2->Sumw2();
181 fhJetEtSecond2          = new TH1F("hJetEtSecond2","E_{T}^{reco}",250,0.,250.); 
182 fhJetEtSecond2->Sumw2();
183 fhJetEtRatio2           = new TH1F("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
184 fhJetEtRatio2->Sumw2();
185 fhEtaPhiDist2           = new TH1F("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
186 fhEtaPhiDist2->Sumw2();
187   
188   fInitialised = kTRUE; 
189   
190 }
191
192 AliEMCALJetFinderPlots::~AliEMCALJetFinderPlots()
193 {
194         // To ensure that all requested memory is returned
195 delete    fhFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
196 delete    fhPartonFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
197 delete   fhPartonJT;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
198 delete    fhPartonPL;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
199 delete    fhJetJT;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
200 delete    fhJetPL;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
201 delete    fhJetEt;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
202 delete    fhJetEta;// = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
203 delete    fhJetPhi;// = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
204 delete    fhPartonEta;// = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
205 delete    fhPartonPhi;// = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
206 delete    fhEtaDiff;// = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
207 delete    fhPhiDiff;//  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
208 delete    fhNJets;// = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
209 delete    fhEtaPhiSpread;           
210
211         delete                          fhFragmFcn2;    // ("hFragmFcn2","Fragmentation Function",100,0,1);
212         delete                          fhPartonFragmFcn2;// ("hFragmFcn2","Parton Fragmentation Function",100,0,1);
213         delete                          fhPartonJT2;    // ("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
214         delete                          fhPartonPL2;    // ("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
215         delete                          fhJetJT2;       // ("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
216         delete                          fhJetPL2;       // ("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
217         delete                          fhJetEt2;       // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
218         delete                          fhJetEta2;      // ("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
219         delete                          fhJetPhi2;      // ("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
220         delete                          fhPartonEta2;   // ("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
221         delete                          fhPartonPhi2;   // ("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
222         delete                          fhEtaDiff2;     // ("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
223         delete                          fhPhiDiff2;     // ("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
224         delete                          fhEtaPhiSpread2;        // ("hEtaPhiSpread2","#eta - #phi Distribution 
225                                                         //of Reconstructed Jets",192,-0.7,0.7,288,pi/3,pi);
226         delete                          fhNJets2;       // ("hNJets2","N Reconstructed jets",11,-0.5,10.5);
227         delete                          fhJetEtSecond2; //("hJetEtSecond2","E_{T}^{reco}",250,0.,250.); 
228         delete                          fhJetEtRatio2;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
229         delete                          fhEtaPhiDist2;  //("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
230
231
232
233 }       
234
235 void AliEMCALJetFinderPlots::FillFromOutput(AliEMCALJetFinderOutput* output)
236 {
237         // Fill histograms from an output object
238 if (!fInitialised) InitPlots(); 
239   fOutput = output;
240   if (!fOutput) return;
241   fhNJets->Fill(fOutput->GetNJets());
242 if (fOutput->GetNJets()>1)
243 {       
244 //========================= CASE 2 ===========================  
245   Int_t nPartons = fOutput->GetNPartons();
246   fhNJets2->Fill(fOutput->GetNJets());
247   AliEMCALParton* parton;
248   AliEMCALJet* jethighest=0; 
249   AliEMCALJet* jetsecond=0; 
250   //  Find Highest and Second Highest Jet
251   for (Int_t counter = 0; counter<fOutput->GetNJets();counter++)
252   {
253           if (counter==0){ 
254                   jethighest = fOutput->GetJet(0);
255                   jetsecond  = fOutput->GetJet(1);
256           }
257           if (counter>0)
258           {
259                   Float_t energyhighest = jethighest->Energy();
260                   Float_t energysecond = jetsecond->Energy();
261                   
262                   if ((fOutput->GetJet(counter))->Energy()>energyhighest) 
263                   {
264                           jetsecond=jethighest;
265                           jethighest=fOutput->GetJet(counter);
266                   }else if ((fOutput->GetJet(counter))->Energy()>energysecond) 
267                   {
268                           jetsecond=fOutput->GetJet(counter);
269                   }
270           }
271   }
272
273   // End finding highest and second highest and continue
274   fhJetEt2->Fill(jethighest->Energy());
275   fhJetEta2->Fill(jethighest->Eta()  );
276   fhJetPhi2->Fill(jethighest->Phi() );
277   if (nPartons ==0) return;
278   parton = fOutput->GetParton(0);
279   
280   fhPartonEta2->Fill( parton->Eta() );
281   fhPartonPhi2->Fill( parton->Phi() );
282
283   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
284   fhEtaDiff2->Fill( jethighest->Eta() - parton->Eta() );
285   fhPhiDiff2->Fill( jethighest->Phi() - parton->Phi() );
286   fhEtaPhiSpread2->Fill(jethighest->Eta()-parton->Eta(),jethighest->Phi() - parton->Phi());
287   fhJetEtSecond2->Fill(jetsecond->Energy()); 
288   fhJetEtRatio2->Fill(jetsecond->Energy()/jethighest->Energy()); 
289   fhEtaPhiDist2->Fill( TMath::Sqrt((jethighest->Eta() - jetsecond->Eta())*(jethighest->Eta() - jetsecond->Eta())
290                       + (jethighest->Phi() - jetsecond->Phi())*(jethighest->Phi() - jetsecond->Phi())     ));  
291   /* 
292   Float_t *pt,*phi,*eta;
293   Int_t *pdg;
294   pt  = new Float_t[parton->GetNTracks()];
295   eta = new Float_t[parton->GetNTracks()];
296   phi = new Float_t[parton->GetNTracks()];
297   pdg = new Int_t[parton->GetNTracks()];*/
298
299
300
301  Float_t pt[2000];
302  Float_t eta[2000];
303  Float_t phi[2000];
304  Int_t pdg[2000];
305   
306   parton->GetTrackList(pt,eta,phi,pdg);
307   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ ) 
308   {
309       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
310            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue; 
311       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
312       Double_t rt = 2.0*atan(exp(-parton->Eta()));
313       Double_t ctt = cos(tt);
314       Double_t crt = cos(rt);
315       Double_t stt = sin(tt);
316       Double_t srt = sin(rt);
317       Double_t ctp = cos(phi[iT]);
318       Double_t crp = cos(parton->Phi());
319       Double_t stp = sin(phi[iT]);
320       Double_t srp = sin(parton->Phi());
321       Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
322       Double_t correctp = pt[iT]/stt;   
323       fhPartonPL2->Fill( correctp*cos(alpha));
324       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +  
325       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 ) 
326               fhPartonJT2->Fill( correctp*sin(alpha));
327       if (fNominalEnergy == 0.0) {
328               fhPartonFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
329       }else 
330       {
331               fhPartonFragmFcn2->Fill(correctp*sin(tt)/fNominalEnergy);
332       }
333   }// loop over tracks
334
335 /*
336   pt  = new Float_t[jet->NTracks()];
337   eta = new Float_t[jet->NTracks()];
338   phi = new Float_t[jet->NTracks()];
339   pdg = new Int_t[jet->NTracks()];*/
340   jethighest->TrackList(pt,eta,phi,pdg);
341   for(Int_t iT=0; iT< jethighest->NTracks() ; iT++ )
342   {
343           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
344           Double_t rt = 2.0*atan(exp(-jethighest->Eta()));
345           Double_t ctt = cos(tt);                   
346           Double_t crt = cos(rt);                         
347           Double_t stt = sin(tt);                               
348           Double_t srt = sin(rt);
349           Double_t ctp = cos(phi[iT]);
350           Double_t crp = cos(jethighest->Phi());
351           Double_t stp = sin(phi[iT]);
352           Double_t srp = sin(jethighest->Phi());
353           Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
354           Double_t correctp = pt[iT]/stt;
355           fhJetPL2->Fill( correctp*cos(alpha));      
356           if ( (jethighest->Eta()-eta[iT])*(jethighest->Eta()-eta[iT]) +  
357                           (jethighest->Phi()-phi[iT])*(jethighest->Phi()-phi[iT]) < 0.2*0.2 )   
358                   fhJetJT2->Fill( correctp*sin(alpha));
359           if (fNominalEnergy==0.0){
360                   fhFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
361           } else
362           {
363                   fhFragmFcn2->Fill(  correctp*sin(tt)/fNominalEnergy );
364           }
365   }// loop over tracks
366   }
367
368   if (fOutput->GetNJets()==1)
369   {
370
371 //========================= CASE 1 ===========================  
372   Int_t nPartons = fOutput->GetNPartons();
373   if (fOutput->GetNJets()!=1) return;
374   AliEMCALParton* parton;
375   AliEMCALJet* jet; 
376   jet = fOutput->GetJet(0);
377   fhJetEt->Fill(jet->Energy());
378   fhJetEta->Fill(jet->Eta()  );
379   fhJetPhi->Fill(jet->Phi() );
380   if (nPartons ==0) return;
381   parton = fOutput->GetParton(0);
382   
383   fhPartonEta->Fill( parton->Eta() );
384   fhPartonPhi->Fill( parton->Phi() );
385
386   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
387   fhEtaDiff->Fill( jet->Eta() - parton->Eta() );
388   fhPhiDiff->Fill( jet->Phi() - parton->Phi() );
389   fhEtaPhiSpread->Fill(jet->Eta()-parton->Eta(),jet->Phi() - parton->Phi());
390  /* 
391   Float_t *pt,*phi,*eta;
392   Int_t *pdg;
393   pt  = new Float_t[parton->GetNTracks()];
394   eta = new Float_t[parton->GetNTracks()];
395   phi = new Float_t[parton->GetNTracks()];
396   pdg = new Int_t[parton->GetNTracks()];*/
397
398
399
400  Float_t pt[2000];
401  Float_t eta[2000];
402  Float_t phi[2000];
403  Int_t pdg[2000];
404   
405   parton->GetTrackList(pt,eta,phi,pdg);
406   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ ) 
407   {
408       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
409            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue; 
410       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
411       Double_t rt = 2.0*atan(exp(-parton->Eta()));
412       Double_t ctt = cos(tt);
413       Double_t crt = cos(rt);
414       Double_t stt = sin(tt);
415       Double_t srt = sin(rt);
416       Double_t ctp = cos(phi[iT]);
417       Double_t crp = cos(parton->Phi());
418       Double_t stp = sin(phi[iT]);
419       Double_t srp = sin(parton->Phi());
420       Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
421       Double_t correctp = pt[iT]/stt;   
422       fhPartonPL->Fill( correctp*cos(alpha));
423       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +  
424       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 ) 
425               fhPartonJT->Fill( correctp*sin(alpha));
426       if (fNominalEnergy == 0.0) {
427               fhPartonFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
428       }else 
429       {
430               fhPartonFragmFcn->Fill(correctp*sin(tt)/fNominalEnergy);
431       }
432   }// loop over tracks
433
434 /*
435   pt  = new Float_t[jet->NTracks()];
436   eta = new Float_t[jet->NTracks()];
437   phi = new Float_t[jet->NTracks()];
438   pdg = new Int_t[jet->NTracks()];*/
439   jet->TrackList(pt,eta,phi,pdg);
440   for(Int_t iT=0; iT< jet->NTracks() ; iT++ )
441   {
442           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
443           Double_t rt = 2.0*atan(exp(-jet->Eta()));
444           Double_t ctt = cos(tt);                   
445           Double_t crt = cos(rt);                         
446           Double_t stt = sin(tt);                               
447           Double_t srt = sin(rt);
448           Double_t ctp = cos(phi[iT]);
449           Double_t crp = cos(jet->Phi());
450           Double_t stp = sin(phi[iT]);
451           Double_t srp = sin(jet->Phi());
452           Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);   
453           Double_t correctp = pt[iT]/stt;
454           fhJetPL->Fill( correctp*cos(alpha));      
455           if ( (jet->Eta()-eta[iT])*(jet->Eta()-eta[iT]) +  
456                           (jet->Phi()-phi[iT])*(jet->Phi()-phi[iT]) < 0.2*0.2 )   
457                   fhJetJT->Fill( correctp*sin(alpha));
458           if (fNominalEnergy==0.0){
459                   fhFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
460           } else
461           {
462                   fhFragmFcn->Fill(  correctp*sin(tt)/fNominalEnergy );
463           }
464   }// loop over tracks
465   }
466 }
467
468