PWGGA/CaloTrackCorrelations/AliAnaParticleJetFinderCorrelation.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 // Class for the analysis of particle (direct gamma) -jet (jet found with finder) correlations
  18 //*-- Author: Gustavo Conesa (LNF-INFN)
  19 //*-- Modified: Adam Matyja (INP PAN Cracow)
  20 //////////////////////////////////////////////////////////////////////////////
  21
  22
  23 // --- ROOT system ---
  24 #include "TH2F.h"
  25 #include "TClonesArray.h"
  26 #include "TClass.h"
  27 //#include "Riostream.h"
  28
  29 //---- AliRoot system ----
  30 #include "AliCaloTrackReader.h"
  31 #include "AliAODJet.h"
  32 #include "AliAnaParticleJetFinderCorrelation.h"
  33 #include "AliAODPWG4ParticleCorrelation.h"
  34 #include "AliVTrack.h"
  35 #include "AliAODCaloCluster.h"
  36 #include "AliAODEvent.h"
  37
  38 //jets
  39 #include "AliAODJetEventBackground.h"
  40 #include "TRandom2.h"
  41 //MC access
  42 #include "AliStack.h"
  43 #include "AliMCAnalysisUtils.h"
  44 #include "AliAODMCParticle.h"
  45 // --- Detectors ---
  46 #include "AliEMCALGeometry.h"
  47
  48 ClassImp(AliAnaParticleJetFinderCorrelation)
  49
  50 //________________________________________________________________________
  51 AliAnaParticleJetFinderCorrelation::AliAnaParticleJetFinderCorrelation() :
  52 AliAnaCaloTrackCorrBaseClass(),
  53   fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fRatioMaxCut(0.),  fRatioMinCut(0.),
  54   fConeSize(0.), fPtThresholdInCone(0.),fUseJetRefTracks(kTRUE),
  55   fMakeCorrelationInHistoMaker(kFALSE), fSelectIsolated(kTRUE),
  56   fJetConeSize(0.4),fJetMinPt(5),fJetAreaFraction(0.6),
  57 //fNonStandardJetFromReader(kTRUE),
  58   fJetBranchName("jets"),
  59   fBackgroundJetFromReader(kTRUE),
  60   fBkgJetBranchName("jets"),
  61   fGammaConeSize(0.3),fUseBackgroundSubtractionGamma(kFALSE),fSaveGJTree(kTRUE),
  62   fMostEnergetic(kFALSE),fMostOpposite(kTRUE), fUseHistogramJetBkg(kTRUE),
  63   fUseHistogramTracks(kTRUE),fUseHistogramJetTracks(kTRUE),fMCStudies(kFALSE),fGenerator(0),
  64   fMomentum(),
  65   fhDeltaEta(0), /*fhDeltaPhi(0),*/fhDeltaPhiCorrect(0),fhDeltaPhi0PiCorrect(0), fhDeltaPt(0), fhPtRatio(0), fhPt(0),
  66   fhFFz(0),fhFFxi(0),fhFFpt(0),fhNTracksInCone(0),
  67   fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetFFzCor(0),fhJetFFxiCor(0),
  68   fhBkgFFz(),fhBkgFFxi(),fhBkgFFpt(),fhBkgNTracksInCone(),fhBkgSumPtInCone(),fhBkgSumPtnTracksInCone(),//<<---new
  69   fhNjetsNgammas(0),fhCuts(0),
  70   fhDeltaEtaBefore(0),fhDeltaPhiBefore(0),fhDeltaPtBefore(0),fhPtRatioBefore(0),
  71   fhPtBefore(0),fhDeltaPhi0PiCorrectBefore(0),
  72   fhJetPtBefore(0),fhJetPt(0),fhJetPtMostEne(0),fhJetPhi(0),fhJetEta(0),fhJetEtaVsPt(0),
  73   fhJetPhiVsEta(0),fhJetEtaVsNpartInJet(0),fhJetEtaVsNpartInJetBkg(0),fhJetChBkgEnergyVsPt(0),fhJetChAreaVsPt(0),/*fhJetNjet(0),*/
  74   fhTrackPhiVsEta(0),fhTrackAveTrackPt(0),fhJetNjetOverPtCut(),
  75 /*fhJetChBkgEnergyVsPtEtaGt05(0),fhJetChBkgEnergyVsPtEtaLe05(0),fhJetChAreaVsPtEtaGt05(0),fhJetChAreaVsPtEtaLe05(0),*/
  76   fhJetChBkgEnergyVsArea(0),fhJetRhoVsPt(0),fhJetRhoVsCentrality(0),//fhJetBkgRho(0),
  77   fhJetNparticlesInJet(0),fhJetDeltaEtaDeltaPhi(0),fhJetDeltaEtaDeltaPhiAllTracks(0),
  78   fhJetAveTrackPt(0),fhJetNtracksInJetAboveThr(),fhJetRatioNTrkAboveToNTrk(),fhJetNtrackRatioMostEne(),
  79   fhJetNtrackRatioJet5GeV(),fhJetNtrackRatioLead5GeV(),
  80   fhBkgJetBackground(),fhBkgJetSigma(),fhBkgJetArea(),fhPhotonPtMostEne(0),
  81   fhPhotonAverageEnergy(0),fhPhotonRatioAveEneToMostEne(0),fhPhotonAverageEnergyMinus1(0),fhPhotonRatioAveEneMinus1ToMostEne(0),
  82   fhPhotonNgammaMoreAverageToNgamma(0),fhPhotonNgammaMoreAverageMinus1ToNgamma(0),fhPhotonNgammaOverPtCut(),
  83   fhPhotonBkgRhoVsNtracks(0),fhPhotonBkgRhoVsNclusters(0),fhPhotonBkgRhoVsCentrality(0),
  84   fhPhotonBkgRhoVsNcells(0),fhPhotonPt(0),fhPhotonPtCorrected(0),fhPhotonPtCorrectedZoom(0),fhPhotonPtDiff(0),
  85   fhPhotonPtDiffVsCentrality(0),fhPhotonPtDiffVsNcells(0),fhPhotonPtDiffVsNtracks(0),fhPhotonPtDiffVsNclusters(0),
  86   fhPhotonSumPtInCone(0),fhPhotonSumPtCorrectInCone(0),fhPhotonSumPtChargedInCone(0),
  87   fhSelectedJetPhiVsEta(0),fhSelectedJetChBkgEnergyVsPtJet(0),fhSelectedJetChAreaVsPtJet(0),fhSelectedJetNjet(0),fhSelectedNtracks(0),
  88   fhSelectedTrackPhiVsEta(0),fhCuts2(0),
  89   fhSelectedPhotonNLMVsPt(0),fhSelectedPhotonLambda0VsPt(0), fhRandomPhiEta(),
  90   fhMCPhotonCuts(0),fhMCPhotonPt(0),fhMCPhotonEtaPhi(0),fhMCJetOrigin(0),
  91   fhMCJetNPartVsPt(0),fhMCJetChNPartVsPt(0),fhMCJetNPart150VsPt(0),fhMCJetChNPart150VsPt(0),fhMCJetChNPart150ConeVsPt(0),
  92   fhMCJetRatioChFull(0),fhMCJetRatioCh150Ch(0),
  93   fhMCJetEtaPhi(0),fhMCJetChEtaPhi(0),fhMCJet150EtaPhi(0),fhMCJetCh150EtaPhi(0),fhMCJetCh150ConeEtaPhi(0),
  94 fTreeGJ     (0),
  95 fGamPt      (0),
  96 fGamLambda0 (0),
  97 fGamNLM     (0),
  98 fGamSumPtCh (0),
  99 fGamTime    (0),
 100 fGamNcells  (0),
 101 fGamEta     (0),
 102 fGamPhi     (0),
 103 fGamSumPtNeu(0),
 104 fGamNtracks (0),
 105 fGamNclusters(0),
 106 fGamAvEne   (0),
 107 fJetPhi     (0),
 108 fJetEta     (0),
 109 fJetPt      (0),
 110 fJetBkgChEne(0),
 111 fJetArea    (0),
 112 fJetNtracks (0),
 113 fJetNtracks1(0),
 114 fJetNtracks2(0),
 115 fJetRho(0),
 116 fEventNumber(0),
 117 fNtracks    (0),
 118 fZvertex    (0),
 119 fCentrality (0),
 120 fIso(0),
 121 fGamRho(0),
 122 fMCGamPt        (0),
 123 fMCGamEta       (0),
 124 fMCGamPhi       (0),
 125 fMCPartonType   (0),
 126 fMCJetPt        (0),
 127 fMCJetChPt      (0),
 128 fMCJet150Pt     (0),
 129 fMCJetCh150Pt   (0),
 130 fMCJetNPart     (0),
 131 fMCJetChNPart   (0),
 132 fMCJet150NPart  (0),
 133 fMCJetCh150NPart(0),
 134 fMCJetEta       (0),
 135 fMCJetPhi       (0),
 136 fMCJetChEta     (0),
 137 fMCJetChPhi     (0),
 138 fMCJet150Eta    (0),
 139 fMCJet150Phi    (0),
 140 fMCJetCh150Eta  (0),
 141 fMCJetCh150Phi  (0),
 142 fMCJetCh150ConePt(0),
 143 fMCJetCh150ConeNPart(0),
 144 fMCJetCh150ConeEta(0),
 145 fMCJetCh150ConePhi(0)
 146
 147 {
 148   //Default Ctor
 149   //printf("constructor\n");
 150
 151   //Initialize parameters
 152   InitParameters();
 153   for(Int_t i=0;i<10;i++){
 154     fhJetNjetOverPtCut[i] = 0;
 155     fhPhotonNgammaOverPtCut[i] = 0;
 156   }
 157   fGenerator = new TRandom2();
 158   fGenerator->SetSeed(0);
 159 }
 160
 161 //___________________________________________________________________
 162 AliAnaParticleJetFinderCorrelation::~AliAnaParticleJetFinderCorrelation(){
 163   delete fGenerator;
 164 }
 165
 166
 167 //___________________________________________________________________
 168 TList *  AliAnaParticleJetFinderCorrelation::GetCreateOutputObjects()
 169 {
 170   // Create histograms to be saved in output file and
 171   // store them in fOutputContainer
 172   //printf("GetCreateOutputObjects\n");
 173
 174   TList * outputContainer = new TList() ;
 175   outputContainer->SetName("ParticleJetFinderHistos") ;
 176
 177   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();
 178   //    Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
 179   //    Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
 180   Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();
 181   //    Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
 182   //    Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
 183   Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
 184   //    Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
 185 //      Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 186
 187 //  fhDeltaPhi  = new TH2F("DeltaPhi","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-6,4);
 188 //  fhDeltaPhi->SetYTitle("#Delta #phi");
 189 //  fhDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 190 //  outputContainer->Add(fhDeltaPhi);
 191
 192   fhDeltaPhiCorrect  = new TH2F("DeltaPhiCorrect","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5);
 193   fhDeltaPhiCorrect->SetYTitle("#Delta #phi");
 194   fhDeltaPhiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
 195   outputContainer->Add(fhDeltaPhiCorrect);
 196
 197   fhDeltaPhi0PiCorrect  = new TH2F("DeltaPhi0PiCorrect","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5);
 198   fhDeltaPhi0PiCorrect->SetYTitle("#Delta #phi");
 199   fhDeltaPhi0PiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
 200   outputContainer->Add(fhDeltaPhi0PiCorrect);
 201
 202
 203   fhDeltaEta  = new TH2F("DeltaEta","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2);
 204   fhDeltaEta->SetYTitle("#Delta #eta");
 205   fhDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 206   outputContainer->Add(fhDeltaEta);
 207
 208   fhDeltaPt  = new TH2F("DeltaPt","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,150,-50,100);
 209   fhDeltaPt->SetYTitle("#Delta p_{T}");
 210   fhDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
 211   outputContainer->Add(fhDeltaPt);
 212
 213   fhPtRatio  = new TH2F("PtRatio","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.);
 214   fhPtRatio->SetYTitle("ratio");
 215   fhPtRatio->SetXTitle("p_{T trigger} (GeV/c)");
 216   outputContainer->Add(fhPtRatio);
 217
 218   fhPt  = new TH2F("Pt","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 219   fhPt->SetYTitle("p_{T jet}(GeV/c)");
 220   fhPt->SetXTitle("p_{T trigger} (GeV/c)");
 221   outputContainer->Add(fhPt);
 222
 223   fhFFz  = new TH2F("FFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,2);
 224   fhFFz->SetYTitle("z");
 225   fhFFz->SetXTitle("p_{T trigger}");
 226   outputContainer->Add(fhFFz) ;
 227
 228   fhFFxi  = new TH2F("FFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,10.);
 229   fhFFxi->SetYTitle("#xi");
 230   fhFFxi->SetXTitle("p_{T trigger}");
 231   outputContainer->Add(fhFFxi) ;
 232
 233   fhFFpt  = new TH2F("FFpt","p_{T i charged} vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,50.);
 234   fhFFpt->SetYTitle("p_{T charged hadron}");
 235   fhFFpt->SetXTitle("p_{T trigger}");
 236   outputContainer->Add(fhFFpt) ;
 237
 238   fhNTracksInCone  = new TH2F("NTracksInCone","Number of tracks in cone vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,150.);
 239   fhNTracksInCone->SetYTitle("p_{T charged hadron}");
 240   fhNTracksInCone->SetXTitle("p_{T trigger}");
 241   outputContainer->Add(fhNTracksInCone) ;
 242
 243   //FF according to jet axis
 244   fhJetFFz  = new TH2F("JetFFz","z = p_{T i charged}/p_{T jet} vs p_{T jet}",nptbins,ptmin,ptmax,200,0.,2);
 245   fhJetFFz->SetYTitle("z");
 246   fhJetFFz->SetXTitle("p_{T jet}");
 247   outputContainer->Add(fhJetFFz) ;
 248
 249   fhJetFFxi  = new TH2F("JetFFxi","#xi = ln(p_{T jet}/p_{T i charged}) vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,10.);
 250   fhJetFFxi->SetYTitle("#xi");
 251   fhJetFFxi->SetXTitle("p_{T jet}");
 252   outputContainer->Add(fhJetFFxi) ;
 253
 254   fhJetFFpt  = new TH2F("JetFFpt","p_{T i charged} vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,50.);
 255   fhJetFFpt->SetYTitle("p_{T charged hadron}");
 256   fhJetFFpt->SetXTitle("p_{T jet}");
 257   outputContainer->Add(fhJetFFpt) ;
 258
 259   fhJetFFzCor  = new TH2F("JetFFzCor","z = -cos(#alpha(jet,trig))*p_{T i charged}/p_{T jet} vs p_{T jet}",nptbins,ptmin,ptmax,200,0.,2);
 260   fhJetFFzCor->SetYTitle("z");
 261   fhJetFFzCor->SetXTitle("p_{T jet}");
 262   outputContainer->Add(fhJetFFzCor) ;
 263
 264   fhJetFFxiCor  = new TH2F("JetFFxiCor","#xi = ln(p_{T jet}/(-cos(#alpha(jet,trig))*p_{T i charged})) vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,10.);
 265   fhJetFFxiCor->SetYTitle("#xi");
 266   fhJetFFxiCor->SetXTitle("p_{T jet}");
 267   outputContainer->Add(fhJetFFxiCor) ;
 268
 269
 270   //background FF
 271   fhBkgFFz[0]  = new TH2F("BkgFFzRC",  "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg RC"  ,nptbins,ptmin,ptmax,200,0.,2);
 272   fhBkgFFz[1]  = new TH2F("BkgFFzPCG", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg PCG" ,nptbins,ptmin,ptmax,200,0.,2);
 273   fhBkgFFz[2]  = new TH2F("BkgFFzPCJ", "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg PCJ" ,nptbins,ptmin,ptmax,200,0.,2);
 274   fhBkgFFz[3]  = new TH2F("BkgFFzMP",  "z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg MP"  ,nptbins,ptmin,ptmax,200,0.,2);
 275   fhBkgFFz[4]  = new TH2F("BkgFFzTest","z = p_{T i charged}/p_{T trigger} vs p_{T trigger} Bkg Test",nptbins,ptmin,ptmax,200,0.,2);
 276   for(Int_t i=0;i<5;i++){
 277     fhBkgFFz[i]->SetYTitle("z");
 278     fhBkgFFz[i]->SetXTitle("p_{T trigger}");
 279     outputContainer->Add(fhBkgFFz[i]) ;
 280   }
 281
 282   fhBkgFFxi[0]  = new TH2F("BkgFFxiRC",  "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg RC",  nptbins,ptmin,ptmax,100,0.,10.);
 283   fhBkgFFxi[1]  = new TH2F("BkgFFxiPCG", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg PCG", nptbins,ptmin,ptmax,100,0.,10.);
 284   fhBkgFFxi[2]  = new TH2F("BkgFFxiPCJ", "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg PCJ", nptbins,ptmin,ptmax,100,0.,10.);
 285   fhBkgFFxi[3]  = new TH2F("BkgFFxiMP",  "#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg MP",  nptbins,ptmin,ptmax,100,0.,10.);
 286   fhBkgFFxi[4]  = new TH2F("BkgFFxiTest","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger} Bkg Test",nptbins,ptmin,ptmax,100,0.,10.);
 287   for(Int_t i=0;i<5;i++){
 288     fhBkgFFxi[i]->SetYTitle("#xi");
 289     fhBkgFFxi[i]->SetXTitle("p_{T trigger}");
 290     outputContainer->Add(fhBkgFFxi[i]) ;
 291   }
 292
 293   fhBkgFFpt[0]  = new TH2F("BkgFFptRC",  "p_{T i charged} vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,50.);
 294   fhBkgFFpt[1]  = new TH2F("BkgFFptPCG", "p_{T i charged} vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,50.);
 295   fhBkgFFpt[2]  = new TH2F("BkgFFptPCJ", "p_{T i charged} vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,50.);
 296   fhBkgFFpt[3]  = new TH2F("BkgFFptMP",  "p_{T i charged} vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,50.);
 297   fhBkgFFpt[4]  = new TH2F("BkgFFptTest","p_{T i charged} vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,50.);
 298   for(Int_t i=0;i<5;i++){
 299     fhBkgFFpt[i]->SetYTitle("p_{T charged hadron}");
 300     fhBkgFFpt[i]->SetXTitle("p_{T trigger}");
 301     outputContainer->Add(fhBkgFFpt[i]) ;
 302   }
 303
 304   fhBkgNTracksInCone[0]  = new TH2F("BkgNTracksInConeRC",  "Number of tracks in cone vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,150.);
 305   fhBkgNTracksInCone[1]  = new TH2F("BkgNTracksInConePCG", "Number of tracks in cone vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,150.);
 306   fhBkgNTracksInCone[2]  = new TH2F("BkgNTracksInConePCJ", "Number of tracks in cone vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,150.);
 307   fhBkgNTracksInCone[3]  = new TH2F("BkgNTracksInConeMP",  "Number of tracks in cone vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,150.);
 308   fhBkgNTracksInCone[4]  = new TH2F("BkgNTracksInConeTest","Number of tracks in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,150.);
 309   for(Int_t i=0;i<5;i++){
 310     fhBkgNTracksInCone[i]->SetYTitle("Number of tracks");
 311     fhBkgNTracksInCone[i]->SetXTitle("p_{T trigger}");
 312     outputContainer->Add(fhBkgNTracksInCone[i]) ;
 313   }
 314
 315   fhBkgSumPtInCone[0]  = new TH2F("BkgSumPtInConeRC",  "Sum P_{T} in cone vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,100.);
 316   fhBkgSumPtInCone[1]  = new TH2F("BkgSumPtInConePCG", "Sum P_{T} in cone vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,100.);
 317   fhBkgSumPtInCone[2]  = new TH2F("BkgSumPtInConePCJ", "Sum P_{T} in cone vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,100.);
 318   fhBkgSumPtInCone[3]  = new TH2F("BkgSumPtInConeMP",  "Sum P_{T} in cone vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,100.);
 319   fhBkgSumPtInCone[4]  = new TH2F("BkgSumPtInConeTest","Sum P_{T} in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,100.);
 320   for(Int_t i=0;i<5;i++){
 321     fhBkgSumPtInCone[i]->SetYTitle("Sum P_{T}");
 322     fhBkgSumPtInCone[i]->SetXTitle("p_{T trigger}");
 323     outputContainer->Add(fhBkgSumPtInCone[i]) ;
 324   }
 325
 326   fhBkgSumPtnTracksInCone[0]  = new TH2F("BkgSumPtnTracksInConeRC",  "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,20.);
 327   fhBkgSumPtnTracksInCone[1]  = new TH2F("BkgSumPtnTracksInConePCG", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,20.);
 328   fhBkgSumPtnTracksInCone[2]  = new TH2F("BkgSumPtnTracksInConePCJ", "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,20.);
 329   fhBkgSumPtnTracksInCone[3]  = new TH2F("BkgSumPtnTracksInConeMP",  "Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,20.);
 330   fhBkgSumPtnTracksInCone[4]  = new TH2F("BkgSumPtnTracksInConeTest","Sum p_{T} / Number of tracks in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,20.);
 331   for(Int_t i=0;i<5;i++){
 332     fhBkgSumPtnTracksInCone[i]->SetYTitle("Sum p_{T}/Number of tracks");
 333     fhBkgSumPtnTracksInCone[i]->SetXTitle("p_{T trigger}");
 334     outputContainer->Add(fhBkgSumPtnTracksInCone[i]) ;
 335   }
 336
 337
 338   //temporary histograms
 339   fhNjetsNgammas  = new TH2F("NjetsNgammas"," Number of jets vs number of gammas in event",20,0.,100.,10,0.,80.);
 340   fhNjetsNgammas->SetYTitle("Number of gammas");
 341   fhNjetsNgammas->SetXTitle("Number of jets");
 342   outputContainer->Add(fhNjetsNgammas) ;
 343
 344   fhCuts  = new TH1F("Cuts"," Cuts",10,0.,10.);
 345   fhCuts->SetYTitle("Counts");
 346   fhCuts->SetXTitle("Cut number");
 347   outputContainer->Add(fhCuts) ;
 348
 349   fhDeltaPhiBefore  = new TH2F("DeltaPhiBefore","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5);
 350   fhDeltaPhiBefore->SetYTitle("#Delta #phi");
 351   fhDeltaPhiBefore->SetXTitle("p_{T trigger} (GeV/c)");
 352   outputContainer->Add(fhDeltaPhiBefore);
 353
 354   fhDeltaEtaBefore  = new TH2F("DeltaEtaBefore","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2);
 355   fhDeltaEtaBefore->SetYTitle("#Delta #eta");
 356   fhDeltaEtaBefore->SetXTitle("p_{T trigger} (GeV/c)");
 357   outputContainer->Add(fhDeltaEtaBefore);
 358
 359   fhDeltaPtBefore  = new TH2F("DeltaPtBefore","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-50,50);
 360   fhDeltaPtBefore->SetYTitle("#Delta p_{T}");
 361   fhDeltaPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
 362   outputContainer->Add(fhDeltaPtBefore);
 363
 364   fhPtRatioBefore  = new TH2F("PtRatioBefore","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.);
 365   fhPtRatioBefore->SetYTitle("ratio");
 366   fhPtRatioBefore->SetXTitle("p_{T trigger} (GeV/c)");
 367   outputContainer->Add(fhPtRatioBefore);
 368
 369   fhPtBefore  = new TH2F("PtBefore","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 370   fhPtBefore->SetYTitle("p_{T jet}(GeV/c)");
 371   fhPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
 372   outputContainer->Add(fhPtBefore);
 373
 374   fhDeltaPhi0PiCorrectBefore  = new TH2F("DeltaPhi0PiCorrectBefore","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5);
 375   fhDeltaPhi0PiCorrectBefore->SetYTitle("#Delta #phi");
 376   fhDeltaPhi0PiCorrectBefore->SetXTitle("p_{T trigger} (GeV/c)");
 377   outputContainer->Add(fhDeltaPhi0PiCorrectBefore);
 378
 379   //temporary jet histograms
 380   fhJetPtBefore            = new TH1F("JetPtBefore","JetPtBefore",150,-50,100);
 381   fhJetPtBefore->SetYTitle("Counts");
 382   fhJetPtBefore->SetXTitle("p_{T jet}(GeV/c)");
 383   outputContainer->Add(fhJetPtBefore) ;
 384
 385   fhJetPt            = new TH1F("JetPt","JetPt",150,-50,100);
 386   fhJetPt->SetYTitle("Counts");
 387   fhJetPt->SetXTitle("p_{T jet}(GeV/c)");
 388   outputContainer->Add(fhJetPt) ;
 389
 390   fhJetPtMostEne            = new TH1F("JetPtMostEne","JetPtMostEne",150,0,150);
 391   fhJetPtMostEne->SetYTitle("Counts");
 392   fhJetPtMostEne->SetXTitle("p_{T jet}(GeV/c)");
 393   outputContainer->Add(fhJetPtMostEne) ;
 394
 395   fhJetPhi           = new TH1F("JetPhi","JetPhi",130,0,6.5);
 396   fhJetPhi->SetYTitle("Counts");
 397   fhJetPhi->SetXTitle("#phi_{jet}");
 398   outputContainer->Add(fhJetPhi) ;
 399
 400   fhJetEta           = new TH1F("JetEta","JetEta",100,-1,1);
 401   fhJetEta->SetYTitle("Counts");
 402   fhJetEta->SetXTitle("#eta_{jet}");
 403   outputContainer->Add(fhJetEta) ;
 404
 405   fhJetEtaVsPt      = new TH2F("JetEtaVsPt","JetEtaVsPt",100,0,100,50,-1,1);
 406   fhJetEtaVsPt->SetYTitle("#eta_{jet}");
 407   fhJetEtaVsPt->SetXTitle("p_{T,jet}(GeV/c)");
 408   outputContainer->Add(fhJetEtaVsPt) ;
 409
 410   fhJetPhiVsEta      = new TH2F("JetPhiVsEta","JetPhiVsEta",65,0,6.5,50,-1,1);
 411   fhJetPhiVsEta->SetYTitle("#eta_{jet}");
 412   fhJetPhiVsEta->SetXTitle("#phi_{jet}");
 413   outputContainer->Add(fhJetPhiVsEta) ;
 414
 415   fhJetEtaVsNpartInJet= new TH2F("JetEtaVsNpartInJet","JetEtaVsNpartInJet",50,-1,1,100,0.,200.);
 416   fhJetEtaVsNpartInJet->SetYTitle("N_{tracks-in-jet}");
 417   fhJetEtaVsNpartInJet->SetXTitle("#eta_{jet}");
 418   outputContainer->Add(fhJetEtaVsNpartInJet) ;
 419
 420   fhJetEtaVsNpartInJetBkg= new TH2F("JetEtaVsNpartInJetBkg","JetEtaVsNpartInJetBkg",50,-1,1,100,0.,200.);
 421   fhJetEtaVsNpartInJetBkg->SetYTitle("N_{tracks-in-jet}");
 422   fhJetEtaVsNpartInJetBkg->SetXTitle("#eta_{jet}");
 423   outputContainer->Add(fhJetEtaVsNpartInJetBkg) ;
 424
 425   fhJetChBkgEnergyVsPt = new TH2F("JetBkgChEnergyVsPt","JetBkgChEnergyVsPt",100,0,100,90,0,90);
 426   fhJetChBkgEnergyVsPt->SetYTitle("Jet Bkg Energy (GeV)");
 427   fhJetChBkgEnergyVsPt->SetXTitle("p_{T,jet} (GeV/c)");
 428   outputContainer->Add(fhJetChBkgEnergyVsPt);
 429
 430   fhJetChAreaVsPt      = new TH2F("JetChAreaVsPt","JetChAreaVsPt",100,0,100,50,0,1);
 431   fhJetChAreaVsPt->SetYTitle("Jet Area");
 432   fhJetChAreaVsPt->SetXTitle("p_{T,jet} (GeV/c)");
 433   outputContainer->Add(fhJetChAreaVsPt);
 434
 435   if(IsHistogramTracks()){
 436     fhTrackPhiVsEta      = new TH2F("TrackPhiVsEta","TrackPhiVsEta",65,0,6.5,50,-1,1);
 437     fhTrackPhiVsEta->SetYTitle("#eta_{track}");
 438     fhTrackPhiVsEta->SetXTitle("#phi_{track}");
 439     outputContainer->Add(fhTrackPhiVsEta) ;
 440
 441     fhTrackAveTrackPt      = new TH1F("TrackAveTrackPt","TrackAveTrackPt",45,0,1.5);
 442     fhTrackAveTrackPt->SetYTitle("Counts");
 443     fhTrackAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
 444     outputContainer->Add(fhTrackAveTrackPt);
 445
 446   }//end of IsHistogramTracks()
 447
 448   for(Int_t i=0;i<10;i++){
 449     fhJetNjetOverPtCut[i]      = new TH1F(Form("JetNjetOverPtCut%d", i),Form("JetNjetOverPtCut%d", i),100,0,100);
 450     fhJetNjetOverPtCut[i]->SetYTitle("Counts");
 451     fhJetNjetOverPtCut[i]->SetXTitle("N_{jets} over threshold");
 452     outputContainer->Add(fhJetNjetOverPtCut[i]);
 453   }
 454
 455   fhJetChBkgEnergyVsArea = new TH2F("JetBkgChEnergyVsArea","JetBkgChEnergyVsArea",100,0,100,70,0,0.7);
 456   fhJetChBkgEnergyVsArea->SetXTitle("Jet Bkg Energy (GeV)");
 457   fhJetChBkgEnergyVsArea->SetYTitle("Area");
 458   outputContainer->Add(fhJetChBkgEnergyVsArea);
 459
 460   fhJetRhoVsPt           = new TH2F("JetRhoVsPt","JetRhoVsPt",100,0,100,100,0,150);
 461   fhJetRhoVsPt->SetYTitle("Rho");
 462   fhJetRhoVsPt->SetXTitle("p_{T,jet} (GeV/c)");
 463   outputContainer->Add(fhJetRhoVsPt);
 464
 465   if(IsHistogramJetBkg()){
 466     fhJetRhoVsCentrality           = new TH2F("JetRhoVsCentrality","JetRhoVsCentrality",100,0,100,100,0,200);
 467     fhJetRhoVsCentrality->SetYTitle("Rho");
 468     fhJetRhoVsCentrality->SetXTitle("Centrality");
 469     outputContainer->Add(fhJetRhoVsCentrality);
 470   }
 471
 472   fhJetNparticlesInJet           = new TH1F("JetNparticlesInJet","JetNparticlesInJet",100,0,200);
 473   fhJetNparticlesInJet->SetXTitle("N^{particles}");
 474   fhJetNparticlesInJet->SetYTitle("N^{jets}");
 475   outputContainer->Add(fhJetNparticlesInJet);
 476
 477   fhJetDeltaEtaDeltaPhi      = new TH2F("JetDeltaEtaDeltaPhi","#Delta #eta^{jet-track} vs. #Delta #phi^{jet-track} for jet tracks",100,-0.8,0.8,100,-0.8,0.8);
 478   fhJetDeltaEtaDeltaPhi->SetXTitle("#Delta #eta^{jet-track}");
 479   fhJetDeltaEtaDeltaPhi->SetYTitle("#Delta #phi^{jet-track}");
 480   outputContainer->Add(fhJetDeltaEtaDeltaPhi );
 481
 482
 483   fhJetDeltaEtaDeltaPhiAllTracks      = new TH2F("JetDeltaEtaDeltaPhiAllTracks","#Delta #eta^{jet-track} vs. #Delta #phi^{jet-track} for all tracks",100,-3.2,3.2,100,-3.2,3.2);
 484   fhJetDeltaEtaDeltaPhiAllTracks->SetXTitle("#Delta #eta^{jet-track}");
 485   fhJetDeltaEtaDeltaPhiAllTracks->SetYTitle("#Delta #phi^{jet-track}");
 486   outputContainer->Add(fhJetDeltaEtaDeltaPhiAllTracks);
 487
 488
 489   if(IsHistogramJetTracks()){
 490     fhJetAveTrackPt           = new TH1F("JetAveTrackPt","JetAveTrackPt",45,0.,1.5);
 491     fhJetAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
 492     fhJetAveTrackPt->SetYTitle("Counts");
 493     outputContainer->Add(fhJetAveTrackPt);
 494
 495     for(Int_t i=0;i<6;i++){
 496       if(i==0) fhJetNtracksInJetAboveThr[i]      = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,200);
 497       else fhJetNtracksInJetAboveThr[i]      = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,100);
 498       fhJetNtracksInJetAboveThr[i]->SetXTitle("p_{T,jet} (GeV/c)");
 499       fhJetNtracksInJetAboveThr[i]->SetYTitle("N_{tracks} over threshold");
 500       outputContainer->Add(fhJetNtracksInJetAboveThr[i]);
 501     }
 502
 503     for(Int_t i=0;i<5;i++){
 504       fhJetRatioNTrkAboveToNTrk[i]      = new TH2F(Form("JetRatioNTrkAboveToNTrk%d", i),Form("JetRatioNTrkAboveToNTrk%d", i),100,0,100,40,0,1);
 505       fhJetRatioNTrkAboveToNTrk[i]->SetXTitle("p_{T,jet} (GeV/c)");
 506       fhJetRatioNTrkAboveToNTrk[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
 507       outputContainer->Add(fhJetRatioNTrkAboveToNTrk[i]);
 508
 509       fhJetNtrackRatioMostEne[i]      = new TH2F(Form("JetNtrackRatioMostEne%d", i),Form("JetNtrackRatioMostEne%d", i),100,0,100,40,0,1);
 510       fhJetNtrackRatioMostEne[i]->SetXTitle("p_{T,jet} (GeV/c)");
 511       fhJetNtrackRatioMostEne[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
 512       outputContainer->Add(fhJetNtrackRatioMostEne[i]);
 513
 514       fhJetNtrackRatioJet5GeV[i]      = new TH2F(Form("JetNtrackRatioJet5GeV%d", i),Form("JetNtrackRatioJet5GeV%d", i),100,0,100,40,0,1);
 515       fhJetNtrackRatioJet5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
 516       fhJetNtrackRatioJet5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
 517       outputContainer->Add(fhJetNtrackRatioJet5GeV[i]);
 518
 519       fhJetNtrackRatioLead5GeV[i]      = new TH2F(Form("JetNtrackRatioLead5GeV%d", i),Form("JetNtrackRatioLead5GeV%d", i),100,0,100,40,0,1);
 520       fhJetNtrackRatioLead5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
 521       fhJetNtrackRatioLead5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
 522       outputContainer->Add(fhJetNtrackRatioLead5GeV[i]);
 523     }
 524   }//end of if IsHistogramJetTracks
 525
 526   //temporary background jets histograms
 527   if(IsHistogramJetBkg()){
 528     for(Int_t i=0;i<4;i++){
 529       fhBkgJetBackground[i]      = new TH1F(Form("BkgJetBackground%d", i),Form("BkgJetBackground%d", i),100,0,200);
 530       fhBkgJetBackground[i]->SetXTitle("<#rho> (GeV/c)");
 531       fhBkgJetBackground[i]->SetYTitle("Counts");
 532       outputContainer->Add(fhBkgJetBackground[i]);
 533
 534       fhBkgJetSigma[i]      = new TH1F(Form("BkgJetSigma%d", i),Form("BkgJetSigma%d", i),100,0,50);
 535       fhBkgJetSigma[i]->SetXTitle("#sigma (GeV/c)");
 536       fhBkgJetSigma[i]->SetYTitle("Counts");
 537       outputContainer->Add(fhBkgJetSigma[i]);
 538
 539       fhBkgJetArea[i]      = new TH1F(Form("BkgJetArea%d", i),Form("BkgJetArea%d", i),100,0,1);
 540       fhBkgJetArea[i]->SetXTitle("<A>");
 541       fhBkgJetArea[i]->SetYTitle("Counts");
 542       outputContainer->Add(fhBkgJetArea[i]);
 543     }
 544   }
 545
 546   //temporary photon histograms
 547   fhPhotonPtMostEne = new TH1F("PhotonPtMostEne","PhotonPtMostEne",100,0,100);
 548   fhPhotonPtMostEne->SetYTitle("Counts");
 549   fhPhotonPtMostEne->SetXTitle("p_{T,#gamma} (GeV/c)");
 550   outputContainer->Add(fhPhotonPtMostEne);
 551
 552 //  fhPhotonIndexMostEne = new TH1F("PhotonIndexMostEne","PhotonIndexMostEne",100,0,100);
 553 //  fhPhotonIndexMostEne->SetYTitle("Counts");
 554 //  fhPhotonIndexMostEne->SetXTitle("Index");
 555 //  outputContainer->Add(fhPhotonIndexMostEne);
 556
 557   fhPhotonAverageEnergy = new TH1F("PhotonAverageEnergy","PhotonAverageEnergy",100,0,10);
 558   fhPhotonAverageEnergy->SetYTitle("Counts");
 559   fhPhotonAverageEnergy->SetXTitle("p_{T,#gamma} (GeV/c)");
 560   outputContainer->Add(fhPhotonAverageEnergy);
 561
 562   fhPhotonRatioAveEneToMostEne = new TH1F("PhotonRatioAveEneToMostEne","PhotonRatioAveEneToMostEne",100,0,1);
 563   fhPhotonRatioAveEneToMostEne->SetYTitle("Counts");
 564   fhPhotonRatioAveEneToMostEne->SetXTitle("Ratio");
 565   outputContainer->Add(fhPhotonRatioAveEneToMostEne);
 566
 567   fhPhotonAverageEnergyMinus1 = new TH1F("PhotonAverageEnergyMinus1","PhotonAverageEnergyMinus1",100,0,10);
 568   fhPhotonAverageEnergyMinus1->SetYTitle("Counts");
 569   fhPhotonAverageEnergyMinus1->SetXTitle("p_{T,#gamma} (GeV/c)");
 570   outputContainer->Add(fhPhotonAverageEnergyMinus1);
 571
 572   fhPhotonRatioAveEneMinus1ToMostEne = new TH1F("PhotonRatioAveEneMinus1ToMostEne","PhotonRatioAveEneMinus1ToMostEne",100,0,1);
 573   fhPhotonRatioAveEneMinus1ToMostEne->SetYTitle("Counts");
 574   fhPhotonRatioAveEneMinus1ToMostEne->SetXTitle("Ratio");
 575   outputContainer->Add(fhPhotonRatioAveEneMinus1ToMostEne);
 576
 577   fhPhotonNgammaMoreAverageToNgamma = new TH1F("PhotonNgammaMoreAverageToNgamma","PhotonNgammaMoreAverageToNgamma",100,0,1);
 578   fhPhotonNgammaMoreAverageToNgamma->SetYTitle("Counts");
 579   fhPhotonNgammaMoreAverageToNgamma->SetXTitle("Ratio");
 580   outputContainer->Add(fhPhotonNgammaMoreAverageToNgamma);
 581
 582   fhPhotonNgammaMoreAverageMinus1ToNgamma = new TH1F("PhotonNgammaMoreAverageMinus1ToNgamma","PhotonNgammaMoreAverageMinus1ToNgamma",100,0,1);
 583   fhPhotonNgammaMoreAverageMinus1ToNgamma->SetYTitle("Counts");
 584   fhPhotonNgammaMoreAverageMinus1ToNgamma->SetXTitle("Ratio");
 585   outputContainer->Add(fhPhotonNgammaMoreAverageMinus1ToNgamma);
 586
 587   for(Int_t i=0;i<10;i++){
 588     fhPhotonNgammaOverPtCut[i]      = new TH1F(Form("PhotonNgammaOverPtCut%d",i),Form("PhotonNgammaOverPtCut%d",i),100,0,100);
 589     fhPhotonNgammaOverPtCut[i]->SetYTitle("Counts");
 590     fhPhotonNgammaOverPtCut[i]->SetXTitle("N_{#gamma} over threshold");
 591     outputContainer->Add(fhPhotonNgammaOverPtCut[i]);
 592   }
 593
 594   fhPhotonBkgRhoVsNtracks = new TH2F("PhotonBkgRhoVsNtracks","PhotonBkgRhoVsNtracks",200,0,2500,75,0,1.5);
 595   //fhPhotonBkgRhoVsNtracks->SetXTitle("Counts");
 596   fhPhotonBkgRhoVsNtracks->SetXTitle("Ntracks");
 597   fhPhotonBkgRhoVsNtracks->SetYTitle("Rho");
 598   outputContainer->Add(fhPhotonBkgRhoVsNtracks);
 599
 600   fhPhotonBkgRhoVsNclusters = new TH2F("PhotonBkgRhoVsNclusters","PhotonBkgRhoVsNclusters",50,0,100,75,0,1.5);
 601   fhPhotonBkgRhoVsNclusters->SetXTitle("Nclusters");
 602   fhPhotonBkgRhoVsNclusters->SetYTitle("Rho");
 603   outputContainer->Add(fhPhotonBkgRhoVsNclusters);
 604
 605   fhPhotonBkgRhoVsCentrality = new TH2F("PhotonBkgRhoVsCentrality","PhotonBkgRhoVsCentrality",100,0,100,75,0,1.5);
 606   fhPhotonBkgRhoVsCentrality->SetXTitle("Centrality");
 607   fhPhotonBkgRhoVsCentrality->SetYTitle("Rho");
 608   outputContainer->Add(fhPhotonBkgRhoVsCentrality);
 609
 610   fhPhotonBkgRhoVsNcells = new TH2F("PhotonBkgRhoVsNcells","PhotonBkgRhoVsNcells",100,0,200,75,0,1.5);
 611   fhPhotonBkgRhoVsNcells->SetXTitle("N_{cells}");
 612   fhPhotonBkgRhoVsNcells->SetYTitle("Rho");
 613   outputContainer->Add(fhPhotonBkgRhoVsNcells);
 614
 615   fhPhotonPt = new TH1F("PhotonPt","PhotonPt",220,-10,100);
 616   fhPhotonPt->SetXTitle("p_{T,#gamma} (GeV/c)");
 617   fhPhotonPt->SetYTitle("Counts");
 618   outputContainer->Add(fhPhotonPt);
 619
 620   fhPhotonPtCorrected = new TH1F("PhotonPtCorrected","PhotonPtCorrected",220,-10,100);
 621   fhPhotonPtCorrected->SetXTitle("p_{T,#gamma} (GeV/c)");
 622   fhPhotonPtCorrected->SetYTitle("Counts");
 623   outputContainer->Add(fhPhotonPtCorrected);
 624
 625   fhPhotonPtDiff = new TH1F("PhotonPtDiff","PhotonPtDiff",50,0,10);
 626   fhPhotonPtDiff->SetXTitle("p_{T,#gamma} (GeV/c)");
 627   fhPhotonPtDiff->SetYTitle("Counts");
 628   outputContainer->Add(fhPhotonPtDiff);
 629
 630   fhPhotonPtDiffVsNtracks = new TH2F("PhotonPtDiffVsNtracks","PhotonPtDiffVsNtracks",200,0,2500,50,0,5);
 631   fhPhotonPtDiffVsNtracks->SetXTitle("N_{tracks}");
 632   fhPhotonPtDiffVsNtracks->SetYTitle("<#rho^{#gamma}>*N_{cells}");
 633   outputContainer->Add(fhPhotonPtDiffVsNtracks);
 634
 635   fhPhotonPtDiffVsNclusters = new TH2F("PhotonPtDiffVsNclusters","PhotonPtDiffVsNclusters",50,0,100,50,0,5);
 636   fhPhotonPtDiffVsNclusters->SetXTitle("N_{clusters}");
 637   fhPhotonPtDiffVsNclusters->SetYTitle("<#rho^{#gamma}>*N_{cells}");
 638   outputContainer->Add(fhPhotonPtDiffVsNclusters);
 639
 640   fhPhotonPtDiffVsCentrality = new TH2F("PhotonPtDiffVsCentrality","PhotonPtDiffVsCentrality",100,0,100,50,0,5);
 641   fhPhotonPtDiffVsCentrality->SetXTitle("Centrality");
 642   fhPhotonPtDiffVsCentrality->SetYTitle("<#rho^{#gamma}>*N_{cells}");
 643   outputContainer->Add(fhPhotonPtDiffVsCentrality);
 644
 645   fhPhotonPtDiffVsNcells = new TH2F("PhotonPtDiffVsNcells","PhotonPtDiffVsNcells",100,0,200,50,0,5);
 646   fhPhotonPtDiffVsNcells->SetXTitle("N_{cells}");
 647   fhPhotonPtDiffVsNcells->SetYTitle("<#rho^{#gamma}>*N_{cells}");
 648   outputContainer->Add(fhPhotonPtDiffVsNcells);
 649
 650
 651   fhPhotonPtCorrectedZoom = new TH1F("PhotonPtCorrectedZoom","PhotonPtCorrectedZoom",200,-5,5);
 652   fhPhotonPtCorrectedZoom->SetXTitle("p_{T,#gamma} (GeV/c)");
 653   fhPhotonPtCorrectedZoom->SetYTitle("Counts");
 654   outputContainer->Add(fhPhotonPtCorrectedZoom);
 655
 656   fhPhotonSumPtInCone = new TH1F("PhotonSumPtInCone","PhotonSumPtInCone",100,0,100);
 657   fhPhotonSumPtInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
 658   fhPhotonSumPtInCone->SetYTitle("Counts");
 659   outputContainer->Add(fhPhotonSumPtInCone);
 660
 661   fhPhotonSumPtCorrectInCone = new TH1F("PhotonSumPtCorrectInCone","PhotonSumPtCorrectInCone",100,-20,80);
 662   fhPhotonSumPtCorrectInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
 663   fhPhotonSumPtCorrectInCone->SetYTitle("Counts");
 664   outputContainer->Add(fhPhotonSumPtCorrectInCone);
 665
 666   fhPhotonSumPtChargedInCone = new TH1F("PhotonSumPtChargedInCone","PhotonSumPtChargedInCone",100,0,100);
 667   fhPhotonSumPtChargedInCone->SetXTitle("#Sigma p_{T,#gamma}^{ch} (GeV/c)");
 668   fhPhotonSumPtChargedInCone->SetYTitle("Counts");
 669   outputContainer->Add(fhPhotonSumPtChargedInCone);
 670
 671
 672   //temporary jet histograms after selection
 673   fhSelectedJetPhiVsEta      = new TH2F("SelectedJetSelectedPhiVsEta","SelectedJetPhiVsEta",65,0,6.5,50,-1,1);
 674   fhSelectedJetPhiVsEta->SetYTitle("#eta_{jet}");
 675   fhSelectedJetPhiVsEta->SetXTitle("#phi_{jet}");
 676   outputContainer->Add(fhSelectedJetPhiVsEta) ;
 677
 678   fhSelectedJetChBkgEnergyVsPtJet = new TH2F("SelectedJetBkgChEnergyVsPtJet","SelectedJetBkgChEnergyVsPtJet",100,0,100,90,0,90);
 679   fhSelectedJetChBkgEnergyVsPtJet->SetYTitle("Jet Bkg Energy (GeV)");
 680   fhSelectedJetChBkgEnergyVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
 681   outputContainer->Add(fhSelectedJetChBkgEnergyVsPtJet);
 682
 683   fhSelectedJetChAreaVsPtJet      = new TH2F("SelectedJetChAreaVsPtJet","SelectedJetChAreaVsPtJet",100,0,100,50,0,1);
 684   fhSelectedJetChAreaVsPtJet->SetYTitle("Jet Area");
 685   fhSelectedJetChAreaVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
 686   outputContainer->Add(fhSelectedJetChAreaVsPtJet);
 687
 688   fhSelectedJetNjet      = new TH1F("SelectedJetNjet","SelectedJetNjet",100,0,100);
 689   fhSelectedJetNjet->SetYTitle("Counts");
 690   fhSelectedJetNjet->SetXTitle("N_{jets} per event");
 691   outputContainer->Add(fhSelectedJetNjet);
 692
 693   fhSelectedNtracks      = new TH1F("SelectedNtracks","SelectedNtracks",100,0,2000);
 694   fhSelectedNtracks->SetYTitle("Counts");
 695   fhSelectedNtracks->SetXTitle("N_{tracks} per event");
 696   outputContainer->Add(fhSelectedNtracks);
 697
 698   fhSelectedTrackPhiVsEta      = new TH2F("SelectedTrackPhiVsEta","SelectedTrackPhiVsEta",65,0,6.5,50,-1,1);
 699   fhSelectedTrackPhiVsEta->SetYTitle("#eta_{track}");
 700   fhSelectedTrackPhiVsEta->SetXTitle("#phi_{track}");
 701   outputContainer->Add(fhSelectedTrackPhiVsEta) ;
 702
 703   fhCuts2      = new TH1F("Cuts2","Cuts2",10,0.,10.);
 704   fhCuts2->SetYTitle("Counts");
 705   fhCuts2->SetXTitle("Cut number");
 706   outputContainer->Add(fhCuts2);
 707
 708   fhSelectedPhotonNLMVsPt      = new TH2F("SelectedPhotonNLMVsPt","SelectedPhotonNLMVsPt",100,0,100,10,0,10);
 709   fhSelectedPhotonNLMVsPt->SetYTitle("NLM");
 710   fhSelectedPhotonNLMVsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
 711   outputContainer->Add(fhSelectedPhotonNLMVsPt);
 712
 713   fhSelectedPhotonLambda0VsPt      = new TH2F("SelectedPhotonLambda0VsPt","SelectedPhotonLambda0VsPt",100,0,100,50,0,5);
 714   fhSelectedPhotonLambda0VsPt->SetYTitle("#lambda_{0}");
 715   fhSelectedPhotonLambda0VsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
 716   outputContainer->Add(fhSelectedPhotonLambda0VsPt);
 717
 718   //random
 719   fhRandomPhiEta[0]  = new TH2F("RandomPhiEtaRC",  "RandomPhiEtaRC",            100,0,6.5,100,-1.,1.);
 720   fhRandomPhiEta[1]  = new TH2F("RandomPhiEtaPCG", "RandomPhiEtaPerpConePhoton",100,0,6.5,100,-1.,1.);
 721   fhRandomPhiEta[2]  = new TH2F("RandomPhiEtaPCJ", "RandomPhiEtaPerpConeJet",   100,0,6.5,100,-1.,1.);
 722   fhRandomPhiEta[3]  = new TH2F("RandomPhiEtaMP",  "RandomPhiEtaMidPoint",      100,0,6.5,100,-1.,1.);
 723   fhRandomPhiEta[4]  = new TH2F("RandomPhiEtaTest","RandomPhiEtaTest",          100,0,6.5,100,-1.,1.);
 724   for(Int_t i=0;i<5;i++){
 725     fhRandomPhiEta[i]->SetXTitle("#phi");
 726     fhRandomPhiEta[i]->SetYTitle("#eta");
 727     outputContainer->Add(fhRandomPhiEta[i]);
 728   }
 729
 730   //MC generated photons and jets information
 731   if(fMCStudies) {
 732     fhMCPhotonCuts      = new TH1F("MCPhotonCuts","MCPhotonCuts",10,0.,10.);
 733     fhMCPhotonCuts->SetYTitle("Counts");
 734     fhMCPhotonCuts->SetXTitle("Cut number");
 735     outputContainer->Add(fhMCPhotonCuts);
 736
 737     fhMCPhotonPt      = new TH1F("MCPhotonPt","MCPhotonPt",100,0.,100.);
 738     fhMCPhotonPt->SetYTitle("Counts");
 739     fhMCPhotonPt->SetXTitle("p_{T,#gamma}^{gen} (GeV/c)");
 740     outputContainer->Add(fhMCPhotonPt);
 741
 742     fhMCPhotonEtaPhi      = new TH2F("MCPhotonEtaPhi","MCPhotonEtaPhi",65,0,6.5,50,-1,1);
 743     fhMCPhotonEtaPhi->SetYTitle("#eta_{#gamma}");
 744     fhMCPhotonEtaPhi->SetXTitle("#phi_{#gamma}");
 745     outputContainer->Add(fhMCPhotonEtaPhi) ;
 746
 747     fhMCJetOrigin      = new TH1F("MCJetOrigin","MCJetOrigin",35,-10.,25.);
 748     fhMCJetOrigin->SetYTitle("Counts");
 749     fhMCJetOrigin->SetXTitle("PDG number");
 750     outputContainer->Add(fhMCJetOrigin);
 751
 752     fhMCJetNPartVsPt      = new TH2F("MCJetNPartVsPt","MCJetNPartVsPt",100,0.,100.,100,0.,100.);
 753     fhMCJetNPartVsPt->SetYTitle("N_{particles}");
 754     fhMCJetNPartVsPt->SetXTitle("p_{T,full-jet}^{gen} (GeV/c)");
 755     outputContainer->Add(fhMCJetNPartVsPt) ;
 756
 757     fhMCJetChNPartVsPt      = new TH2F("MCJetChNPartVsPt","MCJetChNPartVsPt",100,0.,100.,100,0.,100.);
 758     fhMCJetChNPartVsPt->SetYTitle("N_{particles}");
 759     fhMCJetChNPartVsPt->SetXTitle("p_{T,charged-jet}^{gen} (GeV/c)");
 760     outputContainer->Add(fhMCJetChNPartVsPt) ;
 761
 762     fhMCJetNPart150VsPt      = new TH2F("MCJetNPart150VsPt","MCJetNPart150VsPt",100,0.,100.,100,0.,100.);
 763     fhMCJetNPart150VsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
 764     fhMCJetNPart150VsPt->SetXTitle("p_{T,full-jet}^{gen} (GeV/c)");
 765     outputContainer->Add(fhMCJetNPart150VsPt) ;
 766
 767     fhMCJetChNPart150VsPt      = new TH2F("MCJetChNPart150VsPt","MCJetChNPart150VsPt",100,0.,100.,100,0.,100.);
 768     fhMCJetChNPart150VsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
 769     fhMCJetChNPart150VsPt->SetXTitle("p_{T,charged-jet}^{gen} (GeV/c)");
 770     outputContainer->Add(fhMCJetChNPart150VsPt) ;
 771
 772     fhMCJetChNPart150ConeVsPt      = new TH2F("MCJetChNPart150ConeVsPt","MCJetChNPart150ConeVsPt",100,0.,100.,50,0.,50.);
 773     fhMCJetChNPart150ConeVsPt->SetYTitle("N_{particles (p_{T}>150 MeV/c)}");
 774     fhMCJetChNPart150ConeVsPt->SetXTitle("p_{T,charged-jet,R=0.4}^{gen} (GeV/c)");
 775     outputContainer->Add(fhMCJetChNPart150ConeVsPt) ;
 776
 777     fhMCJetRatioChFull      = new TH1F("MCJetRatioChFull","MCJetRatioChFull",100,0.,1.);
 778     fhMCJetRatioChFull->SetYTitle("Counts");
 779     fhMCJetRatioChFull->SetXTitle("p_{T,charged-jet}^{gen}/p_{T,full-jet}^{gen}");
 780     outputContainer->Add(fhMCJetRatioChFull);
 781
 782     fhMCJetRatioCh150Ch      = new TH1F("MCJetRatioCh150Ch","MCJetRatioCh150Ch",100,0.,1.);
 783     fhMCJetRatioCh150Ch->SetYTitle("Counts");
 784     fhMCJetRatioCh150Ch->SetXTitle("p_{T,charged-jet,pT>150MeV/c}^{gen}/p_{T,charged-jet}^{gen}");
 785     outputContainer->Add(fhMCJetRatioCh150Ch);
 786
 787     fhMCJetEtaPhi      = new TH2F("MCJetEtaPhi","MCJetEtaPhi",65,0,6.5,50,-1,1);
 788     fhMCJetEtaPhi->SetYTitle("#eta_{jet}");
 789     fhMCJetEtaPhi->SetXTitle("#phi_{jet}");
 790     outputContainer->Add(fhMCJetEtaPhi) ;
 791
 792     fhMCJetChEtaPhi      = new TH2F("MCJetChEtaPhi","MCJetChEtaPhi",65,0,6.5,50,-1,1);
 793     fhMCJetChEtaPhi->SetYTitle("#eta_{jet}");
 794     fhMCJetChEtaPhi->SetXTitle("#phi_{jet}");
 795     outputContainer->Add(fhMCJetChEtaPhi) ;
 796
 797     fhMCJet150EtaPhi      = new TH2F("MCJet150EtaPhi","MCJet150EtaPhi",65,0,6.5,50,-1,1);
 798     fhMCJet150EtaPhi->SetYTitle("#eta_{jet}");
 799     fhMCJet150EtaPhi->SetXTitle("#phi_{jet}");
 800     outputContainer->Add(fhMCJet150EtaPhi) ;
 801
 802     fhMCJetCh150EtaPhi      = new TH2F("MCJetCh150EtaPhi","MCJetCh150EtaPhi",65,0,6.5,50,-1,1);
 803     fhMCJetCh150EtaPhi->SetYTitle("#eta_{jet}");
 804     fhMCJetCh150EtaPhi->SetXTitle("#phi_{jet}");
 805     outputContainer->Add(fhMCJetCh150EtaPhi) ;
 806
 807     fhMCJetCh150ConeEtaPhi      = new TH2F("MCJetCh150ConeEtaPhi","MCJetCh150ConeEtaPhi",65,0,6.5,50,-1,1);
 808     fhMCJetCh150ConeEtaPhi->SetYTitle("#eta_{jet}");
 809     fhMCJetCh150ConeEtaPhi->SetXTitle("#phi_{jet}");
 810     outputContainer->Add(fhMCJetCh150ConeEtaPhi) ;
 811   }
 812
 813   //tree with data
 814   if(fSaveGJTree){
 815     fTreeGJ=new TTree("fTreeGJ","fTreeGJ");
 816     fTreeGJ->Branch("fGamPt"      ,&fGamPt    ,"fGamPt/D");//! fGamPt
 817     fTreeGJ->Branch("fGamLambda0" ,&fGamLambda0  ,"fGamLambda0/D");
 818     fTreeGJ->Branch("fGamNLM"     ,&fGamNLM      ,"fGamNLM/I");
 819     fTreeGJ->Branch("fGamSumPtCh" ,&fGamSumPtCh  ,"fGamSumPtCh/D");
 820     fTreeGJ->Branch("fGamNtracks" ,&fGamNtracks  ,"fGamNtracks/I");
 821     fTreeGJ->Branch("fGamTime"    ,&fGamTime     ,"fGamTime/D");
 822     fTreeGJ->Branch("fGamNcells"  ,&fGamNcells   ,"fGamNcells/I");
 823     fTreeGJ->Branch("fGamEta"     ,&fGamEta      ,"fGamEta/D");
 824     fTreeGJ->Branch("fGamPhi"     ,&fGamPhi      ,"fGamPhi/D");
 825     fTreeGJ->Branch("fGamSumPtNeu",&fGamSumPtNeu ,"fGamSumPtNeu/D");
 826     fTreeGJ->Branch("fGamNclusters" ,&fGamNclusters  ,"fGamNclusters/I");
 827     fTreeGJ->Branch("fGamAvEne"   ,&fGamAvEne    ,"fGamAvEne/D");
 828     fTreeGJ->Branch("fJetPhi"     ,&fJetPhi      ,"fJetPhi/D");
 829     fTreeGJ->Branch("fJetEta"     ,&fJetEta      ,"fJetEta/D");
 830     fTreeGJ->Branch("fJetPt"      ,&fJetPt       ,"fJetPt/D");
 831     fTreeGJ->Branch("fJetBkgChEne",&fJetBkgChEne ,"fJetBkgChEne/D");
 832     fTreeGJ->Branch("fJetArea"    ,&fJetArea     ,"fJetArea/D");
 833     fTreeGJ->Branch("fJetNtracks" ,&fJetNtracks  ,"fJetNtracks/I");
 834     fTreeGJ->Branch("fJetNtracks1" ,&fJetNtracks1  ,"fJetNtracks1/I");
 835     fTreeGJ->Branch("fJetNtracks2" ,&fJetNtracks2  ,"fJetNtracks2/I");
 836     fTreeGJ->Branch("fJetRho" ,&fJetRho  ,"fJetRho/D");
 837     fTreeGJ->Branch("fEventNumber",&fEventNumber ,"fEventNumber/I");
 838     fTreeGJ->Branch("fNtracks"    ,&fNtracks     ,"fNtracks/I");
 839     fTreeGJ->Branch("fZvertex"    ,&fZvertex     ,"fZvertex/D");
 840     fTreeGJ->Branch("fCentrality" ,&fCentrality  ,"fCentrality/D");
 841     fTreeGJ->Branch("fIso" ,&fIso  ,"fIso/O");
 842     fTreeGJ->Branch("fGamRho" ,&fGamRho  ,"fGamRho/D");
 843
 844     fTreeGJ->Branch("fMCGamPt"         ,&fMCGamPt         ,"fMCGamPt/D"        );
 845     fTreeGJ->Branch("fMCGamEta"        ,&fMCGamEta        ,"fMCGamEta/D"       );
 846     fTreeGJ->Branch("fMCGamPhi"        ,&fMCGamPhi        ,"fMCGamPhi/D"       );
 847     fTreeGJ->Branch("fMCPartonType"    ,&fMCPartonType    ,"fMCPartonType/I"   );
 848     fTreeGJ->Branch("fMCJetPt"         ,&fMCJetPt         ,"fMCJetPt/D"        );
 849     fTreeGJ->Branch("fMCJetChPt"       ,&fMCJetChPt       ,"fMCJetChPt/D"      );
 850     fTreeGJ->Branch("fMCJet150Pt"      ,&fMCJet150Pt      ,"fMCJet150Pt/D"     );
 851     fTreeGJ->Branch("fMCJetCh150Pt"    ,&fMCJetCh150Pt    ,"fMCJetCh150Pt/D"   );
 852     fTreeGJ->Branch("fMCJetNPart"      ,&fMCJetNPart      ,"fMCJetNPart/I"     );
 853     fTreeGJ->Branch("fMCJetChNPart"    ,&fMCJetChNPart    ,"fMCJetChNPart/I"   );
 854     fTreeGJ->Branch("fMCJet150NPart"   ,&fMCJet150NPart   ,"fMCJet150NPart/I"  );
 855     fTreeGJ->Branch("fMCJetCh150NPart" ,&fMCJetCh150NPart ,"fMCJetCh150NPart/I");
 856     fTreeGJ->Branch("fMCJetEta"        ,&fMCJetEta        ,"fMCJetEta/D"       );
 857     fTreeGJ->Branch("fMCJetPhi"        ,&fMCJetPhi        ,"fMCJetPhi/D"       );
 858     fTreeGJ->Branch("fMCJetChEta"      ,&fMCJetChEta      ,"fMCJetChEta/D"     );
 859     fTreeGJ->Branch("fMCJetChPhi"      ,&fMCJetChPhi      ,"fMCJetChPhi/D"     );
 860     fTreeGJ->Branch("fMCJet150Eta"     ,&fMCJet150Eta     ,"fMCJet150Eta/D"    );
 861     fTreeGJ->Branch("fMCJet150Phi"     ,&fMCJet150Phi     ,"fMCJet150Phi/D"    );
 862     fTreeGJ->Branch("fMCJetCh150Eta"   ,&fMCJetCh150Eta   ,"fMCJetCh150Eta/D"  );
 863     fTreeGJ->Branch("fMCJetCh150Phi"   ,&fMCJetCh150Phi   ,"fMCJetCh150Phi/D"  );
 864
 865     fTreeGJ->Branch("fMCJetCh150ConePt"    ,&fMCJetCh150ConePt    ,"fMCJetCh150ConePt/D"  );
 866     fTreeGJ->Branch("fMCJetCh150ConeNPart" ,&fMCJetCh150ConeNPart ,"fMCJetCh150ConeNPart/I");
 867     fTreeGJ->Branch("fMCJetCh150ConeEta"   ,&fMCJetCh150ConeEta   ,"fMCJetCh150ConeEta/D"  );
 868     fTreeGJ->Branch("fMCJetCh150ConePhi"   ,&fMCJetCh150ConePhi   ,"fMCJetCh150ConePhi/D"  );
 869
 870     outputContainer->Add(fTreeGJ);
 871   }
 872
 873   return outputContainer;
 874
 875 }
 876
 877 //_______________________________________________________
 878 void AliAnaParticleJetFinderCorrelation::InitParameters()
 879 {
 880   //printf("InitParameters\n");
 881   //Initialize the parameters of the analysis.
 882   SetInputAODName("PWG4Particle");
 883   AddToHistogramsName("AnaJetFinderCorr_");
 884
 885   fDeltaPhiMinCut    = 2.6 ;
 886   fDeltaPhiMaxCut    = 3.7 ;
 887   fRatioMaxCut       = 1.2 ;
 888   fRatioMinCut       = 0.3 ;
 889   fConeSize          = 0.4 ;
 890   fPtThresholdInCone = 0.5 ;
 891   fUseJetRefTracks   = kFALSE ;
 892   fMakeCorrelationInHistoMaker = kFALSE ;
 893   fSelectIsolated = kTRUE;
 894   fJetConeSize = 0.4 ;
 895   fJetMinPt = 15. ; //GeV/c
 896   fJetAreaFraction = 0.6 ;
 897   fJetBranchName = "jets";
 898   fBkgJetBranchName = "jets";
 899   fGammaConeSize = 0.4;
 900   fUseBackgroundSubtractionGamma = kFALSE;
 901   fSaveGJTree = kTRUE;
 902   fMostEnergetic = kFALSE;
 903   fMostOpposite = kTRUE;
 904   fUseHistogramJetBkg = kTRUE;
 905   fUseHistogramTracks = kTRUE;
 906   fUseHistogramJetTracks = kTRUE;
 907   fMCStudies = kFALSE;
 908 }
 909
 910 //__________________________________________________________________
 911 Int_t  AliAnaParticleJetFinderCorrelation::SelectJet(AliAODPWG4Particle * particle, TClonesArray *aodRecJets)
 912 {
 913   //Input for jets is TClonesArray *aodRecJets
 914   //Returns the index of the jet that is opposite to the particle
 915   //printf(" SelectJet ");
 916
 917   Double_t particlePt=particle->Pt();
 918   if(fUseBackgroundSubtractionGamma) {
 919     Int_t clusterIDtmp = particle->GetCaloLabel(0) ;
 920     Int_t nCells=0;
 921     AliVCluster *cluster=0;
 922     if(!(clusterIDtmp<0) ){
 923       Int_t iclustmp = -1;
 924       TObjArray* clusters = GetEMCALClusters();
 925       cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
 926       nCells = cluster->GetNCells();
 927     }
 928     particlePt-=(fGamRho*nCells);
 929   }
 930   if(particlePt<=0) {
 931     //printf("Particle with negative  or 0 pt\n");
 932     return -2;
 933   }
 934
 935   Int_t njets = aodRecJets->GetEntriesFast();
 936   AliAODJet * jet = 0 ;
 937   Int_t index = -1;
 938   Double_t dphiprev= 10000;
 939   Double_t particlePhi=particle->Phi();
 940   Double_t deltaPhi=-10000.;// in the range (0; 2*pi)
 941   Double_t jetPt=0.;
 942
 943   for(Int_t ijet = 0; ijet < njets ; ijet++){
 944     jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
 945     if(!jet)
 946     {
 947       AliInfo("Jet not in container");
 948       continue;
 949     }
 950     fhCuts2->Fill(2.,1.);
 951     jetPt=jet->Pt();
 952     if(jetPt<fJetMinPt) continue;
 953     fhCuts2->Fill(3.,1.);
 954     if(fBackgroundJetFromReader ){
 955       jetPt-= (fJetRho * jet->EffectiveAreaCharged() );
 956     }
 957     if(jetPt<0.) continue;
 958     //put jet eta requirement here |eta_jet|<0.9-jet_cone_size
 959     fhCuts2->Fill(4.,1.);
 960     if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
 961     fhCuts2->Fill(5.,1.);
 962     if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
 963     fhCuts2->Fill(6.,1.);
 964     //printf("jet found\n");
 965     Double_t deltaPhi0pi  = TMath::Abs(particle->Phi()-jet->Phi());
 966     Double_t ratio = jetPt/particlePt;
 967
 968     deltaPhi = particlePhi - jet->Phi() ;
 969     if ( deltaPhi0pi > TMath::Pi() ) deltaPhi0pi = 2. * TMath::Pi() - deltaPhi0pi ;
 970     if(deltaPhi<0) deltaPhi +=(TMath::Pi()*2.);
 971
 972     fhDeltaPtBefore ->Fill(particlePt, particlePt - jetPt);
 973     fhDeltaPhiBefore->Fill(particlePt, deltaPhi);
 974     fhDeltaEtaBefore->Fill(particlePt, particle->Eta() - jet->Eta());
 975     fhPtRatioBefore ->Fill(particlePt, jetPt/particlePt);
 976     fhPtBefore      ->Fill(particlePt, jetPt);
 977
 978     fhDeltaPhi0PiCorrectBefore->Fill(particlePt, deltaPhi0pi);//good
 979
 980     AliDebug(5,Form("Jet %d, Ratio pT %2.3f, Delta phi %2.3f",ijet,ratio,deltaPhi));
 981
 982     if((deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) &&
 983        (ratio > fRatioMinCut) && (ratio < fRatioMaxCut)  &&
 984        (TMath::Abs(deltaPhi-TMath::Pi()) < TMath::Abs(dphiprev-TMath::Pi())  )){//In case there is more than one jet select the most opposite.
 985       dphiprev = deltaPhi;
 986       index = ijet ;
 987     }//Selection cuts
 988   }//AOD Jet loop
 989
 990   return index ;
 991
 992 }
 993
 994 //__________________________________________________________________
 995 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD()
 996 {
 997   //Particle-Jet Correlation Analysis, fill AODs
 998   //  printf("I use MakeAnalysisFillAOD\n");
 999   //Get the event, check if there are AODs available, if not, skip this analysis
1000   AliAODEvent * event = NULL;
1001
1002   //  cout<<"GetReader()->GetOutputEvent() "<<GetReader()->GetOutputEvent()<<endl;
1003   //  cout<<"GetReader()->GetDataType() "<<GetReader()->GetDataType() <<endl;
1004   //  cout<<"AliCaloTrackReader::kAOD "<<AliCaloTrackReader::kAOD<<endl;
1005   //  cout<<"GetReader()->GetInputEvent() "<<GetReader()->GetInputEvent()<<endl;
1006
1007   if(GetReader()->GetOutputEvent())
1008   {
1009     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1010   }
1011   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1012   {
1013     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1014   }
1015   else
1016   {
1017     AliDebug(1,"There are no jets available for this analysis");
1018     return;
1019   }
1020
1021   if(!GetInputAODBranch() || !event)
1022   {
1023     AliFatal(Form("No input particles in AOD with name branch < %s > \n",
1024                   GetInputAODName().Data()));
1025     return; // Trick coverity
1026   }
1027
1028   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
1029   {
1030     AliFatal(Form("Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s>",
1031                   GetInputAODBranch()->GetClass()->GetName()));
1032     return; // Trick coverity
1033   }
1034
1035   //
1036   // non-standard jets
1037   //
1038   Int_t nJets=-1;
1039   TClonesArray *aodRecJets = 0;
1040   //if(IsNonStandardJetFromReader()){//jet branch from reader
1041   AliDebug(3,Form("GetNonStandardJets function (from reader) is called"));
1042   aodRecJets = GetNonStandardJets();
1043   AliDebug(3,Form("aodRecJets %p",aodRecJets));
1044   if(aodRecJets==0x0)
1045     {
1046       if(GetDebug() > 3) event->Print();
1047       AliFatal("List of jets is null");
1048       return;
1049     }
1050   nJets=aodRecJets->GetEntries();
1051   AliDebug(3,Form("nJets %d",nJets));
1052   //}
1053   //end of new part
1054
1055   if(nJets==0) {
1056     //printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1057     return;
1058   }
1059
1060   //
1061   //Background jets
1062   //
1063   AliAODJetEventBackground* aodBkgJets = 0;
1064   if(IsBackgroundJetFromReader()){//jet branch from reader
1065    AliDebug(3,"GetBackgroundJets function is called");
1066     aodBkgJets = GetBackgroundJets();
1067     AliDebug(3,Form("aodBkgJets %p",aodBkgJets));
1068     if(aodBkgJets==0x0){
1069       if(GetDebug() > 3) event->Print();
1070       AliFatal("No jet background found\n");
1071       return; // Trick coverity
1072     }
1073     if(GetDebug() > 3) aodBkgJets->Print("c");
1074   }
1075
1076   Double_t rhoEvent=0.;
1077   if(aodBkgJets && IsBackgroundJetFromReader() ) {
1078     rhoEvent = aodBkgJets->GetBackground(2);//hardcoded
1079   }
1080   fJetRho = rhoEvent;
1081
1082
1083   Int_t ntrig =  GetInputAODBranch()->GetEntriesFast() ;
1084
1085   AliDebug(3,"Begin jet finder  correlation analysis, fill AODs");
1086   AliDebug(3,Form("In particle branch aod entries %d\n", ntrig));
1087   AliDebug(3,Form("In standard jet branch aod entries %d\n", event->GetNJets()));
1088   AliDebug(3,Form("In non standard jet branch aod entries %d\n", nJets));
1089
1090   //if(nJets==0)   return;//to speed up
1091   //  cout<<"ntrig po return "<<ntrig<<endl;
1092
1093   //
1094   //calculate average cell energy without most energetic photon
1095   //
1096   Double_t medianPhotonRho=0.;
1097   TObjArray* clusters = GetEMCALClusters();
1098   Int_t clusterIDtmp;
1099   Int_t iclustmp = -1;
1100   AliVCluster *cluster=0;
1101
1102   if(IsBackgroundSubtractionGamma()){
1103     //
1104     // Find most energetic photon without background subtraction
1105     //
1106     Double_t maxPt=0.;
1107     Int_t maxIndex=-1;
1108     AliAODPWG4ParticleCorrelation* particlecorr =0;
1109     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1110       particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1111       if(particlecorr->Pt() > maxPt) {
1112         maxPt = particlecorr->Pt();
1113         maxIndex = iaod;
1114       }
1115     }
1116
1117     //
1118     // calculate background energy per cell
1119     //
1120     Int_t numberOfcells=0;
1121     if(ntrig>1){
1122       Double_t *photonRhoArr=new Double_t[ntrig-1];
1123       Int_t photonRhoArrayIndex=0;
1124       for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1125         particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1126         if(iaod==maxIndex) continue;
1127         clusterIDtmp = particlecorr->GetCaloLabel(0) ;
1128         if(clusterIDtmp < 0) continue;
1129         cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1130         photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1131         numberOfcells+=cluster->GetNCells();
1132         photonRhoArrayIndex++;
1133       }
1134       if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1135       delete [] photonRhoArr;
1136     }
1137   }//end of if background calculation for gamma
1138   fGamRho = medianPhotonRho;
1139
1140
1141   //
1142   //take most energetic photon and most energetic jet and combine
1143   //
1144   if(fMostEnergetic){
1145     //
1146     // most energetic photon with background subtraction
1147     //
1148     Double_t mostEnePhotonPt=0.;
1149     Int_t indexMostEnePhoton=-1;
1150     AliAODPWG4ParticleCorrelation* particle =0;
1151     Double_t ptCorrect=0.;
1152     Int_t nCells=0;
1153     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1154       particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1155       clusterIDtmp = particle->GetCaloLabel(0) ;
1156       if(!(clusterIDtmp<0)){
1157         cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1158         nCells = cluster->GetNCells();
1159       }
1160       ptCorrect = particle->Pt() - medianPhotonRho * nCells;
1161       if( ptCorrect > mostEnePhotonPt ){
1162         mostEnePhotonPt = ptCorrect;
1163         indexMostEnePhoton = iaod ;
1164       }
1165     }
1166     //    printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Photon with index %d selected \n",indexMostEnePhoton);
1167     //
1168     // most energetic jet with background subtraction
1169     //
1170     Double_t mostEneJetPt=0.;
1171     Int_t indexMostEneJet=-1;
1172     AliAODJet * jet = 0 ;
1173     //Double_t ptCorrect=0.;
1174     for(Int_t ijet = 0; ijet < nJets ; ijet++){
1175       jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1176       if(!jet) continue;
1177       if(jet->Pt()<fJetMinPt) continue;
1178       if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
1179       if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
1180       ptCorrect = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();
1181       if(ptCorrect > mostEneJetPt){
1182         mostEneJetPt = ptCorrect;
1183         indexMostEneJet = ijet;
1184       }
1185     }
1186     //    printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",indexMostEneJet);
1187
1188     //
1189     // assign jet to photon
1190     //
1191     if(indexMostEneJet>=0 && indexMostEnePhoton>=0){
1192       particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(indexMostEnePhoton));
1193       jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(indexMostEneJet));
1194       if(jet)particle->SetRefJet(jet);
1195     }
1196   }//end of take most energetic photon and most ene. jet after bkg subtraction
1197
1198   if(fMostOpposite){
1199     //Bool_t isJetFound=kFALSE;//new
1200     //Loop on stored AOD particles, trigger
1201     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1202       AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1203
1204       //Correlate with jets
1205       Int_t ijet = SelectJet(particle,aodRecJets);//input for jets is TClonesArray
1206       if(ijet > -1){
1207         //isJetFound=kTRUE;
1208         AliDebug(2,Form("Jet with index %d selected",ijet));
1209         AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1210         if(jet)particle->SetRefJet(jet);
1211         //printf("Most opposite found\n");
1212       }
1213     } // input aod loop
1214     //  if(GetReader()->WriteDeltaAODToFile() && isJetFound) WriteJetsToOutputBranch(aodRecJets);
1215   }//end of take most opposite photon and jet after bkg subtraction
1216
1217   AliDebug(1," End fill AODs \n");
1218 }
1219
1220 //__________________________________________________________________
1221 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms()
1222 {
1223   //Particle-Jet Correlation Analysis, fill histograms
1224
1225   AliDebug(3,"I use MakeAnalysisFillHistograms");
1226   AliDebug(3,Form("ntracks before iso %d\n",GetCTSTracks()->GetEntriesFast()));
1227
1228
1229   //Get the event, check if there are AODs available, if not, skip this analysis
1230   AliAODEvent * event = NULL;
1231
1232   //printf("GetReader()->GetOutputEvent() %d\n",GetReader()->GetOutputEvent() );
1233   //printf("GetReader()->GetDataType() %d\n",GetReader()->GetDataType() );
1234   //printf("AliCaloTrackReader::kAOD %d\n",AliCaloTrackReader::kAOD );
1235   //printf("GetReader()->GetInputEvent() %d\n",GetReader()->GetInputEvent() );
1236
1237   if(GetReader()->GetOutputEvent())
1238   {
1239     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1240   }
1241   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1242   {
1243     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1244   }
1245   else
1246   {
1247     AliDebug(3,"There are no jets available for this analysis");
1248     return;
1249   }
1250
1251   if(!GetInputAODBranch() || !event)
1252   {
1253     AliFatal(Form("No input particles in AOD with name branch < %s >",
1254                   GetInputAODName().Data()));
1255     return; // Trick coverity
1256   }
1257
1258   Int_t nJets=-1;
1259   TClonesArray *aodRecJets = 0;
1260   //if(IsNonStandardJetFromReader()){//branch read in reader from reader
1261   AliDebug(3,"GetNonStandardJets function (from reader) is called");
1262   aodRecJets = GetNonStandardJets();
1263   if(aodRecJets==0x0)
1264   {
1265     if(GetDebug() > 3) event->Print();
1266     AliFatal("Jets container not found\n");
1267     return; // trick coverity
1268   }
1269   nJets=aodRecJets->GetEntries();
1270   //}
1271   if(nJets==0)
1272   {
1273     //    printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1274     GetReader()->FillInputNonStandardJets();
1275     aodRecJets = GetNonStandardJets();
1276     if(aodRecJets) nJets=aodRecJets->GetEntries();
1277     //    printf("nJets = %d\n",nJets);
1278     return;
1279   }
1280
1281   //
1282   //Background jets
1283   //
1284   AliAODJetEventBackground* aodBkgJets = 0;
1285   if(IsBackgroundJetFromReader()){//jet branch from reader
1286     AliDebug(3,"GetBackgroundJets function is called");
1287     aodBkgJets = GetBackgroundJets();
1288     AliDebug(3,Form("aodBkgJets %p",aodBkgJets));
1289     if(aodBkgJets==0x0)
1290     {
1291       if(GetDebug() > 3) event->Print();
1292       AliFatal("No jet background container found");
1293       return; // trick coverity
1294     }
1295     if(GetDebug() > 3) aodBkgJets->Print("c");
1296   }
1297
1298   //
1299   // only background jets informations
1300   //
1301   //  Float_t pTback = 0;
1302   Double_t rhoEvent=0.;
1303   if(aodBkgJets) {
1304     if(IsBackgroundJetFromReader() ) rhoEvent = aodBkgJets->GetBackground(2);
1305     if(IsHistogramJetBkg()) {
1306       fhJetRhoVsCentrality->Fill(GetEventCentrality(),rhoEvent);
1307       for(Int_t i=0;i<4;i++){
1308         fhBkgJetBackground[i]->Fill(aodBkgJets->GetBackground(i));
1309         fhBkgJetSigma[i]->Fill(aodBkgJets->GetSigma(i));
1310         fhBkgJetArea[i]->Fill(aodBkgJets->GetMeanarea(i));
1311       }
1312     }//end of if fill HistogramJetBkg
1313   }//end if aodBkgJets exists
1314
1315   //
1316   //only track information
1317   //
1318   Int_t nCTSTracks = GetCTSTracks()->GetEntriesFast();
1319   AliAODTrack *aodtrack;
1320   Int_t itrack = 0;
1321   if(IsHistogramTracks()) {
1322     Double_t sumTrackPt=0;
1323     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1324       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1325       if(!aodtrack) continue;
1326       fhTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1327       sumTrackPt+=aodtrack->Pt();
1328     }
1329     if(nCTSTracks)
1330       fhTrackAveTrackPt->Fill(sumTrackPt/nCTSTracks);
1331   }//end if IsHistogramTracks
1332
1333   //
1334   //only jet informations
1335   //
1336   AliAODJet * jettmp = 0 ;
1337   Double_t jetPttmp = 0.;
1338   Double_t var1tmp = 0.;
1339   Double_t var2tmp = 0.;
1340   //  fhJetNjet->Fill(nJets);
1341   Double_t ptMostEne=0;
1342   //  Int_t indexMostEne=-1;
1343   Int_t nJetsOverThreshold[10]={nJets,0,0,0,0,0,0,0,0,0};
1344   Int_t iCounter=0;
1345   Double_t sumJetTrackPt=0.;
1346   Int_t sumNJetTrack=0;
1347   Int_t nTracksInJet=0;
1348   Int_t itrk=0;
1349   for(Int_t ijet = 0; ijet < nJets ; ijet++){
1350     jettmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1351     if(!jettmp) continue;
1352     fhJetPtBefore->Fill(jettmp->Pt());
1353     jetPttmp  = jettmp->Pt() - rhoEvent * jettmp->EffectiveAreaCharged();//<<---changed here
1354
1355     //calculate number of tracks above 1,2,3,4,5 GeV in jet
1356     AliVTrack* jettrack = 0x0 ;
1357
1358     Int_t nTrackThrGeV[5]={0,0,0,0,0};
1359     nTracksInJet=(jettmp->GetRefTracks())->GetEntriesFast();
1360     fhJetNparticlesInJet->Fill(nTracksInJet);
1361     if(nTracksInJet==0) continue;
1362     sumNJetTrack+=nTracksInJet;
1363     for(itrack=0;itrack<nTracksInJet;itrack++){
1364       jettrack=(AliVTrack *) ((jettmp->GetRefTracks())->At(itrack));
1365       if(!jettrack) continue;
1366
1367       fhJetDeltaEtaDeltaPhi->Fill(jettmp->Eta()-jettrack->Eta(),jettmp->Phi()-jettrack->Phi());
1368       sumJetTrackPt+=jettrack->Pt();
1369       if(IsHistogramJetTracks()){
1370         if(jettrack->Pt()>1.) nTrackThrGeV[0]++;
1371         if(jettrack->Pt()>2.) nTrackThrGeV[1]++;
1372         if(jettrack->Pt()>3.) nTrackThrGeV[2]++;
1373         if(jettrack->Pt()>4.) nTrackThrGeV[3]++;
1374         if(jettrack->Pt()>5.) nTrackThrGeV[4]++;
1375       }//end of if IsHistogramJetTracks
1376     }//end of jet track loop
1377     //printf("jetPt %f ntrks %d ntrks>1,2,3,4,5GeV in jet %d, %d, %d, %d, %d\n",jetPttmp,nTracksInJet,nTrackThrGeV[0],nTrackThrGeV[1],nTrackThrGeV[2],nTrackThrGeV[3],nTrackThrGeV[4]);
1378
1379     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1380       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1381       if(aodtrack) fhJetDeltaEtaDeltaPhiAllTracks->Fill(jettmp->Eta()-aodtrack->Eta(),jettmp->Phi()-aodtrack->Phi());
1382     }
1383
1384
1385     if(IsHistogramJetTracks()){
1386       fhJetNtracksInJetAboveThr[0]->Fill(jetPttmp,nTracksInJet);//all jets
1387
1388       for(itrk=0;itrk<5;itrk++) {
1389         fhJetNtracksInJetAboveThr[itrk+1]->Fill(jetPttmp,nTrackThrGeV[itrk]);//all jets
1390         fhJetRatioNTrkAboveToNTrk[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);//all jets
1391       }
1392       if(ijet==0){//most ene jet
1393         for(itrk=0;itrk<5;itrk++)
1394           fhJetNtrackRatioMostEne[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1395       }
1396       if(jetPttmp>5){//jet with pt>5GeV/c
1397         for(itrk=0;itrk<5;itrk++)
1398           fhJetNtrackRatioJet5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1399       }
1400       if(nTrackThrGeV[4]>0){//jet with leading particle pt>5GeV/c
1401         for(itrk=0;itrk<5;itrk++)
1402           fhJetNtrackRatioLead5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1403       }
1404     }//end of if IsHistogramJetTracks
1405
1406
1407     Double_t effectiveChargedBgEnergy=(IsBackgroundJetFromReader()?rhoEvent * jettmp->EffectiveAreaCharged():jettmp->ChargedBgEnergy());
1408
1409
1410     fhJetChBkgEnergyVsArea->Fill(effectiveChargedBgEnergy,jettmp->EffectiveAreaCharged());
1411     //if(jettmp->EffectiveAreaCharged()>0)
1412     fhJetRhoVsPt->Fill(jetPttmp,jettmp->ChargedBgEnergy()*jettmp->EffectiveAreaCharged());
1413
1414     if(jetPttmp>ptMostEne) {
1415       ptMostEne = jetPttmp;
1416       //indexMostEne=ijet;
1417     }
1418     fhJetPt->Fill(jetPttmp);
1419     fhJetChBkgEnergyVsPt->Fill(jetPttmp,effectiveChargedBgEnergy);
1420     fhJetChAreaVsPt->Fill(jetPttmp,jettmp->EffectiveAreaCharged());
1421     AliDebug(5,Form("ChargedBgEnergy %f EffectiveAreaCharged %f\n", jettmp->ChargedBgEnergy(),jettmp->EffectiveAreaCharged()));
1422     for(iCounter=1;iCounter<10;iCounter++)
1423     {
1424       if(jetPttmp>iCounter) nJetsOverThreshold[iCounter]++;
1425     }
1426
1427     var1tmp  = jettmp->Phi();
1428     var2tmp  = jettmp->Eta();
1429     fhJetPhi->Fill(var1tmp);
1430     fhJetEta->Fill(var2tmp);
1431     fhJetPhiVsEta->Fill(var1tmp,var2tmp);
1432     fhJetEtaVsPt->Fill(jetPttmp,var2tmp);
1433     fhJetEtaVsNpartInJet->Fill(var2tmp,nTracksInJet);
1434     if(jetPttmp>0)
1435       fhJetEtaVsNpartInJetBkg->Fill(var2tmp,nTracksInJet);
1436
1437   }//end of jet loop
1438   if(IsHistogramJetTracks()){
1439     if(sumNJetTrack>0){
1440       //printf("average track pt %f\n",sumJetTrackPt/sumNJetTrack);
1441       fhJetAveTrackPt->Fill(sumJetTrackPt/sumNJetTrack);
1442     }
1443   }//end of if IsHistogramJetTracks
1444
1445
1446   fhJetPtMostEne->Fill(ptMostEne);
1447   for(iCounter=0;iCounter<10;iCounter++){
1448     fhJetNjetOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1449     nJetsOverThreshold[iCounter]=0;
1450   }
1451
1452   //end of only jet part
1453
1454   //
1455   // Photons
1456   //
1457   Int_t ntrig   =  GetInputAODBranch()->GetEntriesFast() ;
1458
1459   AliDebug(1,"Begin jet finder  correlation analysis, fill histograms");
1460   AliDebug(1,Form("In particle branch aod entries %d\n", ntrig));
1461   AliDebug(1,Form("In jet output branch aod entries %d\n", event->GetNJets()));
1462
1463   fhNjetsNgammas->Fill(nJets,ntrig);
1464
1465   //if(nJets==0)   return;//to speed up
1466   //printf("ntrig %d, njets %d\n",ntrig,nJets);
1467
1468   //
1469   //Fill only temporary photon histograms
1470   //
1471   Double_t maxPt=0.;
1472   Int_t maxIndex=-1;
1473   Double_t tmpPt=0.;
1474   Double_t sumPt=0.;
1475   nJetsOverThreshold[0]=ntrig;
1476   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1477     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1478     tmpPt = particlecorr->Pt();
1479     if(tmpPt>maxPt) {
1480       maxPt = tmpPt;
1481       maxIndex = iaod;
1482     }
1483     for(iCounter=1;iCounter<10;iCounter++){
1484       if(tmpPt>iCounter) nJetsOverThreshold[iCounter]++;
1485     }
1486     sumPt+=tmpPt;
1487   }
1488   for(iCounter=0;iCounter<10;iCounter++){
1489     fhPhotonNgammaOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1490     //    nJetsOverThreshold[iCounter]=0;
1491   }
1492
1493   fGamAvEne=0;
1494   TObjArray* clusters = GetEMCALClusters();
1495   //printf("calculate median bkg energy for photons ");
1496   Double_t medianPhotonRho=0.;
1497   Int_t clusterID;
1498   Int_t iclustmp = -1;
1499   Int_t numberOfcells=0;
1500   AliVCluster *cluster = 0;
1501   if(ntrig>1){
1502     Double_t *photonRhoArr=new Double_t[ntrig-1];
1503     fhPhotonPtMostEne->Fill(maxPt);
1504     //    fhPhotonIndexMostEne->Fill(indexMaxPt);
1505     fhPhotonAverageEnergy->Fill(sumPt/ntrig);
1506     fhPhotonRatioAveEneToMostEne->Fill(sumPt/(ntrig*maxPt));
1507     fhPhotonAverageEnergyMinus1->Fill((sumPt-maxPt)/(ntrig-1));
1508     fGamAvEne=(sumPt-maxPt)/(ntrig-1);
1509     fhPhotonRatioAveEneMinus1ToMostEne->Fill((sumPt-maxPt)/((ntrig-1)*maxPt));
1510
1511     Int_t counterGamma=0;
1512     Int_t counterGammaMinus1=0;
1513
1514     Int_t photonRhoArrayIndex=0;
1515     //TObjArray* clusterstmp = GetEMCALClusters();
1516     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1517       AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1518       if( particlecorr->Pt() > sumPt/ntrig ) counterGamma++;
1519       if( particlecorr->Pt() > (sumPt-maxPt)/(ntrig-1) ) counterGammaMinus1++;
1520
1521       if(iaod==maxIndex) continue;
1522       clusterID = particlecorr->GetCaloLabel(0) ;
1523       if(clusterID < 0) continue;
1524       cluster = FindCluster(clusters,clusterID,iclustmp);
1525       photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1526       numberOfcells+=cluster->GetNCells();
1527       photonRhoArrayIndex++;
1528     }
1529     if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1530     delete [] photonRhoArr;
1531     fhPhotonNgammaMoreAverageToNgamma->Fill((Double_t)counterGamma / (Double_t)ntrig);
1532     fhPhotonNgammaMoreAverageMinus1ToNgamma->Fill((Double_t)counterGammaMinus1 / (Double_t)ntrig);
1533   }
1534   //printf("median = %f\n",medianPhotonRho);
1535   fhPhotonBkgRhoVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),medianPhotonRho);
1536   fhPhotonBkgRhoVsNclusters->Fill(ntrig,medianPhotonRho);
1537   fhPhotonBkgRhoVsCentrality->Fill(GetEventCentrality(),medianPhotonRho);
1538   fhPhotonBkgRhoVsNcells->Fill(numberOfcells,medianPhotonRho);
1539
1540
1541   AliVCluster *cluster2 = 0;
1542   Double_t photon2Corrected=0;
1543   Double_t sumPtTmp=0.;
1544   Double_t sumPtCorrectTmp=0.;
1545   AliVTrack* trackTmp = 0x0 ;
1546   TVector3 p3Tmp;
1547
1548   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1549     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1550     clusterID = particlecorr->GetCaloLabel(0) ;
1551     if(clusterID < 0) continue;
1552     cluster = FindCluster(clusters,clusterID,iclustmp);
1553     fhPhotonPt->Fill(particlecorr->Pt());
1554     fhPhotonPtCorrected->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1555     fhPhotonPtDiff->Fill(cluster->GetNCells() * medianPhotonRho);
1556     fhPhotonPtDiffVsCentrality->Fill(GetEventCentrality(),cluster->GetNCells() * medianPhotonRho);
1557     fhPhotonPtDiffVsNcells->Fill(numberOfcells,cluster->GetNCells() * medianPhotonRho);
1558     fhPhotonPtDiffVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),cluster->GetNCells() * medianPhotonRho);
1559     fhPhotonPtDiffVsNclusters->Fill(ntrig,cluster->GetNCells() * medianPhotonRho);
1560
1561     fhPhotonPtCorrectedZoom->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1562
1563     //test: sum_pt in the cone 0.3 for each photon
1564     //should be: random fake gamma from MB
1565     //is: each gamma for EMCEGA
1566     sumPtTmp=0.;
1567     sumPtCorrectTmp=0.;
1568
1569     for(Int_t iaod2 = 0; iaod2 < ntrig ; iaod2++){
1570       if(iaod==iaod2) continue;
1571       AliAODPWG4ParticleCorrelation* particlecorr2 =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod2));
1572       clusterID = particlecorr2->GetCaloLabel(0) ;
1573       if(clusterID < 0) continue;
1574       cluster2 = FindCluster(clusters,clusterID,iclustmp);
1575       photon2Corrected = particlecorr2->Pt() - cluster2->GetNCells() * medianPhotonRho;
1576
1577       //if(Pt()<0.5) continue; //<<hardcoded here //FIXME
1578       if( TMath::Sqrt((particlecorr->Eta()-particlecorr2->Eta())*(particlecorr->Eta()-particlecorr2->Eta()) +
1579                       (particlecorr->Phi()-particlecorr2->Phi())*(particlecorr->Phi()-particlecorr2->Phi()) )<fGammaConeSize ){//if(/*cone is correct*/){
1580         sumPtTmp+= particlecorr2->Pt();
1581         sumPtCorrectTmp+=photon2Corrected;
1582       }
1583     }
1584     fhPhotonSumPtInCone->Fill(sumPtTmp);
1585     fhPhotonSumPtCorrectInCone->Fill(sumPtCorrectTmp);
1586
1587     //test: sum_pt in the cone 0.3 for each track
1588     //should be: random fake gamma from MB
1589     //is: each gamma for EMCEGA
1590     sumPtTmp=0.;
1591     for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++){
1592       trackTmp = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1593       p3Tmp.SetXYZ(trackTmp->Px(),trackTmp->Py(),trackTmp->Pz());
1594       if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1595                       (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1596         sumPtTmp+=p3Tmp.Pt();
1597       }
1598     }//end of loop over tracks
1599     fhPhotonSumPtChargedInCone->Fill(sumPtTmp);
1600   }
1601
1602   //End of Fill temporary photon histograms
1603
1604   //
1605   // Apply background subtraction for photons
1606   //
1607   fGamRho = medianPhotonRho;
1608   if(!IsBackgroundSubtractionGamma()) medianPhotonRho=0;
1609
1610
1611   //
1612   //Get vertex for cluster momentum calculation <<----new here
1613   //
1614   Double_t vertex[] = {0,0,0} ; //vertex ;
1615   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1616     GetReader()->GetVertex(vertex);
1617   fZvertex = vertex[2];
1618
1619   //
1620   //Loop on stored AOD particles, trigger
1621   //
1622   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1623     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1624     fhCuts->Fill(0);
1625     fhCuts2->Fill(0.,(Double_t)nJets);
1626     AliDebug(1,Form("OnlyIsolated %d  !particlecorr->IsIsolated() %d \n",OnlyIsolated(), !particlecorr->IsIsolated()));
1627
1628     if(OnlyIsolated() && !particlecorr->IsIsolated()) continue;
1629     fhCuts->Fill(1);
1630     fhCuts2->Fill(1.,nJets);
1631
1632     if(nJets>0) {
1633       fhCuts->Fill(2);
1634     }
1635
1636     //Recover the jet correlated, found previously.
1637     AliAODJet   * jet = particlecorr->GetJet();
1638     //If correlation not made before, do it now.
1639     if(fMakeCorrelationInHistoMaker){
1640       //Correlate with jets
1641       Int_t ijet = SelectJet(particlecorr,aodRecJets);//input for jets is TClonesArray
1642       if(ijet > -1)
1643       {
1644         AliDebug(1,Form("Jet with index %d selected \n",ijet));
1645         //jet = event->GetJet(ijet);
1646         jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1647
1648        if(jet) particlecorr->SetRefJet(jet);
1649
1650       }
1651     }
1652
1653     if (!jet) continue ;
1654     fhCuts->Fill(3);
1655     fhCuts2->Fill(7.,1.);
1656
1657     //check MC genereted information
1658     if(fMCStudies) FindMCgenInfo();
1659
1660     //
1661     //Fill Correlation Histograms
1662     //
1663     clusterID = particlecorr->GetCaloLabel(0) ;
1664     if(!(clusterID<0)){
1665       cluster = FindCluster(clusters,clusterID,iclustmp);
1666       //fill tree variables
1667       fGamNcells = cluster->GetNCells();
1668     }
1669     Double_t ptTrig = particlecorr->Pt() - medianPhotonRho * fGamNcells;//<<---changed here
1670     Double_t ptJet = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();//<<---changed here
1671     Double_t phiTrig = particlecorr->Phi();
1672     Double_t phiJet = jet->Phi();
1673     Double_t etaTrig = particlecorr->Eta();
1674     Double_t etaJet = jet->Eta();
1675     Double_t deltaPhi=phiTrig-phiJet;
1676     if(deltaPhi<0)deltaPhi+=(TMath::Pi()*2.);
1677     //printf("pT trigger %2.3f, pT jet %2.3f, Delta phi %2.3f, Delta eta %2.3f, Delta pT %2.3f, ratio %2.3f \n",
1678     //  ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
1679     fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet);
1680     //    fhDeltaPhi->Fill(ptTrig, phiTrig-phiJet);//need to be shifted by 2pi
1681
1682     fhDeltaPhiCorrect->Fill(ptTrig, deltaPhi);//correct
1683
1684     Double_t deltaPhiCorrect = TMath::Abs( particlecorr->Phi() - jet->Phi() );
1685     if ( deltaPhiCorrect > TMath::Pi() ) deltaPhiCorrect = 2. * TMath::Pi() - deltaPhiCorrect ;
1686     fhDeltaPhi0PiCorrect->Fill(ptTrig, deltaPhiCorrect);
1687
1688     fhDeltaEta->Fill(ptTrig, etaTrig-etaJet);
1689     fhPtRatio ->Fill(ptTrig, ptJet/ptTrig);
1690     fhPt      ->Fill(ptTrig, ptJet);
1691
1692     fhSelectedJetPhiVsEta->Fill(phiJet,etaJet);
1693     fhSelectedJetChBkgEnergyVsPtJet->Fill(ptJet,(IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy()) );
1694     fhSelectedJetChAreaVsPtJet->Fill(ptJet,jet->EffectiveAreaCharged());
1695     fhSelectedJetNjet->Fill(nJets);
1696     fhSelectedNtracks->Fill(GetCTSTracks()->GetEntriesFast());//to be checked
1697     fhSelectedPhotonNLMVsPt->Fill(ptTrig,particlecorr->GetNLM());
1698
1699
1700     if(clusterID < 0 ){
1701       fhSelectedPhotonLambda0VsPt->Fill(ptTrig,-1);
1702       //fill tree variables
1703       fGamLambda0  = -1;
1704       fGamTime = -1;
1705       fGamNcells = 0;
1706       fGamSumPtNeu=0;
1707     }
1708     else{
1709       //Int_t iclus = -1;
1710       //      TObjArray* clusters = GetEMCALClusters();
1711       //cluster = FindCluster(clusters,clusterID,iclustmp);
1712       Double_t lambda0=cluster->GetM02();
1713       fhSelectedPhotonLambda0VsPt->Fill(ptTrig,lambda0);
1714       //fill tree variables
1715       fGamLambda0  = lambda0;
1716       fGamTime = cluster->GetTOF();
1717       //fGamNcells = cluster->GetNCells();
1718
1719       fGamSumPtNeu=0;
1720       fGamNclusters=0;
1721       //TVector3 p3Tmp;
1722       //Double_t scalarProduct=0;
1723       //Double_t vectorLength=particlecorr->P();
1724       for(Int_t icalo=0; icalo <clusters->GetEntriesFast(); icalo++){
1725         AliVCluster* calo = (AliVCluster *) clusters->At(icalo);
1726         if(clusterID==calo->GetID()) continue;//the same cluster as trigger
1727         calo->GetMomentum(fMomentum,vertex) ;//Assume that come from vertex in straight line
1728         //printf("min pt %f\n",GetMinPt());
1729         if(fMomentum.Pt()<GetMinPt()) continue; //<<hardcoded here //FIXME 0.5 check if correct
1730         p3Tmp.SetXYZ(fMomentum.Px(),fMomentum.Py(),fMomentum.Pz());
1731         //calculate sum pt in the cone
1732         if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1733                         (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1734           //scalarProduct = particlecorr->Px()*fMomentum.Px() + particlecorr->Py()*fMomentum.Py() + particlecorr->Pz()*fMomentum.Pz();
1735           //scalarProduct/=fMomentum.P();
1736           //scalarProduct/=vectorLength;
1737           //if(scalarProduct>TMath::Cos(0.3)) {//FIXME photon radius
1738           fGamSumPtNeu+=fMomentum.Pt();
1739           fGamNclusters++;
1740         }
1741       }
1742     }
1743
1744     //sum pt of charged tracks in the gamma isolation cone
1745     //starts here
1746     fGamSumPtCh=0;
1747     fGamNtracks=0;
1748     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1749       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1750       if(!aodtrack) continue;
1751       fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());//fill histogram here
1752       //      if(aodtrack->Pt()<0.15) continue;//hardcoded
1753       if(aodtrack->Pt()<fPtThresholdInCone) continue;
1754       if(!aodtrack->IsHybridGlobalConstrainedGlobal()) continue;
1755       if(TMath::Sqrt((particlecorr->Phi() - aodtrack->Phi())*(particlecorr->Phi() - aodtrack->Phi()) +
1756                      (particlecorr->Eta() - aodtrack->Eta())*(particlecorr->Eta() - aodtrack->Eta()) ) <fGammaConeSize ) {
1757         fGamSumPtCh+=aodtrack->Pt();
1758         fGamNtracks++;
1759       }
1760     }
1761     //ends here
1762
1763     //    for(Int_t itrack = 0; itrack < nCTSTracks ; itrack++){
1764     //      aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1765     //      fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1766     //    }
1767
1768     //
1769     // Background Fragmentation function
1770     //
1771     TVector3 gammaVector,jetVector;
1772     gammaVector.SetXYZ(particlecorr->Px(),particlecorr->Py(),particlecorr->Pz());
1773     jetVector.SetXYZ(jet->Px(),jet->Py(),jet->Pz());
1774     CalculateBkg(gammaVector,jetVector,vertex,1);//jet perp
1775     CalculateBkg(gammaVector,jetVector,vertex,2);//RC
1776     CalculateBkg(gammaVector,jetVector,vertex,3);//mid point
1777     CalculateBkg(gammaVector,jetVector,vertex,4);//gamma perp
1778     //CalculateBkg(gammaVector,jetVector,vertex,5);/test
1779     Double_t angleJetGam = gammaVector.Angle(jetVector);
1780     //printf("angleJetGam %f\n",angleJetGam*180/TMath::Pi());
1781
1782     //
1783     // Fragmentation function
1784     //
1785     Float_t      rad = 0, pt = 0, eta = 0, phi = 0;
1786     Int_t        npartcone = 0;
1787     TVector3 p3;
1788
1789     Int_t ntracks =  0;
1790
1791     AliDebug(1,Form("fUseJetRefTracks %d"   ,fUseJetRefTracks   ));
1792     AliDebug(1,Form("jet->GetRefTracks() %p",jet->GetRefTracks()));
1793     AliDebug(1,Form("GetCTSTracks() %p"     ,GetCTSTracks()     ));
1794
1795     if(!fUseJetRefTracks)
1796       ntracks =GetCTSTracks()->GetEntriesFast();
1797     else //If you want to use jet tracks from JETAN
1798       ntracks =  (jet->GetRefTracks())->GetEntriesFast();
1799
1800     AliDebug(3,Form("ntracks %d\n",ntracks));
1801     AliVTrack* track = 0x0 ;
1802     for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
1803       if(!fUseJetRefTracks)
1804         track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1805       else //If you want to use jet tracks from JETAN
1806         track = (AliVTrack *) ((jet->GetRefTracks())->At(ipr));
1807
1808       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1809       pt    = p3.Pt();
1810       eta  = p3.Eta();
1811       phi  = p3.Phi() ;
1812       if(phi < 0) phi+=TMath::TwoPi();
1813
1814       //Check if there is any particle inside cone with pt larger than  fPtThreshold
1815       rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
1816       if(rad < fConeSize  && pt > fPtThresholdInCone){
1817         //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
1818         npartcone++;
1819         fhFFz ->Fill(ptTrig, pt/ptTrig);
1820         fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt));
1821         fhFFpt->Fill(ptTrig, pt);
1822
1823         //according to jet axis
1824         fhJetFFz ->Fill(ptJet, pt/ptJet);
1825         fhJetFFxi->Fill(ptJet, TMath::Log(ptJet/pt));
1826         fhJetFFpt->Fill(ptJet, pt);
1827
1828
1829         if(TMath::Cos(angleJetGam)<0 && ptJet!=0 && pt!=0 ){
1830           fhJetFFzCor ->Fill(ptJet, -pt*TMath::Cos(angleJetGam)/ptJet);
1831           fhJetFFxiCor->Fill(ptJet, TMath::Log(ptJet/(-pt*TMath::Cos(angleJetGam))));
1832         }
1833       }
1834     }//Tracks
1835     fhNTracksInCone->Fill(ptTrig, npartcone);
1836     //fill tree here for each photon-jet (isolation required usually)
1837
1838     fGamPt      = ptTrig;
1839     //fGamLambda0  = ;//filled earlier
1840     fGamNLM      = particlecorr->GetNLM();
1841     //fGamSumPtCh  = ;//filled earlier
1842     //fGamTime     = particlecorr->GetTOF();//filled earlier
1843     //fGamNcells   = particlecorr->GetNCells();//filled earlier
1844     fGamEta      = etaTrig;
1845     fGamPhi      = phiTrig;
1846     //fGamSumPtNeu = ;//filled earlier
1847     //fGamNtracks  = ;//filled earlier
1848     //fGamNclusters= ;//filled earlier
1849     //fGamAvEne    = ;//filled earlier
1850     fJetPhi      = phiJet;
1851     fJetEta      = etaJet;
1852     fJetPt       = ptJet;
1853     fJetBkgChEne = (IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy());
1854     fJetArea     = jet->EffectiveAreaCharged();
1855     fJetNtracks  = (jet->GetRefTracks())->GetEntriesFast();
1856     fEventNumber = 0;
1857     fNtracks     = GetCTSTracks()->GetEntriesFast();
1858     fCentrality  = GetEventCentrality();
1859     fIso         = particlecorr->IsIsolated();
1860
1861     Int_t nTrk1GeV=0;
1862     Int_t nTrk2GeV=0;
1863     for(itrack=0;itrack < fJetNtracks;itrack++){
1864       track = (AliVTrack *) ((jet->GetRefTracks())->At(itrack));
1865       if(track->Pt()>1.) nTrk1GeV++;
1866       if(track->Pt()>2.) nTrk2GeV++;
1867     }
1868
1869     fJetNtracks1 = nTrk1GeV;
1870     fJetNtracks2 = nTrk2GeV;
1871
1872     if(fSaveGJTree) fTreeGJ->Fill();
1873   }//AOD trigger particle loop
1874   AliDebug(1,"End fill histograms");
1875 }
1876
1877
1878 //__________________________________________________________________
1879 void AliAnaParticleJetFinderCorrelation::Print(const Option_t * opt) const
1880 {
1881
1882   //Print some relevant parameters set for the analysis
1883   if(! opt)
1884     return;
1885
1886   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1887   AliAnaCaloTrackCorrBaseClass::Print(" ");
1888
1889   printf("Phi trigger-jet        <     %3.2f\n", fDeltaPhiMaxCut) ;
1890   printf("Phi trigger-jet        >     %3.2f\n", fDeltaPhiMinCut) ;
1891   printf("pT Ratio trigger/jet   <     %3.2f\n", fRatioMaxCut) ;
1892   printf("pT Ratio trigger/jet   >     %3.2f\n", fRatioMinCut) ;
1893   printf("fConeSize              =     %3.2f\n", fConeSize) ;
1894   printf("fPtThresholdInCone     =     %3.2f\n", fPtThresholdInCone) ;
1895   printf("fUseJetRefTracks         =     %d\n",    fUseJetRefTracks) ;
1896   printf("fMakeCorrelationInHistoMaker     =     %d\n",    fMakeCorrelationInHistoMaker) ;
1897   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
1898   printf("Reconstructed jet cone size = %3.2f\n", fJetConeSize) ;
1899   printf("Reconstructed jet minimum pt before background subtraction = %3.2f\n", fJetMinPt) ;
1900   printf("Reconstructed jet minimum area fraction = %3.2f\n", fJetAreaFraction) ;
1901
1902   //if(!fNonStandardJetFromReader){
1903   printf("fJetBranchName =   %s\n", fJetBranchName.Data()) ;
1904   //}
1905   if(!fBackgroundJetFromReader){
1906     printf("fBkgJetBranchName =   %s\n", fBkgJetBranchName.Data()) ;
1907   }
1908
1909   printf("Isolation cone size = %3.2f\n", fGammaConeSize) ;
1910   printf("fUseBackgroundSubtractionGamma = %d\n",fUseBackgroundSubtractionGamma);
1911   printf("fSaveGJTree = %d\n",fSaveGJTree);
1912   printf("fMostEnergetic = %d\n",fMostEnergetic);
1913   printf("fMostOpposite = %d\n",fMostOpposite);
1914
1915   printf("fUseHistogramJetBkg = %d\n",fUseHistogramJetBkg);
1916   printf("fUseHistogramTracks = %d\n",fUseHistogramTracks);
1917   printf("fUseHistogramJetTracks = %d\n",fUseHistogramJetTracks);
1918   printf("fMCStudies = %d\n",fMCStudies);
1919
1920 }
1921
1922 //__________________________________________________________________
1923 void AliAnaParticleJetFinderCorrelation::CalculateBkg(TVector3 gamma, TVector3 jet,Double_t vertex[3],Int_t type=2){
1924   //
1925   // calculate background for fragmentation function and fill histograms
1926   // 1. 90 degrees from jet axis in random place = perpendicular cone
1927   // 2. Random cone not belonging to jet cone nor photon cone
1928   // 3. In the middle point from jet and photon momentum vectors
1929   // 4. 90 degrees from photon direction in random place = perpendicular cone 2
1930
1931   //
1932   // implementation of 2 works, 1 and 4 works
1933   //
1934   Double_t gammaPt  = gamma.Pt();
1935   Double_t gammaEta = gamma.Eta();
1936   Double_t gammaPhi = gamma.Phi();
1937   Double_t jetEta   = jet.Eta();
1938   Double_t jetPhi   = jet.Phi();
1939
1940   //refference direction of background
1941   Double_t refEta=0.,refPhi=0.;
1942   Double_t rad = 0,rad2 = 0.;
1943   if(type==1){//perpendicular to jet axis
1944     //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
1945
1946     Double_t xVar;
1947     Double_t yVar;
1948     Double_t newX=0.;
1949     Double_t newY=0.;
1950     Double_t newZ=0.;
1951     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1952     Double_t jx=jet.Px();
1953     Double_t jy=jet.Py();
1954     Double_t jz=jet.Pz();
1955     //if(jz==0)printf("problem\n");
1956     //X axis
1957     Double_t Xx=1.0-vertex[0];
1958     Double_t Xy=-1.0*vertex[1];
1959     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1960     //Y axis
1961     Double_t Yx=jy*Xz-jz*Xy;
1962     Double_t Yy=jz*Xx-jx*Xz;
1963     Double_t Yz=jx*Xy-jy*Xx;
1964     //Determinant
1965     Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
1966     if(det==0)AliWarning("problem det==0\n");
1967     Double_t detX = 0.;
1968     Double_t detY = 0.;
1969     Double_t detZ = 0.;
1970
1971 //    Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
1972 //    printf("scalar jet.P o X: %f\n",tmpScalar);
1973 //    tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
1974 //    printf("scalar jet.P o Y: %f\n",tmpScalar);
1975 //    tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
1976 //    printf("scalar X o Y: %f\n",tmpScalar);
1977
1978     TVector3 perp;
1979     //randomise
1980     do{
1981       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1982       //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
1983       xVar=TMath::Cos(refPhi);
1984       yVar=TMath::Sin(refPhi);
1985       //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
1986       //zVar=0 in new surface frame
1987       detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
1988       detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
1989       detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
1990
1991       newX=detX/det;
1992       newY=detY/det;
1993       newZ=detZ/det;
1994
1995       perp.SetXYZ(newX,newY,newZ);
1996       refEta = perp.Eta();
1997       refPhi = perp.Phi();//output in <-pi, pi> range
1998       if(refPhi<0)refPhi+=2*TMath::Pi();
1999       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2000       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2001       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2002     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2003     fhRandomPhiEta[2]->Fill(refPhi,refEta);
2004
2005   }
2006   else if(type==2){//random cone
2007     //randomise
2008     do{
2009       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2010       refEta=fGenerator->Uniform(-(0.9-fJetConeSize),0.9-fJetConeSize);
2011       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2012       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2013       //check if reference is not within jet cone or gamma cone +0.4
2014       //example: FF fConeSize=0.4, fJetConeSize=0.4, fIsoGammaConeSize=0.3
2015     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize);
2016     //photon:0.7=0.4+0.3; jets:0.8=0.4 +0.4 //rad<fConeSize+fJetConeSize rad2<fConeSize+0.3
2017     fhRandomPhiEta[0]->Fill(refPhi,refEta);
2018   }
2019   else if(type==4){//perpendicular to photon axis
2020     Double_t xVar;
2021     Double_t yVar;
2022     Double_t newX=0.;
2023     Double_t newY=0.;
2024     Double_t newZ=0.;
2025     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
2026     Double_t jx=gamma.Px();
2027     Double_t jy=gamma.Py();
2028     Double_t jz=gamma.Pz();
2029     //if(jz==0)printf("problem\n");
2030     //X axis
2031     Double_t Xx=1.0-vertex[0];
2032     Double_t Xy=-1.0*vertex[1];
2033     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2034     //Y axis
2035     Double_t Yx=jy*Xz-jz*Xy;
2036     Double_t Yy=jz*Xx-jx*Xz;
2037     Double_t Yz=jx*Xy-jy*Xx;
2038     //Determinant
2039     Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
2040     if(det==0)AliWarning("problem det==0");
2041     Double_t detX = 0.;
2042     Double_t detY = 0.;
2043     Double_t detZ = 0.;
2044
2045 //    Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
2046 //    printf("scalar jet.P o X: %f\n",tmpScalar);
2047 //    tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
2048 //    printf("scalar jet.P o Y: %f\n",tmpScalar);
2049 //    tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
2050 //    printf("scalar X o Y: %f\n",tmpScalar);
2051
2052     TVector3 perp;
2053     //randomise
2054     do{
2055       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2056       //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
2057       xVar=TMath::Cos(refPhi);
2058       yVar=TMath::Sin(refPhi);
2059       //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
2060       //zVar=0 in new surface frame
2061       detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
2062       detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
2063       detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
2064
2065       newX=detX/det;
2066       newY=detY/det;
2067       newZ=detZ/det;
2068
2069       perp.SetXYZ(newX,newY,newZ);
2070       refEta = perp.Eta();
2071       refPhi = perp.Phi();//output in <-pi, pi> range
2072       if(refPhi<0)refPhi+=2*TMath::Pi();
2073       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2074       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2075       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2076     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2077     fhRandomPhiEta[1]->Fill(refPhi,refEta);
2078
2079   }
2080   else if(type==3){//mid point
2081
2082     Double_t jx=jet.Px();
2083     Double_t jy=jet.Py();
2084     Double_t jz=jet.Pz();
2085     //    if(jz==0)printf("problem\n");
2086     Double_t gx=gamma.Px();
2087     Double_t gy=gamma.Py();
2088     Double_t gz=gamma.Pz();
2089
2090     Double_t cosAlpha=(jx*gx+jy*gy+jz*gz)/(jet.Mag()*gamma.Mag());
2091     Double_t cosinus=TMath::Sqrt((cosAlpha+1.)/2.);
2092     //perpendicular axis
2093     Double_t Zx=gy*jz-gz*jy;
2094     Double_t Zy=gz*jx-gx*jz;
2095     Double_t Zz=gx*jy-gy*jx;
2096
2097     //Determinant
2098     Double_t det = Zx*gy*jz + Zy*gz*jx + Zz*gx*jy - Zz*gy*jx - Zy*gx*jz - Zx*gz*jy;
2099
2100     Double_t newX=0.;
2101     Double_t newY=0.;
2102     Double_t newZ=0.;
2103     if(det!=0) {
2104       Double_t detX =            -Zy*gz*cosinus +Zz*cosinus*jy + Zz*gy*cosinus - Zy*cosinus*jz;
2105       Double_t detY = Zx*cosinus*jz          - Zz*gx*cosinus - Zz*cosinus*jx            + Zx*gz*cosinus;
2106       Double_t detZ = -Zx*gy*cosinus + Zy*cosinus*jx + Zy*gx*cosinus - Zx*cosinus*jy;
2107
2108       newX=detX/det;
2109       newY=detY/det;
2110       newZ=detZ/det;
2111     }
2112
2113     TVector3 perp;
2114     perp.SetXYZ(newX,newY,newZ);
2115     refEta = perp.Eta();
2116     refPhi = perp.Phi();//output in <-pi, pi> range
2117     if(refPhi<0)refPhi+=2*TMath::Pi();
2118     rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2119     rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2120       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2121
2122     if (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize) fhRandomPhiEta[3]->Fill(refPhi,refEta);
2123   }
2124   else if(type==5){//tmp
2125     //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
2126
2127     Double_t xVar;
2128     Double_t newX=0.;
2129     Double_t newY=0.;
2130     Double_t newZ=0.;
2131     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
2132     Double_t jx=jet.Px();
2133     Double_t jy=jet.Py();
2134     Double_t jz=jet.Pz();
2135     //    if(jz==0)printf("problem\n");
2136     //X axis
2137     Double_t Xx=1.0-vertex[0];
2138     Double_t Xy=-1.0*vertex[1];
2139     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2140     //Y axis
2141     Double_t Yx=jy*Xz-jz*Xy;
2142     Double_t Yy=jz*Xx-jx*Xz;
2143     Double_t Yz=jx*Xy-jy*Xx;
2144
2145     // X and Y length
2146     Double_t Xlength=TMath::Sqrt(Xx*Xx+Xy*Xy+Xz*Xz);
2147     Double_t Ylength=TMath::Sqrt(Yx*Yx+Yy*Yy+Yz*Yz);
2148     Double_t ratio=Ylength/Xlength;
2149
2150     TVector3 perp;
2151     //randomise
2152     do{
2153       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2154       xVar=TMath::Tan(refPhi)/ratio;
2155       newX=xVar*Yx+Xx;
2156       newY=xVar*Yy+Xy;
2157       newZ=xVar*Yz+Xz;
2158
2159       perp.SetXYZ(newX,newY,newZ);
2160       refEta = perp.Eta();
2161       refPhi = perp.Phi();//output in <-pi, pi> range
2162       if(refPhi<0)refPhi+=2*TMath::Pi();
2163       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2164       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2165       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2166     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2167     fhRandomPhiEta[4]->Fill(refPhi,refEta);
2168   }
2169
2170
2171
2172   //calculate FF in background
2173   Int_t ntracks =  0;
2174   ntracks =GetCTSTracks()->GetEntriesFast();
2175   AliVTrack* track = 0x0 ;
2176   TVector3 p3;
2177
2178   Double_t pt = 0, eta = 0, phi = 0;
2179   Int_t  npartcone = 0;
2180   Double_t sumPt=0.;
2181   for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
2182     track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
2183     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
2184     pt   = p3.Pt();
2185     if(pt<fPtThresholdInCone) {//0.150
2186       //printf("problem: track pt < %f MeV/c \n",fPtThresholdInCone);
2187       continue;
2188     }
2189     eta  = p3.Eta() ;
2190     phi  = p3.Phi() ;
2191     if(phi < 0) phi+=TMath::TwoPi();
2192     //Check if there is any particle inside cone with pt larger than  fPtThreshold
2193     rad = TMath::Sqrt((eta-refEta)*(eta-refEta) + (phi-refPhi)*(phi-refPhi));
2194     if(rad < fConeSize  && pt > fPtThresholdInCone){
2195       //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
2196         npartcone++;
2197         sumPt+=pt;
2198         if(type==1){//perp jet
2199           fhBkgFFz[1] ->Fill(gammaPt, pt/gammaPt);
2200           fhBkgFFxi[1]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2201           fhBkgFFpt[1]->Fill(gammaPt, pt);
2202         }
2203         else if(type==2){//RC
2204           fhBkgFFz[0] ->Fill(gammaPt, pt/gammaPt);
2205           fhBkgFFxi[0]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2206           fhBkgFFpt[0]->Fill(gammaPt, pt);
2207         }
2208         else if(type==3){//mid point
2209           fhBkgFFz[3] ->Fill(gammaPt, pt/gammaPt);
2210           fhBkgFFxi[3]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2211           fhBkgFFpt[3]->Fill(gammaPt, pt);
2212         }
2213         else if(type==4){//perp jet
2214           fhBkgFFz[2] ->Fill(gammaPt, pt/gammaPt);
2215           fhBkgFFxi[2]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2216           fhBkgFFpt[2]->Fill(gammaPt, pt);
2217         }
2218         else if(type==5){//test
2219           fhBkgFFz[4] ->Fill(gammaPt, pt/gammaPt);
2220           fhBkgFFxi[4]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2221           fhBkgFFpt[4]->Fill(gammaPt, pt);
2222         }
2223
2224
2225     }
2226   }//end of loop over tracks
2227   Double_t sumOverTracks=0.;
2228   if(npartcone!=0) sumOverTracks = sumPt/npartcone;
2229   if(type==1) {
2230     fhBkgNTracksInCone[1]->Fill(gammaPt, npartcone);
2231     fhBkgSumPtInCone[1]->Fill(gammaPt,sumPt);
2232     fhBkgSumPtnTracksInCone[1]->Fill(gammaPt,sumOverTracks);
2233   }
2234   else if(type==2) {
2235     fhBkgNTracksInCone[0]->Fill(gammaPt, npartcone);
2236     fhBkgSumPtInCone[0]->Fill(gammaPt,sumPt);
2237     fhBkgSumPtnTracksInCone[0]->Fill(gammaPt,sumOverTracks);
2238   }
2239   else if(type==3) {
2240     fhBkgNTracksInCone[3]->Fill(gammaPt, npartcone);
2241     fhBkgSumPtInCone[3]->Fill(gammaPt,sumPt);
2242     fhBkgSumPtnTracksInCone[3]->Fill(gammaPt,sumOverTracks);
2243   }
2244   else if(type==4) {
2245     fhBkgNTracksInCone[2]->Fill(gammaPt, npartcone);
2246     fhBkgSumPtInCone[2]->Fill(gammaPt,sumPt);
2247     fhBkgSumPtnTracksInCone[2]->Fill(gammaPt,sumOverTracks);
2248   }
2249   else if(type==5) {
2250     fhBkgNTracksInCone[4]->Fill(gammaPt, npartcone);
2251     fhBkgSumPtInCone[4]->Fill(gammaPt,sumPt);
2252     fhBkgSumPtnTracksInCone[4]->Fill(gammaPt,sumOverTracks);
2253   }
2254 }
2255
2256
2257
2258
2259
2260 //__________________________________________________________________
2261 void AliAnaParticleJetFinderCorrelation::FindMCgenInfo(){
2262   //
2263   // Find information about photon and (quark or gluon) on generated level
2264   //
2265
2266   //frequently used variables
2267   Int_t pdg    = 0 ;
2268   Int_t mother = -1 ;
2269   Int_t absID  = 0 ;
2270
2271   //Double_t photonY   = -100 ;
2272   //Double_t photonE   = -1 ;
2273   Double_t photonPt  = -1 ;
2274   Double_t photonPhi =  100 ;
2275   Double_t photonEta = -1 ;
2276   Bool_t   inacceptance = kFALSE;
2277   AliAODMCParticle * primTmp = NULL;
2278
2279   //jet counters
2280   Int_t nParticlesInJet=0;
2281   Int_t nChargedParticlesInJet=0;
2282   Int_t nParticlesInJet150=0;
2283   Int_t nChargedParticlesInJet150=0;
2284   Int_t nChargedParticlesInJet150Cone=0;
2285
2286   Double_t eneParticlesInJet=0.;
2287   Double_t eneChargedParticlesInJet=0.;
2288   Double_t eneParticlesInJet150=0.;
2289   Double_t eneChargedParticlesInJet150=0.;
2290   Double_t eneChargedParticlesInJet150Cone=0.;
2291
2292   Double_t pxParticlesInJet=0.;
2293   Double_t pxChargedParticlesInJet=0.;
2294   Double_t pxParticlesInJet150=0.;
2295   Double_t pxChargedParticlesInJet150=0.;
2296   Double_t pxChargedParticlesInJet150Cone=0.;
2297
2298   Double_t pyParticlesInJet=0.;
2299   Double_t pyChargedParticlesInJet=0.;
2300   Double_t pyParticlesInJet150=0.;
2301   Double_t pyChargedParticlesInJet150=0.;
2302   Double_t pyChargedParticlesInJet150Cone=0.;
2303
2304   Double_t etaParticlesInJet=0.;
2305   Double_t etaChargedParticlesInJet=0.;
2306   Double_t etaParticlesInJet150=0.;
2307   Double_t etaChargedParticlesInJet150=0.;
2308   Double_t etaChargedParticlesInJet150Cone=0.;
2309
2310   Double_t phiParticlesInJet=0.;
2311   Double_t phiChargedParticlesInJet=0.;
2312   Double_t phiParticlesInJet150=0.;
2313   Double_t phiChargedParticlesInJet150=0.;
2314   Double_t phiChargedParticlesInJet150Cone=0.;
2315
2316   Double_t ptParticlesInJet=0.;
2317   Double_t ptChargedParticlesInJet=0.;
2318   Double_t ptParticlesInJet150=0.;
2319   Double_t ptChargedParticlesInJet150=0.;
2320   Double_t ptChargedParticlesInJet150Cone=0.;
2321
2322   Double_t coneJet=0.;
2323   Double_t coneChargedJet=0.;
2324   Double_t coneJet150=0.;
2325   Double_t coneChargedJet150=0.;
2326
2327   std::vector<Int_t> jetParticleIndex;
2328
2329   if(GetReader()->ReadStack()) {//ESD
2330      AliDebug(3,"I use stack");
2331   }//end Stack
2332   else if(GetReader()->ReadAODMCParticles()) {//AOD
2333     TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
2334     if(mcparticles){
2335       //jet origin
2336       //index =6 and 7 is hard scattering (jet-quark or photon)
2337       primTmp = (AliAODMCParticle *) mcparticles->At(6);
2338       pdg=primTmp->GetPdgCode();
2339        AliDebug(3,Form("id 6 pdg %d, pt %f ",pdg,primTmp->Pt() ));
2340       if(TMath::Abs(pdg)<=6 ||pdg==21) {
2341         fhMCJetOrigin->Fill(pdg);
2342         fMCPartonType=pdg;
2343       }
2344       primTmp = (AliAODMCParticle *) mcparticles->At(7);
2345       pdg=primTmp->GetPdgCode();
2346        AliDebug(3,Form("id 7 pdg %d, pt %f",pdg,primTmp->Pt() ));
2347       if(TMath::Abs(pdg)<=6 ||pdg==21) {
2348         fhMCJetOrigin->Fill(pdg);
2349         fMCPartonType=pdg;
2350       }
2351       //end of jet origin
2352
2353       Int_t nprim = mcparticles->GetEntriesFast();
2354       for(Int_t i=0; i < nprim; i++) {
2355         if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
2356         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
2357         pdg = prim->GetPdgCode();
2358         mother=prim->GetMother();
2359         //photon=22, gluon=21, quarks=(1,...,6), antiquarks=(-1,...,-6)
2360         if(pdg == 22){//photon
2361           fhMCPhotonCuts->Fill(0);
2362           if(prim->GetStatus()!=1) continue;
2363           fhMCPhotonCuts->Fill(1);
2364            AliDebug(5,Form("id %d, prim %d, physPrim %d, status %d\n",i,prim->IsPrimary(),prim->IsPhysicalPrimary(),prim->GetStatus()));
2365           while(mother>7){
2366             primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2367             mother=primTmp->GetMother();
2368           }
2369           if(mother<6)continue;
2370           fhMCPhotonCuts->Fill(2);
2371           primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2372           if(primTmp->GetPdgCode()!=22)continue;
2373           fhMCPhotonCuts->Fill(3);
2374
2375           //Get photon kinematics
2376           photonPt  = prim->Pt() ;
2377           photonPhi = prim->Phi() ;
2378           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
2379           photonEta = prim->Eta() ;
2380           fhMCPhotonPt->Fill(photonPt);
2381           fhMCPhotonEtaPhi->Fill(photonPhi,photonEta);
2382
2383           //Check if photons hit the Calorimeter
2384           fMomentum.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
2385           inacceptance = kFALSE;
2386           if(GetCaloUtils()->IsEMCALGeoMatrixSet()){
2387             fhMCPhotonCuts->Fill(4);
2388             //check if in EMCAL
2389             if(GetEMCALGeometry()) {
2390               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
2391               if(absID >= 0) inacceptance = kTRUE;
2392               AliDebug(3,Form("In EMCAL Real acceptance? %d",inacceptance));
2393             }
2394             else{
2395               if(GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) inacceptance = kTRUE ;
2396               AliDebug(1,Form("In EMCAL fiducial cut acceptance? %d",inacceptance));
2397             }
2398           }else{//no EMCAL nor EMCALGeoMatrixSet
2399             AliWarning("not EMCALGeoMatrix set");
2400           }//end of check if EMCAL
2401           if(inacceptance)fhMCPhotonCuts->Fill(5);
2402           AliDebug(5,Form("Photon Energy %f, Pt %f",prim->E(),prim->Pt()));
2403           fMCGamPt=photonPt;
2404           fMCGamEta=photonEta;
2405           fMCGamPhi=photonPhi;
2406         }//end of check if photon
2407         else
2408   {//not photon
2409           if(prim->GetStatus()!=1) continue;
2410           AliDebug(5,Form("id %d, prim %d, physPrim %d, status %d, pdg %d, E %f",
2411                     i,prim->IsPrimary(),prim->IsPhysicalPrimary(),prim->GetStatus(),prim->GetPdgCode(),prim->E()));
2412
2413           while(mother>7)
2414     {
2415             primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2416             mother=primTmp->GetMother();
2417             AliDebug(5,Form("next mother %d",mother));
2418           }
2419           if(mother<6)continue;//soft part
2420
2421           primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2422           pdg=primTmp->GetPdgCode();
2423           if( !(TMath::Abs(pdg)<=6 || pdg==21) ) continue;//origin not hard q or g
2424
2425           //jetParticleIndex.Add(&i);
2426           jetParticleIndex.push_back(i);
2427
2428         nParticlesInJet++;
2429         eneParticlesInJet+=prim->E();
2430         pxParticlesInJet+=prim->Px();
2431         pyParticlesInJet+=prim->Py();
2432         etaParticlesInJet+=(prim->E()*prim->Eta());
2433         photonPhi = prim->Phi() ;
2434         if(photonPhi < 0) photonPhi+=TMath::TwoPi();
2435         phiParticlesInJet+=(prim->E()*photonPhi);
2436
2437         if(prim->Charge()!=0) {
2438           nChargedParticlesInJet++;
2439           eneChargedParticlesInJet+=prim->E();
2440           pxChargedParticlesInJet+=prim->Px();
2441           pyChargedParticlesInJet+=prim->Py();
2442           etaChargedParticlesInJet+=(prim->E()*prim->Eta());
2443           phiChargedParticlesInJet+=(prim->E()*photonPhi);
2444         }
2445         if(prim->Pt()>0.150 && TMath::Abs(prim->Eta())<0.9 ) {//in TPC acceptance :)
2446           nParticlesInJet150++;
2447           eneParticlesInJet150+=prim->E();
2448           pxParticlesInJet150+=prim->Px();
2449           pyParticlesInJet150+=prim->Py();
2450           etaParticlesInJet150+=(prim->E()*prim->Eta());
2451           phiParticlesInJet150+=(prim->E()*photonPhi);
2452         }
2453         if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(prim->Eta())<0.9 ) {//in TPC acceptance :)
2454           nChargedParticlesInJet150++;
2455           eneChargedParticlesInJet150+=prim->E();
2456           pxChargedParticlesInJet150+=prim->Px();
2457           pyChargedParticlesInJet150+=prim->Py();
2458           etaChargedParticlesInJet150+=(prim->E()*prim->Eta());
2459           phiChargedParticlesInJet150+=(prim->E()*photonPhi);
2460         }
2461
2462         }//end of check pdg
2463       }//end of loop over primaries
2464
2465       if(eneParticlesInJet != 0.) {
2466         etaParticlesInJet/=eneParticlesInJet ;
2467         phiParticlesInJet/=eneParticlesInJet ;
2468       }
2469       if(eneChargedParticlesInJet != 0) {
2470         etaChargedParticlesInJet/=eneChargedParticlesInJet;
2471         phiChargedParticlesInJet/=eneChargedParticlesInJet;
2472       }
2473       if(eneParticlesInJet150 != 0) {
2474         etaParticlesInJet150/=eneParticlesInJet150;
2475         phiParticlesInJet150/=eneParticlesInJet150;
2476       }
2477       if(eneChargedParticlesInJet150 != 0) {
2478         etaChargedParticlesInJet150/=eneChargedParticlesInJet150;
2479         phiChargedParticlesInJet150/=eneChargedParticlesInJet150;
2480       }
2481
2482       ptParticlesInJet=TMath::Sqrt(pxParticlesInJet*pxParticlesInJet+pyParticlesInJet*pyParticlesInJet);
2483       ptChargedParticlesInJet=TMath::Sqrt(pxChargedParticlesInJet*pxChargedParticlesInJet+pyChargedParticlesInJet*pyChargedParticlesInJet);
2484       ptParticlesInJet150=TMath::Sqrt(pxParticlesInJet150*pxParticlesInJet150+pyParticlesInJet150*pyParticlesInJet150);
2485       ptChargedParticlesInJet150=TMath::Sqrt(pxChargedParticlesInJet150*pxChargedParticlesInJet150+pyChargedParticlesInJet150*pyChargedParticlesInJet150);
2486
2487       Double_t distance=0.;
2488       Double_t eta=0.;
2489       Double_t phi=0.;
2490       Double_t mostPtCharged=0.;
2491       Int_t mostmostPtChargedId=-1;
2492       std::vector<Int_t>::iterator it;
2493       for( it=jetParticleIndex.begin(); it!=jetParticleIndex.end(); ++it ){
2494         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(*it);
2495         eta = prim->Eta();
2496         phi = prim->Phi();
2497         if(phi < 0) phi+=TMath::TwoPi();
2498         //full jet
2499         distance=TMath::Sqrt((eta-etaParticlesInJet)*(eta-etaParticlesInJet)+(phi-phiParticlesInJet)*(phi-phiParticlesInJet));
2500         if(distance>coneJet) coneJet=distance;
2501         //charged jet
2502         distance=TMath::Sqrt((eta-etaChargedParticlesInJet)*(eta-etaChargedParticlesInJet)+(phi-phiChargedParticlesInJet)*(phi-phiChargedParticlesInJet));
2503         if(distance>coneChargedJet) coneChargedJet=distance;
2504         //
2505         distance=TMath::Sqrt((eta-etaParticlesInJet150)*(eta-etaParticlesInJet150)+(phi-phiParticlesInJet150)*(phi-phiParticlesInJet150));
2506         if(distance>coneJet150 && TMath::Abs(eta)<0.9 ) coneJet150=distance;
2507         //
2508         distance=TMath::Sqrt((eta-etaChargedParticlesInJet150)*(eta-etaChargedParticlesInJet150)+(phi-phiChargedParticlesInJet150)*(phi-phiChargedParticlesInJet150));
2509         if(distance>coneChargedJet150 && TMath::Abs(eta)<0.9) coneChargedJet150=distance;
2510
2511         if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(eta)<0.9) {
2512           if(prim->Pt()>mostPtCharged) {
2513             mostPtCharged=prim->Pt();
2514             mostmostPtChargedId=(*it);
2515           }
2516         }
2517
2518         if(distance<=0.4){
2519           if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(eta)<0.9) {
2520             nChargedParticlesInJet150Cone++;
2521             eneChargedParticlesInJet150Cone+=prim->E();
2522             pxChargedParticlesInJet150Cone+=prim->Px();
2523             pyChargedParticlesInJet150Cone+=prim->Py();
2524             etaChargedParticlesInJet150Cone+=(prim->E()*eta);
2525             phiChargedParticlesInJet150Cone+=(prim->E()*phi);
2526           }
2527
2528         }
2529
2530       }//end of loop over jet particle indexes
2531       if(eneChargedParticlesInJet150Cone != 0) {
2532         etaChargedParticlesInJet150Cone/=eneChargedParticlesInJet150Cone;
2533         phiChargedParticlesInJet150Cone/=eneChargedParticlesInJet150Cone;
2534       }
2535       ptChargedParticlesInJet150Cone=TMath::Sqrt(pxChargedParticlesInJet150Cone*pxChargedParticlesInJet150Cone+pyChargedParticlesInJet150Cone*pyChargedParticlesInJet150Cone);
2536       if(nChargedParticlesInJet150>0 && nChargedParticlesInJet150Cone<1){//no particles in cone? take the most energetic one
2537         nChargedParticlesInJet150Cone=1;
2538         etaChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Eta();
2539         phiChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Phi();
2540         ptChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Pt();
2541       }
2542
2543
2544     }//mc array exists and data is MC
2545   }// read AOD MC
2546
2547   jetParticleIndex.clear();
2548
2549
2550   //printouts
2551
2552   AliDebug(3,Form("cone full %f, charged %f, full150 %f, charged150 %f",coneJet,coneChargedJet,coneJet150,coneChargedJet150));
2553   AliDebug(3,Form("Npart %d, NchPart %d, Npart(pt>150M) %d, NchPart(pt>150M) %d, NchPart(pt>150M)Cone %d\n",nParticlesInJet,nChargedParticlesInJet,nParticlesInJet150,nChargedParticlesInJet150,nChargedParticlesInJet150Cone));
2554   AliDebug(3,Form("Etot %f, Ech %f, E(pt>150M) %f, Ech(pt>150M) %f\n",eneParticlesInJet,eneChargedParticlesInJet,eneParticlesInJet150,eneChargedParticlesInJet150));
2555   AliDebug(3,Form("pt %f, ptch %f, pt(pt>150M) %f,ptch(pt>150M) %f,ptch(pt>150M)Cone %f\n",ptParticlesInJet,ptChargedParticlesInJet,ptParticlesInJet150,ptChargedParticlesInJet150,ptChargedParticlesInJet150Cone));
2556   AliDebug(3,Form("eta/phi tot %f/%f, ch %f/%f, tot150 %f/%f,  ch150 %f/%f, ch150cone %f/%f\n",etaParticlesInJet,phiParticlesInJet,etaChargedParticlesInJet,phiChargedParticlesInJet,etaParticlesInJet150,phiParticlesInJet150,etaChargedParticlesInJet150,phiChargedParticlesInJet150,etaChargedParticlesInJet150Cone,phiChargedParticlesInJet150Cone));
2557
2558   //fill histograms
2559   if(ptParticlesInJet) fhMCJetRatioChFull->Fill(ptChargedParticlesInJet/ptParticlesInJet);
2560   if(ptChargedParticlesInJet) fhMCJetRatioCh150Ch->Fill(ptChargedParticlesInJet150/ptChargedParticlesInJet);
2561
2562   fhMCJetNPartVsPt     ->Fill(ptParticlesInJet,nParticlesInJet);
2563   fhMCJetChNPartVsPt   ->Fill(ptChargedParticlesInJet,nChargedParticlesInJet);
2564   fhMCJetNPart150VsPt  ->Fill(ptParticlesInJet150,nParticlesInJet150);
2565   fhMCJetChNPart150VsPt->Fill(ptChargedParticlesInJet150,nChargedParticlesInJet150);
2566   fhMCJetChNPart150ConeVsPt->Fill(ptChargedParticlesInJet150Cone,nChargedParticlesInJet150Cone);
2567
2568   fhMCJetEtaPhi->Fill(phiParticlesInJet,etaParticlesInJet);
2569   fhMCJetChEtaPhi->Fill(phiChargedParticlesInJet,etaChargedParticlesInJet);
2570   fhMCJet150EtaPhi->Fill(phiParticlesInJet150,etaParticlesInJet150);
2571   fhMCJetCh150EtaPhi->Fill(phiChargedParticlesInJet150,etaChargedParticlesInJet150);
2572   fhMCJetCh150ConeEtaPhi->Fill(phiChargedParticlesInJet150Cone,etaChargedParticlesInJet150Cone);
2573
2574   //fill tree
2575   fMCJetPt      = ptParticlesInJet;
2576   fMCJetChPt    = ptChargedParticlesInJet;
2577   fMCJet150Pt   = ptParticlesInJet150;
2578   fMCJetCh150Pt = ptChargedParticlesInJet150;
2579   fMCJetNPart      = nParticlesInJet;
2580   fMCJetChNPart    = nChargedParticlesInJet;
2581   fMCJet150NPart   = nParticlesInJet150;
2582   fMCJetCh150NPart = nChargedParticlesInJet150;
2583   fMCJetEta      = etaParticlesInJet          ;
2584   fMCJetPhi      = phiParticlesInJet          ;
2585   fMCJetChEta    = etaChargedParticlesInJet   ;
2586   fMCJetChPhi    = phiChargedParticlesInJet   ;
2587   fMCJet150Eta   = etaParticlesInJet150       ;
2588   fMCJet150Phi   = phiParticlesInJet150       ;
2589   fMCJetCh150Eta = etaChargedParticlesInJet150;
2590   fMCJetCh150Phi = phiChargedParticlesInJet150;
2591
2592   fMCJetCh150ConePt    = ptChargedParticlesInJet150Cone;
2593   fMCJetCh150ConeNPart = nChargedParticlesInJet150Cone;
2594   fMCJetCh150ConeEta   = etaChargedParticlesInJet150Cone;
2595   fMCJetCh150ConePhi   = phiChargedParticlesInJet150Cone;
2596
2597 }//end of method FindMCgenInfo