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