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     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 delete    fhFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
 195 delete    fhPartonFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
 196 delete   fhPartonJT;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 197 delete    fhPartonPL;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 198 delete    fhJetJT;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 199 delete    fhJetPL;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 200 delete    fhJetEt;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
 201 delete    fhJetEta;// = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
 202 delete    fhJetPhi;// = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
 203 delete    fhPartonEta;// = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
 204 delete    fhPartonPhi;// = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
 205 delete    fhEtaDiff;// = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 206 delete    fhPhiDiff;//  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 207 delete    fhNJets;// = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
 208 delete    fhEtaPhiSpread;
 209
 210         delete                          fhFragmFcn2;    // ("hFragmFcn2","Fragmentation Function",100,0,1);
 211         delete                          fhPartonFragmFcn2;// ("hFragmFcn2","Parton Fragmentation Function",100,0,1);
 212         delete                          fhPartonJT2;    // ("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 213         delete                          fhPartonPL2;    // ("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 214         delete                          fhJetJT2;       // ("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 215         delete                          fhJetPL2;       // ("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 216         delete                          fhJetEt2;       // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
 217         delete                          fhJetEta2;      // ("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
 218         delete                          fhJetPhi2;      // ("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
 219         delete                          fhPartonEta2;   // ("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
 220         delete                          fhPartonPhi2;   // ("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
 221         delete                          fhEtaDiff2;     // ("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 222         delete                          fhPhiDiff2;     // ("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 223         delete                          fhEtaPhiSpread2;        // ("hEtaPhiSpread2","#eta - #phi Distribution
 224                                                         //of Reconstructed Jets",192,-0.7,0.7,288,pi/3,pi);
 225         delete                          fhNJets2;       // ("hNJets2","N Reconstructed jets",11,-0.5,10.5);
 226         delete                          fhJetEtSecond2; //("hJetEtSecond2","E_{T}^{reco}",250,0.,250.);
 227         delete                          fhJetEtRatio2;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
 228         delete                          fhEtaPhiDist2;  //("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   fhNJets->Fill(fOutput->GetNJets());
 240 if (fOutput->GetNJets()>1)
 241 {
 242 //========================= CASE 2 ===========================
 243   Int_t nPartons = fOutput->GetNPartons();
 244   fhNJets2->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   fhJetEt2->Fill(jethighest->Energy());
 273   fhJetEta2->Fill(jethighest->Eta()  );
 274   fhJetPhi2->Fill(jethighest->Phi() );
 275   if (nPartons ==0) return;
 276   parton = fOutput->GetParton(0);
 277
 278   fhPartonEta2->Fill( parton->Eta() );
 279   fhPartonPhi2->Fill( parton->Phi() );
 280
 281   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
 282   fhEtaDiff2->Fill( jethighest->Eta() - parton->Eta() );
 283   fhPhiDiff2->Fill( jethighest->Phi() - parton->Phi() );
 284   fhEtaPhiSpread2->Fill(jethighest->Eta()-parton->Eta(),jethighest->Phi() - parton->Phi());
 285   fhJetEtSecond2->Fill(jetsecond->Energy());
 286   fhJetEtRatio2->Fill(jetsecond->Energy()/jethighest->Energy());
 287   fhEtaPhiDist2->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       fhPartonPL2->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               fhPartonJT2->Fill( correctp*sin(alpha));
 325       if (fNominalEnergy == 0.0) {
 326               fhPartonFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
 327       }else
 328       {
 329               fhPartonFragmFcn2->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           fhJetPL2->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                   fhJetJT2->Fill( correctp*sin(alpha));
 357           if (fNominalEnergy==0.0){
 358                   fhFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
 359           } else
 360           {
 361                   fhFragmFcn2->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   fhJetEt->Fill(jet->Energy());
 376   fhJetEta->Fill(jet->Eta()  );
 377   fhJetPhi->Fill(jet->Phi() );
 378   if (nPartons ==0) return;
 379   parton = fOutput->GetParton(0);
 380
 381   fhPartonEta->Fill( parton->Eta() );
 382   fhPartonPhi->Fill( parton->Phi() );
 383
 384   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
 385   fhEtaDiff->Fill( jet->Eta() - parton->Eta() );
 386   fhPhiDiff->Fill( jet->Phi() - parton->Phi() );
 387   fhEtaPhiSpread->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       fhPartonPL->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               fhPartonJT->Fill( correctp*sin(alpha));
 424       if (fNominalEnergy == 0.0) {
 425               fhPartonFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
 426       }else
 427       {
 428               fhPartonFragmFcn->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           fhJetPL->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                   fhJetJT->Fill( correctp*sin(alpha));
 456           if (fNominalEnergy==0.0){
 457                   fhFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
 458           } else
 459           {
 460                   fhFragmFcn->Fill(  correctp*sin(tt)/fNominalEnergy );
 461           }
 462   }// loop over tracks
 463   }
 464 }
 465
 466