PWGGA/CaloTrackCorrelations/AliAnaParticleJetLeadingConeCorrelation.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 that contains the algorithm for the reconstruction of jet, cone around leading particle
  18 // The seed is a backward particle (direct photon)
  19 // 1) Take the trigger particle stored in AliAODPWG4ParticleCorrelation,
  20 // 2) Search for the highest pt leading particle opposite to the photon within a phi, pt window
  21 // 3) Take all particles around leading in a cone R with pt larger than threshold and construct the jet
  22 //
  23 //  Class created from old AliPHOSGammaJet
  24 //  (see AliRoot versions previous Release 4-09)
  25 //
  26 //*-- Author: Gustavo Conesa (LNF-INFN)
  27 //////////////////////////////////////////////////////////////////////////////
  28
  29
  30 // --- ROOT system ---
  31 #include "TH2F.h"
  32 #include "TClonesArray.h"
  33 #include "TClass.h"
  34 //#include "Riostream.h"
  35
  36 //---- Analysis system ----
  37 #include "AliVTrack.h"
  38 #include "AliVCluster.h"
  39 #include "AliCaloTrackReader.h"
  40 #include "AliNeutralMesonSelection.h"
  41 #include "AliAnaParticleJetLeadingConeCorrelation.h"
  42 #include "AliCaloPID.h"
  43 #include "AliAODPWG4ParticleCorrelation.h"
  44 #include "AliFiducialCut.h"
  45
  46 ClassImp(AliAnaParticleJetLeadingConeCorrelation)
  47
  48
  49 //_________________________________________________________________________________
  50 AliAnaParticleJetLeadingConeCorrelation::AliAnaParticleJetLeadingConeCorrelation() :
  51 AliAnaCaloTrackCorrBaseClass(), fJetsOnlyInCTS(kFALSE), fPbPb(kFALSE),
  52 fSeveralConeAndPtCuts(0),  fReMakeJet(0),
  53 fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
  54 fLeadingRatioMaxCut(0.),  fLeadingRatioMinCut(0.),
  55 fJetCTSRatioMaxCut(0.), fJetCTSRatioMinCut(0.),
  56 fJetRatioMaxCut(0.),  fJetRatioMinCut(0.),
  57 fJetNCone(0),fJetNPt(0), fJetCone(0),
  58 fJetPtThreshold(0),fJetPtThresPbPb(0),
  59 fPtTriggerSelectionCut(0.0), fSelect(0),fSelectIsolated(0),
  60 fTrackVector(),fBkgMom(),fJetMom(),fJetConstMom(),
  61 fLeadingMom(),fLeadingPi0Mom(),fLeadingPhoMom1(),fLeadingPhoMom2(),fLeadingChargeMom(),
  62 //Histograms
  63 fOutCont(0x0),
  64 //Leading
  65 fhChargedLeadingPt(0),fhChargedLeadingPhi(0),fhChargedLeadingEta(0),
  66 fhChargedLeadingDeltaPt(0),fhChargedLeadingDeltaPhi(0),fhChargedLeadingDeltaEta(0),
  67 fhChargedLeadingRatioPt(0),
  68 fhNeutralLeadingPt(0),fhNeutralLeadingPhi(0),fhNeutralLeadingEta(0),
  69 fhNeutralLeadingDeltaPt(0),fhNeutralLeadingDeltaPhi(0),fhNeutralLeadingDeltaEta(0),
  70 fhNeutralLeadingRatioPt(0),fhChargedLeadingXi(0), fhNeutralLeadingXi(0),
  71 fhChargedLeadingDeltaPhiRatioPt30(0), fhNeutralLeadingDeltaPhiRatioPt30(0),
  72 fhChargedLeadingDeltaPhiRatioPt50(0), fhNeutralLeadingDeltaPhiRatioPt50(0),
  73 //Jet
  74 fhJetPt(0),fhJetRatioPt(0),fhJetDeltaPhi(0), fhJetDeltaEta(0),
  75 fhJetLeadingRatioPt(0),fhJetLeadingDeltaPhi(0),fhJetLeadingDeltaEta(0),
  76 fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetNTracksInCone(0),
  77 fhBkgPt(0),fhBkgRatioPt(0),fhBkgDeltaPhi(0), fhBkgDeltaEta(0),
  78 fhBkgLeadingRatioPt(0),fhBkgLeadingDeltaPhi(0),fhBkgLeadingDeltaEta(0),
  79 fhBkgFFz(0),fhBkgFFxi(0),fhBkgFFpt(0),fhBkgNTracksInCone(0),
  80 //Several cones and thres histograms
  81 fhJetPts(),fhJetRatioPts(),fhJetDeltaPhis(), fhJetDeltaEtas(),
  82 fhJetLeadingRatioPts(),fhJetLeadingDeltaPhis(),fhJetLeadingDeltaEtas(),
  83 fhJetFFzs(),fhJetFFxis(),fhJetFFpts(),fhJetNTracksInCones(),
  84 fhBkgPts(),fhBkgRatioPts(),fhBkgDeltaPhis(), fhBkgDeltaEtas(),
  85 fhBkgLeadingRatioPts(),fhBkgLeadingDeltaPhis(),fhBkgLeadingDeltaEtas(),
  86 fhBkgFFzs(),fhBkgFFxis(),fhBkgFFpts(),fhBkgNTracksInCones()
  87 {
  88   //Default Ctor
  89
  90   //Initialize parameters
  91
  92   for(Int_t i = 0; i<6; i++){
  93     fJetXMin1[i]     = 0.0 ;
  94     fJetXMin2[i]     = 0.0 ;
  95     fJetXMax1[i]     = 0.0 ;
  96     fJetXMax2[i]     = 0.0 ;
  97     fBkgMean[i]      = 0.0 ;
  98     fBkgRMS[i]       = 0.0 ;
  99     if( i < 2 ){
 100       fJetE1[i]        = 0.0 ;
 101       fJetE2[i]        = 0.0 ;
 102       fJetSigma1[i]    = 0.0 ;
 103       fJetSigma2[i]    = 0.0 ;
 104     }
 105   }
 106
 107   //Several cones and thres histograms
 108   for(Int_t i = 0; i<5; i++){
 109     fJetCones[i]         = 0.0 ;
 110     fJetNameCones[i]     = ""  ;
 111     fJetPtThres[i]      = 0.0 ;
 112     fJetNamePtThres[i]  = ""  ;
 113     for(Int_t j = 0; j<5; j++){
 114       fhJetPts[i][j]=0 ;
 115       fhJetRatioPts[i][j]=0 ;
 116       fhJetDeltaPhis[i][j]=0 ;
 117       fhJetDeltaEtas[i][j]=0 ;
 118       fhJetLeadingRatioPts[i][j]=0 ;
 119       fhJetLeadingDeltaPhis[i][j]=0 ;
 120       fhJetLeadingDeltaEtas[i][j]=0 ;
 121       fhJetFFzs[i][j]=0 ;
 122       fhJetFFxis[i][j]=0 ;
 123       fhJetFFpts[i][j]=0 ;
 124       fhJetNTracksInCones[i][j]=0 ;
 125       fhBkgPts[i][j]=0 ;
 126       fhBkgRatioPts[i][j]=0 ;
 127       fhBkgDeltaPhis[i][j]=0 ;
 128       fhBkgDeltaEtas[i][j]=0 ;
 129       fhBkgLeadingRatioPts[i][j]=0 ;
 130       fhBkgLeadingDeltaPhis[i][j]=0 ;
 131       fhBkgLeadingDeltaEtas[i][j]=0 ;
 132       fhBkgFFzs[i][j]=0 ;
 133       fhBkgFFxis[i][j]=0 ;
 134       fhBkgFFpts[i][j]=0 ;
 135       fhBkgNTracksInCones[i][j]=0 ;
 136     }
 137   }
 138
 139   InitParameters();
 140
 141 }
 142
 143 //____________________________________________________________________________
 144 Double_t AliAnaParticleJetLeadingConeCorrelation::CalculateJetRatioLimit(Double_t ptg, const Double_t *par, const Double_t *x) const {
 145   //Calculate the ratio of the jet and trigger particle limit for the selection
 146   //WARNING: need to check what it does
 147   //printf("CalculateLimit: x1 %2.3f, x2%2.3f\n",x[0],x[1]);
 148   Double_t ePP = par[0] + par[1] * ptg ;
 149   Double_t sPP = par[2] + par[3] * ptg ;
 150   Double_t f   = x[0]   + x[1]   * ptg ;
 151   Double_t ePbPb = ePP + par[4] ;
 152   Double_t sPbPb = TMath::Sqrt(sPP*sPP+ par[5]*par[5]) ;
 153   Double_t rat = (ePbPb - sPbPb * f) / ptg ;
 154   //printf("CalculateLimit: ePP %2.3f, sPP %2.3f, f %2.3f\n", ePP, sPP, f);
 155   //printf("CalculateLimit: ePbPb %2.3f, sPbPb %2.3f, rat %2.3f\n", ePbPb, sPbPb, rat);
 156   return rat ;
 157 }
 158
 159 //___________________________________________________________________________________________________
 160 void AliAnaParticleJetLeadingConeCorrelation::FillJetHistos(AliAODPWG4ParticleCorrelation * particle,
 161                                                             const TLorentzVector jet,
 162                                                             const TString & type, const TString & lastname)
 163 {
 164   //Fill jet and background histograms
 165   Double_t ptTrig  = particle->Pt();
 166   Double_t ptJet   = jet.Pt();
 167   Double_t ptLead  = fLeadingMom.Pt();
 168   Double_t phiTrig = particle->Phi();
 169   Double_t phiJet  = jet.Phi();
 170   if(phiJet < 0) phiJet+=TMath::TwoPi();
 171   Double_t phiLead = fLeadingMom.Phi();
 172   if(phiLead < 0) phiLead+=TMath::TwoPi();
 173   Double_t etaTrig = particle->Eta();
 174   Double_t etaJet  = jet.Eta();
 175   Double_t etaLead = fLeadingMom.Eta();
 176
 177   TH2F *h1 = 0x0;
 178   h1 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 179   if(h1)h1->Fill(ptTrig,ptJet);
 180
 181   TH2F *h2 = 0x0;
 182   h2 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sRatioPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 183   if(h2) h2->Fill(ptTrig,ptJet/ptTrig);
 184
 185   TH2F *h3 = 0x0;
 186   h3 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingRatioPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 187   if(h3)h3->Fill(ptTrig,ptLead/ptJet);
 188
 189   //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))->
 190   //     Fill(ptTrig,phiJet);
 191   TH2F *h4 = 0x0;
 192   h4 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sDeltaPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 193   if(h4) h4->Fill(ptTrig,phiJet-phiTrig);
 194   TH2F *h5 = 0x0;
 195   h5 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingDeltaPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 196   if(h5) h5->Fill(ptTrig,phiJet-phiLead);
 197
 198   //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))->
 199   //     Fill(ptTrig,etaJet);
 200   TH2F *h6 = 0x0;
 201   h6 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sDeltaEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 202   if(h6) h6->Fill(ptTrig,etaJet-etaTrig);
 203   TH2F *h7 = 0x0;
 204   h7 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingDeltaEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 205   if(h7) h7->Fill(ptTrig,etaJet-etaLead);
 206
 207   //Construct fragmentation function
 208   TObjArray * pl = new TObjArray;
 209
 210   if(type == "Jet") pl = particle->GetObjArray(Form("%sTracks",GetAODObjArrayName().Data()));
 211   else if(type == "Bkg") particle->GetObjArray(Form("%sTracksBkg",GetAODObjArrayName().Data()));
 212
 213   if(!pl) return ;
 214
 215   //Different pt cut for jet particles in different collisions systems
 216   //Only needed when jet is recalculated from AODs
 217   //Float_t ptcut = fJetPtThreshold;
 218   //if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
 219
 220   Int_t nTracksInCone = 0;
 221
 222   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ )
 223   {
 224     AliVTrack* track = dynamic_cast<AliVTrack *>(pl->At(ipr)) ;
 225     if(track)fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
 226     else AliDebug(1,"Track not available");
 227
 228     //Recheck if particle is in jet cone
 229     if(fReMakeJet || fSeveralConeAndPtCuts)
 230       if(!IsParticleInJetCone(fTrackVector.Eta(), fTrackVector.Phi(), fLeadingMom.Eta(), fLeadingMom.Phi()) ) continue ;
 231
 232     nTracksInCone++;
 233
 234     TH2F *ha =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFz%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 235     if(ha) ha->Fill(ptTrig,fTrackVector.Pt()/ptTrig);
 236     TH2F *hb  =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFxi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 237     if(hb) hb->Fill(ptTrig,TMath::Log(ptTrig/fTrackVector.Pt()));
 238     TH2F *hc =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFpt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 239     if(hc) hc->Fill(ptTrig,fTrackVector.Pt());
 240
 241   }//track loop
 242
 243   if(nTracksInCone > 0)
 244   {
 245     TH2F *hd = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sNTracksInCone%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 246     if(hd)hd->Fill(ptTrig, nTracksInCone);
 247   }
 248
 249 }
 250
 251 //________________________________________________________________________
 252 TList *  AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects()
 253 {
 254   // Create histograms to be saved in output file and
 255   // store them in fOutCont
 256
 257   fOutCont = new TList() ;
 258   fOutCont->SetName("ParticleJetLeadingInConeCorrelationHistograms") ;
 259
 260   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();
 261   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();
 262   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
 263   Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();
 264   Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();
 265   Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
 266   Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
 267   Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
 268   Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 269
 270   fhChargedLeadingPt  = new TH2F("ChargedLeadingPt","p_{T leading charge} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 271   fhChargedLeadingPt->SetYTitle("p_{T leading charge}");
 272   fhChargedLeadingPt->SetXTitle("p_{T trigger} (GeV/c)");
 273
 274   fhChargedLeadingPhi  = new TH2F("ChargedLeadingPhi","#phi_{h^{#pm}}  vs p_{T trigger}", nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 275   fhChargedLeadingPhi->SetYTitle("#phi_{h^{#pm}} (rad)");
 276   fhChargedLeadingPhi->SetXTitle("p_{T trigger} (GeV/c)");
 277
 278   fhChargedLeadingEta  = new TH2F("ChargedLeadingEta","#eta_{h^{#pm}}  vs p_{T trigger}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 279   fhChargedLeadingEta->SetYTitle("#eta_{h^{#pm}} ");
 280   fhChargedLeadingEta->SetXTitle("p_{T trigger} (GeV/c)");
 281
 282   fhChargedLeadingDeltaPt  = new TH2F("ChargedLeadingDeltaPt","p_{T trigger} - p_{T h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 283   fhChargedLeadingDeltaPt->SetYTitle("#Delta p_{T} (GeV/c)");
 284   fhChargedLeadingDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
 285
 286   fhChargedLeadingDeltaPhi  = new TH2F("ChargedLeadingDeltaPhi","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 287   fhChargedLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
 288   fhChargedLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 289
 290   fhChargedLeadingDeltaEta  = new TH2F("ChargedLeadingDeltaEta","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 291   fhChargedLeadingDeltaEta->SetYTitle("#Delta #eta");
 292   fhChargedLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 293
 294   fhChargedLeadingRatioPt  = new TH2F("ChargedLeadingRatioPt","p_{T leading charge} /p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 295   fhChargedLeadingRatioPt->SetYTitle("p_{T lead charge} /p_{T trigger}");
 296   fhChargedLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 297
 298   fhChargedLeadingXi  = new TH2F("ChargedLeadingXi","ln(p_{T trigger} / p_{T leading charge} ) vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,10);
 299   fhChargedLeadingXi->SetYTitle("#xi");
 300   fhChargedLeadingXi->SetXTitle("p_{T trigger} (GeV/c)");
 301
 302   fOutCont->Add(fhChargedLeadingPt) ;
 303   fOutCont->Add(fhChargedLeadingPhi) ;
 304   fOutCont->Add(fhChargedLeadingEta) ;
 305   fOutCont->Add(fhChargedLeadingDeltaPt) ;
 306   fOutCont->Add(fhChargedLeadingDeltaPhi) ;
 307   fOutCont->Add(fhChargedLeadingDeltaEta) ;
 308   fOutCont->Add(fhChargedLeadingRatioPt) ;
 309   fOutCont->Add(fhChargedLeadingXi) ;
 310
 311   fhChargedLeadingDeltaPhiRatioPt30  = new TH2F("ChargedLeadingDeltaPhiRatioPt30","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, charged leading, p_{T trigger} > 30 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1);
 312   fhChargedLeadingDeltaPhiRatioPt30->SetXTitle("#Delta #phi (rad)");
 313   fhChargedLeadingDeltaPhiRatioPt30->SetYTitle("p_{T leading} / p_{T trigger}");
 314
 315   fhChargedLeadingDeltaPhiRatioPt50  = new TH2F("ChargedLeadingDeltaPhiRatioPt50","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, charged leading, p_{T trigger} > 50 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1);
 316   fhChargedLeadingDeltaPhiRatioPt50->SetXTitle("#Delta #phi (rad)");
 317   fhChargedLeadingDeltaPhiRatioPt50->SetYTitle("p_{T leading} / p_{T trigger}");
 318
 319   fOutCont->Add(fhChargedLeadingDeltaPhiRatioPt30) ;
 320   fOutCont->Add(fhChargedLeadingDeltaPhiRatioPt50) ;
 321
 322   if(!fJetsOnlyInCTS){
 323
 324     fhNeutralLeadingPt  = new TH2F("NeutralLeadingPt","p_{T leading #pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 325     fhNeutralLeadingPt->SetYTitle("p_{T leading #pi^{0}}");
 326     fhNeutralLeadingPt->SetXTitle("p_{T trigger} (GeV/c)");
 327
 328     fhNeutralLeadingPhi  = new TH2F("NeutralLeadingPhi","#phi_{#pi^{0}}  vs p_{T trigger}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 329     fhNeutralLeadingPhi->SetYTitle("#phi_{#pi^{0}} (rad)");
 330     fhNeutralLeadingPhi->SetXTitle("p_{T trigger} (GeV/c)");
 331
 332     fhNeutralLeadingEta  = new TH2F("NeutralLeadingEta","#eta_{#pi^{0}}  vs p_{T trigger}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 333     fhNeutralLeadingEta->SetYTitle("#eta_{#pi^{0}} ");
 334     fhNeutralLeadingEta->SetXTitle("p_{T trigger} (GeV/c)");
 335
 336     fhNeutralLeadingDeltaPt  = new TH2F("NeutralLeadingDeltaPt","p_{T trigger} - p_{T #pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 337     fhNeutralLeadingDeltaPt->SetYTitle("#Delta p_{T} (GeV/c)");
 338     fhNeutralLeadingDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
 339
 340     fhNeutralLeadingDeltaPhi  = new TH2F("NeutralLeadingDeltaPhi","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 341     fhNeutralLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
 342     fhNeutralLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 343
 344     fhNeutralLeadingDeltaEta  = new TH2F("NeutralLeadingDeltaEta","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 345     fhNeutralLeadingDeltaEta->SetYTitle("#Delta #eta");
 346     fhNeutralLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 347
 348     fhNeutralLeadingRatioPt  = new TH2F("NeutralLeadingRatioPt","p_{T leading #pi^{0}} /p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 349     fhNeutralLeadingRatioPt->SetYTitle("p_{T lead #pi^{0}} /p_{T trigger}");
 350     fhNeutralLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 351
 352     fhNeutralLeadingXi  = new TH2F("NeutralLeadingXi","ln(p_{T trigger} / p_{T leading #pi^{0}} ) vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,10);
 353     fhNeutralLeadingXi->SetYTitle("#xi");
 354     fhNeutralLeadingXi->SetXTitle("p_{T trigger} (GeV/c)");
 355
 356     fOutCont->Add(fhNeutralLeadingPt) ;
 357     fOutCont->Add(fhNeutralLeadingPhi) ;
 358     fOutCont->Add(fhNeutralLeadingEta) ;
 359     fOutCont->Add(fhNeutralLeadingDeltaPt) ;
 360     fOutCont->Add(fhNeutralLeadingDeltaPhi) ;
 361     fOutCont->Add(fhNeutralLeadingDeltaEta) ;
 362     fOutCont->Add(fhNeutralLeadingRatioPt) ;
 363     fOutCont->Add(fhNeutralLeadingXi) ;
 364
 365     fhNeutralLeadingDeltaPhiRatioPt30  = new TH2F("NeutralLeadingDeltaPhiRatioPt30","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, neutral leading, p_{T trigger} > 30 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1);
 366     fhNeutralLeadingDeltaPhiRatioPt30->SetXTitle("#Delta #phi (rad)");
 367     fhNeutralLeadingDeltaPhiRatioPt30->SetYTitle("p_{T leading} / p_{T trigger}");
 368
 369     fhNeutralLeadingDeltaPhiRatioPt50  = new TH2F("NeutralLeadingDeltaPhiRatioPt50","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, neutral leading, p_{T trigger} > 50 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1);
 370     fhNeutralLeadingDeltaPhiRatioPt50->SetXTitle("#Delta #phi (rad)");
 371     fhNeutralLeadingDeltaPhiRatioPt50->SetYTitle("p_{T leading} / p_{T trigger}");
 372     fOutCont->Add(fhNeutralLeadingDeltaPhiRatioPt30) ;
 373     fOutCont->Add(fhNeutralLeadingDeltaPhiRatioPt50) ;
 374
 375   }
 376
 377   if(!fSeveralConeAndPtCuts){// not several cones
 378
 379     //Jet Distributions
 380     fhJetPt  = new TH2F("JetPt","p_{T  jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 381     fhJetPt->SetYTitle("p_{T  jet}");
 382     fhJetPt->SetXTitle("p_{T trigger} (GeV/c)");
 383
 384     fhJetRatioPt  = new TH2F("JetRatioPt","p_{T  jet}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 385     fhJetRatioPt->SetYTitle("p_{T  jet}/p_{T trigger}");
 386     fhJetRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 387
 388     fhJetDeltaPhi  = new TH2F("JetDeltaPhi","#phi_{jet} - #phi_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 389     fhJetDeltaPhi->SetYTitle("#Delta #phi (rad)");
 390     fhJetDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 391
 392     fhJetDeltaEta  = new TH2F("JetDeltaEta","#eta_{jet} - #eta_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 393     fhJetDeltaEta->SetYTitle("#Delta #eta");
 394     fhJetDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 395
 396     fhJetLeadingRatioPt  = new TH2F("JetLeadingRatioPt","p_{T  jet} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 397     fhJetLeadingRatioPt->SetYTitle("p_{T  leading}/p_{T jet}");
 398     fhJetLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 399
 400     fhJetLeadingDeltaPhi  = new TH2F("JetLeadingDeltaPhi","#phi_{jet} - #phi_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-TMath::Pi(),TMath::Pi());
 401     fhJetLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
 402     fhJetLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 403
 404     fhJetLeadingDeltaEta  = new TH2F("JetLeadingDeltaEta","#eta_{jet} - #eta_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 405     fhJetLeadingDeltaEta->SetYTitle("#Delta #eta");
 406     fhJetLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 407
 408     fhJetFFz  = new TH2F("JetFFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,2);
 409     fhJetFFz->SetYTitle("z");
 410     fhJetFFz->SetXTitle("p_{T trigger}");
 411
 412     fhJetFFxi  = new TH2F("JetFFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0.,10.);
 413     fhJetFFxi->SetYTitle("#xi");
 414     fhJetFFxi->SetXTitle("p_{T trigger}");
 415
 416     fhJetFFpt  = new TH2F("JetFFpt","#xi = p_{T i charged}) vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,50.);
 417     fhJetFFpt->SetYTitle("p_{T charged hadron}");
 418     fhJetFFpt->SetXTitle("p_{T trigger}");
 419
 420     fhJetNTracksInCone  = new TH2F("JetNTracksInCone","N particles in cone vs p_{T trigger}",nptbins,ptmin,ptmax,5000,0, 5000);
 421     fhJetNTracksInCone->SetYTitle("N tracks in jet cone");
 422     fhJetNTracksInCone->SetXTitle("p_{T trigger} (GeV/c)");
 423
 424     fOutCont->Add(fhJetPt) ;
 425     fOutCont->Add(fhJetRatioPt) ;
 426     fOutCont->Add(fhJetDeltaPhi) ;
 427     fOutCont->Add(fhJetDeltaEta) ;
 428     fOutCont->Add(fhJetLeadingRatioPt) ;
 429     fOutCont->Add(fhJetLeadingDeltaPhi) ;
 430     fOutCont->Add(fhJetLeadingDeltaEta) ;
 431     fOutCont->Add(fhJetFFz) ;
 432     fOutCont->Add(fhJetFFxi) ;
 433     fOutCont->Add(fhJetFFpt) ;
 434     fOutCont->Add(fhJetNTracksInCone) ;
 435
 436     //Bkg Distributions
 437     fhBkgPt  = new TH2F("BkgPt","p_{T  bkg} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 438     fhBkgPt->SetYTitle("p_{T  bkg}");
 439     fhBkgPt->SetXTitle("p_{T trigger} (GeV/c)");
 440
 441     fhBkgRatioPt  = new TH2F("BkgRatioPt","p_{T  bkg}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 442     fhBkgRatioPt->SetYTitle("p_{T  bkg}/p_{T trigger}");
 443     fhBkgRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 444
 445     fhBkgDeltaPhi  = new TH2F("BkgDeltaPhi","#phi_{bkg} - #phi_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 446     fhBkgDeltaPhi->SetYTitle("#Delta #phi (rad)");
 447     fhBkgDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 448
 449     fhBkgDeltaEta  = new TH2F("BkgDeltaEta","#eta_{bkg} - #eta_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 450     fhBkgDeltaEta->SetYTitle("#Delta #eta");
 451     fhBkgDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 452
 453     fhBkgLeadingRatioPt  = new TH2F("BkgLeadingRatioPt","p_{T  bkg} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 454     fhBkgLeadingRatioPt->SetYTitle("p_{T  leading}/p_{T bkg}");
 455     fhBkgLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 456
 457     fhBkgLeadingDeltaPhi  = new TH2F("BkgLeadingDeltaPhi","#phi_{bkg} - #phi_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 458     fhBkgLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
 459     fhBkgLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 460
 461     fhBkgLeadingDeltaEta  = new TH2F("BkgLeadingDeltaEta","#eta_{bkg} - #eta_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 462     fhBkgLeadingDeltaEta->SetYTitle("#Delta #eta");
 463     fhBkgLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 464
 465     fhBkgFFz  = new TH2F("BkgFFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", nptbins,ptmin,ptmax,200,0.,2);
 466     fhBkgFFz->SetYTitle("z");
 467     fhBkgFFz->SetXTitle("p_{T trigger}");
 468
 469     fhBkgFFxi  = new TH2F("BkgFFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,10.);
 470     fhBkgFFxi->SetYTitle("#xi");
 471     fhBkgFFxi->SetXTitle("p_{T trigger}");
 472
 473     fhBkgFFpt  = new TH2F("BkgFFpt","p_{T charged hadron } vs p_{T trigger}", nptbins,ptmin,ptmax,200,0.,50.);
 474     fhBkgFFpt->SetYTitle("p_{T charged} hadron");
 475     fhBkgFFpt->SetXTitle("p_{T trigger}");
 476
 477     fhBkgNTracksInCone  = new TH2F("BkgNTracksInCone","N particles in cone vs p_{T trigger}",nptbins,ptmin,ptmax,5000,0, 5000);
 478     fhBkgNTracksInCone->SetYTitle("N tracks in bkg cone");
 479     fhBkgNTracksInCone->SetXTitle("p_{T trigger} (GeV/c)");
 480
 481     fOutCont->Add(fhBkgPt) ;
 482     fOutCont->Add(fhBkgRatioPt) ;
 483     fOutCont->Add(fhBkgDeltaPhi) ;
 484     fOutCont->Add(fhBkgDeltaEta) ;
 485     fOutCont->Add(fhBkgLeadingRatioPt) ;
 486     fOutCont->Add(fhBkgLeadingDeltaPhi) ;
 487     fOutCont->Add(fhBkgLeadingDeltaEta) ;
 488     fOutCont->Add(fhBkgFFz) ;
 489     fOutCont->Add(fhBkgFFxi) ;
 490     fOutCont->Add(fhBkgFFpt) ;
 491     fOutCont->Add(fhBkgNTracksInCone) ;
 492
 493   }//not several cones
 494   else{ //If we want to study the jet for different cones and pt
 495     for(Int_t icone = 0; icone<fJetNCone; icone++){//icone
 496       for(Int_t ipt = 0; ipt<fJetNPt;ipt++){ //ipt
 497
 498         TString lastnamehist ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt];
 499         TString lastnametitle =", cone ="+fJetNameCones[icone]+", pt > " +fJetNamePtThres[ipt]+" GeV/c";
 500
 501         //Jet Distributions
 502         fhJetPts[icone][ipt] = new TH2F(Form("JetPt%s",lastnamehist.Data()),Form("p_{T  jet} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 503         fhJetPts[icone][ipt]->SetYTitle("p_{T  jet}");
 504         fhJetPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 505
 506         fhJetRatioPts[icone][ipt] = new TH2F(Form("JetRatioPt%s",lastnamehist.Data()),Form("p_{T  jet}/p_{T trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2);
 507         fhJetRatioPts[icone][ipt]->SetYTitle("p_{T  jet}/p_{T trigger}");
 508         fhJetRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 509
 510         fhJetDeltaPhis[icone][ipt] = new TH2F(Form("JetDeltaPhi%s",lastnamehist.Data()),Form("#phi_{jet} - #phi_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 511         fhJetDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
 512         fhJetDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 513
 514         fhJetDeltaEtas[icone][ipt] = new TH2F(Form("JetDeltaEta%s",lastnamehist.Data()),Form("#eta_{jet} - #eta_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2);
 515         fhJetDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
 516         fhJetDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 517
 518         fhJetLeadingRatioPts[icone][ipt] = new TH2F(Form("JetLeadingRatioPt%s",lastnamehist.Data()),Form("p_{T  jet} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2);
 519         fhJetLeadingRatioPts[icone][ipt]->SetYTitle("p_{T  leading}/p_{T jet}");
 520         fhJetLeadingRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 521
 522         fhJetLeadingDeltaPhis[icone][ipt] = new TH2F(Form("JetLeadingDeltaPhi%s",lastnamehist.Data()),Form("#phi_{jet} - #phi_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 523         fhJetLeadingDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
 524         fhJetLeadingDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 525
 526         fhJetLeadingDeltaEtas[icone][ipt] = new TH2F(Form("JetLeadingDeltaEta%s",lastnamehist.Data()),Form("#eta_{jet} - #eta_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2);
 527         fhJetLeadingDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
 528         fhJetLeadingDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 529
 530         fhJetFFzs[icone][ipt] = new TH2F(Form("JetFFz%s",lastnamehist.Data()),"z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2);
 531         fhJetFFzs[icone][ipt]->SetYTitle("z");
 532         fhJetFFzs[icone][ipt]->SetXTitle("p_{T trigger}");
 533
 534         fhJetFFxis[icone][ipt] = new TH2F(Form("JetFFxi%s",lastnamehist.Data()),"#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.);
 535         fhJetFFxis[icone][ipt]->SetYTitle("#xi");
 536         fhJetFFxis[icone][ipt]->SetXTitle("p_{T trigger}");
 537
 538         fhJetFFpts[icone][ipt] = new TH2F(Form("JetFFpt%s",lastnamehist.Data()),"p_{T charged hadron } in jet vs p_{T trigger}", 120,0.,120.,200,0.,50.);
 539         fhJetFFpts[icone][ipt]->SetYTitle("p_{T charged hadron}");
 540         fhJetFFpts[icone][ipt]->SetXTitle("p_{T trigger}");
 541
 542         fhJetNTracksInCones[icone][ipt] = new TH2F(Form("JetNTracksInCone%s",lastnamehist.Data()),Form("N particles in cone vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,5000,0, 5000);
 543         fhJetNTracksInCones[icone][ipt]->SetYTitle("N tracks in jet cone");
 544         fhJetNTracksInCones[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 545
 546         fOutCont->Add(fhJetPts[icone][ipt]) ;
 547         fOutCont->Add(fhJetRatioPts[icone][ipt]) ;
 548         fOutCont->Add(fhJetDeltaPhis[icone][ipt]) ;
 549         fOutCont->Add(fhJetDeltaEtas[icone][ipt]) ;
 550         fOutCont->Add(fhJetLeadingRatioPts[icone][ipt]) ;
 551         fOutCont->Add(fhJetLeadingDeltaPhis[icone][ipt]) ;
 552         fOutCont->Add(fhJetLeadingDeltaEtas[icone][ipt]) ;
 553         fOutCont->Add(fhJetFFzs[icone][ipt]) ;
 554         fOutCont->Add(fhJetFFxis[icone][ipt]) ;
 555         fOutCont->Add(fhJetFFpts[icone][ipt]) ;
 556         fOutCont->Add(fhJetNTracksInCones[icone][ipt]) ;
 557
 558         //Bkg Distributions
 559         fhBkgPts[icone][ipt] = new TH2F(Form("BkgPt%s",lastnamehist.Data()),Form("p_{T  bkg} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 560         fhBkgPts[icone][ipt]->SetYTitle("p_{T  bkg}");
 561         fhBkgPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 562
 563         fhBkgRatioPts[icone][ipt] = new TH2F(Form("BkgRatioPt%s",lastnamehist.Data()),Form("p_{T  bkg}/p_{T trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2);
 564         fhBkgRatioPts[icone][ipt]->SetYTitle("p_{T  bkg}/p_{T trigger}");
 565         fhBkgRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 566
 567         fhBkgDeltaPhis[icone][ipt] = new TH2F(Form("BkgDeltaPhi%s",lastnamehist.Data()),Form("#phi_{bkg} - #phi_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 568         fhBkgDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
 569         fhBkgDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 570
 571         fhBkgDeltaEtas[icone][ipt] = new TH2F(Form("BkgDeltaEta%s",lastnamehist.Data()),Form("#eta_{bkg} - #eta_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2);
 572         fhBkgDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
 573         fhBkgDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 574
 575         fhBkgLeadingRatioPts[icone][ipt] = new TH2F(Form("BkgLeadingRatioPt%s",lastnamehist.Data()),Form("p_{T  bkg} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2);
 576         fhBkgLeadingRatioPts[icone][ipt]->SetYTitle("p_{T  leading}/p_{T bkg}");
 577         fhBkgLeadingRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 578
 579         fhBkgLeadingDeltaPhis[icone][ipt] = new TH2F(Form("BkgLeadingDeltaPhi%s",lastnamehist.Data()),Form("#phi_{bkg} - #phi_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 580         fhBkgLeadingDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
 581         fhBkgLeadingDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 582
 583         fhBkgLeadingDeltaEtas[icone][ipt] = new TH2F(Form("BkgLeadingDeltaEta%s",lastnamehist.Data()),Form("#eta_{bkg} - #eta_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2);
 584         fhBkgLeadingDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
 585         fhBkgLeadingDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 586
 587         fhBkgFFzs[icone][ipt] = new TH2F(Form("BkgFFz%s",lastnamehist.Data()),"z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2);
 588         fhBkgFFzs[icone][ipt]->SetYTitle("z");
 589         fhBkgFFzs[icone][ipt]->SetXTitle("p_{T trigger}");
 590
 591         fhBkgFFxis[icone][ipt] = new TH2F(Form("BkgFFxi%s",lastnamehist.Data()),"#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.);
 592         fhBkgFFxis[icone][ipt]->SetYTitle("#xi");
 593         fhBkgFFxis[icone][ipt]->SetXTitle("p_{T trigger}");
 594
 595         fhBkgFFpts[icone][ipt] = new TH2F(Form("BkgFFpt%s",lastnamehist.Data()),"p_{T charged hadron} in jet vs p_{T trigger}", 120,0.,120.,200,0.,50.);
 596         fhBkgFFpts[icone][ipt]->SetYTitle("p_{T charged hadron}");
 597         fhBkgFFpts[icone][ipt]->SetXTitle("p_{T trigger}");
 598
 599         fhBkgNTracksInCones[icone][ipt] = new TH2F(Form("BkgNTracksInCone%s",lastnamehist.Data()),Form("N particles in cone vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,5000,0, 5000);
 600         fhBkgNTracksInCones[icone][ipt]->SetYTitle("N tracks in bkg cone");
 601         fhBkgNTracksInCones[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 602
 603         fOutCont->Add(fhBkgPts[icone][ipt]) ;
 604         fOutCont->Add(fhBkgRatioPts[icone][ipt]) ;
 605         fOutCont->Add(fhBkgDeltaPhis[icone][ipt]) ;
 606         fOutCont->Add(fhBkgDeltaEtas[icone][ipt]) ;
 607         fOutCont->Add(fhBkgLeadingRatioPts[icone][ipt]) ;
 608         fOutCont->Add(fhBkgLeadingDeltaPhis[icone][ipt]) ;
 609         fOutCont->Add(fhBkgLeadingDeltaEtas[icone][ipt]) ;
 610         fOutCont->Add(fhBkgFFzs[icone][ipt]) ;
 611         fOutCont->Add(fhBkgFFxis[icone][ipt]) ;
 612         fOutCont->Add(fhBkgFFpts[icone][ipt]) ;
 613         fOutCont->Add(fhBkgNTracksInCones[icone][ipt]) ;
 614
 615       }//ipt
 616     } //icone
 617   }//If we want to study any cone or pt threshold
 618
 619   //Keep neutral meson selection histograms if requiered
 620   //Setting done in AliNeutralMesonSelection
 621   if(GetNeutralMesonSelection()){
 622     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
 623     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
 624       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) fOutCont->Add(nmsHistos->At(i)) ;
 625     delete nmsHistos;
 626   }
 627
 628
 629   //  if(GetDebug() > 2)
 630   //  {
 631   //    printf("AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects() - All histograms names : \n");
 632   //    for(Int_t i  = 0 ;  i<  fOutCont->GetEntries(); i++)
 633   //      printf("Histo i %d name %s\n",i,((fOutCont->At(i))->GetName()));
 634   //    //cout<< (fOutCont->At(i))->GetName()<<endl;
 635   //  }
 636
 637   return fOutCont;
 638
 639 }
 640
 641 //__________________________________________________________________________________________________________
 642 Bool_t  AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle(AliAODPWG4ParticleCorrelation * particle)
 643 {
 644   //Search Charged or Neutral leading particle, select the highest one and fill AOD
 645
 646
 647   GetLeadingCharge(particle) ;
 648   if(!fJetsOnlyInCTS) GetLeadingPi0(particle) ;
 649
 650   Double_t ptch = fLeadingChargeMom.Pt();
 651   Double_t ptpi = fLeadingPi0Mom   .Pt();
 652
 653   if (ptch > 0 || ptpi > 0)
 654   {
 655     if((ptch >= ptpi))
 656     {
 657       AliDebug(1,"Leading found in CTS");
 658
 659       fLeadingMom = fLeadingChargeMom;
 660
 661       AliDebug(1,Form("Found Leading: pt %2.3f, phi %2.3f deg, eta %2.3f",fLeadingMom.Pt(),fLeadingMom.Phi()*TMath::RadToDeg(),fLeadingMom.Eta())) ;
 662
 663       //Put leading in AOD
 664       particle->SetLeading(fLeadingChargeMom);
 665       particle->SetLeadingDetector(kCTS);
 666       return kTRUE;
 667     }
 668     else
 669     {
 670       AliDebug(1,"Leading found in EMCAL");
 671
 672       fLeadingMom = fLeadingPi0Mom;
 673
 674       AliDebug(1,Form("Found Leading: pt %2.3f, phi %2.3f, eta %2.3f",fLeadingMom.Pt(),fLeadingMom.Phi()*TMath::RadToDeg(),fLeadingMom.Eta())) ;
 675       //Put leading in AOD
 676       particle->SetLeading(fLeadingPi0Mom);
 677       particle->SetLeadingDetector(kEMCAL);
 678       return kTRUE;
 679     }
 680   }
 681
 682   AliDebug(1,"NO LEADING PARTICLE FOUND");
 683
 684   return kFALSE;
 685
 686 }
 687
 688 //_______________________________________________________________________________________________________
 689 void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge(AliAODPWG4ParticleCorrelation * particle)
 690 {
 691   //Search for the charged particle with highest pt and with
 692   //Phi=Phi_trigger-Pi and pT=0.1E_gamma
 693
 694   if(!GetCTSTracks()) return;
 695
 696   Double_t ptTrig  = particle->Pt();
 697   Double_t phiTrig = particle->Phi();
 698   Double_t rat     = -100 ;
 699   Double_t ptl     = -100 ;
 700   Double_t phil    = -100 ;
 701   Double_t pt      = -100.;
 702   Double_t phi     = -100.;
 703
 704   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
 705   {
 706     AliVTrack* track = (AliVTrack *)(GetCTSTracks()->At(ipr)) ;
 707     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
 708     pt   = fTrackVector.Pt();
 709     phi  = fTrackVector.Phi() ;
 710     if(phi < 0) phi+=TMath::TwoPi();
 711     rat  = pt/ptTrig ;
 712
 713     //printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge() - Tracks: pt %2.3f eta %2.3f phi %2.3f pt/ptTrig %2.3f \n",
 714     //     pt, fTrackVector.Eta(), phi,pt/ptTrig) ;
 715
 716     Float_t deltaphi = TMath::Abs(phiTrig-phi);
 717     if((deltaphi > fDeltaPhiMinCut)     && (deltaphi < fDeltaPhiMaxCut) &&
 718        (rat      > fLeadingRatioMinCut) && (rat      < fLeadingRatioMaxCut)  && (pt  > ptl))
 719     {
 720       phil = phi ;
 721       ptl  = pt ;
 722       fLeadingChargeMom.SetVect(fTrackVector);
 723     }
 724   }// track loop
 725
 726   if( ptl > 0 )AliDebug(1,Form("Leading in CTS: pt %2.3f eta %2.3f phi %2.3f pt/ptTrig %2.3f, |phiTrig-phi| %2.3f",
 727                                ptl, fLeadingChargeMom.Eta(), phil,ptl/ptTrig, TMath::Abs(phiTrig-phil))) ;
 728
 729 }
 730
 731 //____________________________________________________________________________________________________
 732 void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0(AliAODPWG4ParticleCorrelation * particle)
 733 {
 734   //Search for the neutral pion with highest pt and with
 735   //Phi=Phi_trigger-Pi and pT=0.1E_gamma
 736
 737   if(!GetEMCALClusters()) return ;
 738
 739   Double_t ptTrig  = particle->Pt();
 740   Double_t phiTrig = particle->Phi();
 741   Double_t rat     = -100 ;
 742   Double_t ptl     = -100 ;
 743   Double_t phil    = -100 ;
 744   Double_t pt      = -100.;
 745   Double_t phi     = -100.;
 746
 747   //Get vertex for photon momentum calculation
 748   Double_t vertex [] = {0,0,0} ; //vertex
 749   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 750   {
 751     GetVertex(vertex);
 752   }
 753
 754   //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
 755   for(Int_t iclus = 0;iclus < GetEMCALClusters()->GetEntriesFast() ; iclus ++ )
 756   {
 757     AliVCluster * calo = (AliVCluster *)(GetEMCALClusters()->At(iclus)) ;
 758
 759     //Cluster selection, not charged, with photon or pi0 id and in fiducial cut
 760     Int_t pdgi=0;
 761     if(!SelectCluster(calo, vertex,  fLeadingPhoMom1, pdgi))  continue ;
 762
 763     AliDebug(1,Form("Neutral cluster: pt %2.3f, phi %2.3f",fLeadingPhoMom1.Pt(),fLeadingPhoMom1.Phi()));
 764
 765     //2 gamma overlapped, found with PID
 766     if(pdgi == AliCaloPID::kPi0)
 767     {
 768       AliDebug(1,"Neutral cluster ID as Pi0");
 769
 770       pt  = fLeadingPhoMom1.Pt();
 771       rat = pt/ptTrig;
 772       phi = fLeadingPhoMom1.Phi();
 773       if(phi < 0) phi+=TMath::TwoPi();
 774
 775       //Selection within angular and energy limits
 776       Float_t deltaphi = TMath::Abs(phiTrig-phi);
 777       if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  &&
 778          deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut )
 779       {
 780         phil = phi ;
 781         ptl  = pt ;
 782         fLeadingPi0Mom.SetPxPyPzE(fLeadingPhoMom1.Px(),fLeadingPhoMom1.Py(),fLeadingPhoMom1.Pz(),fLeadingPhoMom1.E());
 783       }// cuts
 784     }// pdg = AliCaloPID::kPi0
 785     //Make invariant mass analysis
 786     else if(pdgi == AliCaloPID::kPhoton)
 787     {
 788       //Search the photon companion in case it comes from  a Pi0 decay
 789       //Apply several cuts to select the good pair
 790       for(Int_t jclus = iclus+1; jclus < GetEMCALClusters()->GetEntriesFast() ; jclus++ )
 791       {
 792         AliVCluster * calo2 = (AliVCluster *) (GetEMCALClusters()->At(jclus)) ;
 793
 794         //Cluster selection, not charged with photon or pi0 id and in fiducial cut
 795         Int_t pdgj=0;
 796
 797         if     (!SelectCluster(calo2, vertex,  fLeadingPhoMom2, pdgj))  continue ;
 798
 799         if(pdgj == AliCaloPID::kPhoton )
 800         {
 801           pt  = (fLeadingPhoMom1+fLeadingPhoMom2).Pt();
 802           phi = (fLeadingPhoMom1+fLeadingPhoMom2).Phi();
 803           if(phi < 0) phi+=TMath::TwoPi();
 804           rat = pt/ptTrig;
 805
 806           //Selection within angular and energy limits
 807           Float_t deltaphi = TMath::Abs(phiTrig-phi);
 808
 809           AliDebug(1,Form("Neutral Hadron Correlation: gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, |phiTrig-phi| %2.3f, pt/ptTrig %2.3f, M %2.3f",
 810                           pt,phi,(fLeadingPhoMom1+fLeadingPhoMom2).Eta(), deltaphi, rat, (fLeadingPhoMom1+fLeadingPhoMom2).M()));
 811
 812           if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  &&
 813              deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut ){
 814             //Select good pair (aperture and invariant mass)
 815             if(GetNeutralMesonSelection()->SelectPair(fLeadingPhoMom1, fLeadingPhoMom2,kEMCAL)){
 816               phil = phi ;
 817               ptl  = pt ;
 818               fLeadingPi0Mom=(fLeadingPhoMom1+fLeadingPhoMom2);
 819
 820               AliDebug(1,Form("Neutral Hadron Correlation: Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f",
 821                               ptl,phil,(fLeadingPhoMom1+fLeadingPhoMom2).Eta(), (fLeadingPhoMom1+fLeadingPhoMom2).M()));
 822             }//pi0 selection
 823
 824
 825           }//Pair selected as leading
 826         }//if pair of gammas
 827       }//2nd loop
 828     }// if pdg = 22
 829   }// 1st Loop
 830
 831   if(fLeadingPi0Mom.Pt() > 0 )
 832     AliDebug(1,Form("Leading EMCAL: pt %2.3f eta %2.3f phi %2.3f pt/Eg %2.3f",
 833                     fLeadingPi0Mom.Pt(), fLeadingPi0Mom.Eta(), phil,  fLeadingPi0Mom.Pt()/ptTrig)) ;
 834
 835 }
 836
 837 //____________________________________________________________________________
 838 void AliAnaParticleJetLeadingConeCorrelation::InitParameters()
 839 {
 840   //Initialize the parameters of the analysis.
 841   SetInputAODName("PWG4Particle");
 842   SetAODObjArrayName("JetLeadingCone");
 843   AddToHistogramsName("AnaJetLCCorr_");
 844
 845   fJetsOnlyInCTS      = kFALSE ;
 846   fPbPb               = kFALSE ;
 847   fReMakeJet          = kFALSE ;
 848
 849   //Leading selection parameters
 850   fDeltaPhiMinCut     = 2.9 ;
 851   fDeltaPhiMaxCut     = 3.4 ;
 852   fLeadingRatioMinCut = 0.1;
 853   fLeadingRatioMaxCut = 1.5;
 854
 855   //Jet selection parameters
 856   //Fixed cut
 857   fJetRatioMaxCut     = 1.2 ;
 858   fJetRatioMinCut     = 0.3 ;
 859   fJetCTSRatioMaxCut  = 1.2 ;
 860   fJetCTSRatioMinCut  = 0.3 ;
 861   fSelect               = 0  ; //0, Accept all jets, 1, selection depends on energy, 2 fixed selection
 862
 863   fSelectIsolated = kFALSE;
 864
 865   //Cut depending on gamma energy
 866   fPtTriggerSelectionCut = 10.; //For Low pt jets+BKG, another limits applied
 867   //Reconstructed jet energy dependence parameters
 868   //e_jet = a1+e_gamma b2.
 869   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
 870   fJetE1[0] = -5.75; fJetE1[1] = -4.1;
 871   fJetE2[0] = 1.005; fJetE2[1] = 1.05;
 872
 873   //Reconstructed sigma of jet energy dependence parameters
 874   //s_jet = a1+e_gamma b2.
 875   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
 876   fJetSigma1[0] = 2.65;   fJetSigma1[1] = 2.75;
 877   fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
 878
 879   //Background mean energy and RMS
 880   //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
 881   //Index 2-> (low pt jets)BKG > 0.5 GeV;
 882   //Index > 2, same for CTS conf
 883   fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
 884   fBkgMean[3] = 0.; fBkgMean[4] = 6.4;  fBkgMean[5] = 48.6;
 885   fBkgRMS[0]  = 0.; fBkgRMS[1]  = 7.5;  fBkgRMS[2]  = 22.0;
 886   fBkgRMS[3]  = 0.; fBkgRMS[4]  = 5.4;  fBkgRMS[5]  = 13.2;
 887
 888   //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
 889   //limits for monoenergetic jets.
 890   //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
 891   //Index 2-> (low pt jets) BKG > 0.5 GeV;
 892   //Index > 2, same for CTS conf
 893
 894   fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
 895   fJetXMin1[3] =-2.0  ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1  ;
 896   fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
 897   fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
 898   fJetXMax1[0] =-3.8  ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6  ;
 899   fJetXMax1[3] =-2.7  ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7  ;
 900   fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
 901   fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
 902
 903
 904   //Different cones and pt thresholds to construct the jet
 905
 906   fJetCone        = 0.3  ;
 907   fJetPtThreshold = 0.5   ;
 908   fJetPtThresPbPb = 2.   ;
 909   fJetNCone       = 4    ;
 910   fJetNPt         = 4    ;
 911   fJetCones[0]    = 0.2  ; fJetNameCones[0]   = "02" ;
 912   fJetCones[1]    = 0.3  ; fJetNameCones[1]   = "03" ;
 913   fJetCones[2]    = 0.4  ; fJetNameCones[2]   = "04" ;
 914   fJetCones[2]    = 0.5  ; fJetNameCones[2]   = "05" ;
 915
 916   fJetPtThres[0]  = 0.0  ; fJetNamePtThres[0] = "00" ;
 917   fJetPtThres[1]  = 0.5  ; fJetNamePtThres[1] = "05" ;
 918   fJetPtThres[2]  = 1.0  ; fJetNamePtThres[2] = "10" ;
 919   fJetPtThres[3]  = 2.0  ; fJetNamePtThres[3] = "20" ;
 920 }
 921
 922 //__________________________________________________________________________________________________
 923 Bool_t AliAnaParticleJetLeadingConeCorrelation::IsJetSelected(Double_t ptTrig, Double_t ptjet) const
 924 {
 925   //Given the pt of the jet and the trigger particle, select the jet or not
 926   //3 options, fSelect=0 accepts all, fSelect=1 selects jets depending on a
 927   //function energy dependent and fSelect=2 selects on simple fixed cuts
 928
 929   if(ptjet == 0) return kFALSE;
 930
 931   Double_t rat = ptTrig / ptjet ;
 932
 933   //###############################################################
 934   if(fSelect == 0)
 935     return kTRUE; //Accept all jets, no restriction
 936   //###############################################################
 937   else if(fSelect == 1){
 938     //Check if the energy of the reconstructed jet is within an energy window
 939     //WARNING: to be rechecked, don't remember what all the steps mean
 940     Double_t par[6];
 941     Double_t xmax[2];
 942     Double_t xmin[2];
 943
 944     Int_t iCTS = 0;
 945     if(fJetsOnlyInCTS)
 946       iCTS = 3 ;
 947
 948     if(!fPbPb){
 949       //Phythia alone, jets with pt_th > 0.2, r = 0.3
 950       par[0] = fJetE1[0]; par[1] = fJetE2[0];
 951       //Energy of the jet peak
 952       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
 953       par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
 954       //Sigma  of the jet peak
 955       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
 956       par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
 957       //Parameters reserved for PbPb bkg.
 958       xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
 959       xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
 960       //Factor that multiplies sigma to obtain the best limits,
 961       //by observation, of mono jet ratios (ptjet/ptTrig)
 962       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
 963
 964     }
 965     else{
 966       if(ptTrig > fPtTriggerSelectionCut){
 967         //Phythia +PbPb with  pt_th > 2 GeV/c, r = 0.3
 968         par[0] = fJetE1[0]; par[1] = fJetE2[0];
 969         //Energy of the jet peak, same as in pp
 970         //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
 971         par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
 972         //Sigma  of the jet peak, same as in pp
 973         //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
 974         par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
 975         //Mean value and RMS of PbPb Bkg
 976         xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
 977         xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
 978         //Factor that multiplies sigma to obtain the best limits,
 979         //by observation, of mono jet ratios (ptjet/ptTrig) mixed with PbPb Bkg,
 980         //pt_th > 2 GeV, r = 0.3
 981         //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
 982
 983       }
 984       else{
 985         //Phythia + PbPb with  pt_th > 0.5 GeV/c, r = 0.3
 986         par[0] = fJetE1[1]; par[1] = fJetE2[1];
 987         //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
 988         //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
 989         par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
 990         //Sigma  of the jet peak, pt_th > 2 GeV/c, r = 0.3
 991         //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
 992         par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
 993         //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
 994         xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
 995         xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
 996         //Factor that multiplies sigma to obtain the best limits,
 997         //by observation, of mono jet ratios (ptjet/ptTrig) mixed with PbPb Bkg,
 998         //pt_th > 2 GeV, r = 0.3
 999         //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1000
1001       }//If low pt jet in bkg
1002     }//if Bkg
1003
1004     //Calculate minimum and maximum limits of the jet ratio.
1005     Double_t min = CalculateJetRatioLimit(ptTrig, par, xmin);
1006     Double_t max = CalculateJetRatioLimit(ptTrig, par, xmax);
1007
1008     AliDebug(3,Form("Jet selection? : Limits min %2.3f, max %2.3f,  pt_jet %2.3f,  pt_gamma %2.3f, pt_jet / pt_gamma %2.3f",min,max,ptjet,ptTrig,rat));
1009
1010     if(( min < rat ) && ( max > ptjet/rat))
1011       return kTRUE;
1012     else
1013       return kFALSE;
1014   }//fSelect = 1
1015   //###############################################################
1016   else if(fSelect == 2)
1017   {
1018     //Simple selection
1019     if(!fJetsOnlyInCTS)
1020     {
1021       if((rat <  fJetRatioMaxCut) && (rat > fJetRatioMinCut )) return kTRUE;
1022     }
1023     else
1024     {
1025       if((rat <  fJetCTSRatioMaxCut) && (rat > fJetCTSRatioMinCut )) return kTRUE;
1026     }
1027   }// fSelect = 2
1028   //###############################################################
1029   else
1030   {
1031     AliWarning(")Jet selection option larger than 2, DON'T SELECT JETS");
1032     return kFALSE ;
1033   }
1034
1035   return kFALSE;
1036
1037 }
1038
1039 //___________________________________________________________________
1040 Bool_t AliAnaParticleJetLeadingConeCorrelation::IsParticleInJetCone(Double_t eta, Double_t phi, Double_t etal,  Double_t phil)
1041 const {
1042   //Check if the particle is inside the cone defined by the leading particle
1043   //WARNING: To be rechecked
1044
1045   if(phi < 0) phi+=TMath::TwoPi();
1046   if(phil < 0) phil+=TMath::TwoPi();
1047   Double_t  rad = 10000 + fJetCone;
1048
1049   if(TMath::Abs(phi-phil) <= (TMath::TwoPi() - fJetCone))
1050     rad = TMath::Sqrt(TMath::Power(eta-etal,2)+TMath::Power(phi-phil,2));
1051   else{
1052     if(phi-phil > TMath::TwoPi() - fJetCone)
1053       rad = TMath::Sqrt(TMath::Power(eta-etal,2)+TMath::Power((phi-TMath::TwoPi())-phil,2));
1054     if(phi-phil < -(TMath::TwoPi() - fJetCone))
1055       rad = TMath::Sqrt(TMath::Power(eta-etal,2)+TMath::Power((phi+TMath::TwoPi())-phil,2));
1056   }
1057
1058   if(rad < fJetCone) return kTRUE ;
1059   else return kFALSE ;
1060
1061 }
1062
1063 //__________________________________________________________________
1064 void  AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD()
1065 {
1066   //Particle-Hadron Correlation Analysis, fill AODs
1067
1068   if(!GetInputAODBranch())
1069   {
1070     AliFatal(Form("No input particles in AOD with name branch < %s >",GetInputAODName().Data()));
1071     return;// coverity
1072   }
1073
1074   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
1075     AliFatal(Form("Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s>",
1076                   GetInputAODBranch()->GetClass()->GetName()));
1077
1078
1079   AliDebug(1,"Begin jet leading cone  correlation analysis, fill AODs");
1080   AliDebug(1,Form("In particle branch aod entries %d", GetInputAODBranch()->GetEntriesFast()));
1081   AliDebug(1,Form("In CTS aod entries %d",             GetCTSTracks()     ->GetEntriesFast()));
1082   AliDebug(1,Form("In EMCAL aod entries %d",           GetEMCALClusters() ->GetEntriesFast()));
1083
1084   //Loop on stored AOD particles, trigger
1085   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1086   for(Int_t iaod = 0; iaod < naod ; iaod++)
1087   {
1088     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1089
1090     //  printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - Trigger : pt %3.2f, phi %2.2f, eta %2.2f\n",particle->Pt(), particle->Phi(), particle->Eta());
1091
1092     //Search leading particles in CTS and EMCAL
1093     if(GetLeadingParticle(particle))
1094     {
1095       //Construct the jet around the leading, Fill AOD jet particle list, select jet
1096       //and fill AOD with jet and background
1097       MakeAODJet(particle);
1098
1099     }//Leading found
1100   }//AOD trigger particle loop
1101
1102   AliDebug(1,"End of jet leading cone analysis, fill AODs");
1103
1104 }
1105
1106 //_________________________________________________________________________
1107 void  AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms()
1108 {
1109
1110   //Particle-Hadron Correlation Analysis, fill histograms
1111
1112   if(!GetInputAODBranch())
1113   {
1114     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >",
1115            GetInputAODName().Data());
1116     return;
1117   }
1118
1119   AliDebug(1,"Begin jet leading cone  correlation analysis, fill histograms");
1120   AliDebug(1,Form("In particle branch aod entries %d", GetInputAODBranch()->GetEntriesFast()));
1121   AliDebug(1,Form("In CTS aod entries %d",             GetCTSTracks()     ->GetEntriesFast()));
1122   AliDebug(1,Form("In EMCAL aod entries %d",           GetEMCALClusters() ->GetEntriesFast()));
1123
1124   //Loop on stored AOD particles, trigger
1125   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1126   for(Int_t iaod = 0; iaod < naod ; iaod++)
1127   {
1128     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1129
1130     if(OnlyIsolated() && !particle->IsIsolated()) continue;
1131
1132     Double_t pt  = particle->Pt();
1133     Double_t phi = particle->Phi();
1134     Double_t eta = particle->Eta();
1135
1136     //Get leading particle, fill histograms
1137     fLeadingMom = particle->GetLeading();
1138     Int_t det   = particle->GetLeadingDetector();
1139
1140     if(det > 0 && fLeadingMom.Pt() > 0)
1141     {
1142       Double_t ptL = fLeadingMom.Pt();
1143       Double_t phiL = fLeadingMom.Phi();
1144       if(phiL < 0 ) phiL+=TMath::TwoPi();
1145       Double_t etaL = fLeadingMom.Eta();
1146
1147       AliDebug(1,Form("Trigger with pt %3.2f, phi %2.2f, eta %2.2f", pt, phi, eta));
1148
1149       if(det == kCTS)
1150       {
1151         fhChargedLeadingPt->Fill(pt,ptL);
1152         fhChargedLeadingPhi->Fill(pt,phiL);
1153         fhChargedLeadingEta->Fill(pt,etaL);
1154         fhChargedLeadingDeltaPt->Fill(pt,pt-ptL);
1155         fhChargedLeadingDeltaPhi->Fill(pt,TMath::Abs(phi-phiL));
1156         fhChargedLeadingDeltaEta->Fill(pt,eta-etaL);
1157         fhChargedLeadingRatioPt->Fill(pt,ptL/pt);
1158         fhChargedLeadingXi->Fill(pt,TMath::Log(pt/ptL));
1159         if(pt > 30) fhChargedLeadingDeltaPhiRatioPt30->Fill(TMath::Abs(phi-phiL),ptL/pt);
1160         if(pt > 50) fhChargedLeadingDeltaPhiRatioPt50->Fill(TMath::Abs(phi-phiL),ptL/pt);
1161       }
1162       else if(det== kEMCAL)
1163       {
1164         fhNeutralLeadingPt->Fill(pt,ptL);
1165         fhNeutralLeadingPhi->Fill(pt,phiL);
1166         fhNeutralLeadingEta->Fill(pt,etaL);
1167         fhNeutralLeadingDeltaPt->Fill(pt,pt-ptL);
1168         fhNeutralLeadingDeltaPhi->Fill(pt,TMath::Abs(phi-phiL));
1169         fhNeutralLeadingDeltaEta->Fill(pt,eta-etaL);
1170         fhNeutralLeadingRatioPt->Fill(pt,ptL/pt);
1171         fhNeutralLeadingXi->Fill(pt,TMath::Log(pt/ptL));
1172         if(pt > 30) fhNeutralLeadingDeltaPhiRatioPt30->Fill(TMath::Abs(phi-phiL),ptL/pt);
1173         if(pt > 50) fhNeutralLeadingDeltaPhiRatioPt50->Fill(TMath::Abs(phi-phiL),ptL/pt);
1174       }
1175
1176       //Fill Jet histograms
1177
1178       if(!fSeveralConeAndPtCuts)
1179       {//just fill histograms
1180         if(!fReMakeJet)
1181         {
1182           fJetMom=particle->GetCorrelatedJet();
1183           fBkgMom=particle->GetCorrelatedBackground();
1184         }
1185         else  MakeJetFromAOD(particle);
1186
1187         if(fJetMom.Pt() > 0)
1188         {//Jet was found
1189           FillJetHistos(particle, fJetMom,"Jet","");
1190           FillJetHistos(particle, fBkgMom,"Bkg","");
1191         }
1192       }
1193       else if(fSeveralConeAndPtCuts)
1194       {
1195         for(Int_t icone = 0; icone<fJetNCone; icone++)
1196         {
1197           fJetCone=fJetCones[icone];
1198           for(Int_t ipt = 0; ipt<fJetNPt;ipt++)
1199           {
1200             TString lastname ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt];
1201             fJetPtThreshold=fJetPtThres[ipt];
1202
1203             MakeJetFromAOD(particle);
1204
1205             if(fJetMom.Pt() > 0)
1206             {//Jet was found
1207               FillJetHistos(particle, fJetMom,"Jet",lastname);
1208               FillJetHistos(particle, fBkgMom,"Bkg",lastname);
1209             }
1210           }//icone
1211         }//ipt
1212       }//fSeveralConeAndPtCuts
1213     }//Leading
1214   }//AOD trigger particle loop
1215
1216   AliDebug(1,"End of jet leading cone analysis, fill histograms");
1217
1218 }
1219
1220 //_______________________________________________________________________________________________
1221 void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorrelation *particle)
1222 {
1223   //Fill the jet with the particles around the leading particle with
1224   //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and
1225   //fill aod with found information
1226
1227   Double_t ptTrig   = particle->Pt();
1228   Double_t phiTrig  = particle->Phi();
1229   Double_t phil     = fLeadingMom.Phi();
1230   if(phil<0) phil+=TMath::TwoPi();
1231   Double_t etal     = fLeadingMom.Eta();
1232
1233   //Different pt cut for jet particles in different collisions systems
1234   Float_t ptcut = fJetPtThreshold;
1235   if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
1236
1237   //Add charged particles to jet if they are in cone around the leading particle
1238   if(!GetCTSTracks())
1239   {
1240     AliFatal("Cannot construct jets without tracks, STOP analysis");
1241     return;
1242   }
1243
1244   //Fill jet with tracks
1245   //Initialize reference arrays that will contain jet and background tracks
1246   TObjArray * reftracks     = new TObjArray;
1247   TObjArray * reftracksbkg  = new TObjArray;
1248
1249   for(Int_t ipr = 0;ipr < (GetCTSTracks())->GetEntriesFast() ; ipr ++ ){
1250     AliVTrack* track = (AliVTrack *)((GetCTSTracks())->At(ipr)) ;
1251     fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1252
1253     //Particles in jet
1254     if(IsParticleInJetCone(fTrackVector.Eta(), fTrackVector.Phi(), etal, phil)){
1255
1256       reftracks->Add(track);
1257
1258       if(fTrackVector.Pt() > ptcut )
1259       {
1260         fJetConstMom.SetVect(fTrackVector);
1261         fJetMom+=fJetConstMom;
1262       }
1263     }
1264
1265     //Background around (phi_gamma-pi, eta_leading)
1266     else if(IsParticleInJetCone(fTrackVector.Eta(),fTrackVector.Phi(),etal, phiTrig)) {
1267
1268       reftracksbkg->Add(track);
1269
1270       if(fTrackVector.Pt() > ptcut ){
1271         fJetConstMom.SetVect(fTrackVector);
1272         fBkgMom+=fJetConstMom;
1273       }
1274     }
1275   }//Track loop
1276
1277   //Add referenced tracks to AOD
1278   if(reftracks->GetEntriesFast() > 0 )
1279   {
1280     reftracks->SetName(Form("%sTracks",GetAODObjArrayName().Data()));
1281     particle->AddObjArray(reftracks);
1282   }
1283   else  AliDebug(2,"No tracks in jet cone");
1284
1285   if(reftracksbkg->GetEntriesFast() > 0 )
1286   {
1287     reftracksbkg->SetName(Form("%sTracksBkg",GetAODObjArrayName().Data()));
1288     particle->AddObjArray(reftracksbkg);
1289   }
1290   else AliDebug(1,"No background tracks in jet cone");
1291
1292   //Add neutral particles to jet
1293   //Initialize reference arrays that will contain jet and background tracks
1294   TObjArray * refclusters     = new TObjArray;
1295   TObjArray * refclustersbkg  = new TObjArray;
1296   if(!fJetsOnlyInCTS && GetEMCALClusters()){
1297
1298     //Get vertex for photon momentum calculation
1299     Double_t vertex[]  = {0,0,0} ; //vertex
1300     if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
1301     {
1302       GetReader()->GetVertex(vertex);
1303       //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
1304     }
1305
1306     for(Int_t iclus = 0;iclus < (GetEMCALClusters())->GetEntriesFast() ; iclus ++ )
1307     {
1308       AliVCluster * calo = (AliVCluster *) (GetEMCALClusters()->At(iclus)) ;
1309
1310       //Cluster selection, not charged
1311       if(IsTrackMatched(calo,GetReader()->GetInputEvent())) continue ;
1312
1313       //Get Momentum vector,
1314       calo->GetMomentum(fJetConstMom,vertex) ;//Assume that come from vertex in straight line
1315
1316       //Particles in jet
1317       if(IsParticleInJetCone(fJetConstMom.Eta(),fJetConstMom.Phi(), etal, phil)){
1318
1319         refclusters->Add(calo);
1320
1321         if(fJetConstMom.Pt() > ptcut )  fJetMom+=fJetConstMom;
1322       }
1323       //Background around (phi_gamma-pi, eta_leading)
1324       else if(IsParticleInJetCone(fJetConstMom.Eta(),fJetConstMom.Phi(),etal, phiTrig)){
1325
1326
1327         refclustersbkg->Add(calo);
1328
1329         if(fJetConstMom.Pt() > ptcut )  fBkgMom+=fJetConstMom;
1330       }
1331     }//cluster loop
1332   }//jets with neutral particles
1333
1334   //Add referenced clusters to AOD
1335   if(refclusters->GetEntriesFast() > 0 )
1336   {
1337     refclusters->SetName(Form("%sClusters",GetAODObjArrayName().Data()));
1338     particle->AddObjArray(refclusters);
1339   }
1340   else  AliDebug(2,"No clusters in jet cone");
1341
1342   if(refclustersbkg->GetEntriesFast() > 0 )
1343   {
1344     refclustersbkg->SetName(Form("%sClustersBkg",GetAODObjArrayName().Data()));
1345     particle->AddObjArray(refclustersbkg);
1346   }
1347   else AliDebug(1,"No background clusters in jet cone");
1348
1349   //If there is any jet found, select after some criteria and
1350   //and fill AOD with corresponding TLorentzVector kinematics
1351   if(IsJetSelected(particle->Pt(), fJetMom.Pt()))
1352   {
1353     particle->SetCorrelatedJet(fJetMom);
1354     particle->SetCorrelatedBackground(fBkgMom);
1355     AliDebug(1,Form("Found jet: Trigger  pt %2.3f, Jet pt %2.3f, Bkg pt %2.3f",ptTrig,fJetMom.Pt(),fBkgMom.Pt()));
1356   }
1357
1358 }
1359
1360 //______________________________________________________________________________________________________
1361 void AliAnaParticleJetLeadingConeCorrelation::MakeJetFromAOD(AliAODPWG4ParticleCorrelation *particle)
1362 {
1363   //Fill the jet with the particles around the leading particle with
1364   //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and
1365   //fill aod tlorentzvectors with jet and bakcground found
1366
1367   Double_t ptTrig  = particle->Pt();
1368   Double_t phiTrig = particle->Phi();
1369   Double_t etal    = fLeadingMom.Eta();
1370   Double_t phil    = fLeadingMom.Phi();
1371   if(phil < 0) phil+=TMath::TwoPi();
1372
1373   TObjArray * refclusters    = particle->GetObjArray(Form("Clusters%s"   ,GetAODObjArrayName().Data()));
1374   TObjArray * reftracks      = particle->GetObjArray(Form("Tracks%s"     ,GetAODObjArrayName().Data()));
1375   TObjArray * refclustersbkg = particle->GetObjArray(Form("ClustersBkg%s",GetAODObjArrayName().Data()));
1376   TObjArray * reftracksbkg   = particle->GetObjArray(Form("TracksBkg%s"  ,GetAODObjArrayName().Data()));
1377
1378   //Different pt cut for jet particles in different collisions systems
1379   Float_t ptcut = fJetPtThreshold;
1380   if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
1381
1382   //Fill jet with tracks
1383   //Particles in jet
1384   if(reftracks)
1385   {
1386     for(Int_t ipr = 0;ipr < reftracks->GetEntriesFast() ; ipr ++ )
1387     {
1388       AliVTrack* track = (AliVTrack *) reftracks->At(ipr) ;
1389       fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1390       Float_t phi = fTrackVector.Phi();
1391       if(phi < 0) phi+=TMath::TwoPi();
1392       if(fTrackVector.Pt() > ptcut && IsParticleInJetCone(fTrackVector.Eta(), phi, etal, phil) )
1393       {
1394         fJetConstMom.SetVect(fTrackVector);
1395         fJetMom+=fJetConstMom;
1396       }
1397     }//jet Track loop
1398   }
1399   //Particles in background
1400   if(reftracksbkg){
1401     for(Int_t ipr = 0;ipr < reftracksbkg->GetEntriesFast() ; ipr ++ )
1402     {
1403       AliVTrack* track = (AliVTrack *) reftracksbkg->At(ipr) ;
1404       fTrackVector.SetXYZ(track->Px(),track->Py(),track->Pz());
1405       if(fTrackVector.Pt() > ptcut && IsParticleInJetCone(fTrackVector.Eta(),fTrackVector.Phi(),etal, phiTrig) )
1406       {
1407         fJetConstMom.SetVect(fTrackVector);
1408         fBkgMom+=fJetConstMom;
1409       }
1410     }//background Track loop
1411   }
1412
1413   //Add neutral particles to jet
1414   if(!fJetsOnlyInCTS && refclusters)
1415   {
1416     //Get vertex for photon momentum calculation
1417     Double_t vertex[]  = {0,0,0} ; //vertex
1418     if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
1419     {
1420       GetReader()->GetVertex(vertex);
1421     }
1422
1423     //Loop on jet particles
1424     if(refclusters){
1425       for(Int_t iclus = 0;iclus < refclusters->GetEntriesFast() ; iclus ++ )
1426       {
1427         AliVCluster * calo = (AliVCluster *) refclusters->At(iclus) ;
1428
1429         calo->GetMomentum(fJetConstMom,vertex) ;//Assume that come from vertex in straight line
1430
1431         if(fJetConstMom.Pt() > ptcut && IsParticleInJetCone(fJetConstMom.Eta(),fJetConstMom.Phi(), etal, phil)) fJetMom+=fJetConstMom;
1432       }//jet cluster loop
1433     }
1434
1435     //Loop on background particles
1436     if(refclustersbkg)
1437     {
1438       for(Int_t iclus = 0;iclus < refclustersbkg->GetEntriesFast() ; iclus ++ )
1439       {
1440         AliVCluster * calo = (AliVCluster *) refclustersbkg->At(iclus) ;
1441
1442         calo->GetMomentum(fJetConstMom,vertex) ;//Assume that come from vertex in straight line
1443
1444         if( fJetConstMom.Pt() > ptcut && IsParticleInJetCone(fJetConstMom.Eta(),fJetConstMom.Phi(),etal, phiTrig)) fBkgMom+=fJetConstMom;
1445       }//background cluster loop
1446     }
1447   }//clusters in jet
1448
1449   //If there is any jet found, leave jet and bkg as they are,
1450   //if not set them to 0.
1451   if(!IsJetSelected(particle->Pt(), fJetMom.Pt()))
1452   {
1453     fJetMom.SetPxPyPzE(0.,0.,0.,0.);
1454     fBkgMom.SetPxPyPzE(0.,0.,0.,0.);
1455   }
1456   else AliDebug(1,Form("Found jet: Trigger  pt %2.3f, Jet pt %2.3f, Bkg pt %2.3f",ptTrig,fJetMom.Pt(),fBkgMom.Pt()));
1457
1458 }
1459
1460 //__________________________________________________________________
1461 void AliAnaParticleJetLeadingConeCorrelation::Print(const Option_t * opt) const
1462 {
1463
1464   //Print some relevant parameters set for the analysis
1465   if(! opt)
1466     return;
1467
1468   printf("**** Print  %s %s ****\n", GetName(), GetTitle() ) ;
1469   AliAnaCaloTrackCorrBaseClass::Print(" ");
1470
1471   if(fJetsOnlyInCTS)printf("Jets reconstructed in CTS \n");
1472   else printf("Jets reconstructed in CTS+EMCAL \n");
1473
1474   if(fPbPb) printf("PbPb events, pT cut in jet cone energy reconstruction %2.1f \n", fJetPtThreshold);
1475   else printf("pp events, pT cut in jet cone energy reconstruction %2.1f \n", fJetPtThresPbPb);
1476
1477   printf("If pT of trigger < %2.3f, select jets as in pp? \n", fPtTriggerSelectionCut);
1478
1479   printf("Phi gamma-Leading        <     %3.2f\n", fDeltaPhiMaxCut) ;
1480   printf("Phi gamma-Leading        >     %3.2f\n", fDeltaPhiMinCut) ;
1481   printf("pT Leading / pT Trigger  <     %3.2f\n", fLeadingRatioMaxCut) ;
1482   printf("pT Leading / pT Trigger  >     %3.2f\n", fLeadingRatioMinCut) ;
1483
1484   if(fSelect == 2){
1485     printf("pT Jet / pT Gamma                     <    %3.2f\n", fJetRatioMaxCut) ;
1486     printf("pT Jet / pT Gamma                     >    %3.2f\n", fJetRatioMinCut) ;
1487     printf("pT Jet (Only CTS)/ pT Trigger   <    %3.2f\n", fJetCTSRatioMaxCut) ;
1488     printf("pT Jet (Only CTS)/ pT Trigger   >    %3.2f\n", fJetCTSRatioMinCut) ;
1489   }
1490   else if(fSelect == 0)
1491     printf("Accept all reconstructed jets \n") ;
1492   else   if(fSelect == 1)
1493     printf("Accept jets depending on trigger energy \n") ;
1494   else
1495     printf("Wrong jet selection option:   %d \n", fSelect) ;
1496
1497   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
1498
1499 }