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(!jet) continue;
1028       if(TMath::Abs(jet->Eta()) > (0.9 - fJetConeSize) ) continue;
1029       if(jet->EffectiveAreaCharged()<fJetAreaFraction*TMath::Pi()*fJetConeSize*fJetConeSize) continue;
1030       if(jet->Pt()<fJetMinPt) continue;
1031       ptCorrect = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();
1032       if(ptCorrect > mostEneJetPt){
1033         mostEneJetPt = ptCorrect;
1034         indexMostEneJet = ijet;
1035       }
1036     }
1037     //    printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",indexMostEneJet);
1038
1039     //
1040     // assign jet to photon
1041     //
1042     if(indexMostEneJet>=0 && indexMostEnePhoton>=0){
1043       particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(indexMostEnePhoton));
1044       jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(indexMostEneJet));
1045       particle->SetRefJet(jet);
1046     }
1047   }//end of take most energetic photon and most ene. jet after bkg subtraction
1048
1049   if(fMostOpposite){
1050     //Bool_t isJetFound=kFALSE;//new
1051     //Loop on stored AOD particles, trigger
1052     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1053       AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1054
1055       //Correlate with jets
1056       Int_t ijet = SelectJet(particle,aodRecJets);//input for jets is TClonesArray
1057       if(ijet > -1){
1058         //isJetFound=kTRUE;
1059         if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - Jet with index %d selected \n",ijet);
1060         AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1061         particle->SetRefJet(jet);
1062         //printf("Most opposite found\n");
1063       }
1064     } // input aod loop
1065     //  if(GetReader()->WriteDeltaAODToFile() && isJetFound) WriteJetsToOutputBranch(aodRecJets);
1066   }//end of take most opposite photon and jet after bkg subtraction
1067
1068
1069   if(GetDebug() > 1 ) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
1070 }
1071
1072 //__________________________________________________________________
1073 void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms()
1074 {
1075   //Particle-Jet Correlation Analysis, fill histograms
1076   if(GetDebug() > 3 ) {
1077     printf("I use MakeAnalysisFillHistograms\n");
1078     printf("ntracks before iso %d\n",GetCTSTracks()->GetEntriesFast() );
1079   }
1080
1081   //Get the event, check if there are AODs available, if not, skip this analysis
1082   AliAODEvent * event = NULL;
1083
1084   //printf("GetReader()->GetOutputEvent() %d\n",GetReader()->GetOutputEvent() );
1085   //printf("GetReader()->GetDataType() %d\n",GetReader()->GetDataType() );
1086   //printf("AliCaloTrackReader::kAOD %d\n",AliCaloTrackReader::kAOD );
1087   //printf("GetReader()->GetInputEvent() %d\n",GetReader()->GetInputEvent() );
1088
1089   if(GetReader()->GetOutputEvent())
1090   {
1091     event = dynamic_cast<AliAODEvent*>(GetReader()->GetOutputEvent());
1092   }
1093   else if(GetReader()->GetDataType() == AliCaloTrackReader::kAOD)
1094   {
1095     event = dynamic_cast<AliAODEvent*>(GetReader()->GetInputEvent());
1096   }
1097   else {
1098     if(GetDebug() > 3) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - There are no jets available for this analysis\n");
1099     return;
1100   }
1101
1102   if(!GetInputAODBranch() || !event){
1103     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s > \n",GetInputAODName().Data());
1104     abort();
1105   }
1106
1107   Int_t nJets=-1;
1108   TClonesArray *aodRecJets = 0;
1109   if(IsNonStandardJetFromReader()){//branch read in reader from reader
1110     if (GetDebug() > 3) printf("GetNonStandardJets function (from reader) is called\n");
1111     aodRecJets = GetNonStandardJets();
1112     if(aodRecJets==0x0){
1113       if(GetDebug() > 3) event->Print();
1114       AliFatal("Jets container not found");
1115       return; // trick coverity
1116     }
1117     nJets=aodRecJets->GetEntries();
1118   }
1119   if(nJets==0) {
1120     //    printf("Why number of jets = 0? Check what type of collision it is. If PbPb -problem.\n");
1121     GetReader()->FillInputNonStandardJets();
1122     aodRecJets = GetNonStandardJets();
1123     if(aodRecJets) nJets=aodRecJets->GetEntries();
1124     //    printf("nJets = %d\n",nJets);
1125     return;
1126   }
1127
1128   //
1129   //Background jets
1130   //
1131   AliAODJetEventBackground* aodBkgJets = 0;
1132   if(IsBackgroundJetFromReader()){//jet branch from reader
1133     if(GetDebug() > 3) printf("GetBackgroundJets function is called\n");
1134     aodBkgJets = GetBackgroundJets();
1135     if(GetDebug() > 3) printf("aodBkgJets %p\n",aodBkgJets);
1136     if(aodBkgJets==0x0){
1137       if(GetDebug() > 3) event->Print();
1138       abort();
1139     }
1140     if(GetDebug() > 3) aodBkgJets->Print("c");
1141   }
1142
1143
1144   //
1145   // only background jets informations
1146   //
1147   //  Float_t pTback = 0;
1148   Double_t rhoEvent=0.;
1149   if(aodBkgJets) {
1150     if(IsBackgroundJetFromReader() ) rhoEvent = aodBkgJets->GetBackground(2);
1151     if(IsHistogramJetBkg()) {
1152       fhJetRhoVsCentrality->Fill(GetEventCentrality(),rhoEvent);
1153       for(Int_t i=0;i<4;i++){
1154         fhBkgJetBackground[i]->Fill(aodBkgJets->GetBackground(i));
1155         fhBkgJetSigma[i]->Fill(aodBkgJets->GetSigma(i));
1156         fhBkgJetArea[i]->Fill(aodBkgJets->GetMeanarea(i));
1157       }
1158     }//end of if fill HistogramJetBkg
1159   }//end if aodBkgJets exists
1160
1161   //
1162   //only track information
1163   //
1164   Int_t nCTSTracks = GetCTSTracks()->GetEntriesFast();
1165   AliAODTrack *aodtrack;
1166   Int_t itrack = 0;
1167   if(IsHistogramTracks()) {
1168     Double_t sumTrackPt=0;
1169     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1170       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1171       if(!aodtrack) continue;
1172       fhTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1173       sumTrackPt+=aodtrack->Pt();
1174     }
1175     if(nCTSTracks)
1176       fhTrackAveTrackPt->Fill(sumTrackPt/nCTSTracks);
1177   }//end if IsHistogramTracks
1178
1179   //
1180   //only jet informations
1181   //
1182   AliAODJet * jettmp = 0 ;
1183   Double_t jetPttmp = 0.;
1184   Double_t var1tmp = 0.;
1185   Double_t var2tmp = 0.;
1186   //  fhJetNjet->Fill(nJets);
1187   Double_t ptMostEne=0;
1188   //  Int_t indexMostEne=-1;
1189   Int_t nJetsOverThreshold[10]={nJets,0,0,0,0,0,0,0,0,0};
1190   Int_t iCounter=0;
1191   Double_t sumJetTrackPt=0.;
1192   Int_t sumNJetTrack=0;
1193   Int_t nTracksInJet=0;
1194   Int_t itrk=0;
1195   for(Int_t ijet = 0; ijet < nJets ; ijet++){
1196     jettmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ijet));
1197     if(!jettmp) continue;
1198     fhJetPtBefore->Fill(jettmp->Pt());
1199     jetPttmp  = jettmp->Pt() - rhoEvent * jettmp->EffectiveAreaCharged();//<<---changed here
1200
1201     //calculate number of tracks above 1,2,3,4,5 GeV in jet
1202     AliVTrack* jettrack = 0x0 ;
1203
1204     Int_t nTrackThrGeV[5]={0,0,0,0,0};
1205     nTracksInJet=(jettmp->GetRefTracks())->GetEntriesFast();
1206     fhJetNparticlesInJet->Fill(nTracksInJet);
1207     if(nTracksInJet==0) continue;
1208     sumNJetTrack+=nTracksInJet;
1209     for(itrack=0;itrack<nTracksInJet;itrack++){
1210       jettrack=(AliVTrack *) ((jettmp->GetRefTracks())->At(itrack));
1211       if(!jettrack) continue;
1212
1213       fhJetDeltaEtaDeltaPhi->Fill(jettmp->Eta()-jettrack->Eta(),jettmp->Phi()-jettrack->Phi());
1214       sumJetTrackPt+=jettrack->Pt();
1215       if(IsHistogramJetTracks()){
1216         if(jettrack->Pt()>1.) nTrackThrGeV[0]++;
1217         if(jettrack->Pt()>2.) nTrackThrGeV[1]++;
1218         if(jettrack->Pt()>3.) nTrackThrGeV[2]++;
1219         if(jettrack->Pt()>4.) nTrackThrGeV[3]++;
1220         if(jettrack->Pt()>5.) nTrackThrGeV[4]++;
1221       }//end of if IsHistogramJetTracks
1222     }//end of jet track loop
1223     //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]);
1224
1225     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1226       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1227       if(aodtrack) fhJetDeltaEtaDeltaPhiAllTracks->Fill(jettmp->Eta()-aodtrack->Eta(),jettmp->Phi()-aodtrack->Phi());
1228     }
1229
1230
1231     if(IsHistogramJetTracks()){
1232       fhJetNtracksInJetAboveThr[0]->Fill(jetPttmp,nTracksInJet);//all jets
1233
1234       for(itrk=0;itrk<5;itrk++) {
1235         fhJetNtracksInJetAboveThr[itrk+1]->Fill(jetPttmp,nTrackThrGeV[itrk]);//all jets
1236         fhJetRatioNTrkAboveToNTrk[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);//all jets
1237       }
1238       if(ijet==0){//most ene jet
1239         for(itrk=0;itrk<5;itrk++)
1240           fhJetNtrackRatioMostEne[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1241       }
1242       if(jetPttmp>5){//jet with pt>5GeV/c
1243         for(itrk=0;itrk<5;itrk++)
1244           fhJetNtrackRatioJet5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1245       }
1246       if(nTrackThrGeV[4]>0){//jet with leading particle pt>5GeV/c
1247         for(itrk=0;itrk<5;itrk++)
1248           fhJetNtrackRatioLead5GeV[itrk]->Fill(jetPttmp,(Double_t)nTrackThrGeV[itrk]/(Double_t)nTracksInJet);
1249       }
1250     }//end of if IsHistogramJetTracks
1251
1252
1253     Double_t effectiveChargedBgEnergy=(IsBackgroundJetFromReader()?rhoEvent * jettmp->EffectiveAreaCharged():jettmp->ChargedBgEnergy());
1254
1255
1256     fhJetChBkgEnergyVsArea->Fill(effectiveChargedBgEnergy,jettmp->EffectiveAreaCharged());
1257     //if(jettmp->EffectiveAreaCharged()>0)
1258     fhJetRhoVsPt->Fill(jetPttmp,jettmp->ChargedBgEnergy()*jettmp->EffectiveAreaCharged());
1259
1260     if(jetPttmp>ptMostEne) {
1261       ptMostEne = jetPttmp;
1262       //indexMostEne=ijet;
1263     }
1264     fhJetPt->Fill(jetPttmp);
1265     fhJetChBkgEnergyVsPt->Fill(jetPttmp,effectiveChargedBgEnergy);
1266     fhJetChAreaVsPt->Fill(jetPttmp,jettmp->EffectiveAreaCharged());
1267     if(GetDebug()>5) printf("ChargedBgEnergy %f EffectiveAreaCharged %f\n", jettmp->ChargedBgEnergy(),jettmp->EffectiveAreaCharged());
1268     for(iCounter=1;iCounter<10;iCounter++){
1269       if(jetPttmp>iCounter) nJetsOverThreshold[iCounter]++;
1270     }
1271
1272     var1tmp  = jettmp->Phi();
1273     var2tmp  = jettmp->Eta();
1274     fhJetPhi->Fill(var1tmp);
1275     fhJetEta->Fill(var2tmp);
1276     fhJetPhiVsEta->Fill(var1tmp,var2tmp);
1277     fhJetEtaVsPt->Fill(jetPttmp,var2tmp);
1278     fhJetEtaVsNpartInJet->Fill(var2tmp,nTracksInJet);
1279     if(jetPttmp>0)
1280       fhJetEtaVsNpartInJetBkg->Fill(var2tmp,nTracksInJet);
1281
1282   }//end of jet loop
1283   if(IsHistogramJetTracks()){
1284     if(sumNJetTrack>0){
1285       //printf("average track pt %f\n",sumJetTrackPt/sumNJetTrack);
1286       fhJetAveTrackPt->Fill(sumJetTrackPt/sumNJetTrack);
1287     }
1288   }//end of if IsHistogramJetTracks
1289
1290
1291   fhJetPtMostEne->Fill(ptMostEne);
1292   for(iCounter=0;iCounter<10;iCounter++){
1293     fhJetNjetOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1294     nJetsOverThreshold[iCounter]=0;
1295   }
1296
1297   //end of only jet part
1298
1299   //
1300   // Photons
1301   //
1302   Int_t ntrig   =  GetInputAODBranch()->GetEntriesFast() ;
1303   if(GetDebug() > 1){
1304     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Begin jet finder  correlation analysis, fill histograms \n");
1305     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", ntrig);
1306     printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - In jet output branch aod entries %d\n", event->GetNJets());
1307   }
1308   fhNjetsNgammas->Fill(nJets,ntrig);
1309
1310   //if(nJets==0)   return;//to speed up
1311   //printf("ntrig %d, njets %d\n",ntrig,nJets);
1312
1313   //
1314   //Fill only temporary photon histograms
1315   //
1316   Double_t maxPt=0.;
1317   Int_t maxIndex=-1;
1318   Double_t tmpPt=0.;
1319   Double_t sumPt=0.;
1320   nJetsOverThreshold[0]=ntrig;
1321   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1322     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1323     tmpPt = particlecorr->Pt();
1324     if(tmpPt>maxPt) {
1325       maxPt = tmpPt;
1326       maxIndex = iaod;
1327     }
1328     for(iCounter=1;iCounter<10;iCounter++){
1329       if(tmpPt>iCounter) nJetsOverThreshold[iCounter]++;
1330     }
1331     sumPt+=tmpPt;
1332   }
1333   for(iCounter=0;iCounter<10;iCounter++){
1334     fhPhotonNgammaOverPtCut[iCounter]->Fill(nJetsOverThreshold[iCounter]);
1335     //    nJetsOverThreshold[iCounter]=0;
1336   }
1337
1338   fGamAvEne=0;
1339   TObjArray* clusters = GetEMCALClusters();
1340   //printf("calculate median bkg energy for photons ");
1341   Double_t medianPhotonRho=0.;
1342   Int_t clusterID;
1343   Int_t iclustmp = -1;
1344   Int_t numberOfcells=0;
1345   AliVCluster *cluster = 0;
1346   if(ntrig>1){
1347     Double_t *photonRhoArr=new Double_t[ntrig-1];
1348     fhPhotonPtMostEne->Fill(maxPt);
1349     //    fhPhotonIndexMostEne->Fill(indexMaxPt);
1350     fhPhotonAverageEnergy->Fill(sumPt/ntrig);
1351     fhPhotonRatioAveEneToMostEne->Fill(sumPt/(ntrig*maxPt));
1352     fhPhotonAverageEnergyMinus1->Fill((sumPt-maxPt)/(ntrig-1));
1353     fGamAvEne=(sumPt-maxPt)/(ntrig-1);
1354     fhPhotonRatioAveEneMinus1ToMostEne->Fill((sumPt-maxPt)/((ntrig-1)*maxPt));
1355
1356     Int_t counterGamma=0;
1357     Int_t counterGammaMinus1=0;
1358
1359     Int_t photonRhoArrayIndex=0;
1360     //TObjArray* clusterstmp = GetEMCALClusters();
1361     for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1362       AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1363       if( particlecorr->Pt() > sumPt/ntrig ) counterGamma++;
1364       if( particlecorr->Pt() > (sumPt-maxPt)/(ntrig-1) ) counterGammaMinus1++;
1365
1366       if(iaod==maxIndex) continue;
1367       clusterID = particlecorr->GetCaloLabel(0) ;
1368       if(clusterID < 0) continue;
1369       cluster = FindCluster(clusters,clusterID,iclustmp);
1370       photonRhoArr[photonRhoArrayIndex]=particlecorr->Pt()/ cluster->GetNCells();
1371       numberOfcells+=cluster->GetNCells();
1372       photonRhoArrayIndex++;
1373     }
1374     if(photonRhoArrayIndex>0) medianPhotonRho=TMath::Median(photonRhoArrayIndex,photonRhoArr);
1375     delete [] photonRhoArr;
1376     fhPhotonNgammaMoreAverageToNgamma->Fill((Double_t)counterGamma / (Double_t)ntrig);
1377     fhPhotonNgammaMoreAverageMinus1ToNgamma->Fill((Double_t)counterGammaMinus1 / (Double_t)ntrig);
1378   }
1379   //printf("median = %f\n",medianPhotonRho);
1380   fhPhotonBkgRhoVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),medianPhotonRho);
1381   fhPhotonBkgRhoVsNclusters->Fill(ntrig,medianPhotonRho);
1382   fhPhotonBkgRhoVsCentrality->Fill(GetEventCentrality(),medianPhotonRho);
1383   fhPhotonBkgRhoVsNcells->Fill(numberOfcells,medianPhotonRho);
1384
1385
1386   AliVCluster *cluster2 = 0;
1387   Double_t photon2Corrected=0;
1388   Double_t sumPtTmp=0.;
1389   Double_t sumPtCorrectTmp=0.;
1390   AliVTrack* trackTmp = 0x0 ;
1391   TVector3 p3Tmp;
1392
1393   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1394     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1395     clusterID = particlecorr->GetCaloLabel(0) ;
1396     if(clusterID < 0) continue;
1397     cluster = FindCluster(clusters,clusterID,iclustmp);
1398     fhPhotonPt->Fill(particlecorr->Pt());
1399     fhPhotonPtCorrected->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1400     fhPhotonPtDiff->Fill(cluster->GetNCells() * medianPhotonRho);
1401     fhPhotonPtDiffVsCentrality->Fill(GetEventCentrality(),cluster->GetNCells() * medianPhotonRho);
1402     fhPhotonPtDiffVsNcells->Fill(numberOfcells,cluster->GetNCells() * medianPhotonRho);
1403     fhPhotonPtDiffVsNtracks->Fill(GetCTSTracks()->GetEntriesFast(),cluster->GetNCells() * medianPhotonRho);
1404     fhPhotonPtDiffVsNclusters->Fill(ntrig,cluster->GetNCells() * medianPhotonRho);
1405
1406     fhPhotonPtCorrectedZoom->Fill(particlecorr->Pt() - cluster->GetNCells() * medianPhotonRho);
1407
1408     //test: sum_pt in the cone 0.3 for each photon
1409     //should be: random fake gamma from MB
1410     //is: each gamma for EMCEGA
1411     sumPtTmp=0.;
1412     sumPtCorrectTmp=0.;
1413
1414     for(Int_t iaod2 = 0; iaod2 < ntrig ; iaod2++){
1415       if(iaod==iaod2) continue;
1416       AliAODPWG4ParticleCorrelation* particlecorr2 =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod2));
1417       clusterID = particlecorr2->GetCaloLabel(0) ;
1418       if(clusterID < 0) continue;
1419       cluster2 = FindCluster(clusters,clusterID,iclustmp);
1420       photon2Corrected = particlecorr2->Pt() - cluster2->GetNCells() * medianPhotonRho;
1421
1422       //if(Pt()<0.5) continue; //<<hardcoded here //FIXME
1423       if( TMath::Sqrt((particlecorr->Eta()-particlecorr2->Eta())*(particlecorr->Eta()-particlecorr2->Eta()) +
1424                       (particlecorr->Phi()-particlecorr2->Phi())*(particlecorr->Phi()-particlecorr2->Phi()) )<fGammaConeSize ){//if(/*cone is correct*/){
1425         sumPtTmp+= particlecorr2->Pt();
1426         sumPtCorrectTmp+=photon2Corrected;
1427       }
1428     }
1429     fhPhotonSumPtInCone->Fill(sumPtTmp);
1430     fhPhotonSumPtCorrectInCone->Fill(sumPtCorrectTmp);
1431
1432     //test: sum_pt in the cone 0.3 for each track
1433     //should be: random fake gamma from MB
1434     //is: each gamma for EMCEGA
1435     sumPtTmp=0.;
1436     for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++){
1437       trackTmp = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1438       p3Tmp.SetXYZ(trackTmp->Px(),trackTmp->Py(),trackTmp->Pz());
1439       if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1440                       (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1441         sumPtTmp+=p3Tmp.Pt();
1442       }
1443     }//end of loop over tracks
1444     fhPhotonSumPtChargedInCone->Fill(sumPtTmp);
1445   }
1446
1447   //End of Fill temporary photon histograms
1448
1449   //
1450   // Apply background subtraction for photons
1451   //
1452   fGamRho = medianPhotonRho;
1453   if(!IsBackgroundSubtractionGamma()) medianPhotonRho=0;
1454
1455
1456   //
1457   //Get vertex for cluster momentum calculation <<----new here
1458   //
1459   Double_t vertex[] = {0,0,0} ; //vertex ;
1460   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1461     GetReader()->GetVertex(vertex);
1462   fZvertex = vertex[2];
1463
1464   //
1465   //Loop on stored AOD particles, trigger
1466   //
1467   for(Int_t iaod = 0; iaod < ntrig ; iaod++){
1468     AliAODPWG4ParticleCorrelation* particlecorr =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1469     fhCuts->Fill(0);
1470     fhCuts2->Fill(0.,(Double_t)nJets);
1471     if(GetDebug() > 5) printf("OnlyIsolated %d  !particlecorr->IsIsolated() %d \n",OnlyIsolated(), !particlecorr->IsIsolated());
1472
1473     if(OnlyIsolated() && !particlecorr->IsIsolated()) continue;
1474     fhCuts->Fill(1);
1475     fhCuts2->Fill(1.,nJets);
1476
1477     if(nJets>0) {
1478       fhCuts->Fill(2);
1479     }
1480
1481     //Recover the jet correlated, found previously.
1482     AliAODJet   * jet = particlecorr->GetJet();
1483     //If correlation not made before, do it now.
1484     if(fMakeCorrelationInHistoMaker){
1485       //Correlate with jets
1486       Int_t ijet = SelectJet(particlecorr,aodRecJets);//input for jets is TClonesArray
1487       if(ijet > -1){
1488         if(GetDebug() > 2) printf ("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - Jet with index %d selected \n",ijet);
1489         //jet = event->GetJet(ijet);
1490         jet = dynamic_cast<AliAODJet*>(aodRecJets-> At(ijet));
1491
1492         particlecorr->SetRefJet(jet);
1493
1494       }
1495     }
1496
1497     if (!jet) continue ;
1498     fhCuts->Fill(3);
1499     fhCuts2->Fill(7.,1.);
1500
1501     //
1502     //Fill Correlation Histograms
1503     //
1504     clusterID = particlecorr->GetCaloLabel(0) ;
1505     if(!(clusterID<0)){
1506       cluster = FindCluster(clusters,clusterID,iclustmp);
1507       //fill tree variables
1508       fGamNcells = cluster->GetNCells();
1509     }
1510     Double_t ptTrig = particlecorr->Pt() - medianPhotonRho * fGamNcells;//<<---changed here
1511     Double_t ptJet = jet->Pt() - rhoEvent * jet->EffectiveAreaCharged();//<<---changed here
1512     Double_t phiTrig = particlecorr->Phi();
1513     Double_t phiJet = jet->Phi();
1514     Double_t etaTrig = particlecorr->Eta();
1515     Double_t etaJet = jet->Eta();
1516     Double_t deltaPhi=phiTrig-phiJet;
1517     if(deltaPhi<0)deltaPhi+=(TMath::Pi()*2.);
1518     //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",
1519     //  ptTrig,ptJet, phiJet-phiTrig, etaJet-etaTrig, ptTrig-ptJet, ptJet/ptTrig);
1520     fhDeltaPt ->Fill(ptTrig, ptTrig-ptJet);
1521     //    fhDeltaPhi->Fill(ptTrig, phiTrig-phiJet);//need to be shifted by 2pi
1522
1523     fhDeltaPhiCorrect->Fill(ptTrig, deltaPhi);//correct
1524
1525     Double_t deltaPhiCorrect = TMath::Abs( particlecorr->Phi() - jet->Phi() );
1526     if ( deltaPhiCorrect > TMath::Pi() ) deltaPhiCorrect = 2. * TMath::Pi() - deltaPhiCorrect ;
1527     fhDeltaPhi0PiCorrect->Fill(ptTrig, deltaPhiCorrect);
1528
1529     fhDeltaEta->Fill(ptTrig, etaTrig-etaJet);
1530     fhPtRatio ->Fill(ptTrig, ptJet/ptTrig);
1531     fhPt      ->Fill(ptTrig, ptJet);
1532
1533     fhSelectedJetPhiVsEta->Fill(phiJet,etaJet);
1534     fhSelectedJetChBkgEnergyVsPtJet->Fill(ptJet,(IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy()) );
1535     fhSelectedJetChAreaVsPtJet->Fill(ptJet,jet->EffectiveAreaCharged());
1536     fhSelectedJetNjet->Fill(nJets);
1537     fhSelectedNtracks->Fill(GetCTSTracks()->GetEntriesFast());//to be checked
1538     fhSelectedPhotonNLMVsPt->Fill(ptTrig,particlecorr->GetFiducialArea());
1539
1540
1541     if(clusterID < 0 ){
1542       fhSelectedPhotonLambda0VsPt->Fill(ptTrig,-1);
1543       //fill tree variables
1544       fGamLambda0  = -1;
1545       fGamTime = -1;
1546       fGamNcells = 0;
1547       fGamSumPtNeu=0;
1548     }
1549     else{
1550       //Int_t iclus = -1;
1551       //      TObjArray* clusters = GetEMCALClusters();
1552       //cluster = FindCluster(clusters,clusterID,iclustmp);
1553       Double_t lambda0=cluster->GetM02();
1554       fhSelectedPhotonLambda0VsPt->Fill(ptTrig,lambda0);
1555       //fill tree variables
1556       fGamLambda0  = lambda0;
1557       fGamTime = cluster->GetTOF();
1558       //fGamNcells = cluster->GetNCells();
1559
1560       fGamSumPtNeu=0;
1561       fGamNclusters=0;
1562       TLorentzVector mom ;
1563       //TVector3 p3Tmp;
1564       //Double_t scalarProduct=0;
1565       //Double_t vectorLength=particlecorr->P();
1566       for(Int_t icalo=0; icalo <clusters->GetEntriesFast(); icalo++){
1567         AliVCluster* calo = (AliVCluster *) clusters->At(icalo);
1568         if(clusterID==calo->GetID()) continue;//the same cluster as trigger
1569         calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
1570         //printf("min pt %f\n",GetMinPt());
1571         if(mom.Pt()<GetMinPt()) continue; //<<hardcoded here //FIXME 0.5 check if correct
1572         p3Tmp.SetXYZ(mom.Px(),mom.Py(),mom.Pz());
1573         //calculate sum pt in the cone
1574         if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
1575                         (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
1576           //scalarProduct = particlecorr->Px()*mom.Px() + particlecorr->Py()*mom.Py() + particlecorr->Pz()*mom.Pz();
1577           //scalarProduct/=mom.P();
1578           //scalarProduct/=vectorLength;
1579           //if(scalarProduct>TMath::Cos(0.3)) {//FIXME photon radius
1580           fGamSumPtNeu+=mom.Pt();
1581           fGamNclusters++;
1582         }
1583       }
1584     }
1585
1586     //sum pt of charged tracks in the gamma isolation cone
1587     //starts here
1588     fGamSumPtCh=0;
1589     fGamNtracks=0;
1590     for(itrack = 0; itrack < nCTSTracks ; itrack++){
1591       aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1592       fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());//fill histogram here
1593       //      if(aodtrack->Pt()<0.15) continue;//hardcoded
1594       if(aodtrack->Pt()<fPtThresholdInCone) continue;
1595       if(!aodtrack->IsHybridGlobalConstrainedGlobal()) continue;
1596       if(TMath::Sqrt((particlecorr->Phi() - aodtrack->Phi())*(particlecorr->Phi() - aodtrack->Phi()) +
1597                      (particlecorr->Eta() - aodtrack->Eta())*(particlecorr->Eta() - aodtrack->Eta()) ) <fGammaConeSize ) {
1598         fGamSumPtCh+=aodtrack->Pt();
1599         fGamNtracks++;
1600       }
1601     }
1602     //ends here
1603
1604     //    for(Int_t itrack = 0; itrack < nCTSTracks ; itrack++){
1605     //      aodtrack = dynamic_cast <AliAODTrack*>(GetCTSTracks()->At(itrack));
1606     //      fhSelectedTrackPhiVsEta->Fill(aodtrack->Phi(),aodtrack->Eta());
1607     //    }
1608
1609     //
1610     // Background Fragmentation function
1611     //
1612     TVector3 gammaVector,jetVector;
1613     gammaVector.SetXYZ(particlecorr->Px(),particlecorr->Py(),particlecorr->Pz());
1614     jetVector.SetXYZ(jet->Px(),jet->Py(),jet->Pz());
1615     CalculateBkg(gammaVector,jetVector,vertex,1);//jet perp
1616     CalculateBkg(gammaVector,jetVector,vertex,2);//RC
1617     CalculateBkg(gammaVector,jetVector,vertex,3);//mid point
1618     CalculateBkg(gammaVector,jetVector,vertex,4);//gamma perp
1619     //CalculateBkg(gammaVector,jetVector,vertex,5);/test
1620     Double_t angleJetGam = gammaVector.Angle(jetVector);
1621     //printf("angleJetGam %f\n",angleJetGam*180/TMath::Pi());
1622
1623     //
1624     // Fragmentation function
1625     //
1626     Float_t      rad = 0, pt = 0, eta = 0, phi = 0;
1627     Int_t        npartcone = 0;
1628     TVector3 p3;
1629
1630     Int_t ntracks =  0;
1631     if(GetDebug()>3){
1632       printf ("fUseJetRefTracks %d\n",fUseJetRefTracks );
1633       printf ("jet->GetRefTracks() %p\n",jet->GetRefTracks());
1634       printf ("GetCTSTracks() %p\n",GetCTSTracks() );
1635     }
1636
1637     if(!fUseJetRefTracks)
1638       ntracks =GetCTSTracks()->GetEntriesFast();
1639     else //If you want to use jet tracks from JETAN
1640       ntracks =  (jet->GetRefTracks())->GetEntriesFast();
1641
1642     if(GetDebug()>3)    printf ("ntracks %d\n",ntracks);
1643     AliVTrack* track = 0x0 ;
1644     for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
1645       if(!fUseJetRefTracks)
1646         track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1647       else //If you want to use jet tracks from JETAN
1648         track = (AliVTrack *) ((jet->GetRefTracks())->At(ipr));
1649
1650       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1651       pt    = p3.Pt();
1652       eta  = p3.Eta();
1653       phi  = p3.Phi() ;
1654       if(phi < 0) phi+=TMath::TwoPi();
1655
1656       //Check if there is any particle inside cone with pt larger than  fPtThreshold
1657       rad = TMath::Sqrt((eta-etaJet)*(eta-etaJet)+ (phi-phiJet)*(phi-phiJet));
1658       if(rad < fConeSize  && pt > fPtThresholdInCone){
1659         //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
1660         npartcone++;
1661         fhFFz ->Fill(ptTrig, pt/ptTrig);
1662         fhFFxi->Fill(ptTrig, TMath::Log(ptTrig/pt));
1663         fhFFpt->Fill(ptTrig, pt);
1664
1665         //according to jet axis
1666         fhJetFFz ->Fill(ptJet, pt/ptJet);
1667         fhJetFFxi->Fill(ptJet, TMath::Log(ptJet/pt));
1668         fhJetFFpt->Fill(ptJet, pt);
1669
1670
1671         if(TMath::Cos(angleJetGam)<0 && ptJet!=0 && pt!=0 ){
1672           fhJetFFzCor ->Fill(ptJet, -pt*TMath::Cos(angleJetGam)/ptJet);
1673           fhJetFFxiCor->Fill(ptJet, TMath::Log(ptJet/(-pt*TMath::Cos(angleJetGam))));
1674         }
1675       }
1676     }//Tracks
1677     fhNTracksInCone->Fill(ptTrig, npartcone);
1678     //fill tree here for each photon-jet (isolation required usually)
1679
1680     fGamPt      = ptTrig;
1681     //fGamLambda0  = ;//filled earlier
1682     fGamNLM      = particlecorr->GetFiducialArea();
1683     //fGamSumPtCh  = ;//filled earlier
1684     //fGamTime     = particlecorr->GetTOF();//filled earlier
1685     //fGamNcells   = particlecorr->GetNCells();//filled earlier
1686     fGamEta      = etaTrig;
1687     fGamPhi      = phiTrig;
1688     //fGamSumPtNeu = ;//filled earlier
1689     //fGamNtracks  = ;//filled earlier
1690     //fGamNclusters= ;//filled earlier
1691     //fGamAvEne    = ;//filled earlier
1692     fJetPhi      = phiJet;
1693     fJetEta      = etaJet;
1694     fJetPt       = ptJet;
1695     fJetBkgChEne = (IsBackgroundJetFromReader()?rhoEvent * jet->EffectiveAreaCharged():jet->ChargedBgEnergy());
1696     fJetArea     = jet->EffectiveAreaCharged();
1697     fJetNtracks  = (jet->GetRefTracks())->GetEntriesFast();
1698     fEventNumber = 0;
1699     fNtracks     = GetCTSTracks()->GetEntriesFast();
1700     fCentrality  = GetEventCentrality();
1701     fIso         = particlecorr->IsIsolated();
1702
1703     Int_t nTrk1GeV=0;
1704     Int_t nTrk2GeV=0;
1705     for(itrack=0;itrack < fJetNtracks;itrack++){
1706       track = (AliVTrack *) ((jet->GetRefTracks())->At(itrack));
1707       if(track->Pt()>1.) nTrk1GeV++;
1708       if(track->Pt()>2.) nTrk2GeV++;
1709     }
1710
1711     fJetNtracks1 = nTrk1GeV;
1712     fJetNtracks2 = nTrk2GeV;
1713
1714     if(fSaveGJTree) fTreeGJ->Fill();
1715   }//AOD trigger particle loop
1716   if(GetDebug() > 1) printf("AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
1717 }
1718
1719
1720 //__________________________________________________________________
1721 void AliAnaParticleJetFinderCorrelation::Print(const Option_t * opt) const
1722 {
1723
1724   //Print some relevant parameters set for the analysis
1725   if(! opt)
1726     return;
1727
1728   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1729   AliAnaCaloTrackCorrBaseClass::Print(" ");
1730
1731   printf("Phi trigger-jet        <     %3.2f\n", fDeltaPhiMaxCut) ;
1732   printf("Phi trigger-jet        >     %3.2f\n", fDeltaPhiMinCut) ;
1733   printf("pT Ratio trigger/jet   <     %3.2f\n", fRatioMaxCut) ;
1734   printf("pT Ratio trigger/jet   >     %3.2f\n", fRatioMinCut) ;
1735   printf("fConeSize              =     %3.2f\n", fConeSize) ;
1736   printf("fPtThresholdInCone     =     %3.2f\n", fPtThresholdInCone) ;
1737   printf("fUseJetRefTracks         =     %d\n",    fUseJetRefTracks) ;
1738   printf("fMakeCorrelationInHistoMaker     =     %d\n",    fMakeCorrelationInHistoMaker) ;
1739   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
1740   printf("Reconstructed jet cone size = %3.2f\n", fJetConeSize) ;
1741   printf("Reconstructed jet minimum pt = %3.2f\n", fJetMinPt) ;
1742   printf("Reconstructed jet minimum area fraction = %3.2f\n", fJetAreaFraction) ;
1743
1744   if(!fNonStandardJetFromReader){
1745     printf("fJetBranchName =   %s\n", fJetBranchName.Data()) ;
1746   }
1747   if(!fBackgroundJetFromReader){
1748     printf("fBkgJetBranchName =   %s\n", fBkgJetBranchName.Data()) ;
1749   }
1750
1751   printf("Isolation cone size = %3.2f\n", fGammaConeSize) ;
1752   printf("fUseBackgroundSubtractionGamma = %d\n",fUseBackgroundSubtractionGamma);
1753   printf("fSaveGJTree = %d\n",fSaveGJTree);
1754   printf("fMostEnergetic = %d\n",fMostEnergetic);
1755   printf("fMostOpposite = %d\n",fMostOpposite);
1756
1757   printf("fUseHistogramJetBkg = %d\n",fUseHistogramJetBkg);
1758   printf("fUseHistogramTracks = %d\n",fUseHistogramTracks);
1759   printf("fUseHistogramJetTracks = %d\n",fUseHistogramJetTracks);
1760
1761 }
1762
1763 //__________________________________________________________________
1764 void AliAnaParticleJetFinderCorrelation::CalculateBkg(TVector3 gamma, TVector3 jet,Double_t vertex[3],Int_t type=2){
1765   //
1766   // calculate background for fragmentation function and fill histograms
1767   // 1. 90 degrees from jet axis in random place = perpendicular cone
1768   // 2. Random cone not belonging to jet cone nor photon cone
1769   // 3. In the middle point from jet and photon momentum vectors
1770   // 4. 90 degrees from photon direction in random place = perpendicular cone 2
1771
1772   //
1773   // implementation of 2 works, 1 and 4 works
1774   //
1775   Double_t gammaPt  = gamma.Pt();
1776   Double_t gammaEta = gamma.Eta();
1777   Double_t gammaPhi = gamma.Phi();
1778   Double_t jetEta   = jet.Eta();
1779   Double_t jetPhi   = jet.Phi();
1780
1781   //refference direction of background
1782   Double_t refEta=0.,refPhi=0.;
1783   Double_t rad = 0,rad2 = 0.;
1784   if(type==1){//perpendicular to jet axis
1785     //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
1786
1787     Double_t xVar;
1788     Double_t yVar;
1789     Double_t newX=0.;
1790     Double_t newY=0.;
1791     Double_t newZ=0.;
1792     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1793     Double_t jx=jet.Px();
1794     Double_t jy=jet.Py();
1795     Double_t jz=jet.Pz();
1796     //if(jz==0)printf("problem\n");
1797     //X axis
1798     Double_t Xx=1.0-vertex[0];
1799     Double_t Xy=-1.0*vertex[1];
1800     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1801     //Y axis
1802     Double_t Yx=jy*Xz-jz*Xy;
1803     Double_t Yy=jz*Xx-jx*Xz;
1804     Double_t Yz=jx*Xy-jy*Xx;
1805     //Determinant
1806     Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
1807     if(det==0)printf("problem det==0\n");
1808     Double_t detX = 0.;
1809     Double_t detY = 0.;
1810     Double_t detZ = 0.;
1811
1812 //    Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
1813 //    printf("scalar jet.P o X: %f\n",tmpScalar);
1814 //    tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
1815 //    printf("scalar jet.P o Y: %f\n",tmpScalar);
1816 //    tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
1817 //    printf("scalar X o Y: %f\n",tmpScalar);
1818
1819     TVector3 perp;
1820     //randomise
1821     do{
1822       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1823       //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
1824       xVar=TMath::Cos(refPhi);
1825       yVar=TMath::Sin(refPhi);
1826       //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
1827       //zVar=0 in new surface frame
1828       detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
1829       detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
1830       detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
1831
1832       newX=detX/det;
1833       newY=detY/det;
1834       newZ=detZ/det;
1835
1836       perp.SetXYZ(newX,newY,newZ);
1837       refEta = perp.Eta();
1838       refPhi = perp.Phi();//output in <-pi, pi> range
1839       if(refPhi<0)refPhi+=2*TMath::Pi();
1840       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1841       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1842       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
1843     } while (rad<fJetConeSize+fGammaConeSize || rad2<fJetConeSize+fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
1844     fhRandomPhiEta[2]->Fill(refPhi,refEta);
1845
1846   }
1847   else if(type==2){//random cone
1848     //randomise
1849     do{
1850       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1851       refEta=fGenerator->Uniform(-(0.9-fJetConeSize),0.9-fJetConeSize);
1852       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1853       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1854       //check if reference is not within jet cone or gamma cone +0.4
1855       //example: FF fConeSize=0.4, fJetConeSize=0.4, fIsoGammaConeSize=0.3
1856     } while (rad<fJetConeSize+fGammaConeSize || rad2<fJetConeSize+fJetConeSize);
1857     //photon:0.7=0.4+0.3; jets:0.8=0.4 +0.4 //rad<fConeSize+fJetConeSize rad2<fConeSize+0.3
1858     fhRandomPhiEta[0]->Fill(refPhi,refEta);
1859   }
1860   else if(type==4){//perpendicular to photon axis
1861     Double_t xVar;
1862     Double_t yVar;
1863     Double_t newX=0.;
1864     Double_t newY=0.;
1865     Double_t newZ=0.;
1866     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1867     Double_t jx=gamma.Px();
1868     Double_t jy=gamma.Py();
1869     Double_t jz=gamma.Pz();
1870     //if(jz==0)printf("problem\n");
1871     //X axis
1872     Double_t Xx=1.0-vertex[0];
1873     Double_t Xy=-1.0*vertex[1];
1874     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1875     //Y axis
1876     Double_t Yx=jy*Xz-jz*Xy;
1877     Double_t Yy=jz*Xx-jx*Xz;
1878     Double_t Yz=jx*Xy-jy*Xx;
1879     //Determinant
1880     Double_t det = Xx*Yy*jz + Xy*Yz*jx + Xz*Yx*jy - Xx*Yz*jy - Xy*Yx*jz - Xz*Yy*jx;
1881     if(det==0)printf("problem det==0\n");
1882     Double_t detX = 0.;
1883     Double_t detY = 0.;
1884     Double_t detZ = 0.;
1885
1886 //    Double_t tmpScalar = jx*Xx+jy*Xy+jz*Xz;
1887 //    printf("scalar jet.P o X: %f\n",tmpScalar);
1888 //    tmpScalar = jet.Px()*Yx+jet.Py()*Yy+jet.Pz()*Yz;
1889 //    printf("scalar jet.P o Y: %f\n",tmpScalar);
1890 //    tmpScalar = Xx*Yx+Xy*Yy+Xz*Yz;
1891 //    printf("scalar X o Y: %f\n",tmpScalar);
1892
1893     TVector3 perp;
1894     //randomise
1895     do{
1896       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1897       //refPhi=fGenerator->Uniform(-1,1)*TMath::Pi();
1898       xVar=TMath::Cos(refPhi);
1899       yVar=TMath::Sin(refPhi);
1900       //yVar=TMath::Cos(-TMath::Pi()/2.+refPhi);
1901       //zVar=0 in new surface frame
1902       detX = xVar*Yy*jz + Xz*yVar*jy - xVar*Yz*jy - Xy*yVar*jz;
1903       detY = Xx*yVar*jz + xVar*Yz*jx - xVar*Yx*jz - Xz*yVar*jx;
1904       detZ = Xy*yVar*jx + xVar*Yx*jy - Xx*yVar*jy - xVar*Yy*jx;
1905
1906       newX=detX/det;
1907       newY=detY/det;
1908       newZ=detZ/det;
1909
1910       perp.SetXYZ(newX,newY,newZ);
1911       refEta = perp.Eta();
1912       refPhi = perp.Phi();//output in <-pi, pi> range
1913       if(refPhi<0)refPhi+=2*TMath::Pi();
1914       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1915       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1916       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
1917     } while (rad<fJetConeSize+fGammaConeSize || rad2<fJetConeSize+fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
1918     fhRandomPhiEta[1]->Fill(refPhi,refEta);
1919
1920   }
1921   else if(type==3){//mid point
1922
1923     Double_t jx=jet.Px();
1924     Double_t jy=jet.Py();
1925     Double_t jz=jet.Pz();
1926     //    if(jz==0)printf("problem\n");
1927     Double_t gx=gamma.Px();
1928     Double_t gy=gamma.Py();
1929     Double_t gz=gamma.Pz();
1930
1931     Double_t cosAlpha=(jx*gx+jy*gy+jz*gz)/(jet.Mag()*gamma.Mag());
1932     Double_t cosinus=TMath::Sqrt((cosAlpha+1.)/2.);
1933     //perpendicular axis
1934     Double_t Zx=gy*jz-gz*jy;
1935     Double_t Zy=gz*jx-gx*jz;
1936     Double_t Zz=gx*jy-gy*jx;
1937
1938     //Determinant
1939     Double_t det = Zx*gy*jz + Zy*gz*jx + Zz*gx*jy - Zz*gy*jx - Zy*gx*jz - Zx*gz*jy;
1940
1941     Double_t newX=0.;
1942     Double_t newY=0.;
1943     Double_t newZ=0.;
1944     if(det!=0) {
1945       Double_t detX =            -Zy*gz*cosinus +Zz*cosinus*jy + Zz*gy*cosinus - Zy*cosinus*jz;
1946       Double_t detY = Zx*cosinus*jz          - Zz*gx*cosinus - Zz*cosinus*jx            + Zx*gz*cosinus;
1947       Double_t detZ = -Zx*gy*cosinus + Zy*cosinus*jx + Zy*gx*cosinus - Zx*cosinus*jy;
1948
1949       newX=detX/det;
1950       newY=detY/det;
1951       newZ=detZ/det;
1952     }
1953
1954     TVector3 perp;
1955     perp.SetXYZ(newX,newY,newZ);
1956     refEta = perp.Eta();
1957     refPhi = perp.Phi();//output in <-pi, pi> range
1958     if(refPhi<0)refPhi+=2*TMath::Pi();
1959     rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
1960     rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
1961       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
1962
1963     if (rad<fJetConeSize+fGammaConeSize || rad2<fJetConeSize+fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize) fhRandomPhiEta[3]->Fill(refPhi,refEta);
1964   }
1965   else if(type==5){//tmp
1966     //printf("vertex: %f %f %f \n",vertex[0],vertex[1],vertex[2]);
1967
1968     Double_t xVar;
1969     Double_t newX=0.;
1970     Double_t newY=0.;
1971     Double_t newZ=0.;
1972     //Z axis vector: [jet.Px(),jet.Py(),jet.Pz()]
1973     Double_t jx=jet.Px();
1974     Double_t jy=jet.Py();
1975     Double_t jz=jet.Pz();
1976     //    if(jz==0)printf("problem\n");
1977     //X axis
1978     Double_t Xx=1.0-vertex[0];
1979     Double_t Xy=-1.0*vertex[1];
1980     Double_t Xz=jx/jz*(vertex[0]-1.)+vertex[1]*jy/jz;
1981     //Y axis
1982     Double_t Yx=jy*Xz-jz*Xy;
1983     Double_t Yy=jz*Xx-jx*Xz;
1984     Double_t Yz=jx*Xy-jy*Xx;
1985
1986     // X and Y length
1987     Double_t Xlength=TMath::Sqrt(Xx*Xx+Xy*Xy+Xz*Xz);
1988     Double_t Ylength=TMath::Sqrt(Yx*Yx+Yy*Yy+Yz*Yz);
1989     Double_t ratio=Ylength/Xlength;
1990
1991     TVector3 perp;
1992     //randomise
1993     do{
1994       refPhi=fGenerator->Rndm()*TMath::Pi()*2.;
1995       xVar=TMath::Tan(refPhi)/ratio;
1996       newX=xVar*Yx+Xx;
1997       newY=xVar*Yy+Xy;
1998       newZ=xVar*Yz+Xz;
1999
2000       perp.SetXYZ(newX,newY,newZ);
2001       refEta = perp.Eta();
2002       refPhi = perp.Phi();//output in <-pi, pi> range
2003       if(refPhi<0)refPhi+=2*TMath::Pi();
2004       rad = TMath::Sqrt((gammaEta-refEta)*(gammaEta-refEta) + (gammaPhi-refPhi)*(gammaPhi-refPhi));
2005       rad2 = TMath::Sqrt((jetEta-refEta)*(jetEta-refEta) + (jetPhi-refPhi)*(jetPhi-refPhi));
2006       //printf("refEta,refPhi,rad,rad2: %f %f %f %f\n",refEta,refPhi,rad,rad2);
2007     } while (rad<fJetConeSize+fGammaConeSize || rad2<fJetConeSize+fJetConeSize || TMath::Abs(refEta)>0.9-fJetConeSize);
2008     fhRandomPhiEta[4]->Fill(refPhi,refEta);
2009   }
2010
2011
2012
2013   //calculate FF in background
2014   Int_t ntracks =  0;
2015   ntracks =GetCTSTracks()->GetEntriesFast();
2016   AliVTrack* track = 0x0 ;
2017   TVector3 p3;
2018
2019   Double_t pt = 0, eta = 0, phi = 0;
2020   Int_t  npartcone = 0;
2021   Double_t sumPt=0.;
2022   for(Int_t ipr = 0;ipr < ntracks ; ipr ++ ){
2023     track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
2024     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
2025     pt   = p3.Pt();
2026     if(pt<fPtThresholdInCone) {//0.150
2027       //printf("problem: track pt < %f MeV/c \n",fPtThresholdInCone);
2028       continue;
2029     }
2030     eta  = p3.Eta() ;
2031     phi  = p3.Phi() ;
2032     if(phi < 0) phi+=TMath::TwoPi();
2033     //Check if there is any particle inside cone with pt larger than  fPtThreshold
2034     rad = TMath::Sqrt((eta-refEta)*(eta-refEta) + (phi-refPhi)*(phi-refPhi));
2035     if(rad < fConeSize  && pt > fPtThresholdInCone){
2036       //printf("charged in jet cone pt %f, phi %f, eta %f, R %f \n",pt,phi,eta,rad);
2037         npartcone++;
2038         sumPt+=pt;
2039         if(type==1){//perp jet
2040           fhBkgFFz[1] ->Fill(gammaPt, pt/gammaPt);
2041           fhBkgFFxi[1]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2042           fhBkgFFpt[1]->Fill(gammaPt, pt);
2043         }
2044         else if(type==2){//RC
2045           fhBkgFFz[0] ->Fill(gammaPt, pt/gammaPt);
2046           fhBkgFFxi[0]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2047           fhBkgFFpt[0]->Fill(gammaPt, pt);
2048         }
2049         else if(type==3){//mid point
2050           fhBkgFFz[3] ->Fill(gammaPt, pt/gammaPt);
2051           fhBkgFFxi[3]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2052           fhBkgFFpt[3]->Fill(gammaPt, pt);
2053         }
2054         else if(type==4){//perp jet
2055           fhBkgFFz[2] ->Fill(gammaPt, pt/gammaPt);
2056           fhBkgFFxi[2]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2057           fhBkgFFpt[2]->Fill(gammaPt, pt);
2058         }
2059         else if(type==5){//test
2060           fhBkgFFz[4] ->Fill(gammaPt, pt/gammaPt);
2061           fhBkgFFxi[4]->Fill(gammaPt, TMath::Log(gammaPt/pt));
2062           fhBkgFFpt[4]->Fill(gammaPt, pt);
2063         }
2064
2065
2066     }
2067   }//end of loop over tracks
2068   Double_t sumOverTracks=0.;
2069   if(npartcone!=0) sumOverTracks = sumPt/npartcone;
2070   if(type==1) {
2071     fhBkgNTracksInCone[1]->Fill(gammaPt, npartcone);
2072     fhBkgSumPtInCone[1]->Fill(gammaPt,sumPt);
2073     fhBkgSumPtnTracksInCone[1]->Fill(gammaPt,sumOverTracks);
2074   }
2075   else if(type==2) {
2076     fhBkgNTracksInCone[0]->Fill(gammaPt, npartcone);
2077     fhBkgSumPtInCone[0]->Fill(gammaPt,sumPt);
2078     fhBkgSumPtnTracksInCone[0]->Fill(gammaPt,sumOverTracks);
2079   }
2080   else if(type==3) {
2081     fhBkgNTracksInCone[3]->Fill(gammaPt, npartcone);
2082     fhBkgSumPtInCone[3]->Fill(gammaPt,sumPt);
2083     fhBkgSumPtnTracksInCone[3]->Fill(gammaPt,sumOverTracks);
2084   }
2085   else if(type==4) {
2086     fhBkgNTracksInCone[2]->Fill(gammaPt, npartcone);
2087     fhBkgSumPtInCone[2]->Fill(gammaPt,sumPt);
2088     fhBkgSumPtnTracksInCone[2]->Fill(gammaPt,sumOverTracks);
2089   }
2090   else if(type==5) {
2091     fhBkgNTracksInCone[4]->Fill(gammaPt, npartcone);
2092     fhBkgSumPtInCone[4]->Fill(gammaPt,sumPt);
2093     fhBkgSumPtnTracksInCone[4]->Fill(gammaPt,sumOverTracks);
2094   }
2095 }
2096
2097