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