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 //////////////////////////////////////////////////////////////////////////////
  20
  21
  22 // --- ROOT system ---
  23 #include "TH2F.h"
  24 #include "TClonesArray.h"
  25 #include "TClass.h"
  26 //#include "Riostream.h"
  27
  28 //---- AliRoot system ----
  29 #include "AliCaloTrackReader.h"
  30 #include "AliAODJet.h"
  31 #include "AliAnaParticleJetFinderCorrelation.h"
  32 #include "AliAODPWG4ParticleCorrelation.h"
  33 #include "AliVTrack.h"
  34 #include "AliAODCaloCluster.h"
  35 #include "AliAODEvent.h"
  36
  37 //jets
  38 #include "AliAODJetEventBackground.h"
  39 #include "TRandom2.h"
  40
  41 ClassImp(AliAnaParticleJetFinderCorrelation)
  42
  43
  44 //________________________________________________________________________
  45 AliAnaParticleJetFinderCorrelation::AliAnaParticleJetFinderCorrelation() :
  46 AliAnaCaloTrackCorrBaseClass(),
  47   fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fRatioMaxCut(0.),  fRatioMinCut(0.),
  48   fConeSize(0), fPtThresholdInCone(0),fUseJetRefTracks(0),
  49   fMakeCorrelationInHistoMaker(0), fSelectIsolated(0),
  50   fJetConeSize(0.4),fJetMinPt(5),fJetAreaFraction(0.6),
  51   fNonStandardJetFromReader(kTRUE), fJetBranchName("jets"),
  52   fBackgroundJetFromReader(kTRUE),fBkgJetBranchName("jets"),
  53   fGammaConeSize(0.3),fUseBackgroundSubtractionGamma(0),fSaveGJTree(0),
  54   fMostEnergetic(kFALSE),fMostOpposite(kTRUE), fUseHistogramJetBkg(kTRUE),
  55   fUseHistogramTracks(kTRUE),fUseHistogramJetTracks(kTRUE),fGenerator(0),
  56   fhDeltaEta(0), /*fhDeltaPhi(0),*/fhDeltaPhiCorrect(0),fhDeltaPhi0PiCorrect(0), fhDeltaPt(0), fhPtRatio(0), fhPt(0),
  57   fhFFz(0),fhFFxi(0),fhFFpt(0),fhNTracksInCone(0),
  58   fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetFFzCor(0),fhJetFFxiCor(0),
  59   fhBkgFFz(),fhBkgFFxi(),fhBkgFFpt(),fhBkgNTracksInCone(),fhBkgSumPtInCone(),fhBkgSumPtnTracksInCone(),//<<---new
  60   fhNjetsNgammas(0),fhCuts(0),
  61   fhDeltaEtaBefore(0),fhDeltaPhiBefore(0),fhDeltaPtBefore(0),fhPtRatioBefore(0),
  62   fhPtBefore(0),fhDeltaPhi0PiCorrectBefore(0),
  63   fhJetPtBefore(0),fhJetPt(0),fhJetPtMostEne(0),fhJetPhi(0),fhJetEta(0),fhJetEtaVsPt(0),
  64   fhJetPhiVsEta(0),fhJetEtaVsNpartInJet(0),fhJetEtaVsNpartInJetBkg(0),fhJetChBkgEnergyVsPt(0),fhJetChAreaVsPt(0),/*fhJetNjet(0),*/
  65   fhTrackPhiVsEta(0),fhTrackAveTrackPt(0),fhJetNjetOverPtCut(),
  66 /*fhJetChBkgEnergyVsPtEtaGt05(0),fhJetChBkgEnergyVsPtEtaLe05(0),fhJetChAreaVsPtEtaGt05(0),fhJetChAreaVsPtEtaLe05(0),*/
  67   fhJetChBkgEnergyVsArea(0),fhJetRhoVsPt(0),fhJetRhoVsCentrality(0),//fhJetBkgRho(0),
  68   fhJetNparticlesInJet(0),fhJetDeltaEtaDeltaPhi(0),fhJetDeltaEtaDeltaPhiAllTracks(0),
  69   fhJetAveTrackPt(0),fhJetNtracksInJetAboveThr(),fhJetRatioNTrkAboveToNTrk(),fhJetNtrackRatioMostEne(),
  70   fhJetNtrackRatioJet5GeV(),fhJetNtrackRatioLead5GeV(),
  71   fhBkgJetBackground(),fhBkgJetSigma(),fhBkgJetArea(),fhPhotonPtMostEne(0),
  72   fhPhotonAverageEnergy(0),fhPhotonRatioAveEneToMostEne(0),fhPhotonAverageEnergyMinus1(0),fhPhotonRatioAveEneMinus1ToMostEne(0),
  73   fhPhotonNgammaMoreAverageToNgamma(0),fhPhotonNgammaMoreAverageMinus1ToNgamma(0),fhPhotonNgammaOverPtCut(),
  74   fhPhotonBkgRhoVsNtracks(0),fhPhotonBkgRhoVsNclusters(0),fhPhotonBkgRhoVsCentrality(0),
  75   fhPhotonBkgRhoVsNcells(0),fhPhotonPt(0),fhPhotonPtCorrected(0),fhPhotonPtCorrectedZoom(0),fhPhotonPtDiff(0),
  76   fhPhotonPtDiffVsCentrality(0),fhPhotonPtDiffVsNcells(0),fhPhotonPtDiffVsNtracks(0),fhPhotonPtDiffVsNclusters(0),
  77   fhPhotonSumPtInCone(0),fhPhotonSumPtCorrectInCone(0),fhPhotonSumPtChargedInCone(0),
  78   fhSelectedJetPhiVsEta(0),fhSelectedJetChBkgEnergyVsPtJet(0),fhSelectedJetChAreaVsPtJet(0),fhSelectedJetNjet(0),fhSelectedNtracks(0),
  79   fhSelectedTrackPhiVsEta(0),fhCuts2(0),
  80   fhSelectedPhotonNLMVsPt(0),fhSelectedPhotonLambda0VsPt(0), fhRandomPhiEta(),
  81 fTreeGJ     (0),
  82 fGamPt      (0),
  83 fGamLambda0 (0),
  84 fGamNLM     (0),
  85 fGamSumPtCh (0),
  86 fGamTime    (0),
  87 fGamNcells  (0),
  88 fGamEta     (0),
  89 fGamPhi     (0),
  90 fGamSumPtNeu(0),
  91 fGamNtracks (0),
  92 fGamNclusters(0),
  93 fGamAvEne   (0),
  94 fJetPhi     (0),
  95 fJetEta     (0),
  96 fJetPt      (0),
  97 fJetBkgChEne(0),
  98 fJetArea    (0),
  99 fJetNtracks (0),
 100 fJetNtracks1(0),
 101 fJetNtracks2(0),
 102 fJetRho(0),
 103 fEventNumber(0),
 104 fNtracks    (0),
 105 fZvertex    (0),
 106 fCentrality (0),
 107 fIso(0),
 108 fGamRho(0)
 109
 110 {
 111   //Default Ctor
 112   //  printf("constructor\n");
 113
 114   //Initialize parameters
 115   InitParameters();
 116   for(Int_t i=0;i<10;i++){
 117     fhJetNjetOverPtCut[i] = 0;
 118     fhPhotonNgammaOverPtCut[i] = 0;
 119   }
 120   fGenerator = new TRandom2();
 121   fGenerator->SetSeed(0);
 122 }
 123
 124 //___________________________________________________________________
 125 AliAnaParticleJetFinderCorrelation::~AliAnaParticleJetFinderCorrelation(){
 126   delete fGenerator;
 127 }
 128
 129
 130 //___________________________________________________________________
 131 TList *  AliAnaParticleJetFinderCorrelation::GetCreateOutputObjects()
 132 {
 133   // Create histograms to be saved in output file and
 134   // store them in fOutputContainer
 135   //  printf("GetCreateOutputObjects\n");
 136
 137   TList * outputContainer = new TList() ;
 138   outputContainer->SetName("ParticleJetFinderHistos") ;
 139
 140   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();
 141   //    Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
 142   //    Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
 143   Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();
 144   //    Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
 145   //    Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
 146   Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
 147   //    Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
 148 //      Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 149
 150 //  fhDeltaPhi  = new TH2F("DeltaPhi","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-6,4);
 151 //  fhDeltaPhi->SetYTitle("#Delta #phi");
 152 //  fhDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 153 //  outputContainer->Add(fhDeltaPhi);
 154
 155   fhDeltaPhiCorrect  = new TH2F("DeltaPhiCorrect","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5);
 156   fhDeltaPhiCorrect->SetYTitle("#Delta #phi");
 157   fhDeltaPhiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
 158   outputContainer->Add(fhDeltaPhiCorrect);
 159
 160   fhDeltaPhi0PiCorrect  = new TH2F("DeltaPhi0PiCorrect","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5);
 161   fhDeltaPhi0PiCorrect->SetYTitle("#Delta #phi");
 162   fhDeltaPhi0PiCorrect->SetXTitle("p_{T trigger} (GeV/c)");
 163   outputContainer->Add(fhDeltaPhi0PiCorrect);
 164
 165
 166   fhDeltaEta  = new TH2F("DeltaEta","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2);
 167   fhDeltaEta->SetYTitle("#Delta #eta");
 168   fhDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 169   outputContainer->Add(fhDeltaEta);
 170
 171   fhDeltaPt  = new TH2F("DeltaPt","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,150,-50,100);
 172   fhDeltaPt->SetYTitle("#Delta p_{T}");
 173   fhDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
 174   outputContainer->Add(fhDeltaPt);
 175
 176   fhPtRatio  = new TH2F("PtRatio","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.);
 177   fhPtRatio->SetYTitle("ratio");
 178   fhPtRatio->SetXTitle("p_{T trigger} (GeV/c)");
 179   outputContainer->Add(fhPtRatio);
 180
 181   fhPt  = new TH2F("Pt","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 182   fhPt->SetYTitle("p_{T jet}(GeV/c)");
 183   fhPt->SetXTitle("p_{T trigger} (GeV/c)");
 184   outputContainer->Add(fhPt);
 185
 186   fhFFz  = new TH2F("FFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,2);
 187   fhFFz->SetYTitle("z");
 188   fhFFz->SetXTitle("p_{T trigger}");
 189   outputContainer->Add(fhFFz) ;
 190
 191   fhFFxi  = new TH2F("FFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,10.);
 192   fhFFxi->SetYTitle("#xi");
 193   fhFFxi->SetXTitle("p_{T trigger}");
 194   outputContainer->Add(fhFFxi) ;
 195
 196   fhFFpt  = new TH2F("FFpt","p_{T i charged} vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,50.);
 197   fhFFpt->SetYTitle("p_{T charged hadron}");
 198   fhFFpt->SetXTitle("p_{T trigger}");
 199   outputContainer->Add(fhFFpt) ;
 200
 201   fhNTracksInCone  = new TH2F("NTracksInCone","Number of tracks in cone vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,150.);
 202   fhNTracksInCone->SetYTitle("p_{T charged hadron}");
 203   fhNTracksInCone->SetXTitle("p_{T trigger}");
 204   outputContainer->Add(fhNTracksInCone) ;
 205
 206   //FF according to jet axis
 207   fhJetFFz  = new TH2F("JetFFz","z = p_{T i charged}/p_{T jet} vs p_{T jet}",nptbins,ptmin,ptmax,200,0.,2);
 208   fhJetFFz->SetYTitle("z");
 209   fhJetFFz->SetXTitle("p_{T jet}");
 210   outputContainer->Add(fhJetFFz) ;
 211
 212   fhJetFFxi  = new TH2F("JetFFxi","#xi = ln(p_{T jet}/p_{T i charged}) vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,10.);
 213   fhJetFFxi->SetYTitle("#xi");
 214   fhJetFFxi->SetXTitle("p_{T jet}");
 215   outputContainer->Add(fhJetFFxi) ;
 216
 217   fhJetFFpt  = new TH2F("JetFFpt","p_{T i charged} vs p_{T jet}", nptbins,ptmin,ptmax,100,0.,50.);
 218   fhJetFFpt->SetYTitle("p_{T charged hadron}");
 219   fhJetFFpt->SetXTitle("p_{T jet}");
 220   outputContainer->Add(fhJetFFpt) ;
 221
 222   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);
 223   fhJetFFzCor->SetYTitle("z");
 224   fhJetFFzCor->SetXTitle("p_{T jet}");
 225   outputContainer->Add(fhJetFFzCor) ;
 226
 227   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.);
 228   fhJetFFxiCor->SetYTitle("#xi");
 229   fhJetFFxiCor->SetXTitle("p_{T jet}");
 230   outputContainer->Add(fhJetFFxiCor) ;
 231
 232
 233   //background FF
 234   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);
 235   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);
 236   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);
 237   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);
 238   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);
 239   for(Int_t i=0;i<5;i++){
 240     fhBkgFFz[i]->SetYTitle("z");
 241     fhBkgFFz[i]->SetXTitle("p_{T trigger}");
 242     outputContainer->Add(fhBkgFFz[i]) ;
 243   }
 244
 245   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.);
 246   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.);
 247   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.);
 248   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.);
 249   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.);
 250   for(Int_t i=0;i<5;i++){
 251     fhBkgFFxi[i]->SetYTitle("#xi");
 252     fhBkgFFxi[i]->SetXTitle("p_{T trigger}");
 253     outputContainer->Add(fhBkgFFxi[i]) ;
 254   }
 255
 256   fhBkgFFpt[0]  = new TH2F("BkgFFptRC",  "p_{T i charged} vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,50.);
 257   fhBkgFFpt[1]  = new TH2F("BkgFFptPCG", "p_{T i charged} vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,50.);
 258   fhBkgFFpt[2]  = new TH2F("BkgFFptPCJ", "p_{T i charged} vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,50.);
 259   fhBkgFFpt[3]  = new TH2F("BkgFFptMP",  "p_{T i charged} vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,50.);
 260   fhBkgFFpt[4]  = new TH2F("BkgFFptTest","p_{T i charged} vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,50.);
 261   for(Int_t i=0;i<5;i++){
 262     fhBkgFFpt[i]->SetYTitle("p_{T charged hadron}");
 263     fhBkgFFpt[i]->SetXTitle("p_{T trigger}");
 264     outputContainer->Add(fhBkgFFpt[i]) ;
 265   }
 266
 267   fhBkgNTracksInCone[0]  = new TH2F("BkgNTracksInConeRC",  "Number of tracks in cone vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,150.);
 268   fhBkgNTracksInCone[1]  = new TH2F("BkgNTracksInConePCG", "Number of tracks in cone vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,150.);
 269   fhBkgNTracksInCone[2]  = new TH2F("BkgNTracksInConePCJ", "Number of tracks in cone vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,150.);
 270   fhBkgNTracksInCone[3]  = new TH2F("BkgNTracksInConeMP",  "Number of tracks in cone vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,150.);
 271   fhBkgNTracksInCone[4]  = new TH2F("BkgNTracksInConeTest","Number of tracks in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,150.);
 272   for(Int_t i=0;i<5;i++){
 273     fhBkgNTracksInCone[i]->SetYTitle("Number of tracks");
 274     fhBkgNTracksInCone[i]->SetXTitle("p_{T trigger}");
 275     outputContainer->Add(fhBkgNTracksInCone[i]) ;
 276   }
 277
 278   fhBkgSumPtInCone[0]  = new TH2F("BkgSumPtInConeRC",  "Sum P_{T} in cone vs p_{T trigger} Bkg RC",   nptbins,ptmin,ptmax,100,0.,100.);
 279   fhBkgSumPtInCone[1]  = new TH2F("BkgSumPtInConePCG", "Sum P_{T} in cone vs p_{T trigger} Bkg PCG",  nptbins,ptmin,ptmax,100,0.,100.);
 280   fhBkgSumPtInCone[2]  = new TH2F("BkgSumPtInConePCJ", "Sum P_{T} in cone vs p_{T trigger} Bkg PCJ",  nptbins,ptmin,ptmax,100,0.,100.);
 281   fhBkgSumPtInCone[3]  = new TH2F("BkgSumPtInConeMP",  "Sum P_{T} in cone vs p_{T trigger} Bkg MP",   nptbins,ptmin,ptmax,100,0.,100.);
 282   fhBkgSumPtInCone[4]  = new TH2F("BkgSumPtInConeTest","Sum P_{T} in cone vs p_{T trigger} Bkg Test", nptbins,ptmin,ptmax,100,0.,100.);
 283   for(Int_t i=0;i<5;i++){
 284     fhBkgSumPtInCone[i]->SetYTitle("Sum P_{T}");
 285     fhBkgSumPtInCone[i]->SetXTitle("p_{T trigger}");
 286     outputContainer->Add(fhBkgSumPtInCone[i]) ;
 287   }
 288
 289   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.);
 290   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.);
 291   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.);
 292   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.);
 293   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.);
 294   for(Int_t i=0;i<5;i++){
 295     fhBkgSumPtnTracksInCone[i]->SetYTitle("Sum p_{T}/Number of tracks");
 296     fhBkgSumPtnTracksInCone[i]->SetXTitle("p_{T trigger}");
 297     outputContainer->Add(fhBkgSumPtnTracksInCone[i]) ;
 298   }
 299
 300
 301   //temporary histograms
 302   fhNjetsNgammas  = new TH2F("NjetsNgammas"," Number of jets vs number of gammas in event",20,0.,100.,10,0.,80.);
 303   fhNjetsNgammas->SetYTitle("Number of gammas");
 304   fhNjetsNgammas->SetXTitle("Number of jets");
 305   outputContainer->Add(fhNjetsNgammas) ;
 306
 307   fhCuts  = new TH1F("Cuts"," Cuts",10,0.,10.);
 308   fhCuts->SetYTitle("Counts");
 309   fhCuts->SetXTitle("Cut number");
 310   outputContainer->Add(fhCuts) ;
 311
 312   fhDeltaPhiBefore  = new TH2F("DeltaPhiBefore","#phi_{trigger} - #phi_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,6.5);
 313   fhDeltaPhiBefore->SetYTitle("#Delta #phi");
 314   fhDeltaPhiBefore->SetXTitle("p_{T trigger} (GeV/c)");
 315   outputContainer->Add(fhDeltaPhiBefore);
 316
 317   fhDeltaEtaBefore  = new TH2F("DeltaEtaBefore","#eta_{trigger} - #eta_{jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-2,2);
 318   fhDeltaEtaBefore->SetYTitle("#Delta #eta");
 319   fhDeltaEtaBefore->SetXTitle("p_{T trigger} (GeV/c)");
 320   outputContainer->Add(fhDeltaEtaBefore);
 321
 322   fhDeltaPtBefore  = new TH2F("DeltaPtBefore","p_{T trigger} - p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,100,-50,50);
 323   fhDeltaPtBefore->SetYTitle("#Delta p_{T}");
 324   fhDeltaPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
 325   outputContainer->Add(fhDeltaPtBefore);
 326
 327   fhPtRatioBefore  = new TH2F("PtRatioBefore","p_{T jet} / p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0,2.);
 328   fhPtRatioBefore->SetYTitle("ratio");
 329   fhPtRatioBefore->SetXTitle("p_{T trigger} (GeV/c)");
 330   outputContainer->Add(fhPtRatioBefore);
 331
 332   fhPtBefore  = new TH2F("PtBefore","p_{T jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 333   fhPtBefore->SetYTitle("p_{T jet}(GeV/c)");
 334   fhPtBefore->SetXTitle("p_{T trigger} (GeV/c)");
 335   outputContainer->Add(fhPtBefore);
 336
 337   fhDeltaPhi0PiCorrectBefore  = new TH2F("DeltaPhi0PiCorrectBefore","#phi_{trigger} - #phi_{jet} (0,#pi) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0,3.5);
 338   fhDeltaPhi0PiCorrectBefore->SetYTitle("#Delta #phi");
 339   fhDeltaPhi0PiCorrectBefore->SetXTitle("p_{T trigger} (GeV/c)");
 340   outputContainer->Add(fhDeltaPhi0PiCorrectBefore);
 341
 342   //temporary jet histograms
 343   fhJetPtBefore            = new TH1F("JetPtBefore","JetPtBefore",150,-50,100);
 344   fhJetPtBefore->SetYTitle("Counts");
 345   fhJetPtBefore->SetXTitle("p_{T jet}(GeV/c)");
 346   outputContainer->Add(fhJetPtBefore) ;
 347
 348   fhJetPt            = new TH1F("JetPt","JetPt",150,-50,100);
 349   fhJetPt->SetYTitle("Counts");
 350   fhJetPt->SetXTitle("p_{T jet}(GeV/c)");
 351   outputContainer->Add(fhJetPt) ;
 352
 353   fhJetPtMostEne            = new TH1F("JetPtMostEne","JetPtMostEne",150,0,150);
 354   fhJetPtMostEne->SetYTitle("Counts");
 355   fhJetPtMostEne->SetXTitle("p_{T jet}(GeV/c)");
 356   outputContainer->Add(fhJetPtMostEne) ;
 357
 358   fhJetPhi           = new TH1F("JetPhi","JetPhi",130,0,6.5);
 359   fhJetPhi->SetYTitle("Counts");
 360   fhJetPhi->SetXTitle("#phi_{jet}");
 361   outputContainer->Add(fhJetPhi) ;
 362
 363   fhJetEta           = new TH1F("JetEta","JetEta",100,-1,1);
 364   fhJetEta->SetYTitle("Counts");
 365   fhJetEta->SetXTitle("#eta_{jet}");
 366   outputContainer->Add(fhJetEta) ;
 367
 368   fhJetEtaVsPt      = new TH2F("JetEtaVsPt","JetEtaVsPt",100,0,100,50,-1,1);
 369   fhJetEtaVsPt->SetYTitle("#eta_{jet}");
 370   fhJetEtaVsPt->SetXTitle("p_{T,jet}(GeV/c)");
 371   outputContainer->Add(fhJetEtaVsPt) ;
 372
 373   fhJetPhiVsEta      = new TH2F("JetPhiVsEta","JetPhiVsEta",65,0,6.5,50,-1,1);
 374   fhJetPhiVsEta->SetYTitle("#eta_{jet}");
 375   fhJetPhiVsEta->SetXTitle("#phi_{jet}");
 376   outputContainer->Add(fhJetPhiVsEta) ;
 377
 378   fhJetEtaVsNpartInJet= new TH2F("JetEtaVsNpartInJet","JetEtaVsNpartInJet",50,-1,1,100,0.,200.);
 379   fhJetEtaVsNpartInJet->SetYTitle("N_{tracks-in-jet}");
 380   fhJetEtaVsNpartInJet->SetXTitle("#eta_{jet}");
 381   outputContainer->Add(fhJetEtaVsNpartInJet) ;
 382
 383   fhJetEtaVsNpartInJetBkg= new TH2F("JetEtaVsNpartInJetBkg","JetEtaVsNpartInJetBkg",50,-1,1,100,0.,200.);
 384   fhJetEtaVsNpartInJetBkg->SetYTitle("N_{tracks-in-jet}");
 385   fhJetEtaVsNpartInJetBkg->SetXTitle("#eta_{jet}");
 386   outputContainer->Add(fhJetEtaVsNpartInJetBkg) ;
 387
 388   fhJetChBkgEnergyVsPt = new TH2F("JetBkgChEnergyVsPt","JetBkgChEnergyVsPt",100,0,100,90,0,90);
 389   fhJetChBkgEnergyVsPt->SetYTitle("Jet Bkg Energy (GeV)");
 390   fhJetChBkgEnergyVsPt->SetXTitle("p_{T,jet} (GeV/c)");
 391   outputContainer->Add(fhJetChBkgEnergyVsPt);
 392
 393   fhJetChAreaVsPt      = new TH2F("JetChAreaVsPt","JetChAreaVsPt",100,0,100,50,0,1);
 394   fhJetChAreaVsPt->SetYTitle("Jet Area");
 395   fhJetChAreaVsPt->SetXTitle("p_{T,jet} (GeV/c)");
 396   outputContainer->Add(fhJetChAreaVsPt);
 397
 398   if(IsHistogramTracks()){
 399     fhTrackPhiVsEta      = new TH2F("TrackPhiVsEta","TrackPhiVsEta",65,0,6.5,50,-1,1);
 400     fhTrackPhiVsEta->SetYTitle("#eta_{track}");
 401     fhTrackPhiVsEta->SetXTitle("#phi_{track}");
 402     outputContainer->Add(fhTrackPhiVsEta) ;
 403
 404     fhTrackAveTrackPt      = new TH1F("TrackAveTrackPt","TrackAveTrackPt",45,0,1.5);
 405     fhTrackAveTrackPt->SetYTitle("Counts");
 406     fhTrackAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
 407     outputContainer->Add(fhTrackAveTrackPt);
 408
 409   }//end of IsHistogramTracks()
 410
 411   for(Int_t i=0;i<10;i++){
 412     fhJetNjetOverPtCut[i]      = new TH1F(Form("JetNjetOverPtCut%d", i),Form("JetNjetOverPtCut%d", i),100,0,100);
 413     fhJetNjetOverPtCut[i]->SetYTitle("Counts");
 414     fhJetNjetOverPtCut[i]->SetXTitle("N_{jets} over threshold");
 415     outputContainer->Add(fhJetNjetOverPtCut[i]);
 416   }
 417
 418   fhJetChBkgEnergyVsArea = new TH2F("JetBkgChEnergyVsArea","JetBkgChEnergyVsArea",100,0,100,70,0,0.7);
 419   fhJetChBkgEnergyVsArea->SetXTitle("Jet Bkg Energy (GeV)");
 420   fhJetChBkgEnergyVsArea->SetYTitle("Area");
 421   outputContainer->Add(fhJetChBkgEnergyVsArea);
 422
 423   fhJetRhoVsPt           = new TH2F("JetRhoVsPt","JetRhoVsPt",100,0,100,100,0,150);
 424   fhJetRhoVsPt->SetYTitle("Rho");
 425   fhJetRhoVsPt->SetXTitle("p_{T,jet} (GeV/c)");
 426   outputContainer->Add(fhJetRhoVsPt);
 427
 428   if(IsHistogramJetBkg()){
 429     fhJetRhoVsCentrality           = new TH2F("JetRhoVsCentrality","JetRhoVsCentrality",100,0,100,100,0,200);
 430     fhJetRhoVsCentrality->SetYTitle("Rho");
 431     fhJetRhoVsCentrality->SetXTitle("Centrality");
 432     outputContainer->Add(fhJetRhoVsCentrality);
 433   }
 434
 435   fhJetNparticlesInJet           = new TH1F("JetNparticlesInJet","JetNparticlesInJet",100,0,200);
 436   fhJetNparticlesInJet->SetXTitle("N^{particles}");
 437   fhJetNparticlesInJet->SetYTitle("N^{jets}");
 438   outputContainer->Add(fhJetNparticlesInJet);
 439
 440   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);
 441   fhJetDeltaEtaDeltaPhi->SetXTitle("#Delta #eta^{jet-track}");
 442   fhJetDeltaEtaDeltaPhi->SetYTitle("#Delta #phi^{jet-track}");
 443   outputContainer->Add(fhJetDeltaEtaDeltaPhi );
 444
 445
 446   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);
 447   fhJetDeltaEtaDeltaPhiAllTracks->SetXTitle("#Delta #eta^{jet-track}");
 448   fhJetDeltaEtaDeltaPhiAllTracks->SetYTitle("#Delta #phi^{jet-track}");
 449   outputContainer->Add(fhJetDeltaEtaDeltaPhiAllTracks);
 450
 451
 452   if(IsHistogramJetTracks()){
 453     fhJetAveTrackPt           = new TH1F("JetAveTrackPt","JetAveTrackPt",45,0.,1.5);
 454     fhJetAveTrackPt->SetXTitle("Average p_{T,track} (GeV/c)");
 455     fhJetAveTrackPt->SetYTitle("Counts");
 456     outputContainer->Add(fhJetAveTrackPt);
 457
 458     for(Int_t i=0;i<6;i++){
 459       if(i==0) fhJetNtracksInJetAboveThr[i]      = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,200);
 460       else fhJetNtracksInJetAboveThr[i]      = new TH2F(Form("JetNtracksInJetAboveThr%d", i),Form("JetNtracksInJetAboveThr%d", i),100,0,100,100,0,100);
 461       fhJetNtracksInJetAboveThr[i]->SetXTitle("p_{T,jet} (GeV/c)");
 462       fhJetNtracksInJetAboveThr[i]->SetYTitle("N_{tracks} over threshold");
 463       outputContainer->Add(fhJetNtracksInJetAboveThr[i]);
 464     }
 465
 466     for(Int_t i=0;i<5;i++){
 467       fhJetRatioNTrkAboveToNTrk[i]      = new TH2F(Form("JetRatioNTrkAboveToNTrk%d", i),Form("JetRatioNTrkAboveToNTrk%d", i),100,0,100,40,0,1);
 468       fhJetRatioNTrkAboveToNTrk[i]->SetXTitle("p_{T,jet} (GeV/c)");
 469       fhJetRatioNTrkAboveToNTrk[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
 470       outputContainer->Add(fhJetRatioNTrkAboveToNTrk[i]);
 471
 472       fhJetNtrackRatioMostEne[i]      = new TH2F(Form("JetNtrackRatioMostEne%d", i),Form("JetNtrackRatioMostEne%d", i),100,0,100,40,0,1);
 473       fhJetNtrackRatioMostEne[i]->SetXTitle("p_{T,jet} (GeV/c)");
 474       fhJetNtrackRatioMostEne[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
 475       outputContainer->Add(fhJetNtrackRatioMostEne[i]);
 476
 477       fhJetNtrackRatioJet5GeV[i]      = new TH2F(Form("JetNtrackRatioJet5GeV%d", i),Form("JetNtrackRatioJet5GeV%d", i),100,0,100,40,0,1);
 478       fhJetNtrackRatioJet5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
 479       fhJetNtrackRatioJet5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
 480       outputContainer->Add(fhJetNtrackRatioJet5GeV[i]);
 481
 482       fhJetNtrackRatioLead5GeV[i]      = new TH2F(Form("JetNtrackRatioLead5GeV%d", i),Form("JetNtrackRatioLead5GeV%d", i),100,0,100,40,0,1);
 483       fhJetNtrackRatioLead5GeV[i]->SetXTitle("p_{T,jet} (GeV/c)");
 484       fhJetNtrackRatioLead5GeV[i]->SetYTitle("Ratio N_{tracks} over threshold to N_{tracks}");
 485       outputContainer->Add(fhJetNtrackRatioLead5GeV[i]);
 486     }
 487   }//end of if IsHistogramJetTracks
 488
 489   //temporary background jets histograms
 490   if(IsHistogramJetBkg()){
 491     for(Int_t i=0;i<4;i++){
 492       fhBkgJetBackground[i]      = new TH1F(Form("BkgJetBackground%d", i),Form("BkgJetBackground%d", i),100,0,200);
 493       fhBkgJetBackground[i]->SetXTitle("<#rho> (GeV/c)");
 494       fhBkgJetBackground[i]->SetYTitle("Counts");
 495       outputContainer->Add(fhBkgJetBackground[i]);
 496
 497       fhBkgJetSigma[i]      = new TH1F(Form("BkgJetSigma%d", i),Form("BkgJetSigma%d", i),100,0,50);
 498       fhBkgJetSigma[i]->SetXTitle("#sigma (GeV/c)");
 499       fhBkgJetSigma[i]->SetYTitle("Counts");
 500       outputContainer->Add(fhBkgJetSigma[i]);
 501
 502       fhBkgJetArea[i]      = new TH1F(Form("BkgJetArea%d", i),Form("BkgJetArea%d", i),100,0,1);
 503       fhBkgJetArea[i]->SetXTitle("<A>");
 504       fhBkgJetArea[i]->SetYTitle("Counts");
 505       outputContainer->Add(fhBkgJetArea[i]);
 506     }
 507   }
 508
 509   //temporary photon histograms
 510   fhPhotonPtMostEne = new TH1F("PhotonPtMostEne","PhotonPtMostEne",100,0,100);
 511   fhPhotonPtMostEne->SetYTitle("Counts");
 512   fhPhotonPtMostEne->SetXTitle("p_{T,#gamma} (GeV/c)");
 513   outputContainer->Add(fhPhotonPtMostEne);
 514
 515 //  fhPhotonIndexMostEne = new TH1F("PhotonIndexMostEne","PhotonIndexMostEne",100,0,100);
 516 //  fhPhotonIndexMostEne->SetYTitle("Counts");
 517 //  fhPhotonIndexMostEne->SetXTitle("Index");
 518 //  outputContainer->Add(fhPhotonIndexMostEne);
 519
 520   fhPhotonAverageEnergy = new TH1F("PhotonAverageEnergy","PhotonAverageEnergy",100,0,10);
 521   fhPhotonAverageEnergy->SetYTitle("Counts");
 522   fhPhotonAverageEnergy->SetXTitle("p_{T,#gamma} (GeV/c)");
 523   outputContainer->Add(fhPhotonAverageEnergy);
 524
 525   fhPhotonRatioAveEneToMostEne = new TH1F("PhotonRatioAveEneToMostEne","PhotonRatioAveEneToMostEne",100,0,1);
 526   fhPhotonRatioAveEneToMostEne->SetYTitle("Counts");
 527   fhPhotonRatioAveEneToMostEne->SetXTitle("Ratio");
 528   outputContainer->Add(fhPhotonRatioAveEneToMostEne);
 529
 530   fhPhotonAverageEnergyMinus1 = new TH1F("PhotonAverageEnergyMinus1","PhotonAverageEnergyMinus1",100,0,10);
 531   fhPhotonAverageEnergyMinus1->SetYTitle("Counts");
 532   fhPhotonAverageEnergyMinus1->SetXTitle("p_{T,#gamma} (GeV/c)");
 533   outputContainer->Add(fhPhotonAverageEnergyMinus1);
 534
 535   fhPhotonRatioAveEneMinus1ToMostEne = new TH1F("PhotonRatioAveEneMinus1ToMostEne","PhotonRatioAveEneMinus1ToMostEne",100,0,1);
 536   fhPhotonRatioAveEneMinus1ToMostEne->SetYTitle("Counts");
 537   fhPhotonRatioAveEneMinus1ToMostEne->SetXTitle("Ratio");
 538   outputContainer->Add(fhPhotonRatioAveEneMinus1ToMostEne);
 539
 540   fhPhotonNgammaMoreAverageToNgamma = new TH1F("PhotonNgammaMoreAverageToNgamma","PhotonNgammaMoreAverageToNgamma",100,0,1);
 541   fhPhotonNgammaMoreAverageToNgamma->SetYTitle("Counts");
 542   fhPhotonNgammaMoreAverageToNgamma->SetXTitle("Ratio");
 543   outputContainer->Add(fhPhotonNgammaMoreAverageToNgamma);
 544
 545   fhPhotonNgammaMoreAverageMinus1ToNgamma = new TH1F("PhotonNgammaMoreAverageMinus1ToNgamma","PhotonNgammaMoreAverageMinus1ToNgamma",100,0,1);
 546   fhPhotonNgammaMoreAverageMinus1ToNgamma->SetYTitle("Counts");
 547   fhPhotonNgammaMoreAverageMinus1ToNgamma->SetXTitle("Ratio");
 548   outputContainer->Add(fhPhotonNgammaMoreAverageMinus1ToNgamma);
 549
 550   for(Int_t i=0;i<10;i++){
 551     fhPhotonNgammaOverPtCut[i]      = new TH1F(Form("PhotonNgammaOverPtCut%d",i),Form("PhotonNgammaOverPtCut%d",i),100,0,100);
 552     fhPhotonNgammaOverPtCut[i]->SetYTitle("Counts");
 553     fhPhotonNgammaOverPtCut[i]->SetXTitle("N_{#gamma} over threshold");
 554     outputContainer->Add(fhPhotonNgammaOverPtCut[i]);
 555   }
 556
 557   fhPhotonBkgRhoVsNtracks = new TH2F("PhotonBkgRhoVsNtracks","PhotonBkgRhoVsNtracks",200,0,2500,75,0,1.5);
 558   //fhPhotonBkgRhoVsNtracks->SetXTitle("Counts");
 559   fhPhotonBkgRhoVsNtracks->SetXTitle("Ntracks");
 560   fhPhotonBkgRhoVsNtracks->SetYTitle("Rho");
 561   outputContainer->Add(fhPhotonBkgRhoVsNtracks);
 562
 563   fhPhotonBkgRhoVsNclusters = new TH2F("PhotonBkgRhoVsNclusters","PhotonBkgRhoVsNclusters",50,0,100,75,0,1.5);
 564   fhPhotonBkgRhoVsNclusters->SetXTitle("Nclusters");
 565   fhPhotonBkgRhoVsNclusters->SetYTitle("Rho");
 566   outputContainer->Add(fhPhotonBkgRhoVsNclusters);
 567
 568   fhPhotonBkgRhoVsCentrality = new TH2F("PhotonBkgRhoVsCentrality","PhotonBkgRhoVsCentrality",100,0,100,75,0,1.5);
 569   fhPhotonBkgRhoVsCentrality->SetXTitle("Centrality");
 570   fhPhotonBkgRhoVsCentrality->SetYTitle("Rho");
 571   outputContainer->Add(fhPhotonBkgRhoVsCentrality);
 572
 573   fhPhotonBkgRhoVsNcells = new TH2F("PhotonBkgRhoVsNcells","PhotonBkgRhoVsNcells",100,0,200,75,0,1.5);
 574   fhPhotonBkgRhoVsNcells->SetXTitle("N_{cells}");
 575   fhPhotonBkgRhoVsNcells->SetYTitle("Rho");
 576   outputContainer->Add(fhPhotonBkgRhoVsNcells);
 577
 578   fhPhotonPt = new TH1F("PhotonPt","PhotonPt",220,-10,100);
 579   fhPhotonPt->SetXTitle("p_{T,#gamma} (GeV/c)");
 580   fhPhotonPt->SetYTitle("Counts");
 581   outputContainer->Add(fhPhotonPt);
 582
 583   fhPhotonPtCorrected = new TH1F("PhotonPtCorrected","PhotonPtCorrected",220,-10,100);
 584   fhPhotonPtCorrected->SetXTitle("p_{T,#gamma} (GeV/c)");
 585   fhPhotonPtCorrected->SetYTitle("Counts");
 586   outputContainer->Add(fhPhotonPtCorrected);
 587
 588   fhPhotonPtDiff = new TH1F("PhotonPtDiff","PhotonPtDiff",50,0,10);
 589   fhPhotonPtDiff->SetXTitle("p_{T,#gamma} (GeV/c)");
 590   fhPhotonPtDiff->SetYTitle("Counts");
 591   outputContainer->Add(fhPhotonPtDiff);
 592
 593   fhPhotonPtDiffVsNtracks = new TH2F("PhotonPtDiffVsNtracks","PhotonPtDiffVsNtracks",200,0,2500,50,0,5);
 594   fhPhotonPtDiffVsNtracks->SetXTitle("N_{tracks}");
 595   fhPhotonPtDiffVsNtracks->SetYTitle("<#rho^{#gamma}>*N_{cells}");
 596   outputContainer->Add(fhPhotonPtDiffVsNtracks);
 597
 598   fhPhotonPtDiffVsNclusters = new TH2F("PhotonPtDiffVsNclusters","PhotonPtDiffVsNclusters",50,0,100,50,0,5);
 599   fhPhotonPtDiffVsNclusters->SetXTitle("N_{clusters}");
 600   fhPhotonPtDiffVsNclusters->SetYTitle("<#rho^{#gamma}>*N_{cells}");
 601   outputContainer->Add(fhPhotonPtDiffVsNclusters);
 602
 603   fhPhotonPtDiffVsCentrality = new TH2F("PhotonPtDiffVsCentrality","PhotonPtDiffVsCentrality",100,0,100,50,0,5);
 604   fhPhotonPtDiffVsCentrality->SetXTitle("Centrality");
 605   fhPhotonPtDiffVsCentrality->SetYTitle("<#rho^{#gamma}>*N_{cells}");
 606   outputContainer->Add(fhPhotonPtDiffVsCentrality);
 607
 608   fhPhotonPtDiffVsNcells = new TH2F("PhotonPtDiffVsNcells","PhotonPtDiffVsNcells",100,0,200,50,0,5);
 609   fhPhotonPtDiffVsNcells->SetXTitle("N_{cells}");
 610   fhPhotonPtDiffVsNcells->SetYTitle("<#rho^{#gamma}>*N_{cells}");
 611   outputContainer->Add(fhPhotonPtDiffVsNcells);
 612
 613
 614   fhPhotonPtCorrectedZoom = new TH1F("PhotonPtCorrectedZoom","PhotonPtCorrectedZoom",200,-5,5);
 615   fhPhotonPtCorrectedZoom->SetXTitle("p_{T,#gamma} (GeV/c)");
 616   fhPhotonPtCorrectedZoom->SetYTitle("Counts");
 617   outputContainer->Add(fhPhotonPtCorrectedZoom);
 618
 619   fhPhotonSumPtInCone = new TH1F("PhotonSumPtInCone","PhotonSumPtInCone",100,0,100);
 620   fhPhotonSumPtInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
 621   fhPhotonSumPtInCone->SetYTitle("Counts");
 622   outputContainer->Add(fhPhotonSumPtInCone);
 623
 624   fhPhotonSumPtCorrectInCone = new TH1F("PhotonSumPtCorrectInCone","PhotonSumPtCorrectInCone",100,-20,80);
 625   fhPhotonSumPtCorrectInCone->SetXTitle("#Sigma p_{T,#gamma} (GeV/c)");
 626   fhPhotonSumPtCorrectInCone->SetYTitle("Counts");
 627   outputContainer->Add(fhPhotonSumPtCorrectInCone);
 628
 629   fhPhotonSumPtChargedInCone = new TH1F("PhotonSumPtChargedInCone","PhotonSumPtChargedInCone",100,0,100);
 630   fhPhotonSumPtChargedInCone->SetXTitle("#Sigma p_{T,#gamma}^{ch} (GeV/c)");
 631   fhPhotonSumPtChargedInCone->SetYTitle("Counts");
 632   outputContainer->Add(fhPhotonSumPtChargedInCone);
 633
 634
 635   //temporary jet histograms after selection
 636   fhSelectedJetPhiVsEta      = new TH2F("SelectedJetSelectedPhiVsEta","SelectedJetPhiVsEta",65,0,6.5,50,-1,1);
 637   fhSelectedJetPhiVsEta->SetYTitle("#eta_{jet}");
 638   fhSelectedJetPhiVsEta->SetXTitle("#phi_{jet}");
 639   outputContainer->Add(fhSelectedJetPhiVsEta) ;
 640
 641   fhSelectedJetChBkgEnergyVsPtJet = new TH2F("SelectedJetBkgChEnergyVsPtJet","SelectedJetBkgChEnergyVsPtJet",100,0,100,90,0,90);
 642   fhSelectedJetChBkgEnergyVsPtJet->SetYTitle("Jet Bkg Energy (GeV)");
 643   fhSelectedJetChBkgEnergyVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
 644   outputContainer->Add(fhSelectedJetChBkgEnergyVsPtJet);
 645
 646   fhSelectedJetChAreaVsPtJet      = new TH2F("SelectedJetChAreaVsPtJet","SelectedJetChAreaVsPtJet",100,0,100,50,0,1);
 647   fhSelectedJetChAreaVsPtJet->SetYTitle("Jet Area");
 648   fhSelectedJetChAreaVsPtJet->SetXTitle("p_{T,jet} (GeV/c)");
 649   outputContainer->Add(fhSelectedJetChAreaVsPtJet);
 650
 651   fhSelectedJetNjet      = new TH1F("SelectedJetNjet","SelectedJetNjet",100,0,100);
 652   fhSelectedJetNjet->SetYTitle("Counts");
 653   fhSelectedJetNjet->SetXTitle("N_{jets} per event");
 654   outputContainer->Add(fhSelectedJetNjet);
 655
 656   fhSelectedNtracks      = new TH1F("SelectedNtracks","SelectedNtracks",100,0,2000);
 657   fhSelectedNtracks->SetYTitle("Counts");
 658   fhSelectedNtracks->SetXTitle("N_{tracks} per event");
 659   outputContainer->Add(fhSelectedNtracks);
 660
 661   fhSelectedTrackPhiVsEta      = new TH2F("SelectedTrackPhiVsEta","SelectedTrackPhiVsEta",65,0,6.5,50,-1,1);
 662   fhSelectedTrackPhiVsEta->SetYTitle("#eta_{track}");
 663   fhSelectedTrackPhiVsEta->SetXTitle("#phi_{track}");
 664   outputContainer->Add(fhSelectedTrackPhiVsEta) ;
 665
 666   fhCuts2      = new TH1F("Cuts2","Cuts2",10,0.,10.);
 667   fhCuts2->SetYTitle("Counts");
 668   fhCuts2->SetXTitle("Cut number");
 669   outputContainer->Add(fhCuts2);
 670
 671   fhSelectedPhotonNLMVsPt      = new TH2F("SelectedPhotonNLMVsPt","SelectedPhotonNLMVsPt",100,0,100,10,0,10);
 672   fhSelectedPhotonNLMVsPt->SetYTitle("NLM");
 673   fhSelectedPhotonNLMVsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
 674   outputContainer->Add(fhSelectedPhotonNLMVsPt);
 675
 676   fhSelectedPhotonLambda0VsPt      = new TH2F("SelectedPhotonLambda0VsPt","SelectedPhotonLambda0VsPt",100,0,100,50,0,5);
 677   fhSelectedPhotonLambda0VsPt->SetYTitle("#lambda_{0}");
 678   fhSelectedPhotonLambda0VsPt->SetXTitle("p_{T,#gamma} (GeV/c)");
 679   outputContainer->Add(fhSelectedPhotonLambda0VsPt);
 680
 681   //random
 682   fhRandomPhiEta[0]  = new TH2F("RandomPhiEtaRC",  "RandomPhiEtaRC",            100,0,6.5,100,-1.,1.);
 683   fhRandomPhiEta[1]  = new TH2F("RandomPhiEtaPCG", "RandomPhiEtaPerpConePhoton",100,0,6.5,100,-1.,1.);
 684   fhRandomPhiEta[2]  = new TH2F("RandomPhiEtaPCJ", "RandomPhiEtaPerpConeJet",   100,0,6.5,100,-1.,1.);
 685   fhRandomPhiEta[3]  = new TH2F("RandomPhiEtaMP",  "RandomPhiEtaMidPoint",      100,0,6.5,100,-1.,1.);
 686   fhRandomPhiEta[4]  = new TH2F("RandomPhiEtaTest","RandomPhiEtaTest",          100,0,6.5,100,-1.,1.);
 687   for(Int_t i=0;i<5;i++){
 688     fhRandomPhiEta[i]->SetXTitle("#phi");
 689     fhRandomPhiEta[i]->SetYTitle("#eta");
 690     outputContainer->Add(fhRandomPhiEta[i]);
 691   }
 692
 693
 694   //tree with data
 695   if(fSaveGJTree){
 696     fTreeGJ=new TTree("fTreeGJ","fTreeGJ");
 697     fTreeGJ->Branch("fGamPt"      ,&fGamPt    ,"fGamPt/D");//! fGamPt
 698     fTreeGJ->Branch("fGamLambda0" ,&fGamLambda0  ,"fGamLambda0/D");
 699     fTreeGJ->Branch("fGamNLM"     ,&fGamNLM      ,"fGamNLM/I");
 700     fTreeGJ->Branch("fGamSumPtCh" ,&fGamSumPtCh  ,"fGamSumPtCh/D");
 701     fTreeGJ->Branch("fGamNtracks" ,&fGamNtracks  ,"fGamNtracks/I");
 702     fTreeGJ->Branch("fGamTime"    ,&fGamTime     ,"fGamTime/D");
 703     fTreeGJ->Branch("fGamNcells"  ,&fGamNcells   ,"fGamNcells/I");
 704     fTreeGJ->Branch("fGamEta"     ,&fGamEta      ,"fGamEta/D");
 705     fTreeGJ->Branch("fGamPhi"     ,&fGamPhi      ,"fGamPhi/D");
 706     fTreeGJ->Branch("fGamSumPtNeu",&fGamSumPtNeu ,"fGamSumPtNeu/D");
 707     fTreeGJ->Branch("fGamNclusters" ,&fGamNclusters  ,"fGamNclusters/I");
 708     fTreeGJ->Branch("fGamAvEne"   ,&fGamAvEne    ,"fGamAvEne/D");
 709     fTreeGJ->Branch("fJetPhi"     ,&fJetPhi      ,"fJetPhi/D");
 710     fTreeGJ->Branch("fJetEta"     ,&fJetEta      ,"fJetEta/D");
 711     fTreeGJ->Branch("fJetPt"      ,&fJetPt       ,"fJetPt/D");
 712     fTreeGJ->Branch("fJetBkgChEne",&fJetBkgChEne ,"fJetBkgChEne/D");
 713     fTreeGJ->Branch("fJetArea"    ,&fJetArea     ,"fJetArea/D");
 714     fTreeGJ->Branch("fJetNtracks" ,&fJetNtracks  ,"fJetNtracks/I");
 715     fTreeGJ->Branch("fJetNtracks1" ,&fJetNtracks1  ,"fJetNtracks1/I");
 716     fTreeGJ->Branch("fJetNtracks2" ,&fJetNtracks2  ,"fJetNtracks2/I");
 717     fTreeGJ->Branch("fJetRho" ,&fJetRho  ,"fJetRho/D");
 718     fTreeGJ->Branch("fEventNumber",&fEventNumber ,"fEventNumber/I");
 719     fTreeGJ->Branch("fNtracks"    ,&fNtracks     ,"fNtracks/I");
 720     fTreeGJ->Branch("fZvertex"    ,&fZvertex     ,"fZvertex/D");
 721     fTreeGJ->Branch("fCentrality" ,&fCentrality  ,"fCentrality/D");
 722     fTreeGJ->Branch("fIso" ,&fIso  ,"fIso/O");
 723     fTreeGJ->Branch("fGamRho" ,&fGamRho  ,"fGamRho/D");
 724
 725     outputContainer->Add(fTreeGJ);
 726   }
 727
 728   return outputContainer;
 729
 730 }
 731
 732 //_______________________________________________________
 733 void AliAnaParticleJetFinderCorrelation::InitParameters()
 734 {
 735   //  printf("InitParameters\n");
 736   //Initialize the parameters of the analysis.
 737   SetInputAODName("PWG4Particle");
 738   AddToHistogramsName("AnaJetFinderCorr_");
 739
 740   fDeltaPhiMinCut    = 1.5 ;
 741   fDeltaPhiMaxCut    = 4.5 ;
 742   fRatioMaxCut       = 1.2 ;
 743   fRatioMinCut       = 0.3 ;
 744   fConeSize          = 0.3 ;
 745   fPtThresholdInCone = 0.5 ;
 746   fUseJetRefTracks   = kFALSE ;
 747   fMakeCorrelationInHistoMaker = kFALSE ;
 748   fSelectIsolated = kFALSE;
 749   fJetConeSize = 0.4 ;
 750   fJetMinPt = 5 ; //GeV/c
 751   fJetAreaFraction = 0.6 ;
 752   fJetBranchName = "jets";
 753   fBkgJetBranchName = "jets";
 754   fGammaConeSize = 0.3;
 755   fUseBackgroundSubtractionGamma = kFALSE;
 756   fSaveGJTree = kTRUE;
 757   fMostEnergetic = kFALSE;
 758   fMostOpposite = kTRUE;
 759   fUseHistogramJetBkg = kTRUE;
 760   fUseHistogramTracks = kTRUE;
 761   fUseHistogramJetTracks = kTRUE;
 762
 763 }
 764
 765 //__________________________________________________________________
 766 Int_t  AliAnaParticleJetFinderCorrelation::SelectJet(AliAODPWG4Particle * particle, TClonesArray *aodRecJets)
 767 {
 768   //Input for jets is TClonesArray *aodRecJets
 769   //Returns the index of the jet that is opposite to the particle
 770   //  printf(" SelectJet ");
 771
 772   Double_t particlePt=particle->Pt();
 773   if(fUseBackgroundSubtractionGamma) {
 774     Int_t clusterIDtmp = particle->GetCaloLabel(0) ;
 775     Int_t nCells=0;
 776     AliVCluster *cluster=0;
 777     if(!(clusterIDtmp<0) ){
 778       Int_t iclustmp = -1;
 779       TObjArray* clusters = GetEMCALClusters();
 780       cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
 781       nCells = cluster->GetNCells();
 782     }
 783     particlePt-=(fGamRho*nCells);
 784   }
 785   if(particlePt<=0) {
 786     //printf("Particle with negative  or 0 pt\n");
 787     return -2;
 788   }
 789
 790   Int_t njets = aodRecJets->GetEntriesFast();
 791   AliAODJet * jet = 0 ;
 792   Int_t index = -1;
 793   Double_t dphiprev= 10000;
 794   Double_t particlePhi=particle->Phi();
 795   Double_t deltaPhi=-10000.;// in the range (0; 2*pi)
 796   Double_t jetPt=0.;
 797
 798   for(Int_t ijet = 0; ijet < njets ; ijet++){
 799     jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
 800     if(!jet)
 801     {
 802       AliInfo("Jet not in container");
 803       continue;
 804     }
 805     fhCuts2->Fill(2.,1.);
 806     jetPt=jet->Pt();
 807     if(fBackgroundJetFromReader ){
 808       jetPt-= (fJetRho * jet->EffectiveAreaCharged() );
 809     }
 810     if(jetPt<0.) continue;
 811     //put jet eta requirement here |eta_jet|<0.9-jet_cone_size
 812     fhCuts2->Fill(3.,1.);
 813     if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
 814     fhCuts2->Fill(4.,1.);
 815     //    if(jet->Pt()<5) continue;
 816     if(jetPt<fJetMinPt) continue;
 817     fhCuts2->Fill(5.,1.);
 818     if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
 819     fhCuts2->Fill(6.,1.);
 820     //printf("jet found\n");
 821     Double_t deltaPhi0pi  = TMath::Abs(particle->Phi()-jet->Phi());
 822     Double_t ratio = jetPt/particlePt;
 823
 824     deltaPhi = particlePhi - jet->Phi() ;
 825     if ( deltaPhi0pi > TMath::Pi() ) deltaPhi0pi = 2. * TMath::Pi() - deltaPhi0pi ;
 826     if(deltaPhi<0) deltaPhi +=(TMath::Pi()*2.);
 827
 828     fhDeltaPtBefore ->Fill(particlePt, particlePt - jetPt);
 829     fhDeltaPhiBefore->Fill(particlePt, deltaPhi);
 830     fhDeltaEtaBefore->Fill(particlePt, particle->Eta() - jet->Eta());
 831     fhPtRatioBefore ->Fill(particlePt, jetPt/particlePt);
 832     fhPtBefore      ->Fill(particlePt, jetPt);
 833
 834     fhDeltaPhi0PiCorrectBefore->Fill(particlePt, deltaPhi0pi);//good
 835
 836     if(GetDebug() > 5)
 837       printf("AliAnaParticleJetFinderCorrelation::SelectJet() - Jet %d, Ratio pT %2.3f, Delta phi %2.3f\n",ijet,ratio,deltaPhi);
 838
 839     if((deltaPhi > fDeltaPhiMinCut) && (deltaPhi < fDeltaPhiMaxCut) &&
 840        (ratio > fRatioMinCut) && (ratio < fRatioMaxCut)  &&
 841        (TMath::Abs(deltaPhi-TMath::Pi()) < TMath::Abs(dphiprev-TMath::Pi())  )){//In case there is more than one jet select the most opposite.
 842       dphiprev = deltaPhi;
 843       index = ijet ;
 844     }//Selection cuts
 845   }//AOD Jet loop
 846
 847   return index ;
 848
 849 }
 850
 851 //__________________________________________________________________
 852 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD()
 853 {
 854   //Particle-Jet Correlation Analysis, fill AODs
 855   //  printf("I use MakeAnalysisFillAOD\n");
 856   //Get the event, check if there are AODs available, if not, skip this analysis
 857   AliAODEvent * event = NULL;
 858
 859   //  cout<<"GetReader()->GetOutputEvent() "<<GetReader()->GetOutputEvent()<<endl;
 860   //  cout<<"GetReader()->GetDataType() "<<GetReader()->GetDataType() <<endl;
 861   //  cout<<"AliCaloTrackReader::kAOD "<<AliCaloTrackReader::kAOD<<endl;
 862   //  cout<<"GetReader()->GetInputEvent() "<<GetReader()->GetInputEvent()<<endl;
 863
 864   if(GetReader()->GetOutputEvent())
 865   {
 866     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
 867   }
 868   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
 869   {
 870     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
 871   }
 872   else
 873   {
 874     if(GetDebug() > 3) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - There are no jets available for this analysis\n");
 875     return;
 876   }
 877
 878   if(!GetInputAODBranch() || !event){
 879     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s > \n",GetInputAODName().Data());
 880     abort();
 881   }
 882
 883   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
 884     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
 885     abort();
 886   }
 887
 888   //
 889   // non-standard jets
 890   //
 891   Int_t nJets=-1;
 892   TClonesArray *aodRecJets = 0;
 893   if(IsNonStandardJetFromReader()){//jet branch from reader
 894     if(GetDebug() > 3) printf("GetNonStandardJets function (from reader) is called\n");
 895     aodRecJets = GetNonStandardJets();
 896     if(GetDebug() > 3) printf("aodRecJets %p\n",aodRecJets);
 897     if(aodRecJets==0x0){
 898       if(GetDebug() > 3) event->Print();
 899       abort();
 900     }
 901     nJets=aodRecJets->GetEntries();
 902     if(GetDebug() > 3) printf("nJets %d\n",nJets);
 903   }
 904   //end of new part
 905
 906   if(nJets==0) {
 907     //printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
 908     return;
 909   }
 910
 911   //
 912   //Background jets
 913   //
 914   AliAODJetEventBackground* aodBkgJets = 0;
 915   if(IsBackgroundJetFromReader()){//jet branch from reader
 916     if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
 917     aodBkgJets = GetBackgroundJets();
 918     if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
 919     if(aodBkgJets==0x0){
 920       if(GetDebug() > 3) event->Print();
 921       abort();
 922     }
 923     if(GetDebug() > 3) aodBkgJets->Print("c");
 924   }
 925
 926   Double_t rhoEvent=0.;
 927   if(aodBkgJets && IsBackgroundJetFromReader() ) {
 928     rhoEvent = aodBkgJets->GetBackground(2);//hardcoded
 929   }
 930   fJetRho = rhoEvent;
 931
 932
 933   Int_t ntrig =  GetInputAODBranch()->GetEntriesFast() ;
 934   if(GetDebug() > 3){
 935     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Begin jet finder  correlation analysis, fill AODs \n");
 936     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", ntrig);
 937     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In standard jet branch aod entries %d\n", event->GetNJets());
 938     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - In non standard jet branch aod entries %d\n", nJets);
 939   }
 940
 941   //if(nJets==0)   return;//to speed up
 942   //  cout<<"ntrig po return "<<ntrig<<endl;
 943
 944   //
 945   //calculate average cell energy without most energetic photon
 946   //
 947   Double_t medianPhotonRho=0.;
 948   TObjArray* clusters = GetEMCALClusters();
 949   Int_t clusterIDtmp;
 950   Int_t iclustmp = -1;
 951   AliVCluster *cluster=0;
 952
 953   if(IsBackgroundSubtractionGamma()){
 954     //
 955     // Find most energetic photon without background subtraction
 956     //
 957     Double_t maxPt=0.;
 958     Int_t maxIndex=-1;
 959     AliAODPWG4ParticleCorrelation* particlecorr =0;
 960     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
 961       particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 962       if(particlecorr->Pt() > maxPt) {
 963         maxPt = particlecorr->Pt();
 964         maxIndex = iaod;
 965       }
 966     }
 967
 968     //
 969     // calculate background energy per cell
 970     //
 971     Int_t numberOfcells=0;
 972     if(ntrig>1){
 973       Double_t *photonRhoArr=new Double_t[ntrig-1];
 974       Int_t photonRhoArrayIndex=0;
 975       for(Int_t iaod = 0; iaod < ntrig ; iaod++){
 976         particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 977         if(iaod==maxIndex) continue;
 978         clusterIDtmp = particlecorr->GetCaloLabel(0) ;
 979         if(clusterIDtmp < 0) continue;
 980         cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
 981         photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
 982         numberOfcells+=cluster->GetNCells();
 983         photonRhoArrayIndex++;
 984       }
 985       if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
 986       delete photonRhoArr;
 987     }
 988   }//end of if background calculation for gamma
 989   fGamRho = medianPhotonRho;
 990
 991
 992   //
 993   //take most energetic photon and most energetic jet and combine
 994   //
 995   if(fMostEnergetic){
 996     //
 997     // most energetic photon with background subtraction
 998     //
 999     Double_t mostEnePhotonPt=0.;
1000     Int_t indexMostEnePhoton=-1;
1001     AliAODPWG4ParticleCorrelation* particle =0;
1002     Double_t ptCorrect=0.;
1003     Int_t nCells=0;
1004     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1005       particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1006       clusterIDtmp = particle->GetCaloLabel(0) ;
1007       if(!(clusterIDtmp<0)){
1008         cluster = FindCluster(clusters,clusterIDtmp,iclustmp);
1009         nCells = cluster->GetNCells();
1010       }
1011       ptCorrect = particle->Pt() - medianPhotonRho * nCells;
1012       if( ptCorrect > mostEnePhotonPt ){
1013         mostEnePhotonPt = ptCorrect;
1014         indexMostEnePhoton = iaod ;
1015       }
1016     }
1017     //    printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Photon with index %d selected \n",indexMostEnePhoton);
1018     //
1019     // most energetic jet with background subtraction
1020     //
1021     Double_t mostEneJetPt=0.;
1022     Int_t indexMostEneJet=-1;
1023     AliAODJet * jet = 0 ;
1024     //Double_t ptCorrect=0.;
1025     for(Int_t ijet = 0; ijet < nJets ; ijet++){
1026       jet = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1027       if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
1028       if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
1029       if(jet->Pt()<fJetMinPt) continue;
1030       ptCorrect = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();
1031       if(ptCorrect > mostEneJetPt){
1032         mostEneJetPt = ptCorrect;
1033         indexMostEneJet = ijet;
1034       }
1035     }
1036     //    printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",indexMostEneJet);
1037
1038     //
1039     // assign jet to photon
1040     //
1041     if(indexMostEneJet>=0 && indexMostEnePhoton>=0){
1042       particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(indexMostEnePhoton));
1043       jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(indexMostEneJet));
1044       particle->SetRefJet(jet);
1045     }
1046   }//end of take most energetic photon and most ene. jet after bkg subtraction
1047
1048   if(fMostOpposite){
1049     //Bool_t isJetFound=kFALSE;//new
1050     //Loop on stored AOD particles, trigger
1051     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1052       AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1053
1054       //Correlate with jets
1055       Int_t ijet = SelectJet(particle,aodRecJets);//input for jets is TClonesArray
1056       if(ijet > -1){
1057         //isJetFound=kTRUE;
1058         if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",ijet);
1059         AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1060         particle->SetRefJet(jet);
1061         //printf("Most opposite found\n");
1062       }
1063     } // input aod loop
1064     //  if(GetReader()->WriteDeltaAODToFile() && isJetFound) WriteJetsToOutputBranch(aodRecJets);
1065   }//end of take most opposite photon and jet after bkg subtraction
1066
1067
1068   if(GetDebug() > 1 ) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
1069 }
1070
1071 //__________________________________________________________________
1072 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms()
1073 {
1074   //Particle-Jet Correlation Analysis, fill histograms
1075   if(GetDebug() > 3 ) {
1076     printf("I use MakeAnalysisFillHistograms\n");
1077     printf("ntracks before iso %d\n",GetCTSTracks()->GetEntriesFast() );
1078   }
1079
1080   //Get the event, check if there are AODs available, if not, skip this analysis
1081   AliAODEvent * event = NULL;
1082
1083   //printf("GetReader()->GetOutputEvent() %d\n",GetReader()->GetOutputEvent() );
1084   //printf("GetReader()->GetDataType() %d\n",GetReader()->GetDataType() );
1085   //printf("AliCaloTrackReader::kAOD %d\n",AliCaloTrackReader::kAOD );
1086   //printf("GetReader()->GetInputEvent() %d\n",GetReader()->GetInputEvent() );
1087
1088   if(GetReader()->GetOutputEvent())
1089   {
1090     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1091   }
1092   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1093   {
1094     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1095   }
1096   else {
1097     if(GetDebug() > 3) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - There are no jets available for this analysis\n");
1098     return;
1099   }
1100
1101   if(!GetInputAODBranch() || !event){
1102     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s > \n",GetInputAODName().Data());
1103     abort();
1104   }
1105
1106   Int_t nJets=-1;
1107   TClonesArray *aodRecJets = 0;
1108   if(IsNonStandardJetFromReader()){//branch read in reader from reader
1109     if (GetDebug() > 3) printf("GetNonStandardJets function (from reader) is called\n");
1110     aodRecJets = GetNonStandardJets();
1111     if(aodRecJets==0x0){
1112       if(GetDebug() > 3) event->Print();
1113       AliFatal("Jets container not found");
1114       return; // trick coverity
1115     }
1116     nJets=aodRecJets->GetEntries();
1117   }
1118   if(nJets==0) {
1119     //    printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1120     GetReader()->FillInputNonStandardJets();
1121     aodRecJets = GetNonStandardJets();
1122     if(aodRecJets) nJets=aodRecJets->GetEntries();
1123     //    printf("nJets = %d\n",nJets);
1124     return;
1125   }
1126
1127   //
1128   //Background jets
1129   //
1130   AliAODJetEventBackground* aodBkgJets = 0;
1131   if(IsBackgroundJetFromReader()){//jet branch from reader
1132     if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
1133     aodBkgJets = GetBackgroundJets();
1134     if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
1135     if(aodBkgJets==0x0){
1136       if(GetDebug() > 3) event->Print();
1137       abort();
1138     }
1139     if(GetDebug() > 3) aodBkgJets->Print("c");
1140   }
1141
1142
1143   //
1144   // only background jets informations
1145   //
1146   //  Float_t pTback = 0;
1147   Double_t rhoEvent=0.;
1148   if(aodBkgJets) {
1149     if(IsBackgroundJetFromReader() ) rhoEvent = aodBkgJets->GetBackground(2);
1150     if(IsHistogramJetBkg()) {
1151       fhJetRhoVsCentrality->Fill(GetEventCentrality(),rhoEvent);
1152       for(Int_t i=0;i<4;i++){
1153         fhBkgJetBackground[i]->Fill(aodBkgJets->GetBackground(i));
1154         fhBkgJetSigma[i]->Fill(aodBkgJets->GetSigma(i));
1155         fhBkgJetArea[i]->Fill(aodBkgJets->GetMeanarea(i));
1156       }
1157     }//end of if fill HistogramJetBkg
1158   }//end if aodBkgJets exists
1159
1160   //
1161   //only track information
1162   //
1163   Int_t nCTSTracks = GetCTSTracks()->GetEntriesFast();
1164   AliAODTrack *aodtrack;
1165   Int_t itrack = 0;
1166   if(IsHistogramTracks()) {
1167     Double_t sumTrackPt=0;
1168     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1169       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1170       fhTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1171       sumTrackPt+=aodtrack->Pt();
1172     }
1173     if(nCTSTracks)
1174       fhTrackAveTrackPt->Fill(sumTrackPt/nCTSTracks);
1175   }//end if IsHistogramTracks
1176
1177   //
1178   //only jet informations
1179   //
1180   AliAODJet * jettmp = 0 ;
1181   Double_t jetPttmp = 0.;
1182   Double_t var1tmp = 0.;
1183   Double_t var2tmp = 0.;
1184   //  fhJetNjet->Fill(nJets);
1185   Double_t ptMostEne=0;
1186   //  Int_t indexMostEne=-1;
1187   Int_t nJetsOverThreshold[10]={nJets,0,0,0,0,0,0,0,0,0};
1188   Int_t iCounter=0;
1189   Double_t sumJetTrackPt=0.;
1190   Int_t sumNJetTrack=0;
1191   Int_t nTracksInJet=0;
1192   Int_t itrk=0;
1193   for(Int_t ijet = 0; ijet < nJets ; ijet++){
1194     jettmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1195     fhJetPtBefore->Fill(jettmp->Pt());
1196     jetPttmp  = jettmp->Pt() - rhoEvent * jettmp->EffectiveAreaCharged();//<<---changed here
1197
1198     //calculate number of tracks above 1,2,3,4,5 GeV in jet
1199     AliVTrack* jettrack = 0x0 ;
1200
1201     Int_t nTrackThrGeV[5]={0,0,0,0,0};
1202     nTracksInJet=(jettmp->GetRefTracks())->GetEntriesFast();
1203     fhJetNparticlesInJet->Fill(nTracksInJet);
1204     if(nTracksInJet==0) continue;
1205     sumNJetTrack+=nTracksInJet;
1206     for(itrack=0;itrack<nTracksInJet;itrack++){
1207       jettrack=(AliVTrack *) ((jettmp->GetRefTracks())->At(itrack));
1208       if(!jettrack) continue;
1209
1210       fhJetDeltaEtaDeltaPhi->Fill(jettmp->Eta()-jettrack->Eta(),jettmp->Phi()-jettrack->Phi());
1211       sumJetTrackPt+=jettrack->Pt();
1212       if(IsHistogramJetTracks()){
1213         if(jettrack->Pt()>1.) nTrackThrGeV[0]++;
1214         if(jettrack->Pt()>2.) nTrackThrGeV[1]++;
1215         if(jettrack->Pt()>3.) nTrackThrGeV[2]++;
1216         if(jettrack->Pt()>4.) nTrackThrGeV[3]++;
1217         if(jettrack->Pt()>5.) nTrackThrGeV[4]++;
1218       }//end of if IsHistogramJetTracks
1219     }//end of jet track loop
1220     //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]);
1221
1222     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1223       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1224       if(aodtrack) fhJetDeltaEtaDeltaPhiAllTracks->Fill(jettmp->Eta()-aodtrack->Eta(),jettmp->Phi()-aodtrack->Phi());
1225     }
1226
1227
1228     if(IsHistogramJetTracks()){
1229       fhJetNtracksInJetAboveThr[0]->Fill(jetPttmp,nTracksInJet);//all jets
1230
1231       for(itrk=0;itrk<5;itrk++) {
1232         fhJetNtracksInJetAboveThr[itrk+1]->Fill(jetPttmp,nTrackThrGeV[itrk]);//all jets
1233         fhJetRatioNTrkAboveToNTrk[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);//all jets
1234       }
1235       if(ijet==0){//most ene jet
1236         for(itrk=0;itrk<5;itrk++)
1237           fhJetNtrackRatioMostEne[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1238       }
1239       if(jetPttmp>5){//jet with pt>5GeV/c
1240         for(itrk=0;itrk<5;itrk++)
1241           fhJetNtrackRatioJet5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1242       }
1243       if(nTrackThrGeV[4]>0){//jet with leading particle pt>5GeV/c
1244         for(itrk=0;itrk<5;itrk++)
1245           fhJetNtrackRatioLead5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1246       }
1247     }//end of if IsHistogramJetTracks
1248
1249
1250     Double_t effectiveChargedBgEnergy=(IsBackgroundJetFromReader()?rhoEvent * jettmp->EffectiveAreaCharged():jettmp->ChargedBgEnergy());
1251
1252
1253     fhJetChBkgEnergyVsArea->Fill(effectiveChargedBgEnergy,jettmp->EffectiveAreaCharged());
1254     //if(jettmp->EffectiveAreaCharged()>0)
1255     fhJetRhoVsPt->Fill(jetPttmp,jettmp->ChargedBgEnergy()*jettmp->EffectiveAreaCharged());
1256
1257     if(jetPttmp>ptMostEne) {
1258       ptMostEne = jetPttmp;
1259       //indexMostEne=ijet;
1260     }
1261     fhJetPt->Fill(jetPttmp);
1262     fhJetChBkgEnergyVsPt->Fill(jetPttmp,effectiveChargedBgEnergy);
1263     fhJetChAreaVsPt->Fill(jetPttmp,jettmp->EffectiveAreaCharged());
1264     if(GetDebug()>5) printf("ChargedBgEnergy %f EffectiveAreaCharged %f\n", jettmp->ChargedBgEnergy(),jettmp->EffectiveAreaCharged());
1265     for(iCounter=1;iCounter<10;iCounter++){
1266       if(jetPttmp>iCounter) nJetsOverThreshold[iCounter]++;
1267     }
1268
1269     var1tmp  = jettmp->Phi();
1270     var2tmp  = jettmp->Eta();
1271     fhJetPhi->Fill(var1tmp);
1272     fhJetEta->Fill(var2tmp);
1273     fhJetPhiVsEta->Fill(var1tmp,var2tmp);
1274     fhJetEtaVsPt->Fill(jetPttmp,var2tmp);
1275     fhJetEtaVsNpartInJet->Fill(var2tmp,nTracksInJet);
1276     if(jetPttmp>0)
1277       fhJetEtaVsNpartInJetBkg->Fill(var2tmp,nTracksInJet);
1278
1279   }//end of jet loop
1280   if(IsHistogramJetTracks()){
1281     if(sumNJetTrack>0){
1282       //printf("average track pt %f\n",sumJetTrackPt/sumNJetTrack);
1283       fhJetAveTrackPt->Fill(sumJetTrackPt/sumNJetTrack);
1284     }
1285   }//end of if IsHistogramJetTracks
1286
1287
1288   fhJetPtMostEne->Fill(ptMostEne);
1289   for(iCounter=0;iCounter<10;iCounter++){
1290     fhJetNjetOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1291     nJetsOverThreshold[iCounter]=0;
1292   }
1293
1294   //end of only jet part
1295
1296   //
1297   // Photons
1298   //
1299   Int_t ntrig   =  GetInputAODBranch()->GetEntriesFast() ;
1300   if(GetDebug() > 1){
1301     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Begin jet finder  correlation analysis, fill histograms \n");
1302     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", ntrig);
1303     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In jet output branch aod entries %d\n", event->GetNJets());
1304   }
1305   fhNjetsNgammas->Fill(nJets,ntrig);
1306
1307   //if(nJets==0)   return;//to speed up
1308   //printf("ntrig %d, njets %d\n",ntrig,nJets);
1309
1310   //
1311   //Fill only temporary photon histograms
1312   //
1313   Double_t maxPt=0.;
1314   Int_t maxIndex=-1;
1315   Double_t tmpPt=0.;
1316   Double_t sumPt=0.;
1317   nJetsOverThreshold[0]=ntrig;
1318   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1319     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1320     tmpPt = particlecorr->Pt();
1321     if(tmpPt>maxPt) {
1322       maxPt = tmpPt;
1323       maxIndex = iaod;
1324     }
1325     for(iCounter=1;iCounter<10;iCounter++){
1326       if(tmpPt>iCounter) nJetsOverThreshold[iCounter]++;
1327     }
1328     sumPt+=tmpPt;
1329   }
1330   for(iCounter=0;iCounter<10;iCounter++){
1331     fhPhotonNgammaOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1332     //    nJetsOverThreshold[iCounter]=0;
1333   }
1334
1335   fGamAvEne=0;
1336   TObjArray* clusters = GetEMCALClusters();
1337   //printf("calculate median bkg energy for photons ");
1338   Double_t medianPhotonRho=0.;
1339   Int_t clusterID;
1340   Int_t iclustmp = -1;
1341   Int_t numberOfcells=0;
1342   AliVCluster *cluster = 0;
1343   if(ntrig>1){
1344     Double_t *photonRhoArr=new Double_t[ntrig-1];
1345     fhPhotonPtMostEne->Fill(maxPt);
1346     //    fhPhotonIndexMostEne->Fill(indexMaxPt);
1347     fhPhotonAverageEnergy->Fill(sumPt/ntrig);
1348     fhPhotonRatioAveEneToMostEne->Fill(sumPt/(ntrig*maxPt));
1349     fhPhotonAverageEnergyMinus1->Fill((sumPt-maxPt)/(ntrig-1));
1350     fGamAvEne=(sumPt-maxPt)/(ntrig-1);
1351     fhPhotonRatioAveEneMinus1ToMostEne->Fill((sumPt-maxPt)/((ntrig-1)*maxPt));
1352
1353     Int_t counterGamma=0;
1354     Int_t counterGammaMinus1=0;
1355
1356     Int_t photonRhoArrayIndex=0;
1357     //TObjArray* clusterstmp = GetEMCALClusters();
1358     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1359       AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1360       if( particlecorr->Pt() > sumPt/ntrig ) counterGamma++;
1361       if( particlecorr->Pt() > (sumPt-maxPt)/(ntrig-1) ) counterGammaMinus1++;
1362
1363       if(iaod==maxIndex) continue;
1364       clusterID = particlecorr->GetCaloLabel(0) ;
1365       if(clusterID < 0) continue;
1366       cluster = FindCluster(clusters,clusterID,iclustmp);
1367       photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1368       numberOfcells+=cluster->GetNCells();
1369       photonRhoArrayIndex++;
1370     }
1371     if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1372     delete photonRhoArr;
1373     fhPhotonNgammaMoreAverageToNgamma->Fill((Double_t)counterGamma / (Double_t)ntrig);
1374     fhPhotonNgammaMoreAverageMinus1ToNgamma->Fill((Double_t)counterGammaMinus1 / (Double_t)ntrig);
1375   }
1376   //printf("median = %f\n",medianPhotonRho);
1377   fhPhotonBkgRhoVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),medianPhotonRho);
1378   fhPhotonBkgRhoVsNclusters->Fill(ntrig,medianPhotonRho);
1379   fhPhotonBkgRhoVsCentrality->Fill(GetEventCentrality(),medianPhotonRho);
1380   fhPhotonBkgRhoVsNcells->Fill(numberOfcells,medianPhotonRho);
1381
1382
1383   AliVCluster *cluster2 = 0;
1384   Double_t photon2Corrected=0;
1385   Double_t sumPtTmp=0.;
1386   Double_t sumPtCorrectTmp=0.;
1387   AliVTrack* trackTmp = 0x0 ;
1388   TVector3 p3Tmp;
1389
1390   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1391     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1392     clusterID = particlecorr->GetCaloLabel(0) ;
1393     if(clusterID < 0) continue;
1394     cluster = FindCluster(clusters,clusterID,iclustmp);
1395     fhPhotonPt->Fill(particlecorr->Pt());
1396     fhPhotonPtCorrected->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1397     fhPhotonPtDiff->Fill(cluster->GetNCells() * medianPhotonRho);
1398     fhPhotonPtDiffVsCentrality->Fill(GetEventCentrality(),cluster->GetNCells() * medianPhotonRho);
1399     fhPhotonPtDiffVsNcells->Fill(numberOfcells,cluster->GetNCells() * medianPhotonRho);
1400     fhPhotonPtDiffVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),cluster->GetNCells() * medianPhotonRho);
1401     fhPhotonPtDiffVsNclusters->Fill(ntrig,cluster->GetNCells() * medianPhotonRho);
1402
1403     fhPhotonPtCorrectedZoom->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1404
1405     //test: sum_pt in the cone 0.3 for each photon
1406     //should be: random fake gamma from MB
1407     //is: each gamma for EMCEGA
1408     sumPtTmp=0.;
1409     sumPtCorrectTmp=0.;
1410
1411     for(Int_t iaod2 = 0; iaod2 < ntrig ; iaod2++){
1412       if(iaod==iaod2) continue;
1413       AliAODPWG4ParticleCorrelation* particlecorr2 =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod2));
1414       clusterID = particlecorr2->GetCaloLabel(0) ;
1415       if(clusterID < 0) continue;
1416       cluster2 = FindCluster(clusters,clusterID,iclustmp);
1417       photon2Corrected = particlecorr2->Pt() - cluster2->GetNCells() * medianPhotonRho;
1418
1419       //if(Pt()<0.5) continue; //<<hardcoded here //FIXME
1420       if( TMath::Sqrt((particlecorr->Eta()-particlecorr2->Eta())*(particlecorr->Eta()-particlecorr2->Eta()) +
1421                       (particlecorr->Phi()-particlecorr2->Phi())*(particlecorr->Phi()-particlecorr2->Phi()) )<fGammaConeSize ){//if(/*cone is correct*/){
1422         sumPtTmp+= particlecorr2->Pt();
1423         sumPtCorrectTmp+=photon2Corrected;
1424       }
1425     }
1426     fhPhotonSumPtInCone->Fill(sumPtTmp);
1427     fhPhotonSumPtCorrectInCone->Fill(sumPtCorrectTmp);
1428
1429     //test: sum_pt in the cone 0.3 for each track
1430     //should be: random fake gamma from MB
1431     //is: each gamma for EMCEGA
1432     sumPtTmp=0.;
1433     for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++){
1434       trackTmp = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1435       p3Tmp.SetXYZ(trackTmp->Px(),trackTmp->Py(),trackTmp->Pz());
1436       if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1437                       (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1438         sumPtTmp+=p3Tmp.Pt();
1439       }
1440     }//end of loop over tracks
1441     fhPhotonSumPtChargedInCone->Fill(sumPtTmp);
1442   }
1443
1444   //End of Fill temporary photon histograms
1445
1446   //
1447   // Apply background subtraction for photons
1448   //
1449   fGamRho = medianPhotonRho;
1450   if(!IsBackgroundSubtractionGamma()) medianPhotonRho=0;
1451
1452
1453   //
1454   //Get vertex for cluster momentum calculation <<----new here
1455   //
1456   Double_t vertex[] = {0,0,0} ; //vertex ;
1457   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1458     GetReader()->GetVertex(vertex);
1459   fZvertex = vertex[2];
1460
1461   //
1462   //Loop on stored AOD particles, trigger
1463   //
1464   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1465     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1466     fhCuts->Fill(0);
1467     fhCuts2->Fill(0.,(Double_t)nJets);
1468     if(GetDebug() > 5) printf("OnlyIsolated %d  !particlecorr->IsIsolated() %d \n",OnlyIsolated(), !particlecorr->IsIsolated());
1469
1470     if(OnlyIsolated() && !particlecorr->IsIsolated()) continue;
1471     fhCuts->Fill(1);
1472     fhCuts2->Fill(1.,nJets);
1473
1474     if(nJets>0) {
1475       fhCuts->Fill(2);
1476     }
1477
1478     //Recover the jet correlated, found previously.
1479     AliAODJet   * jet = particlecorr->GetJet();
1480     //If correlation not made before, do it now.
1481     if(fMakeCorrelationInHistoMaker){
1482       //Correlate with jets
1483       Int_t ijet = SelectJet(particlecorr,aodRecJets);//input for jets is TClonesArray
1484       if(ijet > -1){
1485         if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Jet with index %d selected \n",ijet);
1486         //jet = event->GetJet(ijet);
1487         jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1488
1489         particlecorr->SetRefJet(jet);
1490
1491       }
1492     }
1493
1494     if (!jet) continue ;
1495     fhCuts->Fill(3);
1496     fhCuts2->Fill(7.,1.);
1497
1498     //
1499     //Fill Correlation Histograms
1500     //
1501     clusterID = particlecorr->GetCaloLabel(0) ;
1502     if(!(clusterID<0)){
1503       cluster = FindCluster(clusters,clusterID,iclustmp);
1504       //fill tree variables
1505       fGamNcells = cluster->GetNCells();
1506     }
1507     Double_t ptTrig = particlecorr->Pt() - medianPhotonRho * fGamNcells;//<<---changed here
1508     Double_t ptJet = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();//<<---changed here
1509     Double_t phiTrig = particlecorr->Phi();
1510     Double_t phiJet = jet->Phi();
1511     Double_t etaTrig = particlecorr->Eta();
1512     Double_t etaJet = jet->Eta();
1513     Double_t deltaPhi=phiTrig-phiJet;
1514     if(deltaPhi<0)deltaPhi+=(TMath::Pi()*2.);
1515     //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",
1516     //  ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
1517     fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet);
1518     //    fhDeltaPhi->Fill(ptTrig, phiTrig-phiJet);//need to be shifted by 2pi
1519
1520     fhDeltaPhiCorrect->Fill(ptTrig, deltaPhi);//correct
1521
1522     Double_t deltaPhiCorrect = TMath::Abs( particlecorr->Phi() - jet->Phi() );
1523     if ( deltaPhiCorrect > TMath::Pi() ) deltaPhiCorrect = 2. * TMath::Pi() - deltaPhiCorrect ;
1524     fhDeltaPhi0PiCorrect->Fill(ptTrig, deltaPhiCorrect);
1525
1526     fhDeltaEta->Fill(ptTrig, etaTrig-etaJet);
1527     fhPtRatio ->Fill(ptTrig, ptJet/ptTrig);
1528     fhPt      ->Fill(ptTrig, ptJet);
1529
1530     fhSelectedJetPhiVsEta->Fill(phiJet,etaJet);
1531     fhSelectedJetChBkgEnergyVsPtJet->Fill(ptJet,(IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy()) );
1532     fhSelectedJetChAreaVsPtJet->Fill(ptJet,jet->EffectiveAreaCharged());
1533     fhSelectedJetNjet->Fill(nJets);
1534     fhSelectedNtracks->Fill(GetCTSTracks()->GetEntriesFast());//to be checked
1535     fhSelectedPhotonNLMVsPt->Fill(ptTrig,particlecorr->GetFiducialArea());
1536
1537
1538     if(clusterID < 0 ){
1539       fhSelectedPhotonLambda0VsPt->Fill(ptTrig,-1);
1540       //fill tree variables
1541       fGamLambda0  = -1;
1542       fGamTime = -1;
1543       fGamNcells = 0;
1544       fGamSumPtNeu=0;
1545     }
1546     else{
1547       //Int_t iclus = -1;
1548       //      TObjArray* clusters = GetEMCALClusters();
1549       //cluster = FindCluster(clusters,clusterID,iclustmp);
1550       Double_t lambda0=cluster->GetM02();
1551       fhSelectedPhotonLambda0VsPt->Fill(ptTrig,lambda0);
1552       //fill tree variables
1553       fGamLambda0  = lambda0;
1554       fGamTime = cluster->GetTOF();
1555       //fGamNcells = cluster->GetNCells();
1556
1557       fGamSumPtNeu=0;
1558       fGamNclusters=0;
1559       TLorentzVector mom ;
1560       //TVector3 p3Tmp;
1561       //Double_t scalarProduct=0;
1562       //Double_t vectorLength=particlecorr->P();
1563       for(Int_t icalo=0; icalo <clusters->GetEntriesFast(); icalo++){
1564         AliVCluster* calo = (AliVCluster *) clusters->At(icalo);
1565         if(clusterID==calo->GetID()) continue;//the same cluster as trigger
1566         calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
1567         //printf("min pt %f\n",GetMinPt());
1568         if(mom.Pt()<GetMinPt()) continue; //<<hardcoded here //FIXME 0.5 check if correct
1569         p3Tmp.SetXYZ(mom.Px(),mom.Py(),mom.Pz());
1570         //calculate sum pt in the cone
1571         if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1572                         (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1573           //scalarProduct = particlecorr->Px()*mom.Px() + particlecorr->Py()*mom.Py() + particlecorr->Pz()*mom.Pz();
1574           //scalarProduct/=mom.P();
1575           //scalarProduct/=vectorLength;
1576           //if(scalarProduct>TMath::Cos(0.3)) {//FIXME photon radius
1577           fGamSumPtNeu+=mom.Pt();
1578           fGamNclusters++;
1579         }
1580       }
1581     }
1582
1583     //sum pt of charged tracks in the gamma isolation cone
1584     //starts here
1585     fGamSumPtCh=0;
1586     fGamNtracks=0;
1587     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1588       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1589       fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());//fill histogram here
1590       //      if(aodtrack->Pt()<0.15) continue;//hardcoded
1591       if(aodtrack->Pt()<fPtThresholdInCone) continue;
1592       if(!aodtrack->IsHybridGlobalConstrainedGlobal()) continue;
1593       if(TMath::Sqrt((particlecorr->Phi() - aodtrack->Phi())*(particlecorr->Phi() - aodtrack->Phi()) +
1594                      (particlecorr->Eta() - aodtrack->Eta())*(particlecorr->Eta() - aodtrack->Eta()) ) <fGammaConeSize ) {
1595         fGamSumPtCh+=aodtrack->Pt();
1596         fGamNtracks++;
1597       }
1598     }
1599     //ends here
1600
1601     //    for(Int_t itrack = 0; itrack < nCTSTracks ; itrack++){
1602     //      aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1603     //      fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1604     //    }
1605
1606     //
1607     // Background Fragmentation function
1608     //
1609     TVector3 gammaVector,jetVector;
1610     gammaVector.SetXYZ(particlecorr->Px(),particlecorr->Py(),particlecorr->Pz());
1611     jetVector.SetXYZ(jet->Px(),jet->Py(),jet->Pz());
1612     CalculateBkg(gammaVector,jetVector,vertex,1);//jet perp
1613     CalculateBkg(gammaVector,jetVector,vertex,2);//RC
1614     CalculateBkg(gammaVector,jetVector,vertex,3);//mid point
1615     CalculateBkg(gammaVector,jetVector,vertex,4);//gamma perp
1616     //CalculateBkg(gammaVector,jetVector,vertex,5);/test
1617     Double_t angleJetGam = gammaVector.Angle(jetVector);
1618     //printf("angleJetGam %f\n",angleJetGam*180/TMath::Pi());
1619
1620     //
1621     // Fragmentation function
1622     //
1623     Float_t      rad = 0, pt = 0, eta = 0, phi = 0;
1624     Int_t        npartcone = 0;
1625     TVector3 p3;
1626
1627     Int_t ntracks =  0;
1628     if(GetDebug()>3){
1629       printf ("fUseJetRefTracks %d\n",fUseJetRefTracks );
1630       printf ("jet->GetRefTracks() %p\n",jet->GetRefTracks());
1631       printf ("GetCTSTracks() %p\n",GetCTSTracks() );
1632     }
1633
1634     if(!fUseJetRefTracks)
1635       ntracks =GetCTSTracks()->GetEntriesFast();
1636     else //If you want to use jet tracks from JETAN
1637       ntracks =  (jet->GetRefTracks())->GetEntriesFast();
1638
1639     if(GetDebug()>3)    printf ("ntracks %d\n",ntracks);
1640     AliVTrack* track = 0x0 ;
1641     for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
1642       if(!fUseJetRefTracks)
1643         track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1644       else //If you want to use jet tracks from JETAN
1645         track = (AliVTrack *) ((jet->GetRefTracks())->At(ipr));
1646
1647       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1648       pt    = p3.Pt();
1649       eta  = p3.Eta();
1650       phi  = p3.Phi() ;
1651       if(phi < 0) phi+=TMath::TwoPi();
1652
1653       //Check if there is any particle inside cone with pt larger than  fPtThreshold
1654       rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
1655       if(rad < fConeSize  && pt > fPtThresholdInCone){
1656         //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
1657         npartcone++;
1658         fhFFz ->Fill(ptTrig, pt/ptTrig);
1659         fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt));
1660         fhFFpt->Fill(ptTrig, pt);
1661
1662         //according to jet axis
1663         fhJetFFz ->Fill(ptJet, pt/ptJet);
1664         fhJetFFxi->Fill(ptJet, TMath::Log(ptJet/pt));
1665         fhJetFFpt->Fill(ptJet, pt);
1666
1667
1668         if(TMath::Cos(angleJetGam)<0 && ptJet!=0 && pt!=0 ){
1669           fhJetFFzCor ->Fill(ptJet, -pt*TMath::Cos(angleJetGam)/ptJet);
1670           fhJetFFxiCor->Fill(ptJet, TMath::Log(ptJet/(-pt*TMath::Cos(angleJetGam))));
1671         }
1672       }
1673     }//Tracks
1674     fhNTracksInCone->Fill(ptTrig, npartcone);
1675     //fill tree here for each photon-jet (isolation required usually)
1676
1677     fGamPt      = ptTrig;
1678     //fGamLambda0  = ;//filled earlier
1679     fGamNLM      = particlecorr->GetFiducialArea();
1680     //fGamSumPtCh  = ;//filled earlier
1681     //fGamTime     = particlecorr->GetTOF();//filled earlier
1682     //fGamNcells   = particlecorr->GetNCells();//filled earlier
1683     fGamEta      = etaTrig;
1684     fGamPhi      = phiTrig;
1685     //fGamSumPtNeu = ;//filled earlier
1686     //fGamNtracks  = ;//filled earlier
1687     //fGamNclusters= ;//filled earlier
1688     //fGamAvEne    = ;//filled earlier
1689     fJetPhi      = phiJet;
1690     fJetEta      = etaJet;
1691     fJetPt       = ptJet;
1692     fJetBkgChEne = (IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy());
1693     fJetArea     = jet->EffectiveAreaCharged();
1694     fJetNtracks  = (jet->GetRefTracks())->GetEntriesFast();
1695     fEventNumber = 0;
1696     fNtracks     = GetCTSTracks()->GetEntriesFast();
1697     fCentrality  = GetEventCentrality();
1698     fIso         = particlecorr->IsIsolated();
1699
1700     Int_t nTrk1GeV=0;
1701     Int_t nTrk2GeV=0;
1702     for(itrack=0;itrack < fJetNtracks;itrack++){
1703       track = (AliVTrack *) ((jet->GetRefTracks())->At(itrack));
1704       if(track->Pt()>1.) nTrk1GeV++;
1705       if(track->Pt()>2.) nTrk2GeV++;
1706     }
1707
1708     fJetNtracks1 = nTrk1GeV;
1709     fJetNtracks2 = nTrk2GeV;
1710
1711     if(fSaveGJTree) fTreeGJ->Fill();
1712   }//AOD trigger particle loop
1713   if(GetDebug() > 1) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
1714 }
1715
1716
1717 //__________________________________________________________________
1718 void AliAnaParticleJetFinderCorrelation::Print(const Option_t * opt) const
1719 {
1720
1721   //Print some relevant parameters set for the analysis
1722   if(! opt)
1723     return;
1724
1725   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1726   AliAnaCaloTrackCorrBaseClass::Print(" ");
1727
1728   printf("Phi trigger-jet        <     %3.2f\n", fDeltaPhiMaxCut) ;
1729   printf("Phi trigger-jet        >     %3.2f\n", fDeltaPhiMinCut) ;
1730   printf("pT Ratio trigger/jet   <     %3.2f\n", fRatioMaxCut) ;
1731   printf("pT Ratio trigger/jet   >     %3.2f\n", fRatioMinCut) ;
1732   printf("fConeSize              =     %3.2f\n", fConeSize) ;
1733   printf("fPtThresholdInCone     =     %3.2f\n", fPtThresholdInCone) ;
1734   printf("fUseJetRefTracks         =     %d\n",    fUseJetRefTracks) ;
1735   printf("fMakeCorrelationInHistoMaker     =     %d\n",    fMakeCorrelationInHistoMaker) ;
1736   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
1737   printf("Reconstructed jet cone size = %3.2f\n", fJetConeSize) ;
1738   printf("Reconstructed jet minimum pt = %3.2f\n", fJetMinPt) ;
1739   printf("Reconstructed jet minimum area fraction = %3.2f\n", fJetAreaFraction) ;
1740
1741   if(!fNonStandardJetFromReader){
1742     printf("fJetBranchName =   %s\n", fJetBranchName.Data()) ;
1743   }
1744   if(!fBackgroundJetFromReader){
1745     printf("fBkgJetBranchName =   %s\n", fBkgJetBranchName.Data()) ;
1746   }
1747
1748   printf("Isolation cone size = %3.2f\n", fGammaConeSize) ;
1749   printf("fUseBackgroundSubtractionGamma = %d\n",fUseBackgroundSubtractionGamma);
1750   printf("fSaveGJTree = %d\n",fSaveGJTree);
1751   printf("fMostEnergetic = %d\n",fMostEnergetic);
1752   printf("fMostOpposite = %d\n",fMostOpposite);
1753
1754   printf("fUseHistogramJetBkg = %d\n",fUseHistogramJetBkg);
1755   printf("fUseHistogramTracks = %d\n",fUseHistogramTracks);
1756   printf("fUseHistogramJetTracks = %d\n",fUseHistogramJetTracks);
1757
1758 }
1759
1760 //__________________________________________________________________
1761 void AliAnaParticleJetFinderCorrelation::CalculateBkg(TVector3 gamma, TVector3 jet,Double_t vertex[3],Int_t type=2){
1762   //
1763   // calculate background for fragmentation function and fill histograms
1764   // 1. 90 degrees from jet axis in random place = perpendicular cone
1765   // 2. Random cone not belonging to jet cone nor photon cone
1766   // 3. In the middle point from jet and photon momentum vectors
1767   // 4. 90 degrees from photon direction in random place = perpendicular cone 2
1768
1769   //
1770   // implementation of 2 works, 1 and 4 works
1771   //
1772   Double_t gammaPt  = gamma.Pt();
1773   Double_t gammaEta = gamma.Eta();
1774   Double_t gammaPhi = gamma.Phi();
1775   Double_t jetEta   = jet.Eta();
1776   Double_t jetPhi   = jet.Phi();
1777
1778   //refference direction of background
1779   Double_t refEta=0.,refPhi=0.;
1780   Double_t rad = 0,rad2 = 0.;
1781   if(type==1){//perpendicular to jet axis
1782     //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
1783
1784     Double_t xVar;
1785     Double_t yVar;
1786     Double_t newX=0.;
1787     Double_t newY=0.;
1788     Double_t newZ=0.;
1789     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1790     Double_t jx=jet.Px();
1791     Double_t jy=jet.Py();
1792     Double_t jz=jet.Pz();
1793     //if(jz==0)printf("problem\n");
1794     //X axis
1795     Double_t Xx=1.0-vertex[0];
1796     Double_t Xy=-1.0*vertex[1];
1797     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1798     //Y axis
1799     Double_t Yx=jy*Xz-jz*Xy;
1800     Double_t Yy=jz*Xx-jx*Xz;
1801     Double_t Yz=jx*Xy-jy*Xx;
1802     //Determinant
1803     Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
1804     if(det==0)printf("problem det==0\n");
1805     Double_t detX = 0.;
1806     Double_t detY = 0.;
1807     Double_t detZ = 0.;
1808
1809 //    Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
1810 //    printf("scalar jet.P o X: %f\n",tmpScalar);
1811 //    tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
1812 //    printf("scalar jet.P o Y: %f\n",tmpScalar);
1813 //    tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
1814 //    printf("scalar X o Y: %f\n",tmpScalar);
1815
1816     TVector3 perp;
1817     //randomise
1818     do{
1819       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1820       //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
1821       xVar=TMath::Cos(refPhi);
1822       yVar=TMath::Sin(refPhi);
1823       //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
1824       //zVar=0 in new surface frame
1825       detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
1826       detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
1827       detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
1828
1829       newX=detX/det;
1830       newY=detY/det;
1831       newZ=detZ/det;
1832
1833       perp.SetXYZ(newX,newY,newZ);
1834       refEta = perp.Eta();
1835       refPhi = perp.Phi();//output in <-pi, pi> range
1836       if(refPhi<0)refPhi+=2*TMath::Pi();
1837       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1838       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1839       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
1840     } while (rad<fJetConeSize+fGammaConeSize || rad2<fJetConeSize+fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
1841     fhRandomPhiEta[2]->Fill(refPhi,refEta);
1842
1843   }
1844   else if(type==2){//random cone
1845     //randomise
1846     do{
1847       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1848       refEta=fGenerator->Uniform(-(0.9-fJetConeSize),0.9-fJetConeSize);
1849       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1850       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1851       //check if reference is not within jet cone or gamma cone +0.4
1852       //example: FF fConeSize=0.4, fJetConeSize=0.4, fIsoGammaConeSize=0.3
1853     } while (rad<fJetConeSize+fGammaConeSize || rad2<fJetConeSize+fJetConeSize);
1854     //photon:0.7=0.4+0.3; jets:0.8=0.4 +0.4 //rad<fConeSize+fJetConeSize rad2<fConeSize+0.3
1855     fhRandomPhiEta[0]->Fill(refPhi,refEta);
1856   }
1857   else if(type==4){//perpendicular to photon axis
1858     Double_t xVar;
1859     Double_t yVar;
1860     Double_t newX=0.;
1861     Double_t newY=0.;
1862     Double_t newZ=0.;
1863     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1864     Double_t jx=gamma.Px();
1865     Double_t jy=gamma.Py();
1866     Double_t jz=gamma.Pz();
1867     //if(jz==0)printf("problem\n");
1868     //X axis
1869     Double_t Xx=1.0-vertex[0];
1870     Double_t Xy=-1.0*vertex[1];
1871     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1872     //Y axis
1873     Double_t Yx=jy*Xz-jz*Xy;
1874     Double_t Yy=jz*Xx-jx*Xz;
1875     Double_t Yz=jx*Xy-jy*Xx;
1876     //Determinant
1877     Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
1878     if(det==0)printf("problem det==0\n");
1879     Double_t detX = 0.;
1880     Double_t detY = 0.;
1881     Double_t detZ = 0.;
1882
1883 //    Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
1884 //    printf("scalar jet.P o X: %f\n",tmpScalar);
1885 //    tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
1886 //    printf("scalar jet.P o Y: %f\n",tmpScalar);
1887 //    tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
1888 //    printf("scalar X o Y: %f\n",tmpScalar);
1889
1890     TVector3 perp;
1891     //randomise
1892     do{
1893       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1894       //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
1895       xVar=TMath::Cos(refPhi);
1896       yVar=TMath::Sin(refPhi);
1897       //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
1898       //zVar=0 in new surface frame
1899       detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
1900       detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
1901       detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
1902
1903       newX=detX/det;
1904       newY=detY/det;
1905       newZ=detZ/det;
1906
1907       perp.SetXYZ(newX,newY,newZ);
1908       refEta = perp.Eta();
1909       refPhi = perp.Phi();//output in <-pi, pi> range
1910       if(refPhi<0)refPhi+=2*TMath::Pi();
1911       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1912       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1913       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
1914     } while (rad<fJetConeSize+fGammaConeSize || rad2<fJetConeSize+fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
1915     fhRandomPhiEta[1]->Fill(refPhi,refEta);
1916
1917   }
1918   else if(type==3){//mid point
1919
1920     Double_t jx=jet.Px();
1921     Double_t jy=jet.Py();
1922     Double_t jz=jet.Pz();
1923     //    if(jz==0)printf("problem\n");
1924     Double_t gx=gamma.Px();
1925     Double_t gy=gamma.Py();
1926     Double_t gz=gamma.Pz();
1927
1928     Double_t cosAlpha=(jx*gx+jy*gy+jz*gz)/(jet.Mag()*gamma.Mag());
1929     Double_t cosinus=TMath::Sqrt((cosAlpha+1.)/2.);
1930     //perpendicular axis
1931     Double_t Zx=gy*jz-gz*jy;
1932     Double_t Zy=gz*jx-gx*jz;
1933     Double_t Zz=gx*jy-gy*jx;
1934
1935     //Determinant
1936     Double_t det = Zx*gy*jz + Zy*gz*jx + Zz*gx*jy - Zz*gy*jx - Zy*gx*jz - Zx*gz*jy;
1937
1938     Double_t newX=0.;
1939     Double_t newY=0.;
1940     Double_t newZ=0.;
1941     if(det!=0) {
1942       Double_t detX =            -Zy*gz*cosinus +Zz*cosinus*jy + Zz*gy*cosinus - Zy*cosinus*jz;
1943       Double_t detY = Zx*cosinus*jz          - Zz*gx*cosinus - Zz*cosinus*jx            + Zx*gz*cosinus;
1944       Double_t detZ = -Zx*gy*cosinus + Zy*cosinus*jx + Zy*gx*cosinus - Zx*cosinus*jy;
1945
1946       newX=detX/det;
1947       newY=detY/det;
1948       newZ=detZ/det;
1949     }
1950
1951     TVector3 perp;
1952     perp.SetXYZ(newX,newY,newZ);
1953     refEta = perp.Eta();
1954     refPhi = perp.Phi();//output in <-pi, pi> range
1955     if(refPhi<0)refPhi+=2*TMath::Pi();
1956     rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1957     rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1958       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
1959
1960     if (rad<fJetConeSize+fGammaConeSize || rad2<fJetConeSize+fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize) fhRandomPhiEta[3]->Fill(refPhi,refEta);
1961   }
1962   else if(type==5){//tmp
1963     //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
1964
1965     Double_t xVar;
1966     Double_t newX=0.;
1967     Double_t newY=0.;
1968     Double_t newZ=0.;
1969     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1970     Double_t jx=jet.Px();
1971     Double_t jy=jet.Py();
1972     Double_t jz=jet.Pz();
1973     //    if(jz==0)printf("problem\n");
1974     //X axis
1975     Double_t Xx=1.0-vertex[0];
1976     Double_t Xy=-1.0*vertex[1];
1977     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1978     //Y axis
1979     Double_t Yx=jy*Xz-jz*Xy;
1980     Double_t Yy=jz*Xx-jx*Xz;
1981     Double_t Yz=jx*Xy-jy*Xx;
1982
1983     // X and Y length
1984     Double_t Xlength=TMath::Sqrt(Xx*Xx+Xy*Xy+Xz*Xz);
1985     Double_t Ylength=TMath::Sqrt(Yx*Yx+Yy*Yy+Yz*Yz);
1986     Double_t ratio=Ylength/Xlength;
1987
1988     TVector3 perp;
1989     //randomise
1990     do{
1991       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1992       xVar=TMath::Tan(refPhi)/ratio;
1993       newX=xVar*Yx+Xx;
1994       newY=xVar*Yy+Xy;
1995       newZ=xVar*Yz+Xz;
1996
1997       perp.SetXYZ(newX,newY,newZ);
1998       refEta = perp.Eta();
1999       refPhi = perp.Phi();//output in <-pi, pi> range
2000       if(refPhi<0)refPhi+=2*TMath::Pi();
2001       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2002       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2003       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2004     } while (rad<fJetConeSize+fGammaConeSize || rad2<fJetConeSize+fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2005     fhRandomPhiEta[4]->Fill(refPhi,refEta);
2006   }
2007
2008
2009
2010   //calculate FF in background
2011   Int_t ntracks =  0;
2012   ntracks =GetCTSTracks()->GetEntriesFast();
2013   AliVTrack* track = 0x0 ;
2014   TVector3 p3;
2015
2016   Double_t pt = 0, eta = 0, phi = 0;
2017   Int_t  npartcone = 0;
2018   Double_t sumPt=0.;
2019   for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
2020     track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
2021     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
2022     pt   = p3.Pt();
2023     if(pt<fPtThresholdInCone) {//0.150
2024       //printf("problem: track pt < %f MeV/c \n",fPtThresholdInCone);
2025       continue;
2026     }
2027     eta  = p3.Eta() ;
2028     phi  = p3.Phi() ;
2029     if(phi < 0) phi+=TMath::TwoPi();
2030     //Check if there is any particle inside cone with pt larger than  fPtThreshold
2031     rad = TMath::Sqrt((eta-refEta)*(eta-refEta) + (phi-refPhi)*(phi-refPhi));
2032     if(rad < fConeSize  && pt > fPtThresholdInCone){
2033       //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
2034         npartcone++;
2035         sumPt+=pt;
2036         if(type==1){//perp jet
2037           fhBkgFFz[1] ->Fill(gammaPt, pt/gammaPt);
2038           fhBkgFFxi[1]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2039           fhBkgFFpt[1]->Fill(gammaPt, pt);
2040         }
2041         else if(type==2){//RC
2042           fhBkgFFz[0] ->Fill(gammaPt, pt/gammaPt);
2043           fhBkgFFxi[0]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2044           fhBkgFFpt[0]->Fill(gammaPt, pt);
2045         }
2046         else if(type==3){//mid point
2047           fhBkgFFz[3] ->Fill(gammaPt, pt/gammaPt);
2048           fhBkgFFxi[3]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2049           fhBkgFFpt[3]->Fill(gammaPt, pt);
2050         }
2051         else if(type==4){//perp jet
2052           fhBkgFFz[2] ->Fill(gammaPt, pt/gammaPt);
2053           fhBkgFFxi[2]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2054           fhBkgFFpt[2]->Fill(gammaPt, pt);
2055         }
2056         else if(type==5){//test
2057           fhBkgFFz[4] ->Fill(gammaPt, pt/gammaPt);
2058           fhBkgFFxi[4]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2059           fhBkgFFpt[4]->Fill(gammaPt, pt);
2060         }
2061
2062
2063     }
2064   }//end of loop over tracks
2065   Double_t sumOverTracks=0.;
2066   if(npartcone!=0) sumOverTracks = sumPt/npartcone;
2067   if(type==1) {
2068     fhBkgNTracksInCone[1]->Fill(gammaPt, npartcone);
2069     fhBkgSumPtInCone[1]->Fill(gammaPt,sumPt);
2070     fhBkgSumPtnTracksInCone[1]->Fill(gammaPt,sumOverTracks);
2071   }
2072   else if(type==2) {
2073     fhBkgNTracksInCone[0]->Fill(gammaPt, npartcone);
2074     fhBkgSumPtInCone[0]->Fill(gammaPt,sumPt);
2075     fhBkgSumPtnTracksInCone[0]->Fill(gammaPt,sumOverTracks);
2076   }
2077   else if(type==3) {
2078     fhBkgNTracksInCone[3]->Fill(gammaPt, npartcone);
2079     fhBkgSumPtInCone[3]->Fill(gammaPt,sumPt);
2080     fhBkgSumPtnTracksInCone[3]->Fill(gammaPt,sumOverTracks);
2081   }
2082   else if(type==4) {
2083     fhBkgNTracksInCone[2]->Fill(gammaPt, npartcone);
2084     fhBkgSumPtInCone[2]->Fill(gammaPt,sumPt);
2085     fhBkgSumPtnTracksInCone[2]->Fill(gammaPt,sumOverTracks);
2086   }
2087   else if(type==5) {
2088     fhBkgNTracksInCone[4]->Fill(gammaPt, npartcone);
2089     fhBkgSumPtInCone[4]->Fill(gammaPt,sumPt);
2090     fhBkgSumPtnTracksInCone[4]->Fill(gammaPt,sumOverTracks);
2091   }
2092 }
2093
2094