EMCAL/AliEMCALJetFinderPlots.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17
  18 /* $Id$ */
  19
  20
  21 //_________________________________________________________________________
  22 //  Class for Filling JetFinder Plots
  23 //
  24 //*-- Author: Mark Horner (LBL/UCT)
  25 //
  26 //
  27
  28
  29 #include "TMath.h"
  30 #include "AliEMCALJetFinderPlots.h"
  31
  32 ClassImp(AliEMCALJetFinderPlots)
  33
  34 AliEMCALJetFinderPlots::AliEMCALJetFinderPlots()
  35 {
  36   fInitialised = kFALSE;
  37   fNominalEnergy = 0.0;
  38   fConeRadius = 0.3;
  39   fDebug = 0;
  40   fOutput=0;
  41     hFragmFcn=0;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
  42     hPartonFragmFcn=0;// = new TH1F("hPartonFragmFcn","Fragmentation Function",100,0,1);
  43    hPartonJT=0;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
  44     hPartonPL=0;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
  45     hJetJT=0;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
  46     hJetPL=0;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
  47     hJetEt=0;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
  48     hJetEta=0;// = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
  49     hJetPhi=0;// = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
  50     hPartonEta=0;// = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
  51     hPartonPhi=0;// = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
  52     hEtaDiff=0;// = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
  53     hPhiDiff=0;//  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
  54     hNJets=0;// = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
  55   hEtaPhiSpread=0;
  56
  57 hFragmFcn2=0;   // ("hFragmFcn2","Fragmentation Function",100,0,1);
  58 hPartonFragmFcn2=0;// ("hFragmFcn2","Parton Fragmentation Function",100,0,1);
  59 hPartonJT2=0;   // ("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
  60 hPartonPL2=0;   // ("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
  61 hJetJT2=0;      // ("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
  62 hJetPL2=0;      // ("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
  63 hJetEt2=0;      // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
  64 hJetEta2=0;     // ("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
  65 hJetPhi2=0;     // ("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
  66 hPartonEta2=0;  // ("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
  67 hPartonPhi2=0;  // ("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
  68 hEtaDiff2=0;    // ("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
  69 hPhiDiff2=0;    // ("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
  70 hEtaPhiSpread2=0;       // ("hEtaPhiSpread2","#eta - #phi Distribution
  71                                                 //of Reconstructed Jets",192,-0.7,0.7,288,pi/3,pi);
  72 hNJets2=0;      // ("hNJets2","N Reconstructed jets",11,-0.5,10.5);
  73 hJetEtSecond2=0; //("hJetEtSecond2","E_{T}^{reco}",250,0.,250.);
  74 hJetEtRatio2=0;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
  75 hEtaPhiDist2=0;  //("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
  76
  77 }
  78
  79 void AliEMCALJetFinderPlots::InitPlots()
  80 {
  81 //========================= CASE 1 =======================================
  82     hFragmFcn = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
  83     hFragmFcn->Sumw2();
  84     hPartonFragmFcn = new TH1F("hPartonFragmFcn","Parton Fragmentation Function",100,0,1);
  85     hPartonFragmFcn->Sumw2();
  86     hPartonJT = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
  87     hPartonJT->Sumw2();
  88     hPartonPL = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
  89     hPartonPL->Sumw2();
  90     hJetJT = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
  91     hJetJT->Sumw2();
  92     hJetPL = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
  93     hJetPL->Sumw2();
  94     hJetEt = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
  95     hJetEt->Sumw2();
  96     hJetEta = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
  97     hJetEta->Sumw2();
  98     hJetPhi = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
  99     hJetPhi->Sumw2();
 100     hPartonEta = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
 101     hPartonEta->Sumw2();
 102     hPartonPhi = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
 103     hPartonPhi->Sumw2();
 104     hEtaDiff = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 105     hEtaDiff->Sumw2();
 106     hPhiDiff  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 107     hPhiDiff->Sumw2();
 108     hNJets = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
 109     hNJets->Sumw2();
 110     hEtaPhiSpread = new TH2F("hEtaPhiSpread","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
 111     hEtaPhiSpread->Sumw2();
 112   hNJets->SetXTitle("N_{jets}^{reco}/event");
 113   hNJets->SetYTitle("N_{events}");
 114
 115   //Jet properties
 116   hJetEt->SetFillColor(16);
 117   hJetEt->SetXTitle("E_{T}^{reco}");
 118
 119   hJetEta->SetFillColor(16);
 120   hJetEta->SetXTitle("#eta_{jet}^{reco}");
 121
 122   hJetPhi->SetFillColor(16);
 123   hJetPhi->SetXTitle("#phi_{jet}^{reco}");
 124
 125   hPartonEta->SetFillColor(16);
 126   hPartonEta->SetXTitle("#eta_{parton}");
 127
 128   hPartonPhi->SetFillColor(16);
 129   hPartonPhi->SetXTitle("#phi_{parton}");
 130
 131   hPartonPL->SetXTitle("p (GeV/c)");
 132   hPartonJT->SetXTitle("p (GeV/c)");
 133
 134   hPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{parton}");
 135
 136   //Jet component properties
 137
 138   hJetPL->SetXTitle("p (GeV/c)");
 139   hJetJT->SetXTitle("p (GeV/c)");
 140   hFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
 141   hPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
 142
 143   hEtaDiff->SetXTitle("#eta_{jet}^{reco}-#eta_{jet}^{input}");
 144   hPhiDiff->SetXTitle("#phi_{jet}^{reco}-#phi_{jet}^{input}");
 145   hEtaPhiSpread->SetXTitle("#eta");
 146   hEtaPhiSpread->SetYTitle("#phi");
 147
 148 //======================= CASE 2 ======================================
 149
 150
 151 hFragmFcn2              = new TH1F("hFragmFcn2","Fragmentation Function",100,0,1);
 152 hFragmFcn2->Sumw2();
 153 hPartonFragmFcn2        = new TH1F("hPartonFragmFcn2","Parton Fragmentation Function",100,0,1);
 154 hPartonFragmFcn2->Sumw2();
 155 hPartonJT2              = new TH1F("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 156 hPartonJT2->Sumw2();
 157 hPartonPL2              = new TH1F("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 158 hPartonPL2->Sumw2();
 159 hJetJT2                 = new TH1F("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 160 hJetJT2->Sumw2();
 161 hJetPL2                 = new TH1F("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 162 hJetPL2->Sumw2();
 163 hJetEt2                 = new TH1F("hJetEt2","E_{T}^{reco}",250,0.,250.);
 164 hJetEt2->Sumw2();
 165 hJetEta2                = new TH1F("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
 166 hJetEta2->Sumw2();
 167 hJetPhi2                = new TH1F("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
 168 hJetPhi2->Sumw2();
 169 hPartonEta2             = new TH1F("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
 170 hPartonEta2->Sumw2();
 171 hPartonPhi2             = new TH1F("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
 172 hPartonPhi2->Sumw2();
 173 hEtaDiff2               = new TH1F("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 174 hEtaDiff2->Sumw2();
 175 hPhiDiff2               = new TH1F("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 176 hPhiDiff2->Sumw2();
 177 hEtaPhiSpread2          = new TH2F("hEtaPhiSpread2","#eta - #phi Distribution of Reconstructed Jets",192,-0.7,0.7,288,0.0,2.0/3.0*TMath::Pi());
 178 hEtaPhiSpread2->Sumw2();
 179 hNJets2                 = new TH1F("hNJets2","N Reconstructed jets",11,-0.5,10.5);
 180 hNJets2->Sumw2();
 181 hJetEtSecond2           = new TH1F("hJetEtSecond2","E_{T}^{reco}",250,0.,250.);
 182 hJetEtSecond2->Sumw2();
 183 hJetEtRatio2            = new TH1F("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
 184 hJetEtRatio2->Sumw2();
 185 hEtaPhiDist2            = new TH1F("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
 186 hEtaPhiDist2->Sumw2();
 187
 188   fInitialised = kTRUE;
 189
 190 }
 191
 192 AliEMCALJetFinderPlots::~AliEMCALJetFinderPlots(){
 193
 194 delete    hFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
 195 delete    hPartonFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
 196 delete   hPartonJT;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 197 delete    hPartonPL;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 198 delete    hJetJT;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 199 delete    hJetPL;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 200 delete    hJetEt;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
 201 delete    hJetEta;// = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
 202 delete    hJetPhi;// = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
 203 delete    hPartonEta;// = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
 204 delete    hPartonPhi;// = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
 205 delete    hEtaDiff;// = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 206 delete    hPhiDiff;//  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 207 delete    hNJets;// = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
 208 delete    hEtaPhiSpread;
 209
 210         delete                          hFragmFcn2;     // ("hFragmFcn2","Fragmentation Function",100,0,1);
 211         delete                          hPartonFragmFcn2;// ("hFragmFcn2","Parton Fragmentation Function",100,0,1);
 212         delete                          hPartonJT2;     // ("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 213         delete                          hPartonPL2;     // ("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 214         delete                          hJetJT2;        // ("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 215         delete                          hJetPL2;        // ("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 216         delete                          hJetEt2;        // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
 217         delete                          hJetEta2;       // ("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
 218         delete                          hJetPhi2;       // ("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
 219         delete                          hPartonEta2;    // ("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
 220         delete                          hPartonPhi2;    // ("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
 221         delete                          hEtaDiff2;      // ("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 222         delete                          hPhiDiff2;      // ("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 223         delete                          hEtaPhiSpread2; // ("hEtaPhiSpread2","#eta - #phi Distribution
 224                                                         //of Reconstructed Jets",192,-0.7,0.7,288,pi/3,pi);
 225         delete                          hNJets2;        // ("hNJets2","N Reconstructed jets",11,-0.5,10.5);
 226         delete                          hJetEtSecond2; //("hJetEtSecond2","E_{T}^{reco}",250,0.,250.);
 227         delete                          hJetEtRatio2;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
 228         delete                          hEtaPhiDist2;  //("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
 229
 230
 231
 232 }
 233
 234 void AliEMCALJetFinderPlots::FillFromOutput(AliEMCALJetFinderOutput* output)
 235 {
 236 if (!fInitialised) InitPlots();
 237   fOutput = output;
 238   if (!fOutput) return;
 239   hNJets->Fill(fOutput->GetNJets());
 240 if (fOutput->GetNJets()>1)
 241 {
 242 //========================= CASE 2 ===========================
 243   Int_t nPartons = fOutput->GetNPartons();
 244   hNJets2->Fill(fOutput->GetNJets());
 245   AliEMCALParton* parton;
 246   AliEMCALJet* jethighest=0;
 247   AliEMCALJet* jetsecond=0;
 248   //  Find Highest and Second Highest Jet
 249   for (Int_t counter = 0; counter<fOutput->GetNJets();counter++)
 250   {
 251           if (counter==0){
 252                   jethighest = fOutput->GetJet(0);
 253                   jetsecond  = fOutput->GetJet(1);
 254           }
 255           if (counter>0)
 256           {
 257                   Float_t energyhighest = jethighest->Energy();
 258                   Float_t energysecond = jetsecond->Energy();
 259
 260                   if ((fOutput->GetJet(counter))->Energy()>energyhighest)
 261                   {
 262                           jetsecond=jethighest;
 263                           jethighest=fOutput->GetJet(counter);
 264                   }else if ((fOutput->GetJet(counter))->Energy()>energysecond)
 265                   {
 266                           jetsecond=fOutput->GetJet(counter);
 267                   }
 268           }
 269   }
 270
 271   // End finding highest and second highest and continue
 272   hJetEt2->Fill(jethighest->Energy());
 273   hJetEta2->Fill(jethighest->Eta()  );
 274   hJetPhi2->Fill(jethighest->Phi() );
 275   if (nPartons ==0) return;
 276   parton = fOutput->GetParton(0);
 277
 278   hPartonEta2->Fill( parton->Eta() );
 279   hPartonPhi2->Fill( parton->Phi() );
 280
 281   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
 282   hEtaDiff2->Fill( jethighest->Eta() - parton->Eta() );
 283   hPhiDiff2->Fill( jethighest->Phi() - parton->Phi() );
 284   hEtaPhiSpread2->Fill(jethighest->Eta()-parton->Eta(),jethighest->Phi() - parton->Phi());
 285   hJetEtSecond2->Fill(jetsecond->Energy());
 286   hJetEtRatio2->Fill(jetsecond->Energy()/jethighest->Energy());
 287   hEtaPhiDist2->Fill( TMath::Sqrt((jethighest->Eta() - jetsecond->Eta())*(jethighest->Eta() - jetsecond->Eta())
 288                       + (jethighest->Phi() - jetsecond->Phi())*(jethighest->Phi() - jetsecond->Phi())     ));
 289   /*
 290   Float_t *pt,*phi,*eta;
 291   Int_t *pdg;
 292   pt  = new Float_t[parton->GetNTracks()];
 293   eta = new Float_t[parton->GetNTracks()];
 294   phi = new Float_t[parton->GetNTracks()];
 295   pdg = new Int_t[parton->GetNTracks()];*/
 296
 297
 298
 299  Float_t pt[2000];
 300  Float_t eta[2000];
 301  Float_t phi[2000];
 302  Int_t pdg[2000];
 303
 304   parton->GetTrackList(pt,eta,phi,pdg);
 305   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ )
 306   {
 307       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
 308            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue;
 309       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 310       Double_t rt = 2.0*atan(exp(-parton->Eta()));
 311       Double_t ctt = cos(tt);
 312       Double_t crt = cos(rt);
 313       Double_t stt = sin(tt);
 314       Double_t srt = sin(rt);
 315       Double_t ctp = cos(phi[iT]);
 316       Double_t crp = cos(parton->Phi());
 317       Double_t stp = sin(phi[iT]);
 318       Double_t srp = sin(parton->Phi());
 319       Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 320       Double_t correctp = pt[iT]/stt;
 321       hPartonPL2->Fill( correctp*cos(alpha));
 322       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +
 323       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 )
 324               hPartonJT2->Fill( correctp*sin(alpha));
 325       if (fNominalEnergy == 0.0) {
 326               hPartonFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
 327       }else
 328       {
 329               hPartonFragmFcn2->Fill(correctp*sin(tt)/fNominalEnergy);
 330       }
 331   }// loop over tracks
 332
 333 /*
 334   pt  = new Float_t[jet->NTracks()];
 335   eta = new Float_t[jet->NTracks()];
 336   phi = new Float_t[jet->NTracks()];
 337   pdg = new Int_t[jet->NTracks()];*/
 338   jethighest->TrackList(pt,eta,phi,pdg);
 339   for(Int_t iT=0; iT< jethighest->NTracks() ; iT++ )
 340   {
 341           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 342           Double_t rt = 2.0*atan(exp(-jethighest->Eta()));
 343           Double_t ctt = cos(tt);
 344           Double_t crt = cos(rt);
 345           Double_t stt = sin(tt);
 346           Double_t srt = sin(rt);
 347           Double_t ctp = cos(phi[iT]);
 348           Double_t crp = cos(jethighest->Phi());
 349           Double_t stp = sin(phi[iT]);
 350           Double_t srp = sin(jethighest->Phi());
 351           Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 352           Double_t correctp = pt[iT]/stt;
 353           hJetPL2->Fill( correctp*cos(alpha));
 354           if ( (jethighest->Eta()-eta[iT])*(jethighest->Eta()-eta[iT]) +
 355                           (jethighest->Phi()-phi[iT])*(jethighest->Phi()-phi[iT]) < 0.2*0.2 )
 356                   hJetJT2->Fill( correctp*sin(alpha));
 357           if (fNominalEnergy==0.0){
 358                   hFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
 359           } else
 360           {
 361                   hFragmFcn2->Fill(  correctp*sin(tt)/fNominalEnergy );
 362           }
 363   }// loop over tracks
 364   }
 365
 366   if (fOutput->GetNJets()==1)
 367   {
 368
 369 //========================= CASE 1 ===========================
 370   Int_t nPartons = fOutput->GetNPartons();
 371   if (fOutput->GetNJets()!=1) return;
 372   AliEMCALParton* parton;
 373   AliEMCALJet* jet;
 374   jet = fOutput->GetJet(0);
 375   hJetEt->Fill(jet->Energy());
 376   hJetEta->Fill(jet->Eta()  );
 377   hJetPhi->Fill(jet->Phi() );
 378   if (nPartons ==0) return;
 379   parton = fOutput->GetParton(0);
 380
 381   hPartonEta->Fill( parton->Eta() );
 382   hPartonPhi->Fill( parton->Phi() );
 383
 384   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
 385   hEtaDiff->Fill( jet->Eta() - parton->Eta() );
 386   hPhiDiff->Fill( jet->Phi() - parton->Phi() );
 387   hEtaPhiSpread->Fill(jet->Eta()-parton->Eta(),jet->Phi() - parton->Phi());
 388  /*
 389   Float_t *pt,*phi,*eta;
 390   Int_t *pdg;
 391   pt  = new Float_t[parton->GetNTracks()];
 392   eta = new Float_t[parton->GetNTracks()];
 393   phi = new Float_t[parton->GetNTracks()];
 394   pdg = new Int_t[parton->GetNTracks()];*/
 395
 396
 397
 398  Float_t pt[2000];
 399  Float_t eta[2000];
 400  Float_t phi[2000];
 401  Int_t pdg[2000];
 402
 403   parton->GetTrackList(pt,eta,phi,pdg);
 404   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ )
 405   {
 406       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
 407            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue;
 408       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 409       Double_t rt = 2.0*atan(exp(-parton->Eta()));
 410       Double_t ctt = cos(tt);
 411       Double_t crt = cos(rt);
 412       Double_t stt = sin(tt);
 413       Double_t srt = sin(rt);
 414       Double_t ctp = cos(phi[iT]);
 415       Double_t crp = cos(parton->Phi());
 416       Double_t stp = sin(phi[iT]);
 417       Double_t srp = sin(parton->Phi());
 418       Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 419       Double_t correctp = pt[iT]/stt;
 420       hPartonPL->Fill( correctp*cos(alpha));
 421       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +
 422       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 )
 423               hPartonJT->Fill( correctp*sin(alpha));
 424       if (fNominalEnergy == 0.0) {
 425               hPartonFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
 426       }else
 427       {
 428               hPartonFragmFcn->Fill(correctp*sin(tt)/fNominalEnergy);
 429       }
 430   }// loop over tracks
 431
 432 /*
 433   pt  = new Float_t[jet->NTracks()];
 434   eta = new Float_t[jet->NTracks()];
 435   phi = new Float_t[jet->NTracks()];
 436   pdg = new Int_t[jet->NTracks()];*/
 437   jet->TrackList(pt,eta,phi,pdg);
 438   for(Int_t iT=0; iT< jet->NTracks() ; iT++ )
 439   {
 440           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 441           Double_t rt = 2.0*atan(exp(-jet->Eta()));
 442           Double_t ctt = cos(tt);
 443           Double_t crt = cos(rt);
 444           Double_t stt = sin(tt);
 445           Double_t srt = sin(rt);
 446           Double_t ctp = cos(phi[iT]);
 447           Double_t crp = cos(jet->Phi());
 448           Double_t stp = sin(phi[iT]);
 449           Double_t srp = sin(jet->Phi());
 450           Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 451           Double_t correctp = pt[iT]/stt;
 452           hJetPL->Fill( correctp*cos(alpha));
 453           if ( (jet->Eta()-eta[iT])*(jet->Eta()-eta[iT]) +
 454                           (jet->Phi()-phi[iT])*(jet->Phi()-phi[iT]) < 0.2*0.2 )
 455                   hJetJT->Fill( correctp*sin(alpha));
 456           if (fNominalEnergy==0.0){
 457                   hFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
 458           } else
 459           {
 460                   hFragmFcn->Fill(  correctp*sin(tt)/fNominalEnergy );
 461           }
 462   }// loop over tracks
 463   }
 464 }
 465
 466