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