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