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