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