| f9cea31c |
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 | /* $Id$ */ |
| 16 | |
| 17 | /* History of cvs commits: |
| 18 | * |
| 19 | * $Log$ |
| 2a1d8a29 |
20 | * Revision 1.1 2007/01/23 17:17:29 schutz |
| 21 | * New Gamma package |
| 22 | * |
| f9cea31c |
23 | * |
| 24 | */ |
| 25 | |
| 26 | //_________________________________________________________________________ |
| 27 | // Class for the analysis of gamma-jet correlations |
| 28 | // Basically it seaches for a prompt photon in the Calorimeters acceptance, |
| 29 | // if so we construct a jet around the highest pt particle in the opposite |
| 30 | // side in azimuth, inside the Central Tracking System (ITS+TPC) and |
| 31 | // EMCAL acceptances (only when PHOS detects the gamma). First the leading |
| 32 | // particle and then the jet have to fullfill several conditions |
| 33 | // (energy, direction ..) to be accepted. Then the fragmentation function |
| 34 | // of this jet is constructed |
| 35 | // Class created from old AliPHOSGammaJet |
| 36 | // (see AliRoot versions previous Release 4-09) |
| 37 | // |
| 38 | //*-- Author: Gustavo Conesa (LNF-INFN) |
| 39 | ////////////////////////////////////////////////////////////////////////////// |
| 40 | |
| 41 | |
| 42 | // --- ROOT system --- |
| 43 | |
| 44 | #include <TFile.h> |
| 45 | #include <TParticle.h> |
| 46 | #include <TH2.h> |
| 47 | |
| 48 | #include "AliAnaGammaJet.h" |
| 49 | #include "AliESD.h" |
| 50 | #include "AliESDtrack.h" |
| 51 | #include "AliESDCaloCluster.h" |
| 52 | #include "Riostream.h" |
| 53 | #include "AliLog.h" |
| 54 | |
| 55 | ClassImp(AliAnaGammaJet) |
| 56 | |
| 57 | //____________________________________________________________________________ |
| 58 | AliAnaGammaJet::AliAnaGammaJet(const char *name) : |
| 59 | AliAnaGammaDirect(name), |
| 60 | fSeveralConeAndPtCuts(0), |
| 61 | fPbPb(kFALSE), |
| 62 | fJetsOnlyInCTS(0), |
| 63 | fEtaEMCALCut(0.),fPhiMaxCut(0.), |
| 64 | fPhiMinCut(0.), |
| 65 | fInvMassMaxCut(0.), fInvMassMinCut(0.), |
| 66 | fJetCTSRatioMaxCut(0.), |
| 67 | fJetCTSRatioMinCut(0.), fJetRatioMaxCut(0.), |
| 68 | fJetRatioMinCut(0.), fNCone(0), |
| 69 | fNPt(0), fCone(0), fPtThreshold(0), |
| 70 | fPtJetSelectionCut(0.0), |
| 71 | fAngleMaxParam(), fSelect(0) |
| 72 | { |
| 73 | |
| 74 | // ctor |
| 75 | fAngleMaxParam.Set(4) ; |
| 76 | fAngleMaxParam.Reset(0.); |
| 77 | |
| 78 | for(Int_t i = 0; i<10; i++){ |
| 79 | fCones[i] = 0.0 ; |
| 80 | fNameCones[i] = "" ; |
| 81 | fPtThres[i] = 0.0 ; |
| 82 | fNamePtThres[i] = "" ; |
| 83 | if( i < 6 ){ |
| 84 | fJetXMin1[i] = 0.0 ; |
| 85 | fJetXMin2[i] = 0.0 ; |
| 86 | fJetXMax1[i] = 0.0 ; |
| 87 | fJetXMax2[i] = 0.0 ; |
| 88 | fBkgMean[i] = 0.0 ; |
| 89 | fBkgRMS[i] = 0.0 ; |
| 90 | if( i < 2 ){ |
| 91 | fJetE1[i] = 0.0 ; |
| 92 | fJetE2[i] = 0.0 ; |
| 93 | fJetSigma1[i] = 0.0 ; |
| 94 | fJetSigma2[i] = 0.0 ; |
| 95 | fPhiEMCALCut[i] = 0.0 ; |
| 96 | } |
| 97 | } |
| 98 | } |
| 99 | |
| 100 | TList * list = gDirectory->GetListOfKeys() ; |
| 101 | TIter next(list) ; |
| 102 | TH2F * h = 0 ; |
| 103 | Int_t index ; |
| 104 | for (index = 0 ; index < list->GetSize()-1 ; index++) { |
| 105 | //-1 to avoid GammaJet Task |
| 106 | h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ; |
| 107 | fOutputContainer->Add(h) ; |
| 108 | } |
| 109 | |
| 110 | // Input slot #0 works with an Ntuple |
| 111 | DefineInput(0, TChain::Class()); |
| 112 | // Output slot #0 writes into a TH1 container |
| 113 | DefineOutput(0, TObjArray::Class()) ; |
| 114 | |
| 115 | } |
| 116 | |
| 117 | |
| 118 | //____________________________________________________________________________ |
| 119 | AliAnaGammaJet::AliAnaGammaJet(const AliAnaGammaJet & gj) : |
| 120 | AliAnaGammaDirect(gj), |
| 121 | fSeveralConeAndPtCuts(gj.fSeveralConeAndPtCuts), |
| 122 | fPbPb(gj.fPbPb), fJetsOnlyInCTS(gj.fJetsOnlyInCTS), |
| 123 | fEtaEMCALCut(gj.fEtaEMCALCut), |
| 124 | fPhiMaxCut(gj.fPhiMaxCut), fPhiMinCut(gj.fPhiMinCut), |
| 125 | fInvMassMaxCut(gj.fInvMassMaxCut), fInvMassMinCut(gj.fInvMassMinCut), |
| 126 | fRatioMinCut(gj.fRatioMinCut), |
| 127 | fJetCTSRatioMaxCut(gj.fJetCTSRatioMaxCut), |
| 128 | fJetCTSRatioMinCut(gj.fJetCTSRatioMinCut), fJetRatioMaxCut(gj.fJetRatioMaxCut), |
| 129 | fJetRatioMinCut(gj.fJetRatioMinCut), fNCone(gj.fNCone), |
| 130 | fNPt(gj.fNPt), fCone(gj.fCone), fPtThreshold(gj.fPtThreshold), |
| 131 | fPtJetSelectionCut(gj.fPtJetSelectionCut), |
| 132 | fOutputContainer(0), fAngleMaxParam(gj.fAngleMaxParam), |
| 133 | fSelect(gj.fSelect) |
| 134 | { |
| 135 | // cpy ctor |
| 136 | SetName (gj.GetName()) ; |
| 137 | SetTitle(gj.GetTitle()) ; |
| 138 | |
| 139 | for(Int_t i = 0; i<10; i++){ |
| 140 | fCones[i] = gj.fCones[i] ; |
| 141 | fNameCones[i] = gj.fNameCones[i] ; |
| 142 | fPtThres[i] = gj.fPtThres[i] ; |
| 143 | fNamePtThres[i] = gj.fNamePtThres[i] ; |
| 144 | if( i < 6 ){ |
| 145 | fJetXMin1[i] = gj.fJetXMin1[i] ; |
| 146 | fJetXMin2[i] = gj.fJetXMin2[i] ; |
| 147 | fJetXMax1[i] = gj.fJetXMax1[i] ; |
| 148 | fJetXMax2[i] = gj.fJetXMax2[i] ; |
| 149 | fBkgMean[i] = gj.fBkgMean[i] ; |
| 150 | fBkgRMS[i] = gj.fBkgRMS[i] ; |
| 151 | if( i < 2 ){ |
| 152 | fJetE1[i] = gj.fJetE1[i] ; |
| 153 | fJetE2[i] = gj.fJetE2[i] ; |
| 154 | fJetSigma1[i] = gj.fJetSigma1[i] ; |
| 155 | fJetSigma2[i] = gj.fJetSigma2[i] ; |
| 156 | fPhiEMCALCut[i] = gj.fPhiEMCALCut[i] ; |
| 157 | } |
| 158 | } |
| 159 | } |
| 160 | } |
| 161 | |
| 162 | //____________________________________________________________________________ |
| 163 | AliAnaGammaJet::~AliAnaGammaJet() |
| 164 | { |
| 165 | // Remove all pointers |
| 166 | fOutputContainer->Clear() ; |
| 167 | delete fOutputContainer ; |
| 168 | |
| 169 | delete fhChargeRatio ; |
| 170 | delete fhPi0Ratio ; |
| 171 | delete fhDeltaPhiCharge ; |
| 172 | delete fhDeltaPhiPi0 ; |
| 173 | delete fhDeltaEtaCharge ; |
| 174 | delete fhDeltaEtaPi0 ; |
| 175 | delete fhAnglePair ; |
| 176 | delete fhAnglePairAccepted ; |
| 177 | delete fhAnglePairNoCut ; |
| 178 | delete fhAnglePairLeadingCut ; |
| 179 | delete fhAnglePairAngleCut ; |
| 180 | delete fhAnglePairAllCut ; |
| 181 | delete fhAnglePairLeading ; |
| 182 | delete fhInvMassPairNoCut ; |
| 183 | delete fhInvMassPairLeadingCut ; |
| 184 | delete fhInvMassPairAngleCut ; |
| 185 | delete fhInvMassPairAllCut ; |
| 186 | delete fhInvMassPairLeading ; |
| 187 | delete fhNBkg ; |
| 188 | delete fhNLeading ; |
| 189 | delete fhNJet ; |
| 190 | delete fhJetRatio ; |
| 191 | delete fhJetPt ; |
| 192 | delete fhBkgRatio ; |
| 193 | delete fhBkgPt ; |
| 194 | delete fhJetFragment ; |
| 195 | delete fhBkgFragment ; |
| 196 | delete fhJetPtDist ; |
| 197 | delete fhBkgPtDist ; |
| 198 | |
| 199 | delete [] fhJetRatios; |
| 200 | delete [] fhJetPts; |
| 201 | delete [] fhBkgRatios; |
| 202 | delete [] fhBkgPts; |
| 203 | |
| 204 | delete [] fhNLeadings; |
| 205 | delete [] fhNJets; |
| 206 | delete [] fhNBkgs; |
| 207 | |
| 208 | delete [] fhJetFragments; |
| 209 | delete [] fhBkgFragments; |
| 210 | delete [] fhJetPtDists; |
| 211 | delete [] fhBkgPtDists; |
| 212 | |
| 213 | } |
| 214 | |
| 215 | //____________________________________________________________________________ |
| 216 | Double_t AliAnaGammaJet::CalculateJetRatioLimit(const Double_t ptg, |
| 217 | const Double_t *par, |
| 218 | const Double_t *x) { |
| 219 | |
| 220 | //Parametrized cut for the energy of the jet. |
| 221 | |
| 222 | Double_t epp = par[0] + par[1] * ptg ; |
| 223 | Double_t spp = par[2] + par[3] * ptg ; |
| 224 | Double_t f = x[0] + x[1] * ptg ; |
| 225 | Double_t epb = epp + par[4] ; |
| 226 | Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ; |
| 227 | Double_t rat = (epb - spb * f) / ptg ; |
| 228 | //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f); |
| 229 | //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat); |
| 230 | return rat ; |
| 231 | } |
| 232 | |
| 2a1d8a29 |
233 | //______________________________________________________________________________ |
| 234 | void AliAnaGammaJet::ConnectInputData(const Option_t*) |
| f9cea31c |
235 | { |
| 2a1d8a29 |
236 | // Initialisation of branch container and histograms |
| 237 | AliAnaGammaDirect::ConnectInputData(""); |
| f9cea31c |
238 | |
| f9cea31c |
239 | } |
| 240 | |
| 2a1d8a29 |
241 | //________________________________________________________________________ |
| 242 | void AliAnaGammaJet::CreateOutputObjects() |
| f9cea31c |
243 | { |
| f9cea31c |
244 | |
| 2a1d8a29 |
245 | // Init parameteres and create histograms to be saved in output file and |
| 246 | // stores them in fOutputContainer |
| 247 | InitParameters(); |
| 248 | AliAnaGammaDirect::CreateOutputObjects(); |
| f9cea31c |
249 | |
| 2a1d8a29 |
250 | fOutputContainer = new TObjArray(100) ; |
| 251 | |
| 252 | TObjArray * outputContainer =GetOutputContainer(); |
| 253 | for(Int_t i = 0; i < outputContainer->GetEntries(); i++ ) |
| 254 | fOutputContainer->Add(outputContainer->At(i)) ; |
| f9cea31c |
255 | |
| 2a1d8a29 |
256 | // |
| 257 | fhChargeRatio = new TH2F |
| 258 | ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}", |
| 259 | 120,0,120,120,0,1); |
| 260 | fhChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}"); |
| 261 | fhChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 262 | fOutputContainer->Add(fhChargeRatio) ; |
| f9cea31c |
263 | |
| 2a1d8a29 |
264 | fhDeltaPhiCharge = new TH2F |
| 265 | ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}", |
| 266 | 200,0,120,200,0,6.4); |
| 267 | fhDeltaPhiCharge->SetYTitle("#Delta #phi"); |
| 268 | fhDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 269 | fOutputContainer->Add(fhDeltaPhiCharge) ; |
| 270 | |
| 271 | fhDeltaEtaCharge = new TH2F |
| 272 | ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}", |
| 273 | 200,0,120,200,-2,2); |
| 274 | fhDeltaEtaCharge->SetYTitle("#Delta #eta"); |
| 275 | fhDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 276 | fOutputContainer->Add(fhDeltaEtaCharge) ; |
| 277 | |
| 278 | // |
| 279 | if(!fJetsOnlyInCTS){ |
| 280 | fhPi0Ratio = new TH2F |
| 281 | ("Pi0Ratio","p_{T leading #pi^{0}} /p_{T #gamma} vs p_{T #gamma}", |
| 282 | 120,0,120,120,0,1); |
| 283 | fhPi0Ratio->SetYTitle("p_{T lead #pi^{0}} /p_{T #gamma}"); |
| 284 | fhPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 285 | fOutputContainer->Add(fhPi0Ratio) ; |
| f9cea31c |
286 | |
| 2a1d8a29 |
287 | fhDeltaPhiPi0 = new TH2F |
| 288 | ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}", |
| 289 | 200,0,120,200,0,6.4); |
| 290 | fhDeltaPhiPi0->SetYTitle("#Delta #phi"); |
| 291 | fhDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 292 | fOutputContainer->Add(fhDeltaPhiPi0) ; |
| 293 | |
| 294 | fhDeltaEtaPi0 = new TH2F |
| 295 | ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}", |
| 296 | 200,0,120,200,-2,2); |
| 297 | fhDeltaEtaPi0->SetYTitle("#Delta #eta"); |
| 298 | fhDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 299 | fOutputContainer->Add(fhDeltaEtaPi0) ; |
| 300 | |
| 301 | // |
| 302 | fhAnglePair = new TH2F |
| 303 | ("AnglePair", |
| 304 | "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}", |
| 305 | 200,0,50,200,0,0.2); |
| 306 | fhAnglePair->SetYTitle("Angle (rad)"); |
| 307 | fhAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)"); |
| 308 | fOutputContainer->Add(fhAnglePair) ; |
| 309 | |
| 310 | fhAnglePairAccepted = new TH2F |
| 311 | ("AnglePairAccepted", |
| 312 | "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}, both #gamma in eta<0.7, inside window", |
| 313 | 200,0,50,200,0,0.2); |
| 314 | fhAnglePairAccepted->SetYTitle("Angle (rad)"); |
| 315 | fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)"); |
| 316 | fOutputContainer->Add(fhAnglePairAccepted) ; |
| 317 | |
| 318 | fhAnglePairNoCut = new TH2F |
| 319 | ("AnglePairNoCut", |
| 320 | "Angle between all #gamma pair vs p_{T #pi^{0}}",200,0,50,200,0,0.2); |
| 321 | fhAnglePairNoCut->SetYTitle("Angle (rad)"); |
| 322 | fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)"); |
| 323 | fOutputContainer->Add(fhAnglePairNoCut) ; |
| 324 | |
| 325 | fhAnglePairLeadingCut = new TH2F |
| 326 | ("AnglePairLeadingCut", |
| 327 | "Angle between all #gamma pair that have a good phi and pt vs p_{T #pi^{0}}", |
| 328 | 200,0,50,200,0,0.2); |
| 329 | fhAnglePairLeadingCut->SetYTitle("Angle (rad)"); |
| 330 | fhAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)"); |
| 331 | fOutputContainer->Add(fhAnglePairLeadingCut) ; |
| 332 | |
| 333 | fhAnglePairAngleCut = new TH2F |
| 334 | ("AnglePairAngleCut", |
| 335 | "Angle between all #gamma pair (angle + leading cut) vs p_{T #pi^{0}}" |
| 336 | ,200,0,50,200,0,0.2); |
| 337 | fhAnglePairAngleCut->SetYTitle("Angle (rad)"); |
| 338 | fhAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)"); |
| 339 | fOutputContainer->Add(fhAnglePairAngleCut) ; |
| 340 | |
| 341 | fhAnglePairAllCut = new TH2F |
| 342 | ("AnglePairAllCut", |
| 343 | "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T #pi^{0}}" |
| 344 | ,200,0,50,200,0,0.2); |
| 345 | fhAnglePairAllCut->SetYTitle("Angle (rad)"); |
| 346 | fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)"); |
| 347 | fOutputContainer->Add(fhAnglePairAllCut) ; |
| 348 | |
| 349 | fhAnglePairLeading = new TH2F |
| 350 | ("AnglePairLeading", |
| 351 | "Angle between all #gamma pair finally selected vs p_{T #pi^{0}}", |
| 352 | 200,0,50,200,0,0.2); |
| 353 | fhAnglePairLeading->SetYTitle("Angle (rad)"); |
| 354 | fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)"); |
| 355 | fOutputContainer->Add(fhAnglePairLeading) ; |
| 356 | |
| 357 | // |
| 358 | fhInvMassPairNoCut = new TH2F |
| 359 | ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}", |
| 360 | 120,0,120,360,0,0.5); |
| 361 | fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})"); |
| 362 | fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 363 | fOutputContainer->Add(fhInvMassPairNoCut) ; |
| 364 | |
| 365 | fhInvMassPairLeadingCut = new TH2F |
| 366 | ("InvMassPairLeadingCut", |
| 367 | "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}", |
| 368 | 120,0,120,360,0,0.5); |
| 369 | fhInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})"); |
| 370 | fhInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 371 | fOutputContainer->Add(fhInvMassPairLeadingCut) ; |
| 372 | |
| 373 | fhInvMassPairAngleCut = new TH2F |
| 374 | ("InvMassPairAngleCut", |
| 375 | "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}", |
| 376 | 120,0,120,360,0,0.5); |
| 377 | fhInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})"); |
| 378 | fhInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 379 | fOutputContainer->Add(fhInvMassPairAngleCut) ; |
| 380 | |
| 381 | fhInvMassPairAllCut = new TH2F |
| 382 | ("InvMassPairAllCut", |
| 383 | "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}", |
| 384 | 120,0,120,360,0,0.5); |
| 385 | fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})"); |
| 386 | fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 387 | fOutputContainer->Add(fhInvMassPairAllCut) ; |
| 388 | |
| 389 | fhInvMassPairLeading = new TH2F |
| 390 | ("InvMassPairLeading", |
| 391 | "Invariant Mass of #gamma pair selected vs p_{T #gamma}", |
| 392 | 120,0,120,360,0,0.5); |
| 393 | fhInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})"); |
| 394 | fhInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 395 | fOutputContainer->Add(fhInvMassPairLeading) ; |
| f9cea31c |
396 | } |
| 397 | |
| f9cea31c |
398 | |
| 2a1d8a29 |
399 | if(!fSeveralConeAndPtCuts){ |
| 400 | |
| 401 | //Count |
| 402 | fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000); |
| 403 | fhNBkg->SetYTitle("counts"); |
| 404 | fhNBkg->SetXTitle("N"); |
| 405 | fOutputContainer->Add(fhNBkg) ; |
| 406 | |
| 407 | fhNLeading = new TH2F |
| 408 | ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120); |
| 409 | fhNLeading->SetYTitle("p_{T charge} (GeV/c)"); |
| 410 | fhNLeading->SetXTitle("p_{T #gamma}(GeV/c)"); |
| 411 | fOutputContainer->Add(fhNLeading) ; |
| 412 | |
| 413 | fhNJet = new TH1F("NJet","Accepted jets",240,0,120); |
| 414 | fhNJet->SetYTitle("N"); |
| 415 | fhNJet->SetXTitle("p_{T #gamma}(GeV/c)"); |
| 416 | fOutputContainer->Add(fhNJet) ; |
| 417 | |
| 418 | //Ratios and Pt dist of reconstructed (not selected) jets |
| 419 | //Jet |
| 420 | fhJetRatio = new TH2F |
| 421 | ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}", |
| 422 | 240,0,120,200,0,10); |
| 423 | fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); |
| 424 | fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 425 | fOutputContainer->Add(fhJetRatio) ; |
| 426 | |
| 427 | fhJetPt = new TH2F |
| 428 | ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200); |
| 429 | fhJetPt->SetYTitle("p_{T jet}"); |
| 430 | fhJetPt->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 431 | fOutputContainer->Add(fhJetPt) ; |
| 432 | |
| 433 | //Bkg |
| 434 | |
| 435 | fhBkgRatio = new TH2F |
| 436 | ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}", |
| 437 | 240,0,120,200,0,10); |
| 438 | fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}"); |
| 439 | fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 440 | fOutputContainer->Add(fhBkgRatio) ; |
| 441 | |
| 442 | fhBkgPt = new TH2F |
| 443 | ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200); |
| 444 | fhBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}"); |
| 445 | fhBkgPt->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 446 | fOutputContainer->Add(fhBkgPt) ; |
| 447 | |
| 448 | //Jet Distributions |
| 449 | |
| 450 | fhJetFragment = |
| 451 | new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}", |
| 452 | 240,0.,120.,1000,0.,1.2); |
| 453 | fhJetFragment->SetYTitle("x_{T}"); |
| 454 | fhJetFragment->SetXTitle("p_{T #gamma}"); |
| 455 | fOutputContainer->Add(fhJetFragment) ; |
| 456 | |
| 457 | fhBkgFragment = new TH2F |
| 458 | ("BkgFragment","x = p_{T i charged}/p_{T #gamma}", |
| 459 | 240,0.,120.,1000,0.,1.2); |
| 460 | fhBkgFragment->SetYTitle("x_{T}"); |
| 461 | fhBkgFragment->SetXTitle("p_{T #gamma}"); |
| 462 | fOutputContainer->Add(fhBkgFragment) ; |
| 463 | |
| 464 | fhJetPtDist = |
| 465 | new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); |
| 466 | fhJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 467 | fOutputContainer->Add(fhJetPtDist) ; |
| 468 | |
| 469 | fhBkgPtDist = new TH2F |
| 470 | ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); |
| 471 | fhBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 472 | fOutputContainer->Add(fhBkgPtDist) ; |
| f9cea31c |
473 | |
| 2a1d8a29 |
474 | } |
| 475 | else{ |
| 476 | //If we want to study the jet for different cones and pt |
| 477 | |
| 478 | for(Int_t icone = 0; icone<fNCone; icone++){ |
| 479 | for(Int_t ipt = 0; ipt<fNPt;ipt++){ |
| 480 | |
| 481 | //Jet |
| 482 | |
| 483 | fhJetRatios[icone][ipt] = new TH2F |
| 484 | ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 485 | "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" |
| 486 | +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", |
| 487 | 240,0,120,200,0,10); |
| 488 | fhJetRatios[icone][ipt]-> |
| 489 | SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); |
| 490 | fhJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 491 | fOutputContainer->Add(fhJetRatios[icone][ipt]) ; |
| 492 | |
| 493 | |
| 494 | fhJetPts[icone][ipt] = new TH2F |
| 495 | ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 496 | "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" |
| 497 | +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", |
| 498 | 240,0,120,400,0,200); |
| 499 | fhJetPts[icone][ipt]-> |
| 500 | SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); |
| 501 | fhJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 502 | fOutputContainer->Add(fhJetPts[icone][ipt]) ; |
| 503 | |
| 504 | //Bkg |
| 505 | fhBkgRatios[icone][ipt] = new TH2F |
| 506 | ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 507 | "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" |
| 508 | +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", |
| 509 | 240,0,120,200,0,10); |
| 510 | fhBkgRatios[icone][ipt]-> |
| 511 | SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}"); |
| 512 | fhBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 513 | fOutputContainer->Add(fhBkgRatios[icone][ipt]) ; |
| 514 | |
| 515 | fhBkgPts[icone][ipt] = new TH2F |
| 516 | ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 517 | "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone =" |
| 518 | +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c", |
| 519 | 240,0,120,400,0,200); |
| 520 | fhBkgPts[icone][ipt]-> |
| 521 | SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}"); |
| 522 | fhBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 523 | fOutputContainer->Add(fhBkgPts[icone][ipt]) ; |
| 524 | |
| 525 | //Counts |
| 526 | fhNBkgs[icone][ipt] = new TH1F |
| 527 | ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 528 | "bkg multiplicity cone ="+fNameCones[icone]+", pt>" |
| 529 | +fNamePtThres[ipt]+" GeV/c",9000,0,9000); |
| 530 | fhNBkgs[icone][ipt]->SetYTitle("counts"); |
| 531 | fhNBkgs[icone][ipt]->SetXTitle("N"); |
| 532 | fOutputContainer->Add(fhNBkgs[icone][ipt]) ; |
| 533 | |
| 534 | fhNLeadings[icone][ipt] = new TH2F |
| 535 | ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 536 | "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>" |
| 537 | +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120); |
| 538 | fhNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)"); |
| 539 | fhNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)"); |
| 540 | fOutputContainer->Add(fhNLeadings[icone][ipt]) ; |
| 541 | |
| 542 | fhNJets[icone][ipt] = new TH1F |
| 543 | ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 544 | "Number of neutral jets, cone ="+fNameCones[icone]+", pt>" |
| 545 | +fNamePtThres[ipt]+" GeV/c",120,0,120); |
| 546 | fhNJets[icone][ipt]->SetYTitle("N"); |
| 547 | fhNJets[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)"); |
| 548 | fOutputContainer->Add(fhNJets[icone][ipt]) ; |
| 549 | |
| 550 | //Fragmentation Function |
| 551 | fhJetFragments[icone][ipt] = new TH2F |
| 552 | ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 553 | "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" |
| 554 | +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); |
| 555 | fhJetFragments[icone][ipt]->SetYTitle("x_{T}"); |
| 556 | fhJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}"); |
| 557 | fOutputContainer->Add(fhJetFragments[icone][ipt]) ; |
| 558 | |
| 559 | fhBkgFragments[icone][ipt] = new TH2F |
| 560 | ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 561 | "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" |
| 562 | +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); |
| 563 | fhBkgFragments[icone][ipt]->SetYTitle("x_{T}"); |
| 564 | fhBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}"); |
| 565 | fOutputContainer->Add(fhBkgFragments[icone][ipt]) ; |
| 566 | |
| 567 | //Jet particle distribution |
| 568 | |
| 569 | fhJetPtDists[icone][ipt] = new TH2F |
| 570 | ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 571 | "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+ |
| 572 | " GeV/c",120,0.,120.,120,0.,120.); |
| 573 | fhJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 574 | fOutputContainer->Add(fhJetPtDists[icone][ipt]) ; |
| 575 | |
| 576 | fhBkgPtDists[icone][ipt] = new TH2F |
| 577 | ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], |
| 578 | "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+ |
| 579 | " GeV/c",120,0.,120.,120,0.,120.); |
| 580 | fhBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)"); |
| 581 | fOutputContainer->Add(fhBkgPtDists[icone][ipt]) ; |
| 582 | |
| 583 | }//ipt |
| 584 | } //icone |
| 585 | }//If we want to study any cone or pt threshold |
| 586 | } |
| f9cea31c |
587 | |
| 588 | |
| 2a1d8a29 |
589 | //____________________________________________________________________________ |
| 590 | void AliAnaGammaJet::Exec(Option_t *) |
| 591 | { |
| 592 | |
| 593 | // Processing of one event |
| f9cea31c |
594 | |
| 2a1d8a29 |
595 | //Get ESDs |
| 596 | Long64_t entry = GetChain()->GetReadEntry() ; |
| 597 | |
| 598 | if (!GetESD()) { |
| 599 | AliError("fESD is not connected to the input!") ; |
| 600 | return ; |
| 601 | } |
| 602 | |
| 603 | if (GetPrintInfo()) |
| 604 | AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(GetChain()))->GetFile()->GetName(), entry)) ; |
| 605 | |
| 606 | //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis. |
| 607 | |
| 608 | TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet) |
| 609 | TClonesArray * plCTS = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC) |
| 610 | TClonesArray * plNe = new TClonesArray("TParticle",1000); // All particles measured in Jet Calorimeter |
| 611 | TClonesArray * plCalo = new TClonesArray("TParticle",1000); // All particles measured in Prompt Gamma calorimeter |
| 612 | |
| 613 | |
| 614 | TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma |
| 615 | TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle |
| 616 | |
| 617 | Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter |
| 618 | |
| 619 | //Fill lists with photons, neutral particles and charged particles |
| 620 | //look for the highest energy photon in the event inside fCalorimeter |
| 621 | AliDebug(2, "Fill particle lists, get prompt gamma"); |
| 622 | |
| 623 | //Fill particle lists |
| 624 | |
| 625 | if(GetCalorimeter() == "PHOS") |
| 626 | CreateParticleList(particleList, plCTS,plNe,plCalo); |
| 627 | else if(GetCalorimeter() == "EMCAL") |
| 628 | CreateParticleList(particleList, plCTS,plCalo,plNe); |
| 629 | else |
| 630 | AliError("No calorimeter selected"); |
| 631 | |
| 632 | //Search highest energy prompt gamma in calorimeter |
| 633 | GetPromptGamma(plCalo, plCTS, pGamma, iIsInPHOSorEMCAL) ; |
| 634 | |
| 635 | AliDebug(1, Form("Is Gamma in %s? %d",GetCalorimeter().Data(),iIsInPHOSorEMCAL)); |
| 636 | AliDebug(3,Form("Charged list entries %d, Neutral list entries %d, %s list entries %d", |
| 637 | plCTS->GetEntries(),plNe->GetEntries(), GetCalorimeter().Data(), plCalo->GetEntries())); |
| f9cea31c |
638 | |
| 2a1d8a29 |
639 | //If there is any prompt photon in fCalorimeter, |
| 640 | //search jet leading particle |
| 641 | |
| 642 | if(iIsInPHOSorEMCAL){ |
| 643 | if (GetPrintInfo()) |
| 644 | AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ; |
| 645 | |
| 646 | AliDebug(2, "Get Leading Particles, Make jet"); |
| f9cea31c |
647 | |
| 2a1d8a29 |
648 | //Search leading particles in CTS and EMCAL |
| 649 | if(GetLeadingParticle(plCTS, plNe, pGamma, pLeading)){ |
| 650 | if (GetPrintInfo()) |
| 651 | AliInfo(Form("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ; |
| f9cea31c |
652 | |
| 2a1d8a29 |
653 | //Search Jet |
| 654 | if(!fSeveralConeAndPtCuts) |
| 655 | MakeJet(particleList,pGamma,pLeading,""); |
| 656 | else{ |
| 657 | for(Int_t icone = 0; icone<fNCone; icone++) { |
| 658 | for(Int_t ipt = 0; ipt<fNPt;ipt++) { |
| 659 | TString lastname ="Cone"+ fNameCones[icone]+"Pt"+ fNamePtThres[ipt]; |
| 660 | MakeJet(particleList, pGamma, pLeading,lastname); |
| 661 | }//icone |
| 662 | }//ipt |
| 663 | }//fSeveralConeAndPtCuts |
| 664 | }//Leading |
| 665 | }//Gamma in Calo |
| 666 | |
| 667 | AliDebug(2, "End of analysis, delete pointers"); |
| f9cea31c |
668 | |
| 2a1d8a29 |
669 | particleList->Delete() ; |
| 670 | plCTS->Delete() ; |
| 671 | plNe->Delete() ; |
| 672 | plCalo->Delete() ; |
| 673 | pLeading->Delete(); |
| 674 | pGamma->Delete(); |
| f9cea31c |
675 | |
| 2a1d8a29 |
676 | delete plNe ; |
| 677 | delete plCalo ; |
| 678 | delete plCTS ; |
| 679 | delete particleList ; |
| 680 | // delete pLeading; |
| 681 | // delete pGamma; |
| f9cea31c |
682 | |
| 2a1d8a29 |
683 | PostData(0, fOutputContainer); |
| 684 | } |
| f9cea31c |
685 | |
| 2a1d8a29 |
686 | //____________________________________________________________________________ |
| 687 | void AliAnaGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname) |
| 688 | { |
| 689 | //Fill histograms wth jet fragmentation |
| 690 | //and number of selected jets and leading particles |
| 691 | //and the background multiplicity |
| 692 | TParticle * particle = 0 ; |
| 693 | Int_t ipr = 0; |
| 694 | Float_t charge = 0; |
| f9cea31c |
695 | |
| 2a1d8a29 |
696 | TIter next(pl) ; |
| 697 | while ( (particle = dynamic_cast<TParticle*>(next())) ) { |
| 698 | ipr++ ; |
| 699 | Double_t pt = particle->Pt(); |
| f9cea31c |
700 | |
| 2a1d8a29 |
701 | charge = TDatabasePDG::Instance() |
| 702 | ->GetParticle(particle->GetPdgCode())->Charge(); |
| 703 | if(charge != 0){//Only jet Charged particles |
| 704 | dynamic_cast<TH2F*> |
| 705 | (fOutputContainer->FindObject(type+"Fragment"+lastname)) |
| 706 | ->Fill(ptg,pt/ptg); |
| 707 | dynamic_cast<TH2F*> |
| 708 | (fOutputContainer->FindObject(type+"PtDist"+lastname)) |
| 709 | ->Fill(ptg,pt); |
| 710 | }//charged |
| f9cea31c |
711 | |
| 2a1d8a29 |
712 | }//while |
| f9cea31c |
713 | |
| 2a1d8a29 |
714 | if(type == "Bkg") |
| 715 | dynamic_cast<TH1F*> |
| 716 | (fOutputContainer->FindObject("NBkg"+lastname)) |
| 717 | ->Fill(ipr); |
| 718 | else{ |
| 719 | dynamic_cast<TH1F*> |
| 720 | (fOutputContainer->FindObject("NJet"+lastname))-> |
| 721 | Fill(ptg); |
| 722 | dynamic_cast<TH2F*> |
| 723 | (fOutputContainer->FindObject("NLeading"+lastname)) |
| 724 | ->Fill(ptg,ptl); |
| 725 | } |
| f9cea31c |
726 | |
| 2a1d8a29 |
727 | } |
| f9cea31c |
728 | |
| 2a1d8a29 |
729 | //____________________________________________________________________________ |
| 730 | void AliAnaGammaJet::GetLeadingCharge(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) const |
| 731 | { |
| 732 | //Search for the charged particle with highest with |
| 733 | //Phi=Phi_gamma-Pi and pT=0.1E_gamma |
| 734 | Double_t pt = -100.; |
| 735 | Double_t phi = -100.; |
| 736 | |
| 737 | for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){ |
| 738 | |
| 739 | TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ; |
| 740 | |
| 741 | Double_t ptl = particle->Pt(); |
| 742 | Double_t rat = ptl/pGamma->Pt() ; |
| 743 | Double_t phil = particle->Phi() ; |
| 744 | Double_t phig = pGamma->Phi(); |
| 745 | |
| 746 | //Selection within angular and energy limits |
| 747 | if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) && |
| 748 | (rat > fRatioMinCut) && (rat < fRatioMaxCut) && (ptl > pt)) { |
| 749 | phi = phil ; |
| 750 | pt = ptl ; |
| 751 | pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy()); |
| 752 | AliDebug(4,Form("Charge in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, particle->Eta(), phi,rat)) ; |
| 753 | } |
| 754 | } |
| 755 | |
| 756 | AliDebug(3,Form("Leading in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, pLeading->Eta(), phi,pt/pGamma->Pt())) ; |
| f9cea31c |
757 | |
| 758 | } |
| 759 | |
| 2a1d8a29 |
760 | |
| f9cea31c |
761 | //____________________________________________________________________________ |
| 2a1d8a29 |
762 | void AliAnaGammaJet::GetLeadingPi0(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) |
| 763 | { |
| f9cea31c |
764 | |
| 2a1d8a29 |
765 | //Search for the neutral pion with highest with |
| 766 | //Phi=Phi_gamma-Pi and pT=0.1E_gamma |
| 767 | Double_t pt = -100.; |
| 768 | Double_t phi = -100.; |
| 769 | Double_t ptl = -100.; |
| 770 | Double_t rat = -100.; |
| 771 | Double_t phil = -100. ; |
| 772 | Double_t ptg = pGamma->Pt(); |
| 773 | Double_t phig = pGamma->Phi(); |
| 774 | |
| 775 | TIter next(pl); |
| 776 | TParticle * particle = 0; |
| f9cea31c |
777 | |
| 2a1d8a29 |
778 | Int_t iPrimary = -1; |
| 779 | TLorentzVector gammai,gammaj; |
| 780 | Double_t angle = 0., e = 0., invmass = 0.; |
| 781 | Double_t anglef = 0., ef = 0., invmassf = 0.; |
| 782 | Int_t ksPdg = 0; |
| 783 | Int_t jPrimary=-1; |
| 784 | |
| 785 | while ( (particle = (TParticle*)next()) ) { |
| 786 | iPrimary++; |
| f9cea31c |
787 | |
| 2a1d8a29 |
788 | ksPdg = particle->GetPdgCode(); |
| 789 | ptl = particle->Pt(); |
| 790 | if(ksPdg == 111){ //2 gamma overlapped, found with PID |
| 791 | rat = ptl/ptg ; |
| 792 | phil = particle->Phi() ; |
| 793 | //Selection within angular and energy limits |
| 794 | if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) && |
| 795 | ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){ |
| 796 | phi = phil ; |
| 797 | pt = ptl ; |
| 798 | pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy()); |
| 799 | AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ; |
| 800 | }// cuts |
| 801 | }// pdg = 111 |
| 802 | else if(ksPdg == 22){//1 gamma |
| 803 | //Search the photon companion in case it comes from a Pi0 decay |
| 804 | //Apply several cuts to select the good pair |
| 805 | particle->Momentum(gammai); |
| 806 | jPrimary=-1; |
| 807 | TIter next2(pl); |
| 808 | while ( (particle = (TParticle*)next2()) ) { |
| 809 | jPrimary++; |
| 810 | if(jPrimary>iPrimary){ |
| 811 | ksPdg = particle->GetPdgCode(); |
| 812 | |
| 813 | if(ksPdg == 22){ |
| 814 | particle->Momentum(gammaj); |
| 815 | //Info("GetLeadingPi0","Egammai %f, Egammaj %f", |
| 816 | //gammai.Pt(),gammaj.Pt()); |
| 817 | ptl = (gammai+gammaj).Pt(); |
| 818 | phil = (gammai+gammaj).Phi(); |
| 819 | if(phil < 0) |
| 820 | phil+=TMath::TwoPi(); |
| 821 | rat = ptl/ptg ; |
| 822 | invmass = (gammai+gammaj).M(); |
| 823 | angle = gammaj.Angle(gammai.Vect()); |
| 824 | //Info("GetLeadingPi0","Angle %f", angle); |
| 825 | e = (gammai+gammaj).E(); |
| 826 | //Fill histograms with no cuts applied. |
| 827 | fhAnglePairNoCut->Fill(e,angle); |
| 828 | fhInvMassPairNoCut->Fill(ptg,invmass); |
| 829 | //First cut on the energy and azimuth of the pair |
| 830 | if((rat > fRatioMinCut) && (rat < fRatioMaxCut) && |
| 831 | ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){ |
| 832 | |
| 833 | fhAnglePairLeadingCut->Fill(e,angle); |
| 834 | fhInvMassPairLeadingCut->Fill(ptg,invmass); |
| 835 | //Second cut on the aperture of the pair |
| 836 | if(IsAngleInWindow(angle,e)){ |
| 837 | fhAnglePairAngleCut->Fill(e,angle); |
| 838 | fhInvMassPairAngleCut->Fill(ptg,invmass); |
| 839 | |
| 840 | //Info("GetLeadingPi0","InvMass %f", invmass); |
| 841 | //Third cut on the invariant mass of the pair |
| 842 | if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){ |
| 843 | fhInvMassPairAllCut->Fill(ptg,invmass); |
| 844 | fhAnglePairAllCut->Fill(e,angle); |
| 845 | if(ptl > pt ){ |
| 846 | pt = ptl; |
| 847 | phi = phil ; |
| 848 | ef = e ; |
| 849 | anglef = angle ; |
| 850 | invmassf = invmass ; |
| 851 | pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy()); |
| 852 | AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ; |
| 853 | } |
| 854 | }//cuts |
| 855 | }//(invmass>0.125) && (invmass<0.145) |
| 856 | }//gammaj.Angle(gammai.Vect())<0.04 |
| 857 | }//if pdg = 22 |
| 858 | }//iprimary<jprimary |
| 859 | }//while |
| 860 | }// if pdg = 22 |
| 861 | // } |
| 862 | }//while |
| 863 | |
| 864 | if(ef > 0.){//Final pi0 found, highest pair energy, fill histograms |
| 865 | fhInvMassPairLeading->Fill(ptg,invmassf); |
| 866 | fhAnglePairLeading->Fill(ef,anglef); |
| f9cea31c |
867 | } |
| 868 | |
| 2a1d8a29 |
869 | AliDebug(3,Form("Leading EMCAL: pt %f eta %f phi %f pt/Eg %f \n", pLeading->Pt(), pLeading->Eta(), pLeading->Phi(), pLeading->Pt()/pGamma->Pt())) ; |
| 870 | } |
| 871 | |
| 872 | //____________________________________________________________________________ |
| 873 | Bool_t AliAnaGammaJet::GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe, |
| 874 | TParticle * pGamma, TParticle * pLeading) |
| 875 | { |
| 876 | //Search Charged or Neutral leading particle, select the highest one. |
| f9cea31c |
877 | |
| 2a1d8a29 |
878 | TParticle * pLeadingCh = new TParticle(); |
| 879 | TParticle * pLeadingPi0 = new TParticle(); |
| 880 | |
| 881 | Double_t ptg = pGamma->Pt(); |
| 882 | Double_t phig = pGamma->Phi(); |
| 883 | Double_t etag = pGamma->Eta(); |
| 884 | |
| 885 | if(GetCalorimeter() == "PHOS" && !fJetsOnlyInCTS) |
| 886 | { |
| 887 | AliDebug(3, "GetLeadingPi0"); |
| 888 | GetLeadingPi0 (plNe, pGamma, pLeadingPi0) ; |
| 889 | AliDebug(3, "GetLeadingCharge"); |
| 890 | GetLeadingCharge(plCTS, pGamma, pLeadingCh) ; |
| 891 | |
| 892 | Double_t ptch = pLeadingCh->Pt(); |
| 893 | Double_t phich = pLeadingCh->Phi(); |
| 894 | Double_t etach = pLeadingCh->Eta(); |
| 895 | Double_t ptpi = pLeadingPi0->Pt(); |
| 896 | Double_t phipi = pLeadingPi0->Phi(); |
| 897 | Double_t etapi = pLeadingPi0->Eta(); |
| 898 | |
| 899 | //Is leading cone inside EMCAL acceptance? |
| 900 | |
| 901 | Bool_t insidech = kFALSE ; |
| 902 | if((phich - fCone) > fPhiEMCALCut[0] && |
| 903 | (phich + fCone) < fPhiEMCALCut[1] && |
| 904 | (etach-fCone) < fEtaEMCALCut ) |
| 905 | insidech = kTRUE ; |
| 906 | |
| 907 | Bool_t insidepi = kFALSE ; |
| 908 | if((phipi - fCone) > fPhiEMCALCut[0] && |
| 909 | (phipi + fCone) < fPhiEMCALCut[1] && |
| 910 | (etapi-fCone) < fEtaEMCALCut ) |
| 911 | insidepi = kTRUE ; |
| 912 | |
| 913 | AliDebug(2,Form("Leading: charged pt %f, pi0 pt %f",ptch,ptpi)) ; |
| 914 | |
| 915 | if (ptch > 0 || ptpi > 0) |
| 916 | { |
| 917 | if(insidech && (ptch > ptpi)) |
| 918 | { |
| 919 | if (GetPrintInfo()) |
| 920 | AliInfo("Leading found in CTS"); |
| 921 | pLeading->SetMomentum(pLeadingCh->Px(),pLeadingCh->Py(),pLeadingCh->Pz(),pLeadingCh->Energy()); |
| 922 | AliDebug(3,Form("Final leading found in CTS, pt %f, phi %f, eta %f",ptch,phich,etach)) ; |
| 923 | fhChargeRatio->Fill(ptg,ptch/pGamma->Pt()); |
| 924 | fhDeltaPhiCharge->Fill(ptg,pGamma->Phi()-phich); |
| 925 | fhDeltaEtaCharge->Fill(ptg,pGamma->Eta()-etach); |
| 926 | return 1; |
| 927 | } |
| 928 | |
| 929 | else if((ptpi > ptch) && insidepi) |
| 930 | { |
| 931 | if (GetPrintInfo()) |
| 932 | AliInfo("Leading found in EMCAL"); |
| 933 | pLeading->SetMomentum(pLeadingPi0->Px(),pLeadingPi0->Py(),pLeadingPi0->Pz(),pLeadingPi0->Energy()); |
| 934 | AliDebug(3,Form("Final leading found in EMCAL, pt %f, phi %f, eta %f",ptpi,phipi,etapi)) ; |
| 935 | fhPi0Ratio ->Fill(ptg,ptpi/ptg); |
| 936 | fhDeltaPhiPi0->Fill(ptg,phig-phipi); |
| 937 | fhDeltaEtaPi0->Fill(ptg,etag-etapi); |
| 938 | return 1; |
| 939 | } |
| 940 | |
| 941 | else{ |
| 942 | AliDebug(3,"NO LEADING PARTICLE FOUND");} |
| 943 | return 0; |
| 944 | } |
| 945 | else{ |
| 946 | AliDebug(3,"NO LEADING PARTICLE FOUND"); |
| 947 | return 0; |
| 948 | } |
| 949 | } |
| 950 | |
| 951 | else |
| 952 | { |
| 953 | //No calorimeter present for Leading particle detection |
| 954 | GetLeadingCharge(plCTS, pGamma, pLeading) ; |
| 955 | Double_t ptch = pLeading->Pt(); |
| 956 | Double_t phich = pLeading->Phi(); |
| 957 | Double_t etach = pLeading->Eta(); |
| 958 | if(ptch > 0){ |
| 959 | fhChargeRatio->Fill(ptg,ptch/ptg); |
| 960 | fhDeltaPhiCharge->Fill(ptg,phig-phich); |
| 961 | fhDeltaEtaCharge->Fill(ptg,etag-etach); |
| 962 | AliDebug(3,Form("Leading found : pt %f, phi %f, eta %f",ptch,phich,etach)) ; |
| 963 | return 1; |
| 964 | } |
| 965 | else |
| 966 | { |
| 967 | AliDebug(3,"NO LEADING PARTICLE FOUND"); |
| 968 | return 0; |
| 969 | } |
| 970 | } |
| 971 | } |
| 972 | |
| 973 | //____________________________________________________________________________ |
| 974 | void AliAnaGammaJet::InitParameters() |
| 975 | { |
| 976 | // // Initialisation of branch container |
| 977 | //AliAnaGammaDirect::InitParameters(); |
| 978 | |
| 979 | //Initialize the parameters of the analysis. |
| 980 | //fCalorimeter="PHOS"; |
| 981 | fAngleMaxParam.Set(4) ; |
| 982 | fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4}; |
| 983 | fAngleMaxParam.AddAt(-0.25,1) ; |
| 984 | fAngleMaxParam.AddAt(0.025,2) ; |
| 985 | fAngleMaxParam.AddAt(-2e-4,3) ; |
| 986 | fSeveralConeAndPtCuts = kFALSE ; |
| 987 | //fPrintInfo = kTRUE; |
| 988 | fPbPb = kFALSE ; |
| 989 | fInvMassMaxCut = 0.16 ; |
| 990 | fInvMassMinCut = 0.11 ; |
| 991 | //fJetsOnlyInCTS = kTRUE ; |
| 992 | fEtaEMCALCut = 0.7 ; |
| 993 | fPhiEMCALCut[0] = 80. *TMath::Pi()/180.; |
| 994 | fPhiEMCALCut[1] = 190.*TMath::Pi()/180.; |
| 995 | fPhiMaxCut = 3.4 ; |
| 996 | fPhiMinCut = 2.9 ; |
| 997 | |
| 998 | //Jet selection parameters |
| 999 | //Fixed cut (old) |
| 1000 | fRatioMaxCut = 1.0 ; |
| 1001 | fRatioMinCut = 0.1 ; |
| 1002 | fJetRatioMaxCut = 1.2 ; |
| 1003 | fJetRatioMinCut = 0.3 ; |
| 1004 | fJetsOnlyInCTS = kFALSE ; |
| 1005 | fJetCTSRatioMaxCut = 1.2 ; |
| 1006 | fJetCTSRatioMinCut = 0.3 ; |
| 1007 | fSelect = 0 ; |
| 1008 | |
| 1009 | //Cut depending on gamma energy |
| 1010 | |
| 1011 | fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applied |
| 1012 | //Reconstructed jet energy dependence parameters |
| 1013 | //e_jet = a1+e_gamma b2. |
| 1014 | //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3 |
| 1015 | fJetE1[0] = -5.75; fJetE1[1] = -4.1; |
| 1016 | fJetE2[0] = 1.005; fJetE2[1] = 1.05; |
| f9cea31c |
1017 | |
| 2a1d8a29 |
1018 | //Reconstructed sigma of jet energy dependence parameters |
| 1019 | //s_jet = a1+e_gamma b2. |
| 1020 | //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3 |
| 1021 | fJetSigma1[0] = 2.65; fJetSigma1[1] = 2.75; |
| 1022 | fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033; |
| 1023 | |
| 1024 | //Background mean energy and RMS |
| 1025 | //Index 0-> No BKG; Index 1-> BKG > 2 GeV; |
| 1026 | //Index 2-> (low pt jets)BKG > 0.5 GeV; |
| 1027 | //Index > 2, same for CTS conf |
| 1028 | fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5; |
| 1029 | fBkgMean[3] = 0.; fBkgMean[4] = 6.4; fBkgMean[5] = 48.6; |
| 1030 | fBkgRMS[0] = 0.; fBkgRMS[1] = 7.5; fBkgRMS[2] = 22.0; |
| 1031 | fBkgRMS[3] = 0.; fBkgRMS[4] = 5.4; fBkgRMS[5] = 13.2; |
| 1032 | |
| 1033 | //Factor x of min/max = E -+ x * sigma. Obtained after selecting the |
| 1034 | //limits for monoenergetic jets. |
| 1035 | //Index 0-> No BKG; Index 1-> BKG > 2 GeV; |
| 1036 | //Index 2-> (low pt jets) BKG > 0.5 GeV; |
| 1037 | //Index > 2, same for CTS conf |
| 1038 | |
| 1039 | fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ; |
| 1040 | fJetXMin1[3] =-2.0 ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1 ; |
| 1041 | fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034; |
| 1042 | fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ; |
| 1043 | fJetXMax1[0] =-3.8 ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6 ; |
| 1044 | fJetXMax1[3] =-2.7 ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7 ; |
| 1045 | fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016; |
| 1046 | fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027; |
| 1047 | |
| 1048 | |
| 1049 | //Different cones and pt thresholds to construct the jet |
| 1050 | |
| 1051 | fCone = 0.3 ; |
| 1052 | fPtThreshold = 0. ; |
| 1053 | fNCone = 3 ; |
| 1054 | fNPt = 3 ; |
| 1055 | fCones[1] = 0.2 ; fNameCones[1] = "02" ; |
| 1056 | fPtThres[0] = 0.0 ; fNamePtThres[0] = "00" ; |
| 1057 | fCones[0] = 0.3 ; fNameCones[0] = "03" ; |
| 1058 | fPtThres[1] = 0.5 ; fNamePtThres[1] = "05" ; |
| 1059 | fCones[2] = 0.4 ; fNameCones[2] = "04" ; |
| 1060 | fPtThres[2] = 1.0 ; fNamePtThres[2] = "10" ; |
| 1061 | //Fill particle lists when PID is ok |
| 1062 | // fEMCALPID = kFALSE; |
| 1063 | // fPHOSPID = kFALSE; |
| 1064 | |
| 1065 | } |
| 1066 | |
| 1067 | //__________________________________________________________________________- |
| 1068 | Bool_t AliAnaGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) { |
| 1069 | //Check if the opening angle of the candidate pairs is inside |
| 1070 | //our selection windowd |
| 1071 | |
| 1072 | Bool_t result = kFALSE; |
| 1073 | Double_t mpi0 = 0.1349766; |
| 1074 | Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e) |
| 1075 | +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e; |
| 1076 | Double_t arg = (e*e-2*mpi0*mpi0)/(e*e); |
| 1077 | Double_t min = 100. ; |
| 1078 | if(arg>0.) |
| 1079 | min = TMath::ACos(arg); |
| 1080 | |
| 1081 | if((angle<max)&&(angle>=min)) |
| 1082 | result = kTRUE; |
| 1083 | |
| 1084 | return result; |
| 1085 | } |
| 1086 | |
| 1087 | //__________________________________________________________________________- |
| 1088 | Bool_t AliAnaGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj){ |
| 1089 | //Check if the energy of the reconstructed jet is within an energy window |
| 1090 | |
| 1091 | Double_t par[6]; |
| 1092 | Double_t xmax[2]; |
| 1093 | Double_t xmin[2]; |
| 1094 | |
| 1095 | Int_t iCTS = 0; |
| 1096 | if(fJetsOnlyInCTS) |
| 1097 | iCTS = 3 ; |
| 1098 | |
| 1099 | if(!fPbPb){ |
| 1100 | //Phythia alone, jets with pt_th > 0.2, r = 0.3 |
| 1101 | par[0] = fJetE1[0]; par[1] = fJetE2[0]; |
| 1102 | //Energy of the jet peak |
| 1103 | //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit |
| 1104 | par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0]; |
| 1105 | //Sigma of the jet peak |
| 1106 | //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit |
| 1107 | par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS]; |
| 1108 | //Parameters reserved for PbPb bkg. |
| 1109 | xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS]; |
| 1110 | xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS]; |
| 1111 | //Factor that multiplies sigma to obtain the best limits, |
| 1112 | //by observation, of mono jet ratios (ptjet/ptg) |
| 1113 | //X_jet = fJetX1[0]+fJetX2[0]*e_gamma |
| 1114 | |
| f9cea31c |
1115 | } |
| 1116 | else{ |
| 2a1d8a29 |
1117 | if(ptg > fPtJetSelectionCut){ |
| 1118 | //Phythia +PbPb with pt_th > 2 GeV/c, r = 0.3 |
| 1119 | par[0] = fJetE1[0]; par[1] = fJetE2[0]; |
| 1120 | //Energy of the jet peak, same as in pp |
| 1121 | //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit |
| 1122 | par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0]; |
| 1123 | //Sigma of the jet peak, same as in pp |
| 1124 | //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit |
| 1125 | par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS]; |
| 1126 | //Mean value and RMS of PbPb Bkg |
| 1127 | xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS]; |
| 1128 | xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS]; |
| 1129 | //Factor that multiplies sigma to obtain the best limits, |
| 1130 | //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, |
| 1131 | //pt_th > 2 GeV, r = 0.3 |
| 1132 | //X_jet = fJetX1[0]+fJetX2[0]*e_gamma |
| 1133 | |
| 1134 | } |
| 1135 | else{ |
| 1136 | //Phythia + PbPb with pt_th > 0.5 GeV/c, r = 0.3 |
| 1137 | par[0] = fJetE1[1]; par[1] = fJetE2[1]; |
| 1138 | //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3 |
| 1139 | //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit |
| 1140 | par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1]; |
| 1141 | //Sigma of the jet peak, pt_th > 2 GeV/c, r = 0.3 |
| 1142 | //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit |
| 1143 | par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS]; |
| 1144 | //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV. |
| 1145 | xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS]; |
| 1146 | xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS]; |
| 1147 | //Factor that multiplies sigma to obtain the best limits, |
| 1148 | //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, |
| 1149 | //pt_th > 2 GeV, r = 0.3 |
| 1150 | //X_jet = fJetX1[0]+fJetX2[0]*e_gamma |
| 1151 | |
| 1152 | }//If low pt jet in bkg |
| 1153 | }//if Bkg |
| 1154 | |
| 1155 | //Calculate minimum and maximum limits of the jet ratio. |
| 1156 | Double_t min = CalculateJetRatioLimit(ptg, par, xmin); |
| 1157 | Double_t max = CalculateJetRatioLimit(ptg, par, xmax); |
| 1158 | |
| 1159 | AliDebug(3,Form("Jet selection? : Limits min %f, max %f, pt_jet %f, pt_gamma %f, pt_jet / pt_gamma %f",min,max,ptj,ptg,ptj/ptg)); |
| 1160 | |
| 1161 | if(( min < ptj/ptg ) && ( max > ptj/ptg)) |
| 1162 | return kTRUE; |
| 1163 | else |
| 1164 | return kFALSE; |
| 1165 | |
| f9cea31c |
1166 | } |
| 1167 | |
| 1168 | //____________________________________________________________________________ |
| 1169 | void AliAnaGammaJet::MakeJet(TClonesArray * pl, TParticle * pGamma, TParticle* pLeading,TString lastname) |
| 1170 | { |
| 1171 | //Fill the jet with the particles around the leading particle with |
| 1172 | //R=fCone and pt_th = fPtThres. Caclulate the energy of the jet and |
| 1173 | //check if we select it. Fill jet histograms |
| 1174 | |
| 1175 | TClonesArray * jetList = new TClonesArray("TParticle",1000); |
| 1176 | TClonesArray * bkgList = new TClonesArray("TParticle",1000); |
| 1177 | |
| 1178 | TLorentzVector jet (0,0,0,0); |
| 1179 | TLorentzVector bkg(0,0,0,0); |
| 1180 | TLorentzVector lv (0,0,0,0); |
| 1181 | |
| 1182 | Double_t ptjet = 0.0; |
| 1183 | Double_t ptbkg = 0.0; |
| 1184 | Int_t n0 = 0; |
| 1185 | Int_t n1 = 0; |
| 1186 | Bool_t b1 = kFALSE; |
| 1187 | Bool_t b0 = kFALSE; |
| 1188 | |
| 1189 | Double_t ptg = pGamma->Pt(); |
| 1190 | Double_t phig = pGamma->Phi(); |
| 1191 | Double_t ptl = pLeading->Pt(); |
| 1192 | Double_t phil = pLeading->Phi(); |
| 1193 | Double_t etal = pLeading->Eta(); |
| 1194 | |
| 1195 | Float_t ptcut = 0. ; |
| 1196 | if(fPbPb){ |
| 1197 | if(ptg > fPtJetSelectionCut) ptcut = 2. ; |
| 1198 | else ptcut = 0.5; |
| 1199 | } |
| 1200 | |
| 1201 | TIter next(pl) ; |
| 1202 | TParticle * particle = 0 ; |
| 1203 | while ( (particle = dynamic_cast<TParticle*>(next())) ) { |
| 1204 | |
| 1205 | b0 = kFALSE; |
| 1206 | b1 = kFALSE; |
| 1207 | |
| 1208 | //Particles in jet |
| 1209 | SetJet(particle, b0, fCone, etal, phil) ; |
| 1210 | |
| 1211 | if(b0){ |
| 1212 | new((*jetList)[n0++]) TParticle(*particle) ; |
| 1213 | particle->Momentum(lv); |
| 1214 | if(particle->Pt() > ptcut ){ |
| 1215 | jet+=lv; |
| 1216 | ptjet+=particle->Pt(); |
| 1217 | } |
| 1218 | } |
| 1219 | |
| 1220 | //Background around (phi_gamma-pi, eta_leading) |
| 1221 | SetJet(particle, b1, fCone,etal, phig) ; |
| 1222 | |
| 1223 | if(b1) { |
| 1224 | new((*bkgList)[n1++]) TParticle(*particle) ; |
| 1225 | particle->Momentum(lv); |
| 1226 | if(particle->Pt() > ptcut ){ |
| 1227 | bkg+=lv; |
| 1228 | ptbkg+=particle->Pt(); |
| 1229 | } |
| 1230 | } |
| 1231 | } |
| 1232 | |
| 1233 | ptjet = jet.Pt(); |
| 1234 | ptbkg = bkg.Pt(); |
| 1235 | |
| 1236 | if(ptjet > 0.) { |
| 1237 | |
| 1238 | AliDebug(2,Form("Gamma pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg)); |
| 1239 | |
| 1240 | //Fill histograms |
| 1241 | |
| 1242 | Double_t ratjet = ptjet/ptg ; |
| 1243 | Double_t ratbkg = ptbkg/ptg ; |
| 1244 | |
| 1245 | dynamic_cast<TH2F*> |
| 1246 | (fOutputContainer->FindObject("JetRatio"+lastname)) |
| 1247 | ->Fill(ptg,ratjet); |
| 1248 | dynamic_cast<TH2F*> |
| 1249 | (fOutputContainer->FindObject("JetPt"+lastname)) |
| 1250 | ->Fill(ptg,ptjet); |
| 1251 | |
| 1252 | dynamic_cast<TH2F*> |
| 1253 | (fOutputContainer->FindObject("BkgRatio"+lastname)) |
| 1254 | ->Fill(ptg,ratbkg); |
| 1255 | |
| 1256 | dynamic_cast<TH2F*> |
| 1257 | (fOutputContainer->FindObject("BkgPt"+lastname)) |
| 1258 | ->Fill(ptg,ptbkg); |
| 1259 | |
| 1260 | |
| 1261 | //Jet selection |
| 1262 | Bool_t kSelect = kFALSE; |
| 1263 | if(fSelect == 0) |
| 1264 | kSelect = kTRUE; //Accept all jets, no restriction |
| 1265 | else if(fSelect == 1){ |
| 1266 | //Selection with parametrized cuts |
| 1267 | if(IsJetSelected(ptg,ptjet)) kSelect = kTRUE; |
| 1268 | } |
| 1269 | else if(fSelect == 2){ |
| 1270 | //Simple selection |
| 1271 | if(!fJetsOnlyInCTS){ |
| 1272 | if((ratjet < fJetRatioMaxCut) && (ratjet > fJetRatioMinCut )) kSelect = kTRUE; |
| 1273 | } |
| 1274 | else{ |
| 1275 | if((ratjet < fJetCTSRatioMaxCut) && (ratjet > fJetCTSRatioMinCut )) kSelect = kTRUE; |
| 1276 | } |
| 1277 | } |
| 1278 | else |
| 1279 | AliError("Jet selection option larger than 2, DONT SELECT JETS"); |
| 1280 | |
| 1281 | |
| 1282 | if(kSelect){ |
| 1283 | if (GetPrintInfo()) |
| 1284 | AliInfo(Form("Jet Selected: pt %f ", ptjet)) ; |
| 1285 | |
| 1286 | FillJetHistos(jetList, ptg, ptl,"Jet",lastname); |
| 1287 | FillJetHistos(bkgList, ptg, ptl, "Bkg",lastname); |
| 1288 | } |
| 1289 | } //ptjet > 0 |
| 1290 | |
| 1291 | jetList ->Delete(); |
| 1292 | bkgList ->Delete(); |
| 1293 | |
| 1294 | } |
| 1295 | |
| 1296 | //____________________________________________________________________________ |
| 1297 | void AliAnaGammaJet::Print(const Option_t * opt) const |
| 1298 | { |
| 1299 | |
| 1300 | //Print some relevant parameters set for the analysis |
| 1301 | if(! opt) |
| 1302 | return; |
| 1303 | |
| 1304 | Info("Print", "%s %s", GetName(), GetTitle() ) ; |
| 1305 | if(!fJetsOnlyInCTS) |
| 1306 | printf("EMCAL Acceptance cut : phi min %f, phi max %f, eta %f \n", fPhiEMCALCut[0], fPhiEMCALCut[1], fEtaEMCALCut) ; |
| 1307 | |
| 1308 | printf("Phi_Leading < %f\n", fPhiMaxCut) ; |
| 1309 | printf("Phi_Leading > %f\n", fPhiMinCut) ; |
| 1310 | printf("pT Leading / pT Gamma < %f\n", fRatioMaxCut) ; |
| 1311 | printf("pT Leading / pT Gamma > %f\n", fRatioMinCut) ; |
| 1312 | printf("M_pair < %f\n", fInvMassMaxCut) ; |
| 1313 | printf("M_pair > %f\n", fInvMassMinCut) ; |
| 1314 | if(fSelect == 2){ |
| 1315 | printf("pT Jet / pT Gamma < %f\n", fJetRatioMaxCut) ; |
| 1316 | printf("pT Jet / pT Gamma > %f\n", fJetRatioMinCut) ; |
| 1317 | printf("pT Jet (Only CTS)/ pT Gamma < %f\n", fJetCTSRatioMaxCut) ; |
| 1318 | printf("pT Jet (Only CTS)/ pT Gamma > %f\n", fJetCTSRatioMinCut) ; |
| 1319 | } |
| 1320 | |
| 1321 | |
| 1322 | } |
| 1323 | |
| 1324 | //___________________________________________________________________ |
| 1325 | void AliAnaGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone, |
| 1326 | Double_t eta, Double_t phi) |
| 1327 | { |
| 1328 | |
| 1329 | //Check if the particle is inside the cone defined by the leading particle |
| 1330 | b = kFALSE; |
| 1331 | |
| 1332 | if(phi > TMath::TwoPi()) |
| 1333 | phi-=TMath::TwoPi(); |
| 1334 | if(phi < 0.) |
| 1335 | phi+=TMath::TwoPi(); |
| 1336 | |
| 1337 | Double_t rad = 10000 + cone; |
| 1338 | |
| 1339 | if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone)) |
| 1340 | rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+ |
| 1341 | TMath::Power(part->Phi()-phi,2)); |
| 1342 | else{ |
| 1343 | if(part->Phi()-phi > TMath::TwoPi() - cone) |
| 1344 | rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+ |
| 1345 | TMath::Power((part->Phi()-TMath::TwoPi())-phi,2)); |
| 1346 | if(part->Phi()-phi < -(TMath::TwoPi() - cone)) |
| 1347 | rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+ |
| 1348 | TMath::Power((part->Phi()+TMath::TwoPi())-phi,2)); |
| 1349 | } |
| 1350 | |
| 1351 | if(rad < cone ) |
| 1352 | b = kTRUE; |
| 1353 | |
| 1354 | } |
| 1355 | |
| 1356 | |
| 1357 | void AliAnaGammaJet::Terminate(Option_t *) |
| 1358 | { |
| 1359 | // The Terminate() function is the last function to be called during |
| 1360 | // a query. It always runs on the client, it can be used to present |
| 1361 | // the results graphically or save the results to file. |
| 1362 | |
| 1363 | |
| 1364 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame