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 fhRecoBinPtNoBg=0;             // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
  81 fhRecoBinPartonPt=0;    // ("fhRecoBinPartonPt","Input Pt Distribution",100,0,1);
  82 fhRecoBinJetEt=0;       // ("fhRecoJetEt","E_{T}^{reco}",250,0.,250.);
  83 fhRecoBinInputJetEt=0;  // ("fhRecoInputJetEt","E_{T}^{reco}",250,0.,250.);
  84 fhJetPT =0;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
  85 fhPartonPT =0;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
  86 fhJetPT2 =0;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
  87 fhPartonPT2 =0;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
  88 fhRecoBinFragmFcn =0;//new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
  89 fhRecoBinFragmFcnNoBg =0;//new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
  90 fhRecoBinPartonFragmFcn =0;// new TH1F("fhRecoBinPartonFragmFcn","Input Bin Fragm Fcn Distribution",100,0,1);
  91 fhJetInvE=0;// = new TH1F("fhJetInvE","#frac{1}{E_{R}}",100,0,1);
  92 fhJetInvE2=0;// = new TH1F("fhJetInvE","#frac{1}{E_{R}}",100,0,1);
  93 fScaleFactor = 1.0/0.6731;
  94 fhBackHisto=0;
  95
  96 }
  97
  98 void AliEMCALJetFinderPlots::InitPlots()
  99 {
 100 //========================= CASE 1 =======================================
 101     fhFragmFcn = new TH1F("hFragmFcn","Fragmentation Function",200,0,2);
 102     fhFragmFcn->Sumw2();
 103     fhJetPT = new TH1F("hJetPT","P_{T} Distribution",200,0,200);
 104     fhJetPT->Sumw2();
 105     fhPartonPT = new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,200);
 106     fhPartonPT->Sumw2();
 107     fhPartonFragmFcn = new TH1F("hPartonFragmFcn","Parton Fragmentation Function",200,0,2);
 108     fhPartonFragmFcn->Sumw2();
 109     fhPartonJT = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 110     fhPartonJT->Sumw2();
 111     fhPartonPL = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 112     fhPartonPL->Sumw2();
 113     fhJetJT = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 114     fhJetJT->Sumw2();
 115     fhJetPL = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 116     fhJetPL->Sumw2();
 117     fhJetEt = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
 118     fhJetEt->Sumw2();
 119     fhJetEtDiff = new TH1F("hJetEtDiff","E_{T}^{reco}-E_{T}^{Parton}",250,-124.,125.);
 120     fhJetEtDiff->Sumw2();
 121     fhJetEta = new      TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
 122     fhJetEta->Sumw2();
 123     fhJetPhi = new      TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
 124     fhJetPhi->Sumw2();
 125     fhPartonEta = new   TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
 126     fhPartonEta->Sumw2();
 127     fhPartonPhi = new   TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
 128     fhPartonPhi->Sumw2();
 129     fhEtaDiff = new     TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 130     fhEtaDiff->Sumw2();
 131     fhPhiDiff  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 132     fhPhiDiff->Sumw2();
 133     fhNJets = new       TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
 134     fhNJets->Sumw2();
 135     fhEtaPhiSpread = new TH2F("hEtaPhiSpread","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
 136     fhEtaPhiSpread->Sumw2();
 137   fhNJets->SetXTitle("N_{jets}^{reco}/event");
 138   fhNJets->SetYTitle("N_{events}");
 139
 140   //Jet properties
 141   fhJetEt->SetFillColor(16);
 142   fhJetEt->SetXTitle("E_{T}^{reco}");
 143
 144   fhJetEta->SetFillColor(16);
 145   fhJetEta->SetXTitle("#eta_{jet}^{reco}");
 146
 147   fhJetPhi->SetFillColor(16);
 148   fhJetPhi->SetXTitle("#phi_{jet}^{reco}");
 149
 150   fhPartonEta->SetFillColor(16);
 151   fhPartonEta->SetXTitle("#eta_{parton}");
 152
 153   fhPartonPhi->SetFillColor(16);
 154   fhPartonPhi->SetXTitle("#phi_{parton}");
 155
 156   fhPartonPL->SetXTitle("p (GeV/c)");
 157   fhPartonJT->SetXTitle("p (GeV/c)");
 158
 159   fhPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{parton}");
 160
 161   //Jet component properties
 162
 163   fhJetPL->SetXTitle("p (GeV/c)");
 164   fhJetJT->SetXTitle("p (GeV/c)");
 165   fhFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
 166   fhPartonFragmFcn->SetXTitle("Z = p_{T}^{Chg}/E_{T}^{reco}");
 167
 168   fhEtaDiff->SetXTitle("#eta_{jet}^{reco}-#eta_{jet}^{input}");
 169   fhPhiDiff->SetXTitle("#phi_{jet}^{reco}-#phi_{jet}^{input}");
 170   fhEtaPhiSpread->SetXTitle("#eta");
 171   fhEtaPhiSpread->SetYTitle("#phi");
 172
 173 //======================= CASE 2 ======================================
 174
 175
 176 fhFragmFcn2             = new TH1F("hFragmFcn2","Fragmentation Function",200,0,2);
 177 fhFragmFcn2->Sumw2();
 178     fhJetPT2 = new TH1F("hJetPT2","P_{T} Distribution",200,0,200);
 179     fhJetPT2->Sumw2();
 180     fhPartonPT2 = new TH1F("hPartonPT2","Parton P_{T} Distribution",200,0,1);
 181     fhPartonPT2->Sumw2();
 182 fhPartonFragmFcn2       = new TH1F("hPartonFragmFcn2","Parton Fragmentation Function",200,0,2);
 183 fhPartonFragmFcn2->Sumw2();
 184 fhPartonJT2             = new TH1F("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 185 fhPartonJT2->Sumw2();
 186 fhPartonPL2             = new TH1F("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 187 fhPartonPL2->Sumw2();
 188 fhJetJT2                        = new TH1F("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 189 fhJetJT2->Sumw2();
 190 fhJetPL2                        = new TH1F("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 191 fhJetPL2->Sumw2();
 192 fhJetEt2                        = new TH1F("hJetEt2","E_{T}^{reco}",250,0.,250.);
 193 fhJetEt2->Sumw2();
 194     fhJetEtDiff2 = new TH1F("hJetEtDiff2","E_{T}^{reco}-E_{T}^{Parton}",250,-124.,125.);
 195     fhJetEtDiff2->Sumw2();
 196 fhJetEta2               = new TH1F("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
 197 fhJetEta2->Sumw2();
 198 fhJetPhi2               = new TH1F("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
 199 fhJetPhi2->Sumw2();
 200 fhPartonEta2            = new TH1F("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
 201 fhPartonEta2->Sumw2();
 202 fhPartonPhi2            = new TH1F("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
 203 fhPartonPhi2->Sumw2();
 204 fhEtaDiff2              = new TH1F("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 205 fhEtaDiff2->Sumw2();
 206 fhPhiDiff2              = new TH1F("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 207 fhPhiDiff2->Sumw2();
 208 fhEtaPhiSpread2 = new TH2F("hEtaPhiSpread2","#eta - #phi Distribution of Reconstructed Jets",100,-0.5,0.5,100,-0.5,0.5);
 209 fhEtaPhiSpread2->Sumw2();
 210 fhNJets2                        = new TH1F("hNJets2","N Reconstructed jets",11,-0.5,10.5);
 211 fhNJets2->Sumw2();
 212 fhJetEtSecond2          = new TH1F("hJetEtSecond2","E_{T}^{reco}",250,0.,250.);
 213 fhJetEtSecond2->Sumw2();
 214 fhJetEtRatio2           = new TH1F("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
 215 fhJetEtRatio2->Sumw2();
 216 fhEtaPhiDist2           = new TH1F("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
 217 fhEtaPhiDist2->Sumw2();
 218
 219 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);
 220
 221 //============================== Reconstruction Bin Comparison  ============================================
 222
 223 fhRecoBinPt =new TH1F("fhRecoBinPt","Reconstructed Pt Distribution",200,0,200);
 224 fhRecoBinPt->Sumw2();
 225 fhRecoBinPtNoBg =new TH1F("fhRecoBinPtNoBg","Reconstructed Pt Distribution Background Subtracted",200,0,200);
 226 fhRecoBinPtNoBg->Sumw2();
 227 fhRecoBinPartonPt = new TH1F("fhRecoBinPartonPt","Input Pt Distribution",200,0,200);
 228 fhRecoBinPartonPt->Sumw2();
 229 fhRecoBinFragmFcn =new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",200,0,2);
 230 fhRecoBinFragmFcn->Sumw2();
 231 fhRecoBinFragmFcnNoBg =new TH1F("fhRecoBinFragmFcnNoBg","Reconstructed Frag. Fcn With Background Removed",200,0,2);
 232 fhRecoBinFragmFcnNoBg->Sumw2();
 233 fhRecoBinPartonFragmFcn = new TH1F("fhRecoBinPartonFragmFcn","Input Bin Fragm Fcn Distribution",200,0,2);
 234 fhRecoBinPartonFragmFcn->Sumw2();
 235 fhRecoBinJetEt = new TH1F("fhRecoJetEt","E_{T}^{reco}",250,0.,250.);
 236 fhRecoBinJetEt->Sumw2();
 237 fhRecoBinInputJetEt = new TH1F("fhRecoInputJetEt","E_{T}^{reco}",250,0.,250.);
 238 fhRecoBinInputJetEt->Sumw2();
 239
 240
 241 fhJetInvE = new TH1F("fhJetInvE","#frac{1}{E_{R}}",100,0,1);
 242 fhJetInvE->Sumw2();
 243 fhJetInvE2 = new TH1F("fhJetInvE2","#frac{1}{E_{R}}",100,0,1);
 244 fhJetInvE2->Sumw2();
 245
 246
 247
 248   fInitialised = kTRUE;
 249
 250 }
 251
 252 AliEMCALJetFinderPlots::~AliEMCALJetFinderPlots()
 253 {
 254         // To ensure that all requested memory is returned
 255 delete    fhFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
 256 delete    fhPartonFragmFcn;// = new TH1F("hFragmFcn","Fragmentation Function",100,0,1);
 257 delete   fhPartonJT;// = new TH1F("hPartonJT","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 258 delete    fhPartonPL;// = new TH1F("hPartonPL","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 259 delete    fhJetJT;// = new TH1F("hJetJT","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 260 delete    fhJetPL;// = new TH1F("hJetPL","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 261 delete    fhJetEt;// = new TH1F("hJetEt","E_{T}^{reco}",250,0.,250.);
 262 delete  fhJetEtDiff;    // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
 263 delete    fhJetEta;// = new       TH1F("hJetEta","#eta_{jet}^{reco}",180,-0.9,0.9);
 264 delete    fhJetPhi;// = new       TH1F("hJetPhi","#phi_{jet}^{reco}",62,0.,3.1);
 265 delete    fhPartonEta;// = new    TH1F("hPartonEta","#eta_{Parton}",180,-0.9,0.9);
 266 delete    fhPartonPhi;// = new    TH1F("hPartonPhi","#phi_{Parton}",62,0.,3.1);
 267 delete    fhEtaDiff;// = new      TH1F("hEtaDiff","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 268 delete    fhPhiDiff;//  = new TH1F("hPhiDiff","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 269 delete    fhNJets;// = new        TH1F("hNJets","N Reconstructed jets",11,-0.5,10.5);
 270 delete    fhEtaPhiSpread;
 271
 272         delete                          fhFragmFcn2;    // ("hFragmFcn2","Fragmentation Function",100,0,1);
 273         delete                          fhPartonFragmFcn2;// ("hFragmFcn2","Parton Fragmentation Function",100,0,1);
 274         delete                          fhPartonJT2;    // ("hPartonJT2","Track Momentum Perpendicular to Parton Axis",100,0.,10.);
 275         delete                          fhPartonPL2;    // ("hPartonPL2","Track Momentum Parallel to Parton Axis ",100,0.,100.);
 276         delete                          fhJetJT2;       // ("hJetJT2","Track Momentum Perpendicular to Jet Axis",100,0.,10.);
 277         delete                          fhJetPL2;       // ("hJetPL2","Track Momentum Parallel to Jet Axis ",100,0.,100.);
 278         delete                          fhJetEt2;       // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
 279         delete                          fhJetEtDiff2;   // ("hJetEt2","E_{T}^{reco}",250,0.,250.);
 280         delete                          fhJetEta2;      // ("hJetEta2","#eta_{jet}^{reco}",180,-0.9,0.9);
 281         delete                          fhJetPhi2;      // ("hJetPhi2","#phi_{jet}^{reco}",62,0.,3.1);
 282         delete                          fhPartonEta2;   // ("hPartonEta2","#eta_{Parton}",180,-0.9,0.9);
 283         delete                          fhPartonPhi2;   // ("hPartonPhi2","#phi_{Parton}",62,0.,3.1);
 284         delete                          fhEtaDiff2;     // ("hEtaDiff2","#eta_{jet}^{reco}-#eta_{jet}^{input}",100,-0.5,0.5);
 285         delete                          fhPhiDiff2;     // ("hPhiDiff2","#phi_{jet}^{reco}-#phi_{jet}^{input}",100,-0.5,0.5);
 286         delete                          fhEtaPhiSpread2;        // ("hEtaPhiSpread2","#eta - #phi Distribution
 287                                                         //of Reconstructed Jets",192,-0.7,0.7,288,pi/3,pi);
 288         delete                          fhNJets2;       // ("hNJets2","N Reconstructed jets",11,-0.5,10.5);
 289         delete                          fhJetEtSecond2; //("hJetEtSecond2","E_{T}^{reco}",250,0.,250.);
 290         delete                          fhJetEtRatio2;  //("hJetEtRatio2","Ratio of Second Highest to Highest",100,0,1);
 291         delete                          fhEtaPhiDist2;  //("hEtaPhiDist2","Angular Distance Between First and Second",100,0,3);
 292
 293         delete fhRecoBinPt;          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 294         delete fhRecoBinPtNoBg;          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 295         delete fhRecoBinPartonPt;    // ("fhRecoBinPartonPt","Input Pt Distribution",100,0,1);
 296         delete fhRecoBinJetEt;       // ("fhRecoJetEt","E_{T}^{reco}",250,0.,250.);
 297         delete fhRecoBinInputJetEt;  // ("fhRecoInputJetEt","E_{T}^{reco}",250,0.,250.);
 298
 299         delete fhJetPT ;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
 300         delete fhPartonPT ;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
 301         delete fhJetPT2 ;// new TH1F("hJetPT","P_{T} Distribution",200,0,200);
 302         delete fhPartonPT2 ;// new TH1F("hPartonPT","Parton P_{T} Distribution",200,0,1);
 303         delete fhRecoBinFragmFcn;//new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
 304         delete fhRecoBinFragmFcnNoBg;//new TH1F("fhRecoBinFragmFcn","Reconstructed Frag. Fcn",100,0,1);
 305         delete fhRecoBinPartonFragmFcn;// new TH1F("fhRecoBinPartonFragmFcn","Input Bin Fragm Fcn Distribution",100,0,1);
 306     delete fhJetInvE;// = new TH1F("fhJetInvE","#frac{1}{E_{R}}",100,0,1);
 307     delete fhJetInvE2;// = new TH1F("fhJetInvE","#frac{1}{E_{R}}",100,0,1);
 308
 309 }
 310
 311 void AliEMCALJetFinderPlots::FillFromOutput(AliEMCALJetFinderOutput* output, Float_t weight)
 312 {
 313         // Fill histograms from an output object
 314 if (!fInitialised) InitPlots();
 315   fOutput = output;
 316   if (!fOutput) return;
 317 // Make some temporary histograms to make sure we subtract
 318   // background properly
 319 /*
 320 tempFragmFcnNoBg =new TH1F("tempFragmFcnNoBg","Reconstructed Frag. Fcn With Background Removed",200,0,2);
 321 tempPtNoBg =new TH1F("tempPtNoBg","Reconstructed Frag. Fcn With Background Removed",200,0,200);
 322 tempFragmFcnNoBg->Fill(count/(jethighest->Energy()*fScaleFactor),-fhBackHisto->GetBinContent(count));
 323 tempPtNoBg->AddBinContent(count,-fhBackHisto->GetBinContent(count));
 324 */
 325
 326   fhNJets->Fill(fOutput->GetNJets());
 327   Bool_t doesJetMeetBinCriteria = 0;
 328   AliEMCALJet* jethighest=0;
 329   AliEMCALJet* jetsecond=0;
 330   //  Find Highest and Second Highest Jet
 331   // (NB!!!!!!!)  Pointing into the EMCAL!!!!!!
 332
 333 // =========================== All cases ===================================
 334
 335
 336         // I will make a little array of jet indices for which jets are in
 337         // the EMCAL then my counter can loop from 0 below - but it will
 338         // be the index of the array of applicable jets
 339
 340   Int_t appjet[4];
 341   Int_t numappjet=0;
 342
 343   for (Int_t appc=0;appc<fOutput->GetNJets();appc++)
 344   { // Check all jets for applicability
 345           Float_t eta = fOutput->GetJet(appc)->Eta();
 346           Float_t phi = fOutput->GetJet(appc)->Phi();
 347           if (eta > -0.7 && eta < 0.7 && phi > 1./3.*TMath::Pi() && phi < TMath::Pi())
 348           { // Then jet is applicable
 349                   appjet[numappjet]=appc;
 350                   numappjet++;
 351           }
 352   }
 353
 354
 355 Float_t et=0;
 356 if (numappjet >=1)
 357 {
 358         Float_t theta = 2.0*atan(exp(-fOutput->GetParton(0)->Eta()));
 359         et    = fOutput->GetParton(0)->Energy() * TMath::Sin(theta);
 360         if (fOutput->GetNJets()>1)
 361         {
 362                 for (Int_t counter = 0; counter<numappjet;counter++)
 363                 {
 364                         if (counter==0)
 365                         {
 366                                 jethighest = fOutput->GetJet(appjet[0]);
 367                                 jetsecond  = fOutput->GetJet(appjet[1]);
 368                         }
 369                         if (counter>0)
 370                         {
 371                                 Float_t energyhighest = jethighest->Energy();
 372                                 Float_t energysecond = jetsecond->Energy();
 373
 374                                 if ((fOutput->GetJet(appjet[counter]))->Energy()>energyhighest)
 375                                 {
 376                                         jetsecond=jethighest;
 377                                         jethighest=fOutput->GetJet(appjet[counter]);
 378                                 }else if ((fOutput->GetJet(appjet[counter]))->Energy()>energysecond)
 379                                 {
 380                                         jetsecond=fOutput->GetJet(appjet[counter]);
 381                                 }
 382
 383                         }
 384                 }
 385         }else
 386         {
 387                 Float_t eta = fOutput->GetJet(0)->Eta();
 388                 Float_t phi = fOutput->GetJet(0)->Phi();
 389                 if (eta > -0.7 && eta < 0.7 && phi > 1./3.*TMath::Pi() && phi < TMath::Pi())
 390                 { // Then jet is applicable
 391                         jethighest=fOutput->GetJet(0);
 392                         jetsecond=0;
 393                 }else
 394                 {
 395                         Error("FillFromOutput","There is only one jet and it isn't in the area of applicability");
 396                 }
 397
 398         }
 399         if ( 95.0 < jethighest->Energy()*fScaleFactor && jethighest->Energy()*fScaleFactor < 105.0  )
 400         {
 401                 doesJetMeetBinCriteria = 1;
 402                 fhRecoBinJetEt->Fill(jethighest->Energy()*fScaleFactor,weight);
 403                 fhRecoBinInputJetEt->Fill(et,weight);
 404         }
 405         fhInputOutput->Fill(et,jethighest->Energy());
 406
 407 }
 408
 409 if (numappjet > 1)
 410 {
 411 //========================= CASE 2 ===========================
 412   Int_t nPartons = fOutput->GetNPartons();
 413   fhNJets2->Fill(fOutput->GetNJets());
 414   AliEMCALParton* parton;
 415
 416   // End finding highest and second highest and continue
 417   fhJetEt2->Fill(jethighest->Energy()*fScaleFactor,weight);
 418   fhJetEtDiff2->Fill(jethighest->Energy()*fScaleFactor-et,weight);
 419   fhJetInvE2->Fill(1.0/(jethighest->Energy()*fScaleFactor),weight);
 420   fhJetEta2->Fill(jethighest->Eta(),weight  );
 421   fhJetPhi2->Fill(jethighest->Phi(),weight );
 422   if (nPartons ==0) return;
 423   parton = fOutput->GetParton(0);
 424
 425   fhPartonEta2->Fill( parton->Eta(),weight );
 426   fhPartonPhi2->Fill( parton->Phi(),weight );
 427
 428   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
 429   fhEtaDiff2->Fill( jethighest->Eta() - parton->Eta(),weight );
 430   fhPhiDiff2->Fill( jethighest->Phi() - parton->Phi(),weight );
 431   fhEtaPhiSpread2->Fill(jethighest->Eta()-parton->Eta(),jethighest->Phi() - parton->Phi());
 432   fhJetEtSecond2->Fill(jetsecond->Energy()*fScaleFactor,weight);
 433   fhJetEtRatio2->Fill(jetsecond->Energy()/jethighest->Energy(),weight);
 434   fhEtaPhiDist2->Fill( TMath::Sqrt((jethighest->Eta() - jetsecond->Eta())*(jethighest->Eta() - jetsecond->Eta())
 435                       + (jethighest->Phi() - jetsecond->Phi())*(jethighest->Phi() - jetsecond->Phi())     ),weight);
 436   /*
 437   Float_t *pt,*phi,*eta;
 438   Int_t *pdg;
 439   pt  = new Float_t[parton->GetNTracks()];
 440   eta = new Float_t[parton->GetNTracks()];
 441   phi = new Float_t[parton->GetNTracks()];
 442   pdg = new Int_t[parton->GetNTracks()];*/
 443
 444
 445
 446  Float_t pt[2000];
 447  Float_t eta[2000];
 448  Float_t phi[2000];
 449  Int_t pdg[2000];
 450
 451   parton->GetTrackList(pt,eta,phi,pdg);
 452   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ )
 453   {
 454       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
 455            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue;
 456       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 457       Double_t rt = 2.0*atan(exp(-parton->Eta()));
 458       Double_t ctt = cos(tt);
 459       Double_t crt = cos(rt);
 460       Double_t stt = sin(tt);
 461       Double_t srt = sin(rt);
 462       Double_t ctp = cos(phi[iT]);
 463       Double_t crp = cos(parton->Phi());
 464       Double_t stp = sin(phi[iT]);
 465       Double_t srp = sin(parton->Phi());
 466       //Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 467           Double_t alpha;
 468           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
 469           {
 470                   alpha = 0.0;
 471           }else
 472           {
 473                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 474           }
 475       Double_t correctp = pt[iT]/stt;
 476       fhPartonPL2->Fill( correctp*cos(alpha),weight);
 477       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +
 478       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 )
 479               fhPartonJT2->Fill( correctp*sin(alpha),weight);
 480       fhPartonPT2->Fill(correctp*sin(tt),weight);
 481       if (fNominalEnergy == 0.0) {
 482               fhPartonFragmFcn2->Fill(  correctp*sin(tt)/parton->Energy(),weight  );
 483       }else
 484       {
 485               fhPartonFragmFcn2->Fill(correctp*sin(tt)/fNominalEnergy,weight);
 486       }
 487       if (doesJetMeetBinCriteria)
 488       {
 489               fhRecoBinPartonPt->Fill(correctp*sin(tt),weight);          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 490       }
 491   }// loop over tracks
 492
 493 /*
 494   pt  = new Float_t[jet->NTracks()];
 495   eta = new Float_t[jet->NTracks()];
 496   phi = new Float_t[jet->NTracks()];
 497   pdg = new Int_t[jet->NTracks()];*/
 498   jethighest->TrackList(pt,eta,phi,pdg);
 499   for(Int_t iT=0; iT< jethighest->NTracks() ; iT++ )
 500   {
 501           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 502           Double_t rt = 2.0*atan(exp(-jethighest->Eta()));
 503           Double_t ctt = cos(tt);
 504           Double_t crt = cos(rt);
 505           Double_t stt = sin(tt);
 506           Double_t srt = sin(rt);
 507           Double_t ctp = cos(phi[iT]);
 508           Double_t crp = cos(jethighest->Phi());
 509           Double_t stp = sin(phi[iT]);
 510           Double_t srp = sin(jethighest->Phi());
 511          // Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 512           Double_t alpha;
 513           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
 514           {
 515                   alpha = 0.0;
 516           }else
 517           {
 518                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 519           }
 520           Double_t correctp = pt[iT]/stt;
 521           fhJetPL2->Fill( correctp*cos(alpha),weight);
 522           if ( (jethighest->Eta()-eta[iT])*(jethighest->Eta()-eta[iT]) +
 523                           (jethighest->Phi()-phi[iT])*(jethighest->Phi()-phi[iT]) < 0.2*0.2 )
 524                   fhJetJT2->Fill( correctp*sin(alpha),weight);
 525           fhJetPT2->Fill(correctp*sin(tt),weight);
 526           if (fNominalEnergy==0.0){
 527                   fhFragmFcn2->Fill(  correctp*sin(tt)/(jethighest->Energy()*fScaleFactor),weight  );
 528           } else
 529           {
 530                   fhFragmFcn2->Fill(  correctp*sin(tt)/fNominalEnergy,weight );
 531           }
 532           if (doesJetMeetBinCriteria)
 533           {
 534                   fhRecoBinPt->Fill(correctp*sin(tt),weight);          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 535                   fhRecoBinPtNoBg->Fill(correctp*sin(tt),weight);          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 536                   fhRecoBinFragmFcn->Fill(  correctp*sin(tt)/(jethighest->Energy()*fScaleFactor),weight  ); // This is the jet fragmentation function
 537                   fhRecoBinFragmFcnNoBg->Fill(  correctp*sin(tt)/(jethighest->Energy()*fScaleFactor),weight  ); // This is the jet fragmentation function
 538           }
 539   }// loop over tracks
 540   }
 541
 542   if (numappjet == 1)
 543   {
 544
 545 //========================= CASE 1 ===========================
 546   Int_t nPartons = fOutput->GetNPartons();
 547   if (fOutput->GetNJets()!=1) return;
 548   AliEMCALParton* parton;
 549   AliEMCALJet* jet;
 550   jet = jethighest;//fOutput->GetJet(0);
 551   fhJetEt->Fill(jet->Energy()*fScaleFactor,weight);
 552   fhJetEtDiff->Fill(jethighest->Energy()*fScaleFactor-et,weight);
 553   fhJetInvE->Fill(1.0/(jethighest->Energy()*fScaleFactor),weight);
 554   fhJetEta->Fill(jet->Eta(),weight  );
 555   fhJetPhi->Fill(jet->Phi(),weight );
 556   if (nPartons ==0) return;
 557   parton = fOutput->GetParton(0);
 558
 559   fhPartonEta->Fill( parton->Eta(),weight );
 560   fhPartonPhi->Fill( parton->Phi(),weight );
 561
 562   //hJetEtDiff->Fill( jet->Energy() - parton->Energy() );
 563   fhEtaDiff->Fill( jet->Eta() - parton->Eta(),weight );
 564   fhPhiDiff->Fill( jet->Phi() - parton->Phi(),weight );
 565   fhEtaPhiSpread->Fill(jet->Eta()-parton->Eta(),jet->Phi() - parton->Phi());
 566  /*
 567   Float_t *pt,*phi,*eta;
 568   Int_t *pdg;
 569   pt  = new Float_t[parton->GetNTracks()];
 570   eta = new Float_t[parton->GetNTracks()];
 571   phi = new Float_t[parton->GetNTracks()];
 572   pdg = new Int_t[parton->GetNTracks()];*/
 573
 574
 575
 576  Float_t pt[2000];
 577  Float_t eta[2000];
 578  Float_t phi[2000];
 579  Int_t pdg[2000];
 580
 581   parton->GetTrackList(pt,eta,phi,pdg);
 582   for(Int_t iT=0; iT< parton->GetNTracks() ; iT++ )
 583   {
 584       if ( (eta[iT]-parton->Eta())*(eta[iT]-parton->Eta())+
 585            (phi[iT]-parton->Phi())*(phi[iT]-parton->Phi()) >fConeRadius * fConeRadius ) continue;
 586       Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 587       Double_t rt = 2.0*atan(exp(-parton->Eta()));
 588       Double_t ctt = cos(tt);
 589       Double_t crt = cos(rt);
 590       Double_t stt = sin(tt);
 591       Double_t srt = sin(rt);
 592       Double_t ctp = cos(phi[iT]);
 593       Double_t crp = cos(parton->Phi());
 594       Double_t stp = sin(phi[iT]);
 595       Double_t srp = sin(parton->Phi());
 596      // Double_t alpha = acos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 597           Double_t alpha;
 598           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
 599           {
 600                   alpha = 0.0;
 601           }else
 602           {
 603                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);           }
 604       Double_t correctp = pt[iT]/stt;
 605       fhPartonPL->Fill( correctp*cos(alpha),weight);
 606       if ( (parton->Eta()-eta[iT])*(parton->Eta()-eta[iT]) +
 607       (parton->Phi()-phi[iT])*(parton->Phi()-phi[iT]) < 0.2*0.2 )
 608               fhPartonJT->Fill( correctp*sin(alpha),weight);
 609       if (fNominalEnergy == 0.0) {
 610               fhPartonFragmFcn->Fill(  correctp*sin(tt)/parton->Energy(),weight  );
 611           fhPartonPT->Fill(correctp*sin(tt),weight);
 612       }else
 613       {
 614               fhPartonFragmFcn->Fill(correctp*sin(tt)/fNominalEnergy,weight);
 615       }
 616       if (doesJetMeetBinCriteria)
 617       {
 618               fhRecoBinPartonPt->Fill(correctp*sin(tt),weight);          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 619       }
 620   }// loop over tracks
 621
 622 /*
 623   pt  = new Float_t[jet->NTracks()];
 624   eta = new Float_t[jet->NTracks()];
 625   phi = new Float_t[jet->NTracks()];
 626   pdg = new Int_t[jet->NTracks()];*/
 627   jet->TrackList(pt,eta,phi,pdg);
 628   for(Int_t iT=0; iT< jet->NTracks() ; iT++ )
 629   {
 630           Double_t tt = 2.0*atan(exp(-eta[iT])); // These names are short to make the equation manageable
 631           Double_t rt = 2.0*atan(exp(-jet->Eta()));
 632           Double_t ctt = cos(tt);
 633           Double_t crt = cos(rt);
 634           Double_t stt = sin(tt);
 635           Double_t srt = sin(rt);
 636           Double_t ctp = cos(phi[iT]);
 637           Double_t crp = cos(jet->Phi());
 638           Double_t stp = sin(phi[iT]);
 639           Double_t srp = sin(jet->Phi());
 640           //Info("plots","acos(%1.16f)\nstt=%f\npt=%f",crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt,stt,pt[iT]);
 641           //Info("plots","diff to 1 %f",1.0-crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 642           Double_t alpha;
 643           if (TMath::Abs(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt) > 0.9990)
 644           {
 645                   alpha = 0.0;
 646           }else
 647           {
 648                   alpha = TMath::ACos(crp*ctp*srt*stt+srp*stp*srt*stt+crt*ctt);
 649           }
 650           Double_t correctp = pt[iT]/stt;
 651           fhJetPL->Fill( correctp*cos(alpha),weight);
 652           if ( (jet->Eta()-eta[iT])*(jet->Eta()-eta[iT]) +
 653                (jet->Phi()-phi[iT])*(jet->Phi()-phi[iT]) < 0.2*0.2 )
 654                   fhJetJT->Fill( correctp*sin(alpha),weight);
 655           fhJetPT->Fill(correctp*sin(tt),weight);
 656           if (fNominalEnergy==0.0){
 657                   fhFragmFcn->Fill(  correctp*sin(tt)/(jet->Energy()*fScaleFactor),weight  ); // This is the jet fragmentation function
 658           } else
 659           {
 660                   fhFragmFcn->Fill(  correctp*sin(tt)/fNominalEnergy,weight );
 661           }
 662           if (doesJetMeetBinCriteria)
 663           {
 664                   fhRecoBinPt->Fill(correctp*sin(tt),weight);          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 665                   fhRecoBinPtNoBg->Fill(correctp*sin(tt),weight);          // ("fhRecoBinPt","Reconstructed Pt Distribution",100,0,1);
 666                   fhRecoBinFragmFcn->Fill(  correctp*sin(tt)/(jet->Energy()*fScaleFactor),weight  ); // This is the jet fragmentation function
 667                   fhRecoBinFragmFcnNoBg->Fill(  correctp*sin(tt)/(jet->Energy()*fScaleFactor),weight  ); // This is the jet fragmentation function
 668           }
 669   }// loop over tracks
 670   }
 671
 672 if (numappjet>=1 && fhBackHisto != 0 && doesJetMeetBinCriteria)
 673   {
 674     for (Int_t count=1;count<=100;count++)
 675     {
 676         fhRecoBinFragmFcnNoBg->Fill( ((Float_t)count)/(jethighest->Energy()*fScaleFactor),-fhBackHisto->GetBinContent(count)*weight);
 677         fhRecoBinPtNoBg->AddBinContent(count,-fhBackHisto->GetBinContent(count)*weight);
 678     }
 679   }
 680
 681
 682 }
 683
 684