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     if(GetDebug() > 5)
 981       printf("AliAnaParticleJetFinderCorrelation::SelectJet() - Jet %d, Ratio pT %2.3f, Delta phi %2.3f\n",ijet,ratio,deltaPhi);
 982
 983     if((deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) &&
 984        (ratio > fRatioMinCut) && (ratio < fRatioMaxCut)  &&
 985        (TMath::Abs(deltaPhi-TMath::Pi()) < TMath::Abs(dphiprev-TMath::Pi())  )){//In case there is more than one jet select the most opposite.
 986       dphiprev = deltaPhi;
 987       index = ijet ;
 988     }//Selection cuts
 989   }//AOD Jet loop
 990
 991   return index ;
 992
 993 }
 994
 995 //__________________________________________________________________
 996 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD()
 997 {
 998   //Particle-Jet Correlation Analysis, fill AODs
 999   //  printf("I use MakeAnalysisFillAOD\n");
1000   //Get the event, check if there are AODs available, if not, skip this analysis
1001   AliAODEvent * event = NULL;
1002
1003   //  cout<<"GetReader()->GetOutputEvent() "<<GetReader()->GetOutputEvent()<<endl;
1004   //  cout<<"GetReader()->GetDataType() "<<GetReader()->GetDataType() <<endl;
1005   //  cout<<"AliCaloTrackReader::kAOD "<<AliCaloTrackReader::kAOD<<endl;
1006   //  cout<<"GetReader()->GetInputEvent() "<<GetReader()->GetInputEvent()<<endl;
1007
1008   if(GetReader()->GetOutputEvent())
1009   {
1010     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1011   }
1012   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1013   {
1014     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1015   }
1016   else
1017   {
1018     if(GetDebug() > 3) AliInfo("There are no jets available for this analysis");
1019     return;
1020   }
1021
1022   if(!GetInputAODBranch() || !event)
1023   {
1024     AliFatal(Form("No input particles in AOD with name branch < %s > \n",
1025                   GetInputAODName().Data()));
1026     return; // Trick coverity
1027   }
1028
1029   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
1030   {
1031     AliFatal(Form("Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s>",
1032                   GetInputAODBranch()->GetClass()->GetName()));
1033     return; // Trick coverity
1034   }
1035
1036   //
1037   // non-standard jets
1038   //
1039   Int_t nJets=-1;
1040   TClonesArray *aodRecJets = 0;
1041   //if(IsNonStandardJetFromReader()){//jet branch from reader
1042   if(GetDebug() > 3) AliInfo(Form("GetNonStandardJets function (from reader) is called"));
1043   aodRecJets = GetNonStandardJets();
1044   if(GetDebug() > 3) AliInfo(Form("aodRecJets %p",aodRecJets));
1045   if(aodRecJets==0x0)
1046     {
1047       if(GetDebug() > 3) event->Print();
1048       AliFatal("List of jets is null");
1049       return;
1050     }
1051   nJets=aodRecJets->GetEntries();
1052   if(GetDebug() > 3) printf("nJets %d\n",nJets);
1053   //}
1054   //end of new part
1055
1056   if(nJets==0) {
1057     //printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1058     return;
1059   }
1060
1061   //
1062   //Background jets
1063   //
1064   AliAODJetEventBackground* aodBkgJets = 0;
1065   if(IsBackgroundJetFromReader()){//jet branch from reader
1066     if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
1067     aodBkgJets = GetBackgroundJets();
1068     if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
1069     if(aodBkgJets==0x0){
1070       if(GetDebug() > 3) event->Print();
1071       AliFatal("No jet background found\n");
1072       return; // Trick coverity
1073     }
1074     if(GetDebug() > 3) aodBkgJets->Print("c");
1075   }
1076
1077   Double_t rhoEvent=0.;
1078   if(aodBkgJets && IsBackgroundJetFromReader() ) {
1079     rhoEvent = aodBkgJets->GetBackground(2);//hardcoded
1080   }
1081   fJetRho = rhoEvent;
1082
1083
1084   Int_t ntrig =  GetInputAODBranch()->GetEntriesFast() ;
1085   if(GetDebug() > 3){
1086     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Begin jet finder  correlation analysis, fill AODs \n");
1087     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", ntrig);
1088     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In standard jet branch aod entries %d\n", event->GetNJets());
1089     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In non standard jet branch aod entries %d\n", nJets);
1090   }
1091
1092   //if(nJets==0)   return;//to speed up
1093   //  cout<<"ntrig po return "<<ntrig<<endl;
1094
1095   //
1096   //calculate average cell energy without most energetic photon
1097   //
1098   Double_t medianPhotonRho=0.;
1099   TObjArray* clusters = GetEMCALClusters();
1100   Int_t clusterIDtmp;
1101   Int_t iclustmp = -1;
1102   AliVCluster *cluster=0;
1103
1104   if(IsBackgroundSubtractionGamma()){
1105     //
1106     // Find most energetic photon without background subtraction
1107     //
1108     Double_t maxPt=0.;
1109     Int_t maxIndex=-1;
1110     AliAODPWG4ParticleCorrelation* particlecorr =0;
1111     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1112       particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1113       if(particlecorr->Pt() > maxPt) {
1114         maxPt = particlecorr->Pt();
1115         maxIndex = iaod;
1116       }
1117     }
1118
1119     //
1120     // calculate background energy per cell
1121     //
1122     Int_t numberOfcells=0;
1123     if(ntrig>1){
1124       Double_t *photonRhoArr=new Double_t[ntrig-1];
1125       Int_t photonRhoArrayIndex=0;
1126       for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1127         particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1128         if(iaod==maxIndex) continue;
1129         clusterIDtmp = particlecorr->GetCaloLabel(0) ;
1130         if(clusterIDtmp < 0) continue;
1131         cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1132         photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1133         numberOfcells+=cluster->GetNCells();
1134         photonRhoArrayIndex++;
1135       }
1136       if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1137       delete [] photonRhoArr;
1138     }
1139   }//end of if background calculation for gamma
1140   fGamRho = medianPhotonRho;
1141
1142
1143   //
1144   //take most energetic photon and most energetic jet and combine
1145   //
1146   if(fMostEnergetic){
1147     //
1148     // most energetic photon with background subtraction
1149     //
1150     Double_t mostEnePhotonPt=0.;
1151     Int_t indexMostEnePhoton=-1;
1152     AliAODPWG4ParticleCorrelation* particle =0;
1153     Double_t ptCorrect=0.;
1154     Int_t nCells=0;
1155     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1156       particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1157       clusterIDtmp = particle->GetCaloLabel(0) ;
1158       if(!(clusterIDtmp<0)){
1159         cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1160         nCells = cluster->GetNCells();
1161       }
1162       ptCorrect = particle->Pt() - medianPhotonRho * nCells;
1163       if( ptCorrect > mostEnePhotonPt ){
1164         mostEnePhotonPt = ptCorrect;
1165         indexMostEnePhoton = iaod ;
1166       }
1167     }
1168     //    printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Photon with index %d selected \n",indexMostEnePhoton);
1169     //
1170     // most energetic jet with background subtraction
1171     //
1172     Double_t mostEneJetPt=0.;
1173     Int_t indexMostEneJet=-1;
1174     AliAODJet * jet = 0 ;
1175     //Double_t ptCorrect=0.;
1176     for(Int_t ijet = 0; ijet < nJets ; ijet++){
1177       jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1178       if(!jet) continue;
1179       if(jet->Pt()<fJetMinPt) continue;
1180       if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
1181       if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
1182       ptCorrect = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();
1183       if(ptCorrect > mostEneJetPt){
1184         mostEneJetPt = ptCorrect;
1185         indexMostEneJet = ijet;
1186       }
1187     }
1188     //    printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",indexMostEneJet);
1189
1190     //
1191     // assign jet to photon
1192     //
1193     if(indexMostEneJet>=0 && indexMostEnePhoton>=0){
1194       particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(indexMostEnePhoton));
1195       jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(indexMostEneJet));
1196       if(jet)particle->SetRefJet(jet);
1197     }
1198   }//end of take most energetic photon and most ene. jet after bkg subtraction
1199
1200   if(fMostOpposite){
1201     //Bool_t isJetFound=kFALSE;//new
1202     //Loop on stored AOD particles, trigger
1203     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1204       AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1205
1206       //Correlate with jets
1207       Int_t ijet = SelectJet(particle,aodRecJets);//input for jets is TClonesArray
1208       if(ijet > -1){
1209         //isJetFound=kTRUE;
1210         if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",ijet);
1211         AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1212         if(jet)particle->SetRefJet(jet);
1213         //printf("Most opposite found\n");
1214       }
1215     } // input aod loop
1216     //  if(GetReader()->WriteDeltaAODToFile() && isJetFound) WriteJetsToOutputBranch(aodRecJets);
1217   }//end of take most opposite photon and jet after bkg subtraction
1218
1219
1220   if(GetDebug() > 1 ) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
1221 }
1222
1223 //__________________________________________________________________
1224 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms()
1225 {
1226   //Particle-Jet Correlation Analysis, fill histograms
1227   if(GetDebug() > 3 ) {
1228     printf("I use MakeAnalysisFillHistograms\n");
1229     printf("ntracks before iso %d\n",GetCTSTracks()->GetEntriesFast() );
1230   }
1231
1232   //Get the event, check if there are AODs available, if not, skip this analysis
1233   AliAODEvent * event = NULL;
1234
1235   //printf("GetReader()->GetOutputEvent() %d\n",GetReader()->GetOutputEvent() );
1236   //printf("GetReader()->GetDataType() %d\n",GetReader()->GetDataType() );
1237   //printf("AliCaloTrackReader::kAOD %d\n",AliCaloTrackReader::kAOD );
1238   //printf("GetReader()->GetInputEvent() %d\n",GetReader()->GetInputEvent() );
1239
1240   if(GetReader()->GetOutputEvent())
1241   {
1242     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1243   }
1244   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1245   {
1246     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1247   }
1248   else {
1249     if(GetDebug() > 3) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - There are no jets available for this analysis\n");
1250     return;
1251   }
1252
1253   if(!GetInputAODBranch() || !event){
1254
1255     AliFatal(Form("No input particles in AOD with name branch < %s > \n",
1256                   GetInputAODName().Data()));
1257     return; // Trick coverity
1258   }
1259
1260   Int_t nJets=-1;
1261   TClonesArray *aodRecJets = 0;
1262   //if(IsNonStandardJetFromReader()){//branch read in reader from reader
1263   if (GetDebug() > 3) printf("GetNonStandardJets function (from reader) is called\n");
1264   aodRecJets = GetNonStandardJets();
1265   if(aodRecJets==0x0){
1266     if(GetDebug() > 3) event->Print();
1267     AliFatal("Jets container not found\n");
1268     return; // trick coverity
1269   }
1270   nJets=aodRecJets->GetEntries();
1271   //}
1272   if(nJets==0) {
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     if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
1287     aodBkgJets = GetBackgroundJets();
1288     if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
1289     if(aodBkgJets==0x0){
1290       if(GetDebug() > 3) event->Print();
1291       AliFatal("No jet background container found");
1292       return; // trick coverity
1293     }
1294     if(GetDebug() > 3) aodBkgJets->Print("c");
1295   }
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     if(GetDebug()>5) printf("ChargedBgEnergy %f EffectiveAreaCharged %f\n", jettmp->ChargedBgEnergy(),jettmp->EffectiveAreaCharged());
1422     for(iCounter=1;iCounter<10;iCounter++){
1423       if(jetPttmp>iCounter) nJetsOverThreshold[iCounter]++;
1424     }
1425
1426     var1tmp  = jettmp->Phi();
1427     var2tmp  = jettmp->Eta();
1428     fhJetPhi->Fill(var1tmp);
1429     fhJetEta->Fill(var2tmp);
1430     fhJetPhiVsEta->Fill(var1tmp,var2tmp);
1431     fhJetEtaVsPt->Fill(jetPttmp,var2tmp);
1432     fhJetEtaVsNpartInJet->Fill(var2tmp,nTracksInJet);
1433     if(jetPttmp>0)
1434       fhJetEtaVsNpartInJetBkg->Fill(var2tmp,nTracksInJet);
1435
1436   }//end of jet loop
1437   if(IsHistogramJetTracks()){
1438     if(sumNJetTrack>0){
1439       //printf("average track pt %f\n",sumJetTrackPt/sumNJetTrack);
1440       fhJetAveTrackPt->Fill(sumJetTrackPt/sumNJetTrack);
1441     }
1442   }//end of if IsHistogramJetTracks
1443
1444
1445   fhJetPtMostEne->Fill(ptMostEne);
1446   for(iCounter=0;iCounter<10;iCounter++){
1447     fhJetNjetOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1448     nJetsOverThreshold[iCounter]=0;
1449   }
1450
1451   //end of only jet part
1452
1453   //
1454   // Photons
1455   //
1456   Int_t ntrig   =  GetInputAODBranch()->GetEntriesFast() ;
1457   if(GetDebug() > 1){
1458     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Begin jet finder  correlation analysis, fill histograms \n");
1459     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", ntrig);
1460     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In jet output branch aod entries %d\n", event->GetNJets());
1461   }
1462   fhNjetsNgammas->Fill(nJets,ntrig);
1463
1464   //if(nJets==0)   return;//to speed up
1465   //printf("ntrig %d, njets %d\n",ntrig,nJets);
1466
1467   //
1468   //Fill only temporary photon histograms
1469   //
1470   Double_t maxPt=0.;
1471   Int_t maxIndex=-1;
1472   Double_t tmpPt=0.;
1473   Double_t sumPt=0.;
1474   nJetsOverThreshold[0]=ntrig;
1475   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1476     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1477     tmpPt = particlecorr->Pt();
1478     if(tmpPt>maxPt) {
1479       maxPt = tmpPt;
1480       maxIndex = iaod;
1481     }
1482     for(iCounter=1;iCounter<10;iCounter++){
1483       if(tmpPt>iCounter) nJetsOverThreshold[iCounter]++;
1484     }
1485     sumPt+=tmpPt;
1486   }
1487   for(iCounter=0;iCounter<10;iCounter++){
1488     fhPhotonNgammaOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1489     //    nJetsOverThreshold[iCounter]=0;
1490   }
1491
1492   fGamAvEne=0;
1493   TObjArray* clusters = GetEMCALClusters();
1494   //printf("calculate median bkg energy for photons ");
1495   Double_t medianPhotonRho=0.;
1496   Int_t clusterID;
1497   Int_t iclustmp = -1;
1498   Int_t numberOfcells=0;
1499   AliVCluster *cluster = 0;
1500   if(ntrig>1){
1501     Double_t *photonRhoArr=new Double_t[ntrig-1];
1502     fhPhotonPtMostEne->Fill(maxPt);
1503     //    fhPhotonIndexMostEne->Fill(indexMaxPt);
1504     fhPhotonAverageEnergy->Fill(sumPt/ntrig);
1505     fhPhotonRatioAveEneToMostEne->Fill(sumPt/(ntrig*maxPt));
1506     fhPhotonAverageEnergyMinus1->Fill((sumPt-maxPt)/(ntrig-1));
1507     fGamAvEne=(sumPt-maxPt)/(ntrig-1);
1508     fhPhotonRatioAveEneMinus1ToMostEne->Fill((sumPt-maxPt)/((ntrig-1)*maxPt));
1509
1510     Int_t counterGamma=0;
1511     Int_t counterGammaMinus1=0;
1512
1513     Int_t photonRhoArrayIndex=0;
1514     //TObjArray* clusterstmp = GetEMCALClusters();
1515     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1516       AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1517       if( particlecorr->Pt() > sumPt/ntrig ) counterGamma++;
1518       if( particlecorr->Pt() > (sumPt-maxPt)/(ntrig-1) ) counterGammaMinus1++;
1519
1520       if(iaod==maxIndex) continue;
1521       clusterID = particlecorr->GetCaloLabel(0) ;
1522       if(clusterID < 0) continue;
1523       cluster = FindCluster(clusters,clusterID,iclustmp);
1524       photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1525       numberOfcells+=cluster->GetNCells();
1526       photonRhoArrayIndex++;
1527     }
1528     if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1529     delete [] photonRhoArr;
1530     fhPhotonNgammaMoreAverageToNgamma->Fill((Double_t)counterGamma / (Double_t)ntrig);
1531     fhPhotonNgammaMoreAverageMinus1ToNgamma->Fill((Double_t)counterGammaMinus1 / (Double_t)ntrig);
1532   }
1533   //printf("median = %f\n",medianPhotonRho);
1534   fhPhotonBkgRhoVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),medianPhotonRho);
1535   fhPhotonBkgRhoVsNclusters->Fill(ntrig,medianPhotonRho);
1536   fhPhotonBkgRhoVsCentrality->Fill(GetEventCentrality(),medianPhotonRho);
1537   fhPhotonBkgRhoVsNcells->Fill(numberOfcells,medianPhotonRho);
1538
1539
1540   AliVCluster *cluster2 = 0;
1541   Double_t photon2Corrected=0;
1542   Double_t sumPtTmp=0.;
1543   Double_t sumPtCorrectTmp=0.;
1544   AliVTrack* trackTmp = 0x0 ;
1545   TVector3 p3Tmp;
1546
1547   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1548     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1549     clusterID = particlecorr->GetCaloLabel(0) ;
1550     if(clusterID < 0) continue;
1551     cluster = FindCluster(clusters,clusterID,iclustmp);
1552     fhPhotonPt->Fill(particlecorr->Pt());
1553     fhPhotonPtCorrected->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1554     fhPhotonPtDiff->Fill(cluster->GetNCells() * medianPhotonRho);
1555     fhPhotonPtDiffVsCentrality->Fill(GetEventCentrality(),cluster->GetNCells() * medianPhotonRho);
1556     fhPhotonPtDiffVsNcells->Fill(numberOfcells,cluster->GetNCells() * medianPhotonRho);
1557     fhPhotonPtDiffVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),cluster->GetNCells() * medianPhotonRho);
1558     fhPhotonPtDiffVsNclusters->Fill(ntrig,cluster->GetNCells() * medianPhotonRho);
1559
1560     fhPhotonPtCorrectedZoom->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1561
1562     //test: sum_pt in the cone 0.3 for each photon
1563     //should be: random fake gamma from MB
1564     //is: each gamma for EMCEGA
1565     sumPtTmp=0.;
1566     sumPtCorrectTmp=0.;
1567
1568     for(Int_t iaod2 = 0; iaod2 < ntrig ; iaod2++){
1569       if(iaod==iaod2) continue;
1570       AliAODPWG4ParticleCorrelation* particlecorr2 =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod2));
1571       clusterID = particlecorr2->GetCaloLabel(0) ;
1572       if(clusterID < 0) continue;
1573       cluster2 = FindCluster(clusters,clusterID,iclustmp);
1574       photon2Corrected = particlecorr2->Pt() - cluster2->GetNCells() * medianPhotonRho;
1575
1576       //if(Pt()<0.5) continue; //<<hardcoded here //FIXME
1577       if( TMath::Sqrt((particlecorr->Eta()-particlecorr2->Eta())*(particlecorr->Eta()-particlecorr2->Eta()) +
1578                       (particlecorr->Phi()-particlecorr2->Phi())*(particlecorr->Phi()-particlecorr2->Phi()) )<fGammaConeSize ){//if(/*cone is correct*/){
1579         sumPtTmp+= particlecorr2->Pt();
1580         sumPtCorrectTmp+=photon2Corrected;
1581       }
1582     }
1583     fhPhotonSumPtInCone->Fill(sumPtTmp);
1584     fhPhotonSumPtCorrectInCone->Fill(sumPtCorrectTmp);
1585
1586     //test: sum_pt in the cone 0.3 for each track
1587     //should be: random fake gamma from MB
1588     //is: each gamma for EMCEGA
1589     sumPtTmp=0.;
1590     for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++){
1591       trackTmp = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1592       p3Tmp.SetXYZ(trackTmp->Px(),trackTmp->Py(),trackTmp->Pz());
1593       if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1594                       (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1595         sumPtTmp+=p3Tmp.Pt();
1596       }
1597     }//end of loop over tracks
1598     fhPhotonSumPtChargedInCone->Fill(sumPtTmp);
1599   }
1600
1601   //End of Fill temporary photon histograms
1602
1603   //
1604   // Apply background subtraction for photons
1605   //
1606   fGamRho = medianPhotonRho;
1607   if(!IsBackgroundSubtractionGamma()) medianPhotonRho=0;
1608
1609
1610   //
1611   //Get vertex for cluster momentum calculation <<----new here
1612   //
1613   Double_t vertex[] = {0,0,0} ; //vertex ;
1614   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1615     GetReader()->GetVertex(vertex);
1616   fZvertex = vertex[2];
1617
1618   //
1619   //Loop on stored AOD particles, trigger
1620   //
1621   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1622     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1623     fhCuts->Fill(0);
1624     fhCuts2->Fill(0.,(Double_t)nJets);
1625     if(GetDebug() > 5) printf("OnlyIsolated %d  !particlecorr->IsIsolated() %d \n",OnlyIsolated(), !particlecorr->IsIsolated());
1626
1627     if(OnlyIsolated() && !particlecorr->IsIsolated()) continue;
1628     fhCuts->Fill(1);
1629     fhCuts2->Fill(1.,nJets);
1630
1631     if(nJets>0) {
1632       fhCuts->Fill(2);
1633     }
1634
1635     //Recover the jet correlated, found previously.
1636     AliAODJet   * jet = particlecorr->GetJet();
1637     //If correlation not made before, do it now.
1638     if(fMakeCorrelationInHistoMaker){
1639       //Correlate with jets
1640       Int_t ijet = SelectJet(particlecorr,aodRecJets);//input for jets is TClonesArray
1641       if(ijet > -1){
1642         if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Jet with index %d selected \n",ijet);
1643         //jet = event->GetJet(ijet);
1644         jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1645
1646        if(jet) particlecorr->SetRefJet(jet);
1647
1648       }
1649     }
1650
1651     if (!jet) continue ;
1652     fhCuts->Fill(3);
1653     fhCuts2->Fill(7.,1.);
1654
1655     //check MC genereted information
1656     if(fMCStudies) FindMCgenInfo();
1657
1658     //
1659     //Fill Correlation Histograms
1660     //
1661     clusterID = particlecorr->GetCaloLabel(0) ;
1662     if(!(clusterID<0)){
1663       cluster = FindCluster(clusters,clusterID,iclustmp);
1664       //fill tree variables
1665       fGamNcells = cluster->GetNCells();
1666     }
1667     Double_t ptTrig = particlecorr->Pt() - medianPhotonRho * fGamNcells;//<<---changed here
1668     Double_t ptJet = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();//<<---changed here
1669     Double_t phiTrig = particlecorr->Phi();
1670     Double_t phiJet = jet->Phi();
1671     Double_t etaTrig = particlecorr->Eta();
1672     Double_t etaJet = jet->Eta();
1673     Double_t deltaPhi=phiTrig-phiJet;
1674     if(deltaPhi<0)deltaPhi+=(TMath::Pi()*2.);
1675     //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",
1676     //  ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
1677     fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet);
1678     //    fhDeltaPhi->Fill(ptTrig, phiTrig-phiJet);//need to be shifted by 2pi
1679
1680     fhDeltaPhiCorrect->Fill(ptTrig, deltaPhi);//correct
1681
1682     Double_t deltaPhiCorrect = TMath::Abs( particlecorr->Phi() - jet->Phi() );
1683     if ( deltaPhiCorrect > TMath::Pi() ) deltaPhiCorrect = 2. * TMath::Pi() - deltaPhiCorrect ;
1684     fhDeltaPhi0PiCorrect->Fill(ptTrig, deltaPhiCorrect);
1685
1686     fhDeltaEta->Fill(ptTrig, etaTrig-etaJet);
1687     fhPtRatio ->Fill(ptTrig, ptJet/ptTrig);
1688     fhPt      ->Fill(ptTrig, ptJet);
1689
1690     fhSelectedJetPhiVsEta->Fill(phiJet,etaJet);
1691     fhSelectedJetChBkgEnergyVsPtJet->Fill(ptJet,(IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy()) );
1692     fhSelectedJetChAreaVsPtJet->Fill(ptJet,jet->EffectiveAreaCharged());
1693     fhSelectedJetNjet->Fill(nJets);
1694     fhSelectedNtracks->Fill(GetCTSTracks()->GetEntriesFast());//to be checked
1695     fhSelectedPhotonNLMVsPt->Fill(ptTrig,particlecorr->GetNLM());
1696
1697
1698     if(clusterID < 0 ){
1699       fhSelectedPhotonLambda0VsPt->Fill(ptTrig,-1);
1700       //fill tree variables
1701       fGamLambda0  = -1;
1702       fGamTime = -1;
1703       fGamNcells = 0;
1704       fGamSumPtNeu=0;
1705     }
1706     else{
1707       //Int_t iclus = -1;
1708       //      TObjArray* clusters = GetEMCALClusters();
1709       //cluster = FindCluster(clusters,clusterID,iclustmp);
1710       Double_t lambda0=cluster->GetM02();
1711       fhSelectedPhotonLambda0VsPt->Fill(ptTrig,lambda0);
1712       //fill tree variables
1713       fGamLambda0  = lambda0;
1714       fGamTime = cluster->GetTOF();
1715       //fGamNcells = cluster->GetNCells();
1716
1717       fGamSumPtNeu=0;
1718       fGamNclusters=0;
1719       //TVector3 p3Tmp;
1720       //Double_t scalarProduct=0;
1721       //Double_t vectorLength=particlecorr->P();
1722       for(Int_t icalo=0; icalo <clusters->GetEntriesFast(); icalo++){
1723         AliVCluster* calo = (AliVCluster *) clusters->At(icalo);
1724         if(clusterID==calo->GetID()) continue;//the same cluster as trigger
1725         calo->GetMomentum(fMomentum,vertex) ;//Assume that come from vertex in straight line
1726         //printf("min pt %f\n",GetMinPt());
1727         if(fMomentum.Pt()<GetMinPt()) continue; //<<hardcoded here //FIXME 0.5 check if correct
1728         p3Tmp.SetXYZ(fMomentum.Px(),fMomentum.Py(),fMomentum.Pz());
1729         //calculate sum pt in the cone
1730         if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1731                         (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1732           //scalarProduct = particlecorr->Px()*fMomentum.Px() + particlecorr->Py()*fMomentum.Py() + particlecorr->Pz()*fMomentum.Pz();
1733           //scalarProduct/=fMomentum.P();
1734           //scalarProduct/=vectorLength;
1735           //if(scalarProduct>TMath::Cos(0.3)) {//FIXME photon radius
1736           fGamSumPtNeu+=fMomentum.Pt();
1737           fGamNclusters++;
1738         }
1739       }
1740     }
1741
1742     //sum pt of charged tracks in the gamma isolation cone
1743     //starts here
1744     fGamSumPtCh=0;
1745     fGamNtracks=0;
1746     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1747       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1748       if(!aodtrack) continue;
1749       fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());//fill histogram here
1750       //      if(aodtrack->Pt()<0.15) continue;//hardcoded
1751       if(aodtrack->Pt()<fPtThresholdInCone) continue;
1752       if(!aodtrack->IsHybridGlobalConstrainedGlobal()) continue;
1753       if(TMath::Sqrt((particlecorr->Phi() - aodtrack->Phi())*(particlecorr->Phi() - aodtrack->Phi()) +
1754                      (particlecorr->Eta() - aodtrack->Eta())*(particlecorr->Eta() - aodtrack->Eta()) ) <fGammaConeSize ) {
1755         fGamSumPtCh+=aodtrack->Pt();
1756         fGamNtracks++;
1757       }
1758     }
1759     //ends here
1760
1761     //    for(Int_t itrack = 0; itrack < nCTSTracks ; itrack++){
1762     //      aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1763     //      fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1764     //    }
1765
1766     //
1767     // Background Fragmentation function
1768     //
1769     TVector3 gammaVector,jetVector;
1770     gammaVector.SetXYZ(particlecorr->Px(),particlecorr->Py(),particlecorr->Pz());
1771     jetVector.SetXYZ(jet->Px(),jet->Py(),jet->Pz());
1772     CalculateBkg(gammaVector,jetVector,vertex,1);//jet perp
1773     CalculateBkg(gammaVector,jetVector,vertex,2);//RC
1774     CalculateBkg(gammaVector,jetVector,vertex,3);//mid point
1775     CalculateBkg(gammaVector,jetVector,vertex,4);//gamma perp
1776     //CalculateBkg(gammaVector,jetVector,vertex,5);/test
1777     Double_t angleJetGam = gammaVector.Angle(jetVector);
1778     //printf("angleJetGam %f\n",angleJetGam*180/TMath::Pi());
1779
1780     //
1781     // Fragmentation function
1782     //
1783     Float_t      rad = 0, pt = 0, eta = 0, phi = 0;
1784     Int_t        npartcone = 0;
1785     TVector3 p3;
1786
1787     Int_t ntracks =  0;
1788     if(GetDebug()>3){
1789       printf ("fUseJetRefTracks %d\n",fUseJetRefTracks );
1790       printf ("jet->GetRefTracks() %p\n",jet->GetRefTracks());
1791       printf ("GetCTSTracks() %p\n",GetCTSTracks() );
1792     }
1793
1794     if(!fUseJetRefTracks)
1795       ntracks =GetCTSTracks()->GetEntriesFast();
1796     else //If you want to use jet tracks from JETAN
1797       ntracks =  (jet->GetRefTracks())->GetEntriesFast();
1798
1799     if(GetDebug()>3)    printf ("ntracks %d\n",ntracks);
1800     AliVTrack* track = 0x0 ;
1801     for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
1802       if(!fUseJetRefTracks)
1803         track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1804       else //If you want to use jet tracks from JETAN
1805         track = (AliVTrack *) ((jet->GetRefTracks())->At(ipr));
1806
1807       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1808       pt    = p3.Pt();
1809       eta  = p3.Eta();
1810       phi  = p3.Phi() ;
1811       if(phi < 0) phi+=TMath::TwoPi();
1812
1813       //Check if there is any particle inside cone with pt larger than  fPtThreshold
1814       rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
1815       if(rad < fConeSize  && pt > fPtThresholdInCone){
1816         //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
1817         npartcone++;
1818         fhFFz ->Fill(ptTrig, pt/ptTrig);
1819         fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt));
1820         fhFFpt->Fill(ptTrig, pt);
1821
1822         //according to jet axis
1823         fhJetFFz ->Fill(ptJet, pt/ptJet);
1824         fhJetFFxi->Fill(ptJet, TMath::Log(ptJet/pt));
1825         fhJetFFpt->Fill(ptJet, pt);
1826
1827
1828         if(TMath::Cos(angleJetGam)<0 && ptJet!=0 && pt!=0 ){
1829           fhJetFFzCor ->Fill(ptJet, -pt*TMath::Cos(angleJetGam)/ptJet);
1830           fhJetFFxiCor->Fill(ptJet, TMath::Log(ptJet/(-pt*TMath::Cos(angleJetGam))));
1831         }
1832       }
1833     }//Tracks
1834     fhNTracksInCone->Fill(ptTrig, npartcone);
1835     //fill tree here for each photon-jet (isolation required usually)
1836
1837     fGamPt      = ptTrig;
1838     //fGamLambda0  = ;//filled earlier
1839     fGamNLM      = particlecorr->GetNLM();
1840     //fGamSumPtCh  = ;//filled earlier
1841     //fGamTime     = particlecorr->GetTOF();//filled earlier
1842     //fGamNcells   = particlecorr->GetNCells();//filled earlier
1843     fGamEta      = etaTrig;
1844     fGamPhi      = phiTrig;
1845     //fGamSumPtNeu = ;//filled earlier
1846     //fGamNtracks  = ;//filled earlier
1847     //fGamNclusters= ;//filled earlier
1848     //fGamAvEne    = ;//filled earlier
1849     fJetPhi      = phiJet;
1850     fJetEta      = etaJet;
1851     fJetPt       = ptJet;
1852     fJetBkgChEne = (IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy());
1853     fJetArea     = jet->EffectiveAreaCharged();
1854     fJetNtracks  = (jet->GetRefTracks())->GetEntriesFast();
1855     fEventNumber = 0;
1856     fNtracks     = GetCTSTracks()->GetEntriesFast();
1857     fCentrality  = GetEventCentrality();
1858     fIso         = particlecorr->IsIsolated();
1859
1860     Int_t nTrk1GeV=0;
1861     Int_t nTrk2GeV=0;
1862     for(itrack=0;itrack < fJetNtracks;itrack++){
1863       track = (AliVTrack *) ((jet->GetRefTracks())->At(itrack));
1864       if(track->Pt()>1.) nTrk1GeV++;
1865       if(track->Pt()>2.) nTrk2GeV++;
1866     }
1867
1868     fJetNtracks1 = nTrk1GeV;
1869     fJetNtracks2 = nTrk2GeV;
1870
1871     if(fSaveGJTree) fTreeGJ->Fill();
1872   }//AOD trigger particle loop
1873   if(GetDebug() > 1) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
1874 }
1875
1876
1877 //__________________________________________________________________
1878 void AliAnaParticleJetFinderCorrelation::Print(const Option_t * opt) const
1879 {
1880
1881   //Print some relevant parameters set for the analysis
1882   if(! opt)
1883     return;
1884
1885   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1886   AliAnaCaloTrackCorrBaseClass::Print(" ");
1887
1888   printf("Phi trigger-jet        <     %3.2f\n", fDeltaPhiMaxCut) ;
1889   printf("Phi trigger-jet        >     %3.2f\n", fDeltaPhiMinCut) ;
1890   printf("pT Ratio trigger/jet   <     %3.2f\n", fRatioMaxCut) ;
1891   printf("pT Ratio trigger/jet   >     %3.2f\n", fRatioMinCut) ;
1892   printf("fConeSize              =     %3.2f\n", fConeSize) ;
1893   printf("fPtThresholdInCone     =     %3.2f\n", fPtThresholdInCone) ;
1894   printf("fUseJetRefTracks         =     %d\n",    fUseJetRefTracks) ;
1895   printf("fMakeCorrelationInHistoMaker     =     %d\n",    fMakeCorrelationInHistoMaker) ;
1896   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
1897   printf("Reconstructed jet cone size = %3.2f\n", fJetConeSize) ;
1898   printf("Reconstructed jet minimum pt before background subtraction = %3.2f\n", fJetMinPt) ;
1899   printf("Reconstructed jet minimum area fraction = %3.2f\n", fJetAreaFraction) ;
1900
1901   //if(!fNonStandardJetFromReader){
1902   printf("fJetBranchName =   %s\n", fJetBranchName.Data()) ;
1903   //}
1904   if(!fBackgroundJetFromReader){
1905     printf("fBkgJetBranchName =   %s\n", fBkgJetBranchName.Data()) ;
1906   }
1907
1908   printf("Isolation cone size = %3.2f\n", fGammaConeSize) ;
1909   printf("fUseBackgroundSubtractionGamma = %d\n",fUseBackgroundSubtractionGamma);
1910   printf("fSaveGJTree = %d\n",fSaveGJTree);
1911   printf("fMostEnergetic = %d\n",fMostEnergetic);
1912   printf("fMostOpposite = %d\n",fMostOpposite);
1913
1914   printf("fUseHistogramJetBkg = %d\n",fUseHistogramJetBkg);
1915   printf("fUseHistogramTracks = %d\n",fUseHistogramTracks);
1916   printf("fUseHistogramJetTracks = %d\n",fUseHistogramJetTracks);
1917   printf("fMCStudies = %d\n",fMCStudies);
1918
1919 }
1920
1921 //__________________________________________________________________
1922 void AliAnaParticleJetFinderCorrelation::CalculateBkg(TVector3 gamma, TVector3 jet,Double_t vertex[3],Int_t type=2){
1923   //
1924   // calculate background for fragmentation function and fill histograms
1925   // 1. 90 degrees from jet axis in random place = perpendicular cone
1926   // 2. Random cone not belonging to jet cone nor photon cone
1927   // 3. In the middle point from jet and photon momentum vectors
1928   // 4. 90 degrees from photon direction in random place = perpendicular cone 2
1929
1930   //
1931   // implementation of 2 works, 1 and 4 works
1932   //
1933   Double_t gammaPt  = gamma.Pt();
1934   Double_t gammaEta = gamma.Eta();
1935   Double_t gammaPhi = gamma.Phi();
1936   Double_t jetEta   = jet.Eta();
1937   Double_t jetPhi   = jet.Phi();
1938
1939   //refference direction of background
1940   Double_t refEta=0.,refPhi=0.;
1941   Double_t rad = 0,rad2 = 0.;
1942   if(type==1){//perpendicular to jet axis
1943     //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
1944
1945     Double_t xVar;
1946     Double_t yVar;
1947     Double_t newX=0.;
1948     Double_t newY=0.;
1949     Double_t newZ=0.;
1950     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1951     Double_t jx=jet.Px();
1952     Double_t jy=jet.Py();
1953     Double_t jz=jet.Pz();
1954     //if(jz==0)printf("problem\n");
1955     //X axis
1956     Double_t Xx=1.0-vertex[0];
1957     Double_t Xy=-1.0*vertex[1];
1958     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1959     //Y axis
1960     Double_t Yx=jy*Xz-jz*Xy;
1961     Double_t Yy=jz*Xx-jx*Xz;
1962     Double_t Yz=jx*Xy-jy*Xx;
1963     //Determinant
1964     Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
1965     if(det==0)printf("problem det==0\n");
1966     Double_t detX = 0.;
1967     Double_t detY = 0.;
1968     Double_t detZ = 0.;
1969
1970 //    Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
1971 //    printf("scalar jet.P o X: %f\n",tmpScalar);
1972 //    tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
1973 //    printf("scalar jet.P o Y: %f\n",tmpScalar);
1974 //    tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
1975 //    printf("scalar X o Y: %f\n",tmpScalar);
1976
1977     TVector3 perp;
1978     //randomise
1979     do{
1980       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1981       //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
1982       xVar=TMath::Cos(refPhi);
1983       yVar=TMath::Sin(refPhi);
1984       //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
1985       //zVar=0 in new surface frame
1986       detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
1987       detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
1988       detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
1989
1990       newX=detX/det;
1991       newY=detY/det;
1992       newZ=detZ/det;
1993
1994       perp.SetXYZ(newX,newY,newZ);
1995       refEta = perp.Eta();
1996       refPhi = perp.Phi();//output in <-pi, pi> range
1997       if(refPhi<0)refPhi+=2*TMath::Pi();
1998       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1999       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2000       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2001     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2002     fhRandomPhiEta[2]->Fill(refPhi,refEta);
2003
2004   }
2005   else if(type==2){//random cone
2006     //randomise
2007     do{
2008       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2009       refEta=fGenerator->Uniform(-(0.9-fJetConeSize),0.9-fJetConeSize);
2010       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2011       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2012       //check if reference is not within jet cone or gamma cone +0.4
2013       //example: FF fConeSize=0.4, fJetConeSize=0.4, fIsoGammaConeSize=0.3
2014     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize);
2015     //photon:0.7=0.4+0.3; jets:0.8=0.4 +0.4 //rad<fConeSize+fJetConeSize rad2<fConeSize+0.3
2016     fhRandomPhiEta[0]->Fill(refPhi,refEta);
2017   }
2018   else if(type==4){//perpendicular to photon axis
2019     Double_t xVar;
2020     Double_t yVar;
2021     Double_t newX=0.;
2022     Double_t newY=0.;
2023     Double_t newZ=0.;
2024     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
2025     Double_t jx=gamma.Px();
2026     Double_t jy=gamma.Py();
2027     Double_t jz=gamma.Pz();
2028     //if(jz==0)printf("problem\n");
2029     //X axis
2030     Double_t Xx=1.0-vertex[0];
2031     Double_t Xy=-1.0*vertex[1];
2032     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2033     //Y axis
2034     Double_t Yx=jy*Xz-jz*Xy;
2035     Double_t Yy=jz*Xx-jx*Xz;
2036     Double_t Yz=jx*Xy-jy*Xx;
2037     //Determinant
2038     Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
2039     if(det==0)printf("problem det==0\n");
2040     Double_t detX = 0.;
2041     Double_t detY = 0.;
2042     Double_t detZ = 0.;
2043
2044 //    Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
2045 //    printf("scalar jet.P o X: %f\n",tmpScalar);
2046 //    tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
2047 //    printf("scalar jet.P o Y: %f\n",tmpScalar);
2048 //    tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
2049 //    printf("scalar X o Y: %f\n",tmpScalar);
2050
2051     TVector3 perp;
2052     //randomise
2053     do{
2054       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2055       //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
2056       xVar=TMath::Cos(refPhi);
2057       yVar=TMath::Sin(refPhi);
2058       //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
2059       //zVar=0 in new surface frame
2060       detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
2061       detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
2062       detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
2063
2064       newX=detX/det;
2065       newY=detY/det;
2066       newZ=detZ/det;
2067
2068       perp.SetXYZ(newX,newY,newZ);
2069       refEta = perp.Eta();
2070       refPhi = perp.Phi();//output in <-pi, pi> range
2071       if(refPhi<0)refPhi+=2*TMath::Pi();
2072       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2073       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2074       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2075     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2076     fhRandomPhiEta[1]->Fill(refPhi,refEta);
2077
2078   }
2079   else if(type==3){//mid point
2080
2081     Double_t jx=jet.Px();
2082     Double_t jy=jet.Py();
2083     Double_t jz=jet.Pz();
2084     //    if(jz==0)printf("problem\n");
2085     Double_t gx=gamma.Px();
2086     Double_t gy=gamma.Py();
2087     Double_t gz=gamma.Pz();
2088
2089     Double_t cosAlpha=(jx*gx+jy*gy+jz*gz)/(jet.Mag()*gamma.Mag());
2090     Double_t cosinus=TMath::Sqrt((cosAlpha+1.)/2.);
2091     //perpendicular axis
2092     Double_t Zx=gy*jz-gz*jy;
2093     Double_t Zy=gz*jx-gx*jz;
2094     Double_t Zz=gx*jy-gy*jx;
2095
2096     //Determinant
2097     Double_t det = Zx*gy*jz + Zy*gz*jx + Zz*gx*jy - Zz*gy*jx - Zy*gx*jz - Zx*gz*jy;
2098
2099     Double_t newX=0.;
2100     Double_t newY=0.;
2101     Double_t newZ=0.;
2102     if(det!=0) {
2103       Double_t detX =            -Zy*gz*cosinus +Zz*cosinus*jy + Zz*gy*cosinus - Zy*cosinus*jz;
2104       Double_t detY = Zx*cosinus*jz          - Zz*gx*cosinus - Zz*cosinus*jx            + Zx*gz*cosinus;
2105       Double_t detZ = -Zx*gy*cosinus + Zy*cosinus*jx + Zy*gx*cosinus - Zx*cosinus*jy;
2106
2107       newX=detX/det;
2108       newY=detY/det;
2109       newZ=detZ/det;
2110     }
2111
2112     TVector3 perp;
2113     perp.SetXYZ(newX,newY,newZ);
2114     refEta = perp.Eta();
2115     refPhi = perp.Phi();//output in <-pi, pi> range
2116     if(refPhi<0)refPhi+=2*TMath::Pi();
2117     rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2118     rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2119       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2120
2121     if (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize) fhRandomPhiEta[3]->Fill(refPhi,refEta);
2122   }
2123   else if(type==5){//tmp
2124     //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
2125
2126     Double_t xVar;
2127     Double_t newX=0.;
2128     Double_t newY=0.;
2129     Double_t newZ=0.;
2130     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
2131     Double_t jx=jet.Px();
2132     Double_t jy=jet.Py();
2133     Double_t jz=jet.Pz();
2134     //    if(jz==0)printf("problem\n");
2135     //X axis
2136     Double_t Xx=1.0-vertex[0];
2137     Double_t Xy=-1.0*vertex[1];
2138     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
2139     //Y axis
2140     Double_t Yx=jy*Xz-jz*Xy;
2141     Double_t Yy=jz*Xx-jx*Xz;
2142     Double_t Yz=jx*Xy-jy*Xx;
2143
2144     // X and Y length
2145     Double_t Xlength=TMath::Sqrt(Xx*Xx+Xy*Xy+Xz*Xz);
2146     Double_t Ylength=TMath::Sqrt(Yx*Yx+Yy*Yy+Yz*Yz);
2147     Double_t ratio=Ylength/Xlength;
2148
2149     TVector3 perp;
2150     //randomise
2151     do{
2152       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
2153       xVar=TMath::Tan(refPhi)/ratio;
2154       newX=xVar*Yx+Xx;
2155       newY=xVar*Yy+Xy;
2156       newZ=xVar*Yz+Xz;
2157
2158       perp.SetXYZ(newX,newY,newZ);
2159       refEta = perp.Eta();
2160       refPhi = perp.Phi();//output in <-pi, pi> range
2161       if(refPhi<0)refPhi+=2*TMath::Pi();
2162       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2163       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2164       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2165     } while (rad<fJetConeSize+fGammaConeSize || rad2<2.*fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2166     fhRandomPhiEta[4]->Fill(refPhi,refEta);
2167   }
2168
2169
2170
2171   //calculate FF in background
2172   Int_t ntracks =  0;
2173   ntracks =GetCTSTracks()->GetEntriesFast();
2174   AliVTrack* track = 0x0 ;
2175   TVector3 p3;
2176
2177   Double_t pt = 0, eta = 0, phi = 0;
2178   Int_t  npartcone = 0;
2179   Double_t sumPt=0.;
2180   for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
2181     track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
2182     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
2183     pt   = p3.Pt();
2184     if(pt<fPtThresholdInCone) {//0.150
2185       //printf("problem: track pt < %f MeV/c \n",fPtThresholdInCone);
2186       continue;
2187     }
2188     eta  = p3.Eta() ;
2189     phi  = p3.Phi() ;
2190     if(phi < 0) phi+=TMath::TwoPi();
2191     //Check if there is any particle inside cone with pt larger than  fPtThreshold
2192     rad = TMath::Sqrt((eta-refEta)*(eta-refEta) + (phi-refPhi)*(phi-refPhi));
2193     if(rad < fConeSize  && pt > fPtThresholdInCone){
2194       //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
2195         npartcone++;
2196         sumPt+=pt;
2197         if(type==1){//perp jet
2198           fhBkgFFz[1] ->Fill(gammaPt, pt/gammaPt);
2199           fhBkgFFxi[1]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2200           fhBkgFFpt[1]->Fill(gammaPt, pt);
2201         }
2202         else if(type==2){//RC
2203           fhBkgFFz[0] ->Fill(gammaPt, pt/gammaPt);
2204           fhBkgFFxi[0]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2205           fhBkgFFpt[0]->Fill(gammaPt, pt);
2206         }
2207         else if(type==3){//mid point
2208           fhBkgFFz[3] ->Fill(gammaPt, pt/gammaPt);
2209           fhBkgFFxi[3]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2210           fhBkgFFpt[3]->Fill(gammaPt, pt);
2211         }
2212         else if(type==4){//perp jet
2213           fhBkgFFz[2] ->Fill(gammaPt, pt/gammaPt);
2214           fhBkgFFxi[2]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2215           fhBkgFFpt[2]->Fill(gammaPt, pt);
2216         }
2217         else if(type==5){//test
2218           fhBkgFFz[4] ->Fill(gammaPt, pt/gammaPt);
2219           fhBkgFFxi[4]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2220           fhBkgFFpt[4]->Fill(gammaPt, pt);
2221         }
2222
2223
2224     }
2225   }//end of loop over tracks
2226   Double_t sumOverTracks=0.;
2227   if(npartcone!=0) sumOverTracks = sumPt/npartcone;
2228   if(type==1) {
2229     fhBkgNTracksInCone[1]->Fill(gammaPt, npartcone);
2230     fhBkgSumPtInCone[1]->Fill(gammaPt,sumPt);
2231     fhBkgSumPtnTracksInCone[1]->Fill(gammaPt,sumOverTracks);
2232   }
2233   else if(type==2) {
2234     fhBkgNTracksInCone[0]->Fill(gammaPt, npartcone);
2235     fhBkgSumPtInCone[0]->Fill(gammaPt,sumPt);
2236     fhBkgSumPtnTracksInCone[0]->Fill(gammaPt,sumOverTracks);
2237   }
2238   else if(type==3) {
2239     fhBkgNTracksInCone[3]->Fill(gammaPt, npartcone);
2240     fhBkgSumPtInCone[3]->Fill(gammaPt,sumPt);
2241     fhBkgSumPtnTracksInCone[3]->Fill(gammaPt,sumOverTracks);
2242   }
2243   else if(type==4) {
2244     fhBkgNTracksInCone[2]->Fill(gammaPt, npartcone);
2245     fhBkgSumPtInCone[2]->Fill(gammaPt,sumPt);
2246     fhBkgSumPtnTracksInCone[2]->Fill(gammaPt,sumOverTracks);
2247   }
2248   else if(type==5) {
2249     fhBkgNTracksInCone[4]->Fill(gammaPt, npartcone);
2250     fhBkgSumPtInCone[4]->Fill(gammaPt,sumPt);
2251     fhBkgSumPtnTracksInCone[4]->Fill(gammaPt,sumOverTracks);
2252   }
2253 }
2254
2255
2256
2257
2258
2259 //__________________________________________________________________
2260 void AliAnaParticleJetFinderCorrelation::FindMCgenInfo(){
2261   //
2262   // Find information about photon and (quark or gluon) on generated level
2263   //
2264
2265   //frequently used variables
2266   Int_t pdg    = 0 ;
2267   Int_t mother = -1 ;
2268   Int_t absID  = 0 ;
2269
2270   //Double_t photonY   = -100 ;
2271   //Double_t photonE   = -1 ;
2272   Double_t photonPt  = -1 ;
2273   Double_t photonPhi =  100 ;
2274   Double_t photonEta = -1 ;
2275   Bool_t   inacceptance = kFALSE;
2276   AliAODMCParticle * primTmp = NULL;
2277
2278   //jet counters
2279   Int_t nParticlesInJet=0;
2280   Int_t nChargedParticlesInJet=0;
2281   Int_t nParticlesInJet150=0;
2282   Int_t nChargedParticlesInJet150=0;
2283   Int_t nChargedParticlesInJet150Cone=0;
2284
2285   Double_t eneParticlesInJet=0.;
2286   Double_t eneChargedParticlesInJet=0.;
2287   Double_t eneParticlesInJet150=0.;
2288   Double_t eneChargedParticlesInJet150=0.;
2289   Double_t eneChargedParticlesInJet150Cone=0.;
2290
2291   Double_t pxParticlesInJet=0.;
2292   Double_t pxChargedParticlesInJet=0.;
2293   Double_t pxParticlesInJet150=0.;
2294   Double_t pxChargedParticlesInJet150=0.;
2295   Double_t pxChargedParticlesInJet150Cone=0.;
2296
2297   Double_t pyParticlesInJet=0.;
2298   Double_t pyChargedParticlesInJet=0.;
2299   Double_t pyParticlesInJet150=0.;
2300   Double_t pyChargedParticlesInJet150=0.;
2301   Double_t pyChargedParticlesInJet150Cone=0.;
2302
2303   Double_t etaParticlesInJet=0.;
2304   Double_t etaChargedParticlesInJet=0.;
2305   Double_t etaParticlesInJet150=0.;
2306   Double_t etaChargedParticlesInJet150=0.;
2307   Double_t etaChargedParticlesInJet150Cone=0.;
2308
2309   Double_t phiParticlesInJet=0.;
2310   Double_t phiChargedParticlesInJet=0.;
2311   Double_t phiParticlesInJet150=0.;
2312   Double_t phiChargedParticlesInJet150=0.;
2313   Double_t phiChargedParticlesInJet150Cone=0.;
2314
2315   Double_t ptParticlesInJet=0.;
2316   Double_t ptChargedParticlesInJet=0.;
2317   Double_t ptParticlesInJet150=0.;
2318   Double_t ptChargedParticlesInJet150=0.;
2319   Double_t ptChargedParticlesInJet150Cone=0.;
2320
2321   Double_t coneJet=0.;
2322   Double_t coneChargedJet=0.;
2323   Double_t coneJet150=0.;
2324   Double_t coneChargedJet150=0.;
2325
2326   std::vector<Int_t> jetParticleIndex;
2327
2328   if(GetReader()->ReadStack()) {//ESD
2329      if(GetDebug()>3) printf("I use stack\n");
2330   }//end Stack
2331   else if(GetReader()->ReadAODMCParticles()) {//AOD
2332     TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
2333     if(mcparticles){
2334       //jet origin
2335       //index =6 and 7 is hard scattering (jet-quark or photon)
2336       primTmp = (AliAODMCParticle *) mcparticles->At(6);
2337       pdg=primTmp->GetPdgCode();
2338        if(GetDebug()>3) printf("id 6 pdg %d, pt %f ",pdg,primTmp->Pt() );
2339       if(TMath::Abs(pdg)<=6 ||pdg==21) {
2340         fhMCJetOrigin->Fill(pdg);
2341         fMCPartonType=pdg;
2342       }
2343       primTmp = (AliAODMCParticle *) mcparticles->At(7);
2344       pdg=primTmp->GetPdgCode();
2345        if(GetDebug()>3) printf("id 7 pdg %d, pt %f\n",pdg,primTmp->Pt() );
2346       if(TMath::Abs(pdg)<=6 ||pdg==21) {
2347         fhMCJetOrigin->Fill(pdg);
2348         fMCPartonType=pdg;
2349       }
2350       //end of jet origin
2351
2352       Int_t nprim = mcparticles->GetEntriesFast();
2353       for(Int_t i=0; i < nprim; i++) {
2354         if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
2355         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
2356         pdg = prim->GetPdgCode();
2357         mother=prim->GetMother();
2358         //photon=22, gluon=21, quarks=(1,...,6), antiquarks=(-1,...,-6)
2359         if(pdg == 22){//photon
2360           fhMCPhotonCuts->Fill(0);
2361           if(prim->GetStatus()!=1) continue;
2362           fhMCPhotonCuts->Fill(1);
2363            if(GetDebug()>5) printf("id %d, prim %d, physPrim %d, status %d\n",i,prim->IsPrimary(),prim->IsPhysicalPrimary(),prim->GetStatus());
2364           while(mother>7){
2365             primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2366             mother=primTmp->GetMother();
2367           }
2368           if(mother<6)continue;
2369           fhMCPhotonCuts->Fill(2);
2370           primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2371           if(primTmp->GetPdgCode()!=22)continue;
2372           fhMCPhotonCuts->Fill(3);
2373
2374           //Get photon kinematics
2375           photonPt  = prim->Pt() ;
2376           photonPhi = prim->Phi() ;
2377           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
2378           photonEta = prim->Eta() ;
2379           fhMCPhotonPt->Fill(photonPt);
2380           fhMCPhotonEtaPhi->Fill(photonPhi,photonEta);
2381
2382           //Check if photons hit the Calorimeter
2383           fMomentum.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
2384           inacceptance = kFALSE;
2385           if(GetCaloUtils()->IsEMCALGeoMatrixSet()){
2386             fhMCPhotonCuts->Fill(4);
2387             //check if in EMCAL
2388             if(GetEMCALGeometry()) {
2389               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
2390               if(absID >= 0) inacceptance = kTRUE;
2391               if(GetDebug() > 3) printf("In EMCAL Real acceptance? %d\n",inacceptance);
2392             }
2393             else{
2394               if(GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) inacceptance = kTRUE ;
2395               if(GetDebug() > 3) printf("In EMCAL fiducial cut acceptance? %d\n",inacceptance);
2396             }
2397           }else{//no EMCAL nor EMCALGeoMatrixSet
2398             printf("not EMCALGeoMatrix set\n");
2399           }//end of check if EMCAL
2400           if(inacceptance)fhMCPhotonCuts->Fill(5);
2401           if(GetDebug() > 5)
2402             printf("Photon Energy %f, Pt %f\n",prim->E(),prim->Pt());
2403           fMCGamPt=photonPt;
2404           fMCGamEta=photonEta;
2405           fMCGamPhi=photonPhi;
2406         }//end of check if photon
2407         else{//not photon
2408           if(prim->GetStatus()!=1) continue;
2409           if(GetDebug() > 5)
2410             printf("id %d, prim %d, physPrim %d, status %d, pdg %d, E %f ",i,prim->IsPrimary(),prim->IsPhysicalPrimary(),prim->GetStatus(),prim->GetPdgCode(),prim->E());
2411           while(mother>7){
2412             primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2413             mother=primTmp->GetMother();
2414             if(GetDebug() > 5) printf("next mother %d ",mother);
2415           }
2416           if(GetDebug() > 5) printf("\n");
2417           if(mother<6)continue;//soft part
2418           primTmp = (AliAODMCParticle *) mcparticles->At(mother);
2419           pdg=primTmp->GetPdgCode();
2420           if( !(TMath::Abs(pdg)<=6 || pdg==21) ) continue;//origin not hard q or g
2421
2422           //jetParticleIndex.Add(&i);
2423           jetParticleIndex.push_back(i);
2424
2425         nParticlesInJet++;
2426         eneParticlesInJet+=prim->E();
2427         pxParticlesInJet+=prim->Px();
2428         pyParticlesInJet+=prim->Py();
2429         etaParticlesInJet+=(prim->E()*prim->Eta());
2430         photonPhi = prim->Phi() ;
2431         if(photonPhi < 0) photonPhi+=TMath::TwoPi();
2432         phiParticlesInJet+=(prim->E()*photonPhi);
2433
2434         if(prim->Charge()!=0) {
2435           nChargedParticlesInJet++;
2436           eneChargedParticlesInJet+=prim->E();
2437           pxChargedParticlesInJet+=prim->Px();
2438           pyChargedParticlesInJet+=prim->Py();
2439           etaChargedParticlesInJet+=(prim->E()*prim->Eta());
2440           phiChargedParticlesInJet+=(prim->E()*photonPhi);
2441         }
2442         if(prim->Pt()>0.150 && TMath::Abs(prim->Eta())<0.9 ) {//in TPC acceptance :)
2443           nParticlesInJet150++;
2444           eneParticlesInJet150+=prim->E();
2445           pxParticlesInJet150+=prim->Px();
2446           pyParticlesInJet150+=prim->Py();
2447           etaParticlesInJet150+=(prim->E()*prim->Eta());
2448           phiParticlesInJet150+=(prim->E()*photonPhi);
2449         }
2450         if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(prim->Eta())<0.9 ) {//in TPC acceptance :)
2451           nChargedParticlesInJet150++;
2452           eneChargedParticlesInJet150+=prim->E();
2453           pxChargedParticlesInJet150+=prim->Px();
2454           pyChargedParticlesInJet150+=prim->Py();
2455           etaChargedParticlesInJet150+=(prim->E()*prim->Eta());
2456           phiChargedParticlesInJet150+=(prim->E()*photonPhi);
2457         }
2458
2459         }//end of check pdg
2460       }//end of loop over primaries
2461
2462       if(eneParticlesInJet != 0.) {
2463         etaParticlesInJet/=eneParticlesInJet ;
2464         phiParticlesInJet/=eneParticlesInJet ;
2465       }
2466       if(eneChargedParticlesInJet != 0) {
2467         etaChargedParticlesInJet/=eneChargedParticlesInJet;
2468         phiChargedParticlesInJet/=eneChargedParticlesInJet;
2469       }
2470       if(eneParticlesInJet150 != 0) {
2471         etaParticlesInJet150/=eneParticlesInJet150;
2472         phiParticlesInJet150/=eneParticlesInJet150;
2473       }
2474       if(eneChargedParticlesInJet150 != 0) {
2475         etaChargedParticlesInJet150/=eneChargedParticlesInJet150;
2476         phiChargedParticlesInJet150/=eneChargedParticlesInJet150;
2477       }
2478
2479       ptParticlesInJet=TMath::Sqrt(pxParticlesInJet*pxParticlesInJet+pyParticlesInJet*pyParticlesInJet);
2480       ptChargedParticlesInJet=TMath::Sqrt(pxChargedParticlesInJet*pxChargedParticlesInJet+pyChargedParticlesInJet*pyChargedParticlesInJet);
2481       ptParticlesInJet150=TMath::Sqrt(pxParticlesInJet150*pxParticlesInJet150+pyParticlesInJet150*pyParticlesInJet150);
2482       ptChargedParticlesInJet150=TMath::Sqrt(pxChargedParticlesInJet150*pxChargedParticlesInJet150+pyChargedParticlesInJet150*pyChargedParticlesInJet150);
2483
2484       Double_t distance=0.;
2485       Double_t eta=0.;
2486       Double_t phi=0.;
2487       Double_t mostPtCharged=0.;
2488       Int_t mostmostPtChargedId=-1;
2489       std::vector<Int_t>::iterator it;
2490       for( it=jetParticleIndex.begin(); it!=jetParticleIndex.end(); ++it ){
2491         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(*it);
2492         eta = prim->Eta();
2493         phi = prim->Phi();
2494         if(phi < 0) phi+=TMath::TwoPi();
2495         //full jet
2496         distance=TMath::Sqrt((eta-etaParticlesInJet)*(eta-etaParticlesInJet)+(phi-phiParticlesInJet)*(phi-phiParticlesInJet));
2497         if(distance>coneJet) coneJet=distance;
2498         //charged jet
2499         distance=TMath::Sqrt((eta-etaChargedParticlesInJet)*(eta-etaChargedParticlesInJet)+(phi-phiChargedParticlesInJet)*(phi-phiChargedParticlesInJet));
2500         if(distance>coneChargedJet) coneChargedJet=distance;
2501         //
2502         distance=TMath::Sqrt((eta-etaParticlesInJet150)*(eta-etaParticlesInJet150)+(phi-phiParticlesInJet150)*(phi-phiParticlesInJet150));
2503         if(distance>coneJet150 && TMath::Abs(eta)<0.9 ) coneJet150=distance;
2504         //
2505         distance=TMath::Sqrt((eta-etaChargedParticlesInJet150)*(eta-etaChargedParticlesInJet150)+(phi-phiChargedParticlesInJet150)*(phi-phiChargedParticlesInJet150));
2506         if(distance>coneChargedJet150 && TMath::Abs(eta)<0.9) coneChargedJet150=distance;
2507
2508         if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(eta)<0.9) {
2509           if(prim->Pt()>mostPtCharged) {
2510             mostPtCharged=prim->Pt();
2511             mostmostPtChargedId=(*it);
2512           }
2513         }
2514
2515         if(distance<=0.4){
2516           if(prim->Charge()!=0 && prim->Pt()>0.150 && TMath::Abs(eta)<0.9) {
2517             nChargedParticlesInJet150Cone++;
2518             eneChargedParticlesInJet150Cone+=prim->E();
2519             pxChargedParticlesInJet150Cone+=prim->Px();
2520             pyChargedParticlesInJet150Cone+=prim->Py();
2521             etaChargedParticlesInJet150Cone+=(prim->E()*eta);
2522             phiChargedParticlesInJet150Cone+=(prim->E()*phi);
2523           }
2524
2525         }
2526
2527       }//end of loop over jet particle indexes
2528       if(eneChargedParticlesInJet150Cone != 0) {
2529         etaChargedParticlesInJet150Cone/=eneChargedParticlesInJet150Cone;
2530         phiChargedParticlesInJet150Cone/=eneChargedParticlesInJet150Cone;
2531       }
2532       ptChargedParticlesInJet150Cone=TMath::Sqrt(pxChargedParticlesInJet150Cone*pxChargedParticlesInJet150Cone+pyChargedParticlesInJet150Cone*pyChargedParticlesInJet150Cone);
2533       if(nChargedParticlesInJet150>0 && nChargedParticlesInJet150Cone<1){//no particles in cone? take the most energetic one
2534         nChargedParticlesInJet150Cone=1;
2535         etaChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Eta();
2536         phiChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Phi();
2537         ptChargedParticlesInJet150Cone=((AliAODMCParticle *)mcparticles->At(mostmostPtChargedId))->Pt();
2538       }
2539
2540
2541     }//mc array exists and data is MC
2542   }// read AOD MC
2543
2544   jetParticleIndex.clear();
2545
2546
2547   //printouts
2548   if(GetDebug() > 3) {
2549     printf("cone full %f, charged %f, full150 %f, charged150 %f\n",coneJet,coneChargedJet,coneJet150,coneChargedJet150);
2550     printf("Npart %d, NchPart %d, Npart(pt>150M) %d, NchPart(pt>150M) %d, NchPart(pt>150M)Cone %d\n",nParticlesInJet,nChargedParticlesInJet,nParticlesInJet150,nChargedParticlesInJet150,nChargedParticlesInJet150Cone);
2551     printf("Etot %f, Ech %f, E(pt>150M) %f, Ech(pt>150M) %f\n",eneParticlesInJet,eneChargedParticlesInJet,eneParticlesInJet150,eneChargedParticlesInJet150);
2552     printf("pt %f, ptch %f, pt(pt>150M) %f,ptch(pt>150M) %f,ptch(pt>150M)Cone %f\n",ptParticlesInJet,ptChargedParticlesInJet,ptParticlesInJet150,ptChargedParticlesInJet150,ptChargedParticlesInJet150Cone);
2553     printf("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);
2554   }
2555   //fill histograms
2556   if(ptParticlesInJet) fhMCJetRatioChFull->Fill(ptChargedParticlesInJet/ptParticlesInJet);
2557   if(ptChargedParticlesInJet) fhMCJetRatioCh150Ch->Fill(ptChargedParticlesInJet150/ptChargedParticlesInJet);
2558
2559   fhMCJetNPartVsPt     ->Fill(ptParticlesInJet,nParticlesInJet);
2560   fhMCJetChNPartVsPt   ->Fill(ptChargedParticlesInJet,nChargedParticlesInJet);
2561   fhMCJetNPart150VsPt  ->Fill(ptParticlesInJet150,nParticlesInJet150);
2562   fhMCJetChNPart150VsPt->Fill(ptChargedParticlesInJet150,nChargedParticlesInJet150);
2563   fhMCJetChNPart150ConeVsPt->Fill(ptChargedParticlesInJet150Cone,nChargedParticlesInJet150Cone);
2564
2565   fhMCJetEtaPhi->Fill(phiParticlesInJet,etaParticlesInJet);
2566   fhMCJetChEtaPhi->Fill(phiChargedParticlesInJet,etaChargedParticlesInJet);
2567   fhMCJet150EtaPhi->Fill(phiParticlesInJet150,etaParticlesInJet150);
2568   fhMCJetCh150EtaPhi->Fill(phiChargedParticlesInJet150,etaChargedParticlesInJet150);
2569   fhMCJetCh150ConeEtaPhi->Fill(phiChargedParticlesInJet150Cone,etaChargedParticlesInJet150Cone);
2570
2571   //fill tree
2572   fMCJetPt      = ptParticlesInJet;
2573   fMCJetChPt    = ptChargedParticlesInJet;
2574   fMCJet150Pt   = ptParticlesInJet150;
2575   fMCJetCh150Pt = ptChargedParticlesInJet150;
2576   fMCJetNPart      = nParticlesInJet;
2577   fMCJetChNPart    = nChargedParticlesInJet;
2578   fMCJet150NPart   = nParticlesInJet150;
2579   fMCJetCh150NPart = nChargedParticlesInJet150;
2580   fMCJetEta      = etaParticlesInJet          ;
2581   fMCJetPhi      = phiParticlesInJet          ;
2582   fMCJetChEta    = etaChargedParticlesInJet   ;
2583   fMCJetChPhi    = phiChargedParticlesInJet   ;
2584   fMCJet150Eta   = etaParticlesInJet150       ;
2585   fMCJet150Phi   = phiParticlesInJet150       ;
2586   fMCJetCh150Eta = etaChargedParticlesInJet150;
2587   fMCJetCh150Phi = phiChargedParticlesInJet150;
2588
2589   fMCJetCh150ConePt    = ptChargedParticlesInJet150Cone;
2590   fMCJetCh150ConeNPart = nChargedParticlesInJet150Cone;
2591   fMCJetCh150ConeEta   = etaChargedParticlesInJet150Cone;
2592   fMCJetCh150ConePhi   = phiChargedParticlesInJet150Cone;
2593
2594 }//end of method FindMCgenInfo