PWGGA/CaloTrackCorrelations/AliAnaParticleJetFinderCorrelation.cxx

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