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