EMCAL/AliEMCALJetFinderPlots.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 /* $Id$ */
  18
  19
  20 //_________________________________________________________________________
  21 //  Class for Filling JetFinder Plots
  22 // --
  23 //*-- Author: Mark Horner (LBL/UCT)
  24 // --
  25 // --
  26
  27
  28 #include "TMath.h"
  29 #include "AliEMCALJetFinderPlots.h"
  30
  31 ClassImp(AliEMCALJetFinderPlots)
  32
  33 AliEMCALJetFinderPlots::AliEMCALJetFinderPlots()
  34 {
  35         // Constructor to initialise variables
  36   fInitialised = kFALSE;
  37   fNominalEnergy = 0.0;
  38   fConeRadius = 0.3;
  39   fDebug = 0;
  40   fOutput=0;
  41     fhFragmFcn=0;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
  42     fhPartonFragmFcn=0;// = new TH1F("hPartonFragmFcn","Fragmentation Function",100,0,1);
  43    fhPartonJT=0;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
  44     fhPartonPL=0;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
  45     fhJetJT=0;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
  46     fhJetPL=0;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
  47     fhJetEt=0;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
  48     fhJetEta=0;// = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
  49     fhJetPhi=0;// = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
  50     fhPartonEta=0;// = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
  51     fhPartonPhi=0;// = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
  52     fhEtaDiff=0;// = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
  53     fhPhiDiff=0;//  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
  54     fhNJets=0;// = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
  55   fhEtaPhiSpread=0;
  56
  57 fhFragmFcn2=0;  // ("hFragmFcn2","Fragmentation Function",100,0,1);
  58 fhPartonFragmFcn2=0;// ("hFragmFcn2","Parton Fragmentation Function",100,0,1);
  59 fhPartonJT2=0;  // ("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
  60 fhPartonPL2=0;  // ("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
  61 fhJetJT2=0;     // ("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
  62 fhJetPL2=0;     // ("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
  63 fhJetEt2=0;     // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
  64 fhJetEta2=0;    // ("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
  65 fhJetPhi2=0;    // ("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
  66 fhPartonEta2=0; // ("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
  67 fhPartonPhi2=0; // ("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
  68 fhEtaDiff2=0;   // ("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
  69 fhPhiDiff2=0;   // ("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
  70 fhEtaPhiSpread2=0;      // ("hEtaPhiSpread2","#eta - #phi Distribution
  71                                                 //of Reconstructed Jets",192,-0.7,0.7,288,pi/3,pi);
  72 fhNJets2=0;     // ("hNJets2","N Reconstructed jets",11,-0.5,10.5);
  73 fhJetEtSecond2=0; //("hJetEtSecond2","E_{T}^{reco}",250,0.,250.);
  74 fhJetEtRatio2=0;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
  75 fhEtaPhiDist2=0;  //("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
  76 fhInputOutput=0;
  77 //     TH2F                            *fhInputOutput;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
  78
  79 fhRecoBinPt=0;         // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
  80 fhRecoBinPartonPt=0;    // ("fhRecoBinPartonPt","Input Pt Distribution",100,0,1);
  81 fhRecoBinJetEt=0;       // ("fhRecoJetEt","E_{T}^{reco}",250,0.,250.);
  82 fhRecoBinInputJetEt=0;  // ("fhRecoInputJetEt","E_{T}^{reco}",250,0.,250.);
  83 fhJetPT =0;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
  84 fhPartonPT =0;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
  85 fhJetPT2 =0;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
  86 fhPartonPT2 =0;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
  87 fhRecoBinFragmFcn =0;//new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
  88 fhRecoBinPartonFragmFcn =0;// new TH1F("fhRecoBinPartonFragmFcn","Input Bin Fragm Fcn Distribution",100,0,1);
  89 }
  90
  91 void AliEMCALJetFinderPlots::InitPlots()
  92 {
  93 //========================= CASE 1 =======================================
  94     fhFragmFcn = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
  95     fhFragmFcn->Sumw2();
  96     fhJetPT = new TH1F("hJetPT","P_{T} Distribution",200,0,200);
  97     fhJetPT->Sumw2();
  98     fhPartonPT = new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
  99     fhPartonPT->Sumw2();
 100     fhPartonFragmFcn = new TH1F("hPartonFragmFcn","Parton Fragmentation Function",100,0,1);
 101     fhPartonFragmFcn->Sumw2();
 102     fhPartonJT = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 103     fhPartonJT->Sumw2();
 104     fhPartonPL = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 105     fhPartonPL->Sumw2();
 106     fhJetJT = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 107     fhJetJT->Sumw2();
 108     fhJetPL = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 109     fhJetPL->Sumw2();
 110     fhJetEt = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
 111     fhJetEt->Sumw2();
 112     fhJetEta = new      TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
 113     fhJetEta->Sumw2();
 114     fhJetPhi = new      TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
 115     fhJetPhi->Sumw2();
 116     fhPartonEta = new   TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
 117     fhPartonEta->Sumw2();
 118     fhPartonPhi = new   TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
 119     fhPartonPhi->Sumw2();
 120     fhEtaDiff = new     TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 121     fhEtaDiff->Sumw2();
 122     fhPhiDiff  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 123     fhPhiDiff->Sumw2();
 124     fhNJets = new       TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
 125     fhNJets->Sumw2();
 126     fhEtaPhiSpread = new TH2F("hEtaPhiSpread","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
 127     fhEtaPhiSpread->Sumw2();
 128   fhNJets->SetXTitle("N_{jets}^{reco}/event");
 129   fhNJets->SetYTitle("N_{events}");
 130
 131   //Jet properties
 132   fhJetEt->SetFillColor(16);
 133   fhJetEt->SetXTitle("E_{T}^{reco}");
 134
 135   fhJetEta->SetFillColor(16);
 136   fhJetEta->SetXTitle("#eta_{jet}^{reco}");
 137
 138   fhJetPhi->SetFillColor(16);
 139   fhJetPhi->SetXTitle("#phi_{jet}^{reco}");
 140
 141   fhPartonEta->SetFillColor(16);
 142   fhPartonEta->SetXTitle("#eta_{parton}");
 143
 144   fhPartonPhi->SetFillColor(16);
 145   fhPartonPhi->SetXTitle("#phi_{parton}");
 146
 147   fhPartonPL->SetXTitle("p (GeV/c)");
 148   fhPartonJT->SetXTitle("p (GeV/c)");
 149
 150   fhPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{parton}");
 151
 152   //Jet component properties
 153
 154   fhJetPL->SetXTitle("p (GeV/c)");
 155   fhJetJT->SetXTitle("p (GeV/c)");
 156   fhFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
 157   fhPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
 158
 159   fhEtaDiff->SetXTitle("#eta_{jet}^{reco}-#eta_{jet}^{input}");
 160   fhPhiDiff->SetXTitle("#phi_{jet}^{reco}-#phi_{jet}^{input}");
 161   fhEtaPhiSpread->SetXTitle("#eta");
 162   fhEtaPhiSpread->SetYTitle("#phi");
 163
 164 //======================= CASE 2 ======================================
 165
 166
 167 fhFragmFcn2             = new TH1F("hFragmFcn2","Fragmentation Function",100,0,1);
 168 fhFragmFcn2->Sumw2();
 169     fhJetPT2 = new TH1F("hJetPT2","P_{T} Distribution",200,0,200);
 170     fhJetPT2->Sumw2();
 171     fhPartonPT2 = new TH1F("hPartonPT2","Parton P_{T} Distribution",200,0,1);
 172     fhPartonPT2->Sumw2();
 173 fhPartonFragmFcn2       = new TH1F("hPartonFragmFcn2","Parton Fragmentation Function",100,0,1);
 174 fhPartonFragmFcn2->Sumw2();
 175 fhPartonJT2             = new TH1F("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 176 fhPartonJT2->Sumw2();
 177 fhPartonPL2             = new TH1F("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 178 fhPartonPL2->Sumw2();
 179 fhJetJT2                        = new TH1F("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 180 fhJetJT2->Sumw2();
 181 fhJetPL2                        = new TH1F("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 182 fhJetPL2->Sumw2();
 183 fhJetEt2                        = new TH1F("hJetEt2","E_{T}^{reco}",250,0.,250.);
 184 fhJetEt2->Sumw2();
 185 fhJetEta2               = new TH1F("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
 186 fhJetEta2->Sumw2();
 187 fhJetPhi2               = new TH1F("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
 188 fhJetPhi2->Sumw2();
 189 fhPartonEta2            = new TH1F("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
 190 fhPartonEta2->Sumw2();
 191 fhPartonPhi2            = new TH1F("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
 192 fhPartonPhi2->Sumw2();
 193 fhEtaDiff2              = new TH1F("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 194 fhEtaDiff2->Sumw2();
 195 fhPhiDiff2              = new TH1F("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 196 fhPhiDiff2->Sumw2();
 197 fhEtaPhiSpread2 = new TH2F("hEtaPhiSpread2","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
 198 fhEtaPhiSpread2->Sumw2();
 199 fhNJets2                        = new TH1F("hNJets2","N Reconstructed jets",11,-0.5,10.5);
 200 fhNJets2->Sumw2();
 201 fhJetEtSecond2          = new TH1F("hJetEtSecond2","E_{T}^{reco}",250,0.,250.);
 202 fhJetEtSecond2->Sumw2();
 203 fhJetEtRatio2           = new TH1F("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
 204 fhJetEtRatio2->Sumw2();
 205 fhEtaPhiDist2           = new TH1F("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
 206 fhEtaPhiDist2->Sumw2();
 207
 208 fhInputOutput=  new TH2F("hInputOutput","Input and Reconstruction Correlations;Input;Output",200,0,200,200,0,200);  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
 209
 210 //============================== Reconstruction Bin Comparison  ============================================
 211
 212 fhRecoBinPt =new TH1F("fhRecoBinPt","Reconstructed Pt Distribution",200,0,200);
 213 fhRecoBinPt->Sumw2();
 214 fhRecoBinPartonPt = new TH1F("fhRecoBinPartonPt","Input Pt Distribution",200,0,200);
 215 fhRecoBinPartonPt->Sumw2();
 216 fhRecoBinFragmFcn =new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
 217 fhRecoBinFragmFcn->Sumw2();
 218 fhRecoBinPartonFragmFcn = new TH1F("fhRecoBinPartonFragmFcn","Input Bin Fragm Fcn Distribution",100,0,1);
 219 fhRecoBinPartonFragmFcn->Sumw2();
 220 fhRecoBinJetEt = new TH1F("fhRecoJetEt","E_{T}^{reco}",250,0.,250.);
 221 fhRecoBinJetEt->Sumw2();
 222 fhRecoBinInputJetEt = new TH1F("fhRecoInputJetEt","E_{T}^{reco}",250,0.,250.);
 223 fhRecoBinInputJetEt->Sumw2();
 224
 225   fInitialised = kTRUE;
 226
 227 }
 228
 229 AliEMCALJetFinderPlots::~AliEMCALJetFinderPlots()
 230 {
 231         // To ensure that all requested memory is returned
 232 delete    fhFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
 233 delete    fhPartonFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
 234 delete   fhPartonJT;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 235 delete    fhPartonPL;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 236 delete    fhJetJT;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 237 delete    fhJetPL;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 238 delete    fhJetEt;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
 239 delete    fhJetEta;// = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
 240 delete    fhJetPhi;// = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
 241 delete    fhPartonEta;// = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
 242 delete    fhPartonPhi;// = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
 243 delete    fhEtaDiff;// = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 244 delete    fhPhiDiff;//  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 245 delete    fhNJets;// = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
 246 delete    fhEtaPhiSpread;
 247
 248         delete                          fhFragmFcn2;    // ("hFragmFcn2","Fragmentation Function",100,0,1);
 249         delete                          fhPartonFragmFcn2;// ("hFragmFcn2","Parton Fragmentation Function",100,0,1);
 250         delete                          fhPartonJT2;    // ("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 251         delete                          fhPartonPL2;    // ("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 252         delete                          fhJetJT2;       // ("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 253         delete                          fhJetPL2;       // ("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 254         delete                          fhJetEt2;       // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
 255         delete                          fhJetEta2;      // ("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
 256         delete                          fhJetPhi2;      // ("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
 257         delete                          fhPartonEta2;   // ("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
 258         delete                          fhPartonPhi2;   // ("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
 259         delete                          fhEtaDiff2;     // ("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 260         delete                          fhPhiDiff2;     // ("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 261         delete                          fhEtaPhiSpread2;        // ("hEtaPhiSpread2","#eta - #phi Distribution
 262                                                         //of Reconstructed Jets",192,-0.7,0.7,288,pi/3,pi);
 263         delete                          fhNJets2;       // ("hNJets2","N Reconstructed jets",11,-0.5,10.5);
 264         delete                          fhJetEtSecond2; //("hJetEtSecond2","E_{T}^{reco}",250,0.,250.);
 265         delete                          fhJetEtRatio2;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
 266         delete                          fhEtaPhiDist2;  //("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
 267
 268         delete fhRecoBinPt;          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 269         delete fhRecoBinPartonPt;    // ("fhRecoBinPartonPt","Input Pt Distribution",100,0,1);
 270         delete fhRecoBinJetEt;       // ("fhRecoJetEt","E_{T}^{reco}",250,0.,250.);
 271         delete fhRecoBinInputJetEt;  // ("fhRecoInputJetEt","E_{T}^{reco}",250,0.,250.);
 272
 273         delete fhJetPT ;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
 274         delete fhPartonPT ;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
 275         delete fhJetPT2 ;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
 276         delete fhPartonPT2 ;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
 277         delete fhRecoBinFragmFcn;//new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
 278         delete fhRecoBinPartonFragmFcn;// new TH1F("fhRecoBinPartonFragmFcn","Input Bin Fragm Fcn Distribution",100,0,1);
 279
 280 }
 281
 282 void AliEMCALJetFinderPlots::FillFromOutput(AliEMCALJetFinderOutput* output)
 283 {
 284         // Fill histograms from an output object
 285 if (!fInitialised) InitPlots();
 286   fOutput = output;
 287   if (!fOutput) return;
 288   fhNJets->Fill(fOutput->GetNJets());
 289   Bool_t doesJetMeetBinCriteria = 0;
 290   AliEMCALJet* jethighest=0;
 291   AliEMCALJet* jetsecond=0;
 292   //  Find Highest and Second Highest Jet
 293   // (NB!!!!!!!)  Pointing into the EMCAL!!!!!!
 294
 295 // =========================== All cases ===================================
 296
 297
 298         // I will make a little array of jet indices for which jets are in
 299         // the EMCAL then my counter can loop from 0 below - but it will
 300         // be the index of the array of applicable jets
 301
 302   Int_t appjet[4];
 303   Int_t numappjet=0;
 304
 305   for (Int_t appc=0;appc<fOutput->GetNJets();appc++)
 306   { // Check all jets for applicability
 307           Float_t eta = fOutput->GetJet(appc)->Eta();
 308           Float_t phi = fOutput->GetJet(appc)->Phi();
 309           if (eta > -0.7 && eta < 0.7 && phi > 1./3.*TMath::Pi() && phi < TMath::Pi())
 310           { // Then jet is applicable
 311                   appjet[numappjet]=appc;
 312                   numappjet++;
 313           }
 314   }
 315
 316
 317 if (numappjet >=1)
 318 {
 319         Float_t theta = 2.0*atan(exp(-fOutput->GetParton(0)->Eta()));
 320         Float_t et    = fOutput->GetParton(0)->Energy() * TMath::Sin(theta);
 321         if (fOutput->GetNJets()>1)
 322         {
 323                 for (Int_t counter = 0; counter<numappjet;counter++)
 324                 {
 325                         if (counter==0)
 326                         {
 327                                 jethighest = fOutput->GetJet(appjet[0]);
 328                                 jetsecond  = fOutput->GetJet(appjet[1]);
 329                         }
 330                         if (counter>0)
 331                         {
 332                                 Float_t energyhighest = jethighest->Energy();
 333                                 Float_t energysecond = jetsecond->Energy();
 334
 335                                 if ((fOutput->GetJet(appjet[counter]))->Energy()>energyhighest)
 336                                 {
 337                                         jetsecond=jethighest;
 338                                         jethighest=fOutput->GetJet(appjet[counter]);
 339                                 }else if ((fOutput->GetJet(appjet[counter]))->Energy()>energysecond)
 340                                 {
 341                                         jetsecond=fOutput->GetJet(appjet[counter]);
 342                                 }
 343
 344                         }
 345                 }
 346         }else
 347         {
 348                 Float_t eta = fOutput->GetJet(0)->Eta();
 349                 Float_t phi = fOutput->GetJet(0)->Phi();
 350                 if (eta > -0.7 && eta < 0.7 && phi > 1./3.*TMath::Pi() && phi < TMath::Pi())
 351                 { // Then jet is applicable
 352                         jethighest=fOutput->GetJet(0);
 353                         jetsecond=0;
 354                 }else
 355                 {
 356                         Error("FillFromOutput","There is only one jet and it isn't in the area of applicability");
 357                 }
 358
 359         }
 360         if ( 95.0 < jethighest->Energy() && jethighest->Energy() < 105.0  )
 361         {
 362                 doesJetMeetBinCriteria = 1;
 363                 fhRecoBinJetEt->Fill(jethighest->Energy());
 364                 fhRecoBinInputJetEt->Fill(et);
 365         }
 366         fhInputOutput->Fill(et,jethighest->Energy());
 367
 368 }
 369
 370 if (numappjet > 1)
 371 {
 372 //========================= CASE 2 ===========================
 373   Int_t nPartons = fOutput->GetNPartons();
 374   fhNJets2->Fill(fOutput->GetNJets());
 375   AliEMCALParton* parton;
 376
 377   // End finding highest and second highest and continue
 378   fhJetEt2->Fill(jethighest->Energy());
 379   fhJetEta2->Fill(jethighest->Eta()  );
 380   fhJetPhi2->Fill(jethighest->Phi() );
 381   if (nPartons ==0) return;
 382   parton = fOutput->GetParton(0);
 383
 384   fhPartonEta2->Fill( parton->Eta() );
 385   fhPartonPhi2->Fill( parton->Phi() );
 386
 387   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
 388   fhEtaDiff2->Fill( jethighest->Eta() - parton->Eta() );
 389   fhPhiDiff2->Fill( jethighest->Phi() - parton->Phi() );
 390   fhEtaPhiSpread2->Fill(jethighest->Eta()-parton->Eta(),jethighest->Phi() - parton->Phi());
 391   fhJetEtSecond2->Fill(jetsecond->Energy());
 392   fhJetEtRatio2->Fill(jetsecond->Energy()/jethighest->Energy());
 393   fhEtaPhiDist2->Fill( TMath::Sqrt((jethighest->Eta() - jetsecond->Eta())*(jethighest->Eta() - jetsecond->Eta())
 394                       + (jethighest->Phi() - jetsecond->Phi())*(jethighest->Phi() - jetsecond->Phi())     ));
 395   /*
 396   Float_t *pt,*phi,*eta;
 397   Int_t *pdg;
 398   pt  = new Float_t[parton->GetNTracks()];
 399   eta = new Float_t[parton->GetNTracks()];
 400   phi = new Float_t[parton->GetNTracks()];
 401   pdg = new Int_t[parton->GetNTracks()];*/
 402
 403
 404
 405  Float_t pt[2000];
 406  Float_t eta[2000];
 407  Float_t phi[2000];
 408  Int_t pdg[2000];
 409
 410   parton->GetTrackList(pt,eta,phi,pdg);
 411   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ )
 412   {
 413       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
 414            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue;
 415       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 416       Double_t rt = 2.0*atan(exp(-parton->Eta()));
 417       Double_t ctt = cos(tt);
 418       Double_t crt = cos(rt);
 419       Double_t stt = sin(tt);
 420       Double_t srt = sin(rt);
 421       Double_t ctp = cos(phi[iT]);
 422       Double_t crp = cos(parton->Phi());
 423       Double_t stp = sin(phi[iT]);
 424       Double_t srp = sin(parton->Phi());
 425       //Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 426           Double_t alpha;
 427           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
 428           {
 429                   alpha = 0.0;
 430           }else
 431           {
 432                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 433           }
 434       Double_t correctp = pt[iT]/stt;
 435       fhPartonPL2->Fill( correctp*cos(alpha));
 436       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +
 437       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 )
 438               fhPartonJT2->Fill( correctp*sin(alpha));
 439       if (fNominalEnergy == 0.0) {
 440               fhPartonFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
 441       }else
 442       {
 443               fhPartonFragmFcn2->Fill(correctp*sin(tt)/fNominalEnergy);
 444       }
 445       if (doesJetMeetBinCriteria)
 446       {
 447               fhRecoBinPartonPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 448       }
 449   }// loop over tracks
 450
 451 /*
 452   pt  = new Float_t[jet->NTracks()];
 453   eta = new Float_t[jet->NTracks()];
 454   phi = new Float_t[jet->NTracks()];
 455   pdg = new Int_t[jet->NTracks()];*/
 456   jethighest->TrackList(pt,eta,phi,pdg);
 457   for(Int_t iT=0; iT< jethighest->NTracks() ; iT++ )
 458   {
 459           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 460           Double_t rt = 2.0*atan(exp(-jethighest->Eta()));
 461           Double_t ctt = cos(tt);
 462           Double_t crt = cos(rt);
 463           Double_t stt = sin(tt);
 464           Double_t srt = sin(rt);
 465           Double_t ctp = cos(phi[iT]);
 466           Double_t crp = cos(jethighest->Phi());
 467           Double_t stp = sin(phi[iT]);
 468           Double_t srp = sin(jethighest->Phi());
 469          // Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 470           Double_t alpha;
 471           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
 472           {
 473                   alpha = 0.0;
 474           }else
 475           {
 476                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 477           }
 478           Double_t correctp = pt[iT]/stt;
 479           fhJetPL2->Fill( correctp*cos(alpha));
 480           if ( (jethighest->Eta()-eta[iT])*(jethighest->Eta()-eta[iT]) +
 481                           (jethighest->Phi()-phi[iT])*(jethighest->Phi()-phi[iT]) < 0.2*0.2 )
 482                   fhJetJT2->Fill( correctp*sin(alpha));
 483           if (fNominalEnergy==0.0){
 484                   fhFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy()  );
 485           } else
 486           {
 487                   fhFragmFcn2->Fill(  correctp*sin(tt)/fNominalEnergy );
 488           }
 489           if (doesJetMeetBinCriteria)
 490           {
 491                   fhRecoBinPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 492           }
 493   }// loop over tracks
 494   }
 495
 496   if (numappjet == 1)
 497   {
 498
 499 //========================= CASE 1 ===========================
 500   Int_t nPartons = fOutput->GetNPartons();
 501   if (fOutput->GetNJets()!=1) return;
 502   AliEMCALParton* parton;
 503   AliEMCALJet* jet;
 504   jet = jethighest;//fOutput->GetJet(0);
 505   fhJetEt->Fill(jet->Energy());
 506   fhJetEta->Fill(jet->Eta()  );
 507   fhJetPhi->Fill(jet->Phi() );
 508   if (nPartons ==0) return;
 509   parton = fOutput->GetParton(0);
 510
 511   fhPartonEta->Fill( parton->Eta() );
 512   fhPartonPhi->Fill( parton->Phi() );
 513
 514   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
 515   fhEtaDiff->Fill( jet->Eta() - parton->Eta() );
 516   fhPhiDiff->Fill( jet->Phi() - parton->Phi() );
 517   fhEtaPhiSpread->Fill(jet->Eta()-parton->Eta(),jet->Phi() - parton->Phi());
 518  /*
 519   Float_t *pt,*phi,*eta;
 520   Int_t *pdg;
 521   pt  = new Float_t[parton->GetNTracks()];
 522   eta = new Float_t[parton->GetNTracks()];
 523   phi = new Float_t[parton->GetNTracks()];
 524   pdg = new Int_t[parton->GetNTracks()];*/
 525
 526
 527
 528  Float_t pt[2000];
 529  Float_t eta[2000];
 530  Float_t phi[2000];
 531  Int_t pdg[2000];
 532
 533   parton->GetTrackList(pt,eta,phi,pdg);
 534   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ )
 535   {
 536       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
 537            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue;
 538       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 539       Double_t rt = 2.0*atan(exp(-parton->Eta()));
 540       Double_t ctt = cos(tt);
 541       Double_t crt = cos(rt);
 542       Double_t stt = sin(tt);
 543       Double_t srt = sin(rt);
 544       Double_t ctp = cos(phi[iT]);
 545       Double_t crp = cos(parton->Phi());
 546       Double_t stp = sin(phi[iT]);
 547       Double_t srp = sin(parton->Phi());
 548      // Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 549           Double_t alpha;
 550           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
 551           {
 552                   alpha = 0.0;
 553           }else
 554           {
 555                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);           }
 556       Double_t correctp = pt[iT]/stt;
 557       fhPartonPL->Fill( correctp*cos(alpha));
 558       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +
 559       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 )
 560               fhPartonJT->Fill( correctp*sin(alpha));
 561       if (fNominalEnergy == 0.0) {
 562               fhPartonFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
 563       }else
 564       {
 565               fhPartonFragmFcn->Fill(correctp*sin(tt)/fNominalEnergy);
 566       }
 567       if (doesJetMeetBinCriteria)
 568       {
 569               fhRecoBinPartonPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 570       }
 571   }// loop over tracks
 572
 573 /*
 574   pt  = new Float_t[jet->NTracks()];
 575   eta = new Float_t[jet->NTracks()];
 576   phi = new Float_t[jet->NTracks()];
 577   pdg = new Int_t[jet->NTracks()];*/
 578   jet->TrackList(pt,eta,phi,pdg);
 579   for(Int_t iT=0; iT< jet->NTracks() ; iT++ )
 580   {
 581           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 582           Double_t rt = 2.0*atan(exp(-jet->Eta()));
 583           Double_t ctt = cos(tt);
 584           Double_t crt = cos(rt);
 585           Double_t stt = sin(tt);
 586           Double_t srt = sin(rt);
 587           Double_t ctp = cos(phi[iT]);
 588           Double_t crp = cos(jet->Phi());
 589           Double_t stp = sin(phi[iT]);
 590           Double_t srp = sin(jet->Phi());
 591           //Info("plots","acos(%1.16f)\nstt=%f\npt=%f",crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt,stt,pt[iT]);
 592           //Info("plots","diff to 1 %f",1.0-crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 593           Double_t alpha;
 594           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
 595           {
 596                   alpha = 0.0;
 597           }else
 598           {
 599                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 600           }
 601           Double_t correctp = pt[iT]/stt;
 602           fhJetPL->Fill( correctp*cos(alpha));
 603           if ( (jet->Eta()-eta[iT])*(jet->Eta()-eta[iT]) +
 604                           (jet->Phi()-phi[iT])*(jet->Phi()-phi[iT]) < 0.2*0.2 )
 605                   fhJetJT->Fill( correctp*sin(alpha));
 606           if (fNominalEnergy==0.0){
 607                   fhFragmFcn->Fill(  correctp*sin(tt)/parton->Energy()  );
 608           } else
 609           {
 610                   fhFragmFcn->Fill(  correctp*sin(tt)/fNominalEnergy );
 611           }
 612           if (doesJetMeetBinCriteria)
 613           {
 614                   fhRecoBinPt->Fill(correctp*sin(tt));          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 615           }
 616   }// loop over tracks
 617   }
 618 }
 619
 620