]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/AliAnaGammaHadron.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Possibility to compile with AliAnalysisGoodies (Yves)
[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaHadron.cxx
CommitLineData
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$
c005641f 20 * Revision 1.3 2007/03/08 10:24:32 schutz
21 * Coding convention
22 *
463ee300 23 * Revision 1.2 2007/02/09 18:40:40 schutz
24 * New version from Gustavo
25 *
2a1d8a29 26 * Revision 1.1 2007/01/23 17:17:29 schutz
27 * New Gamma package
28 *
f9cea31c 29 *
30 */
31
32//_________________________________________________________________________
33// Class for the analysis of gamma-jet correlations
34// Basically it seaches for a prompt photon in the Calorimeters acceptance,
35// if so we construct a jet around the highest pt particle in the opposite
36// side in azimuth, inside the Central Tracking System (ITS+TPC) and
37// EMCAL acceptances (only when PHOS detects the gamma). First the leading
38// particle and then the jet have to fullfill several conditions
39// (energy, direction ..) to be accepted. Then the fragmentation function
40// of this jet is constructed
41// Class created from old AliPHOSGammaPion
42// (see AliRoot versions previous Release 4-09)
43//
44//*-- Author: Gustavo Conesa (LNF-INFN)
45//////////////////////////////////////////////////////////////////////////////
46
47
48// --- ROOT system ---
49
50#include <TFile.h>
51#include <TParticle.h>
52#include <TH2.h>
53
54#include "AliAnaGammaHadron.h"
55#include "AliESD.h"
56#include "AliESDtrack.h"
57#include "AliESDCaloCluster.h"
58#include "Riostream.h"
59#include "AliLog.h"
60
61ClassImp(AliAnaGammaHadron)
62
63//____________________________________________________________________________
64AliAnaGammaHadron::AliAnaGammaHadron(const char *name) :
65 AliAnaGammaDirect(name),
66 fPhiMaxCut(0.), fPhiMinCut(0.),
67 fInvMassMaxCut(0.), fInvMassMinCut(0.),
463ee300 68 fMinPtPion(0), fOutputContainer(new TObjArray(100)),
69 fAngleMaxParam(),
70 fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0),
71 fhDeltaPhiGammaCharged(0), fhDeltaPhiGammaNeutral(0),
72 fhDeltaEtaGammaCharged(0), fhDeltaEtaGammaNeutral(0),
73 fhCorrelationGammaNeutral(0), fhCorrelationGammaCharged(0),
74 fhAnglePairAccepted(0), fhAnglePairNoCut(0), fhAnglePairAzimuthCut(0),
75 fhAnglePairOpeningAngleCut(0), fhAnglePairAllCut(0),
76 fhInvMassPairNoCut(0), fhInvMassPairAzimuthCut(0), fhInvMassPairOpeningAngleCut(0), fhInvMassPairAllCut(0)
77
f9cea31c 78{
79
80 // ctor
81 fAngleMaxParam.Set(4) ;
82 fAngleMaxParam.Reset(0.);
463ee300 83
84 //Init Parameters
85 InitParameters();
86
f9cea31c 87 // Input slot #0 works with an Ntuple
88 DefineInput(0, TChain::Class());
89 // Output slot #0 writes into a TH1 container
90 DefineOutput(0, TObjArray::Class()) ;
91
92}
93
94
95//____________________________________________________________________________
463ee300 96AliAnaGammaHadron::AliAnaGammaHadron(const AliAnaGammaHadron & gh) :
97 AliAnaGammaDirect(gh),
98 fPhiMaxCut(gh.fPhiMaxCut), fPhiMinCut(gh.fPhiMinCut),
99 fInvMassMaxCut(gh.fInvMassMaxCut), fInvMassMinCut(gh.fInvMassMinCut),
100 fMinPtPion(gh.fMinPtPion),
101 fOutputContainer(gh.fOutputContainer), fAngleMaxParam(gh.fAngleMaxParam),
102 fhPhiCharged(gh.fhPhiCharged), fhPhiNeutral(gh.fhPhiNeutral), fhEtaCharged(gh.fhEtaCharged), fhEtaNeutral(gh.fhEtaNeutral),
103 fhDeltaPhiGammaCharged(gh.fhDeltaPhiGammaCharged), fhDeltaPhiGammaNeutral(gh.fhDeltaPhiGammaNeutral),
104 fhDeltaEtaGammaCharged(gh.fhDeltaEtaGammaCharged), fhDeltaEtaGammaNeutral(gh.fhDeltaEtaGammaNeutral),
105 fhCorrelationGammaNeutral(gh.fhCorrelationGammaNeutral), fhCorrelationGammaCharged(gh.fhCorrelationGammaCharged),
106 fhAnglePairAccepted(gh.fhAnglePairAccepted), fhAnglePairNoCut(gh. fhAnglePairNoCut), fhAnglePairAzimuthCut(gh.fhAnglePairAzimuthCut),
107 fhAnglePairOpeningAngleCut(gh. fhAnglePairOpeningAngleCut), fhAnglePairAllCut(gh. fhAnglePairAllCut),
108 fhInvMassPairNoCut(gh.fhInvMassPairNoCut), fhInvMassPairAzimuthCut(gh.fhInvMassPairAzimuthCut),
109 fhInvMassPairOpeningAngleCut(gh.fhInvMassPairOpeningAngleCut), fhInvMassPairAllCut(gh.fhInvMassPairAllCut)
f9cea31c 110
111{
112 // cpy ctor
463ee300 113 SetName (gh.GetName()) ;
114 SetTitle(gh.GetTitle()) ;
f9cea31c 115
116}
117
463ee300 118//_________________________________________________________________________
119AliAnaGammaHadron & AliAnaGammaHadron::operator = (const AliAnaGammaHadron & source)
120{
121 //assignment operator
122 if(&source == this) return *this;
123
124 fPhiMaxCut = source.fPhiMaxCut ; fPhiMinCut = source.fPhiMinCut ;
125 fInvMassMaxCut = source.fInvMassMaxCut ; fInvMassMinCut = source.fInvMassMinCut ;
126 fMinPtPion = source.fMinPtPion ;
127 fOutputContainer = source.fOutputContainer ; fAngleMaxParam = source.fAngleMaxParam ;
128 fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ; fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ;
129 fhDeltaPhiGammaCharged = source.fhDeltaPhiGammaCharged ; fhDeltaPhiGammaNeutral = source.fhDeltaPhiGammaNeutral ;
130 fhDeltaEtaGammaCharged = source.fhDeltaEtaGammaCharged ; fhDeltaEtaGammaNeutral = source.fhDeltaEtaGammaNeutral ;
131 fhCorrelationGammaNeutral = source.fhCorrelationGammaNeutral ; fhCorrelationGammaCharged = source.fhCorrelationGammaCharged ;
132 fhAnglePairAccepted = source.fhAnglePairAccepted ; fhAnglePairNoCut = source. fhAnglePairNoCut ; fhAnglePairAzimuthCut = source.fhAnglePairAzimuthCut ;
133 fhAnglePairOpeningAngleCut = source. fhAnglePairOpeningAngleCut ; fhAnglePairAllCut = source. fhAnglePairAllCut ;
134 fhInvMassPairNoCut = source.fhInvMassPairNoCut ; fhInvMassPairAzimuthCut = source.fhInvMassPairAzimuthCut ;
135 fhInvMassPairOpeningAngleCut = source.fhInvMassPairOpeningAngleCut ; fhInvMassPairAllCut = source.fhInvMassPairAllCut ;
136
137 return *this;
138}
139
f9cea31c 140//____________________________________________________________________________
141AliAnaGammaHadron::~AliAnaGammaHadron()
142{
143 // Remove all pointers
144 fOutputContainer->Clear() ;
145 delete fOutputContainer ;
146
147 delete fhPhiCharged ;
148 delete fhPhiNeutral ;
149 delete fhEtaCharged ;
150 delete fhEtaNeutral ;
151 delete fhDeltaPhiGammaCharged ;
152 delete fhDeltaPhiGammaNeutral ;
153 delete fhDeltaEtaGammaCharged ;
154 delete fhDeltaEtaGammaNeutral ;
155
156 delete fhCorrelationGammaNeutral ;
157 delete fhCorrelationGammaCharged ;
158
159 delete fhAnglePairNoCut ;
160 delete fhAnglePairAzimuthCut ;
161 delete fhAnglePairOpeningAngleCut ;
162 delete fhAnglePairAllCut ;
163 delete fhInvMassPairNoCut ;
164 delete fhInvMassPairAzimuthCut ;
165 delete fhInvMassPairOpeningAngleCut ;
166 delete fhInvMassPairAllCut ;
167
168}
169
2a1d8a29 170//______________________________________________________________________________
171void AliAnaGammaHadron::ConnectInputData(const Option_t*)
172{
173 // Initialisation of branch container and histograms
174 AliAnaGammaDirect::ConnectInputData("");
175}
176
177//________________________________________________________________________
178void AliAnaGammaHadron::CreateOutputObjects()
179{
180
463ee300 181 // Create histograms to be saved in output file and
2a1d8a29 182 // stores them in fOutputContainer
2a1d8a29 183 AliAnaGammaDirect::CreateOutputObjects();
184
185 fOutputContainer = new TObjArray(100) ;
186
187 //Use histograms in AliAnaGammaDirect
188 TObjArray * outputContainer =GetOutputContainer();
189 for(Int_t i = 0; i < outputContainer->GetEntries(); i++ )
190 fOutputContainer->Add(outputContainer->At(i)) ;
191
192 fhPhiCharged = new TH2F
193 ("PhiCharged","#phi_{#pi^{#pm}} vs p_{T #gamma}",
194 120,0,120,120,0,7);
195 fhPhiCharged->SetYTitle("#phi_{#pi^{#pm}} (rad)");
196 fhPhiCharged->SetXTitle("p_{T #gamma} (GeV/c)");
197 fOutputContainer->Add(fhPhiCharged) ;
198
199 fhPhiNeutral = new TH2F
200 ("PhiNeutral","#phi_{#pi^{0}} vs p_{T #gamma}",
201 120,0,120,120,0,7);
202 fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
203 fhPhiNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
204 fOutputContainer->Add(fhPhiNeutral) ;
205
206 fhEtaCharged = new TH2F
207 ("EtaCharged","#eta_{#pi^{#pm}} vs p_{T #gamma}",
208 120,0,120,120,-1,1);
209 fhEtaCharged->SetYTitle("#eta_{#pi^{#pm}} (rad)");
210 fhEtaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
211 fOutputContainer->Add(fhEtaCharged) ;
212
213 fhEtaNeutral = new TH2F
214 ("EtaNeutral","#eta_{#pi^{0}} vs p_{T #gamma}",
215 120,0,120,120,-1,1);
216 fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
217 fhEtaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
218 fOutputContainer->Add(fhEtaNeutral) ;
f9cea31c 219
2a1d8a29 220 fhDeltaPhiGammaCharged = new TH2F
221 ("DeltaPhiGammaCharged","#phi_{#gamma} - #phi_{charged #pi} vs p_{T #gamma}",
222 200,0,120,200,0,6.4);
223 fhDeltaPhiGammaCharged->SetYTitle("#Delta #phi");
224 fhDeltaPhiGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
225 fOutputContainer->Add(fhDeltaPhiGammaCharged) ;
226
227 fhDeltaEtaGammaCharged = new TH2F
228 ("DeltaEtaGammaCharged","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
229 200,0,120,200,-2,2);
230 fhDeltaEtaGammaCharged->SetYTitle("#Delta #eta");
231 fhDeltaEtaGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
232 fOutputContainer->Add(fhDeltaEtaGammaCharged) ;
233
234 fhDeltaPhiGammaNeutral = new TH2F
235 ("DeltaPhiGammaNeutral","#phi_{#gamma} - #phi_{#pi^{0}} vs p_{T #gamma}",
236 200,0,120,200,0,6.4);
237 fhDeltaPhiGammaNeutral->SetYTitle("#Delta #phi");
238 fhDeltaPhiGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
239 fOutputContainer->Add(fhDeltaPhiGammaNeutral) ;
240
241 fhDeltaEtaGammaNeutral = new TH2F
242 ("DeltaEtaGammaNeutral","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
243 200,0,120,200,-2,2);
244 fhDeltaEtaGammaNeutral->SetYTitle("#Delta #eta");
245 fhDeltaEtaGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
246 fOutputContainer->Add(fhDeltaEtaGammaNeutral) ;
247
248 //
249 fhAnglePairAccepted = new TH2F
250 ("AnglePairAccepted",
251 "Angle between #pi^{0} #gamma pair vs p_{T #pi^{0}}, both #gamma in eta<0.7, inside window",
252 200,0,50,200,0,0.2);
253 fhAnglePairAccepted->SetYTitle("Angle (rad)");
254 fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
255 fOutputContainer->Add(fhAnglePairAccepted) ;
256
257 fhAnglePairNoCut = new TH2F
258 ("AnglePairNoCut",
259 "Angle between all #gamma pair vs p_{T #pi^{0}}",200,0,50,200,0,0.2);
260 fhAnglePairNoCut->SetYTitle("Angle (rad)");
261 fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
262 fOutputContainer->Add(fhAnglePairNoCut) ;
263
264 fhAnglePairAzimuthCut = new TH2F
265 ("AnglePairAzimuthCut",
266 "Angle between all #gamma pair that have a good phi and pt vs p_{T #pi^{0}}",
267 200,0,50,200,0,0.2);
268 fhAnglePairAzimuthCut->SetYTitle("Angle (rad)");
269 fhAnglePairAzimuthCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
270 fOutputContainer->Add(fhAnglePairAzimuthCut) ;
271
272 fhAnglePairOpeningAngleCut = new TH2F
273 ("AnglePairOpeningAngleCut",
274 "Angle between all #gamma pair (opening angle + azimuth cut) vs p_{T #pi^{0}}"
275 ,200,0,50,200,0,0.2);
276 fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
277 fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
278 fOutputContainer->Add(fhAnglePairOpeningAngleCut) ;
279
280 fhAnglePairAllCut = new TH2F
281 ("AnglePairAllCut",
282 "Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs p_{T #pi^{0}}"
283 ,200,0,50,200,0,0.2);
284 fhAnglePairAllCut->SetYTitle("Angle (rad)");
285 fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
286 fOutputContainer->Add(fhAnglePairAllCut) ;
287
288
289 //
290 fhInvMassPairNoCut = new TH2F
291 ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
292 120,0,120,360,0,0.5);
293 fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
294 fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
295 fOutputContainer->Add(fhInvMassPairNoCut) ;
296
297 fhInvMassPairAzimuthCut = new TH2F
298 ("InvMassPairAzimuthCut",
299 "Invariant Mass of #gamma pair (azimuth cuts) vs p_{T #gamma}",
300 120,0,120,360,0,0.5);
301 fhInvMassPairAzimuthCut->SetYTitle("Invariant Mass (GeV/c^{2})");
302 fhInvMassPairAzimuthCut->SetXTitle("p_{T #gamma} (GeV/c)");
303 fOutputContainer->Add(fhInvMassPairAzimuthCut) ;
304
305 fhInvMassPairOpeningAngleCut = new TH2F
306 ("InvMassPairOpeningAngleCut",
307 "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
308 120,0,120,360,0,0.5);
309 fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
310 fhInvMassPairOpeningAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
311 fOutputContainer->Add(fhInvMassPairOpeningAngleCut) ;
312
313 fhInvMassPairAllCut = new TH2F
314 ("InvMassPairAllCut",
315 "Invariant Mass of #gamma pair (opening angle+invmass cut+azimuth) vs p_{T #gamma}",
316 120,0,120,360,0,0.5);
317 fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
318 fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
319 fOutputContainer->Add(fhInvMassPairAllCut) ;
320
321 //
322 fhCorrelationGammaCharged =
323 new TH2F("CorrelationGammaCharged","z_{#gamma #pi} = p_{T #pi^{#pm}} / p_{T #gamma}",
324 240,0.,120.,1000,0.,1.2);
325 fhCorrelationGammaCharged->SetYTitle("z_{#gamma #pi}");
326 fhCorrelationGammaCharged->SetXTitle("p_{T #gamma}");
327 fOutputContainer->Add(fhCorrelationGammaCharged) ;
328
329 fhCorrelationGammaNeutral =
330 new TH2F("CorrelationGammaNeutral","z_{#gamma #pi} = p_{T #pi^{0}} / p_{T #gamma}",
331 240,0.,120.,1000,0.,1.2);
332 fhCorrelationGammaNeutral->SetYTitle("z_{#gamma #pi}");
333 fhCorrelationGammaNeutral->SetXTitle("p_{T #gamma}");
334 fOutputContainer->Add(fhCorrelationGammaNeutral) ;
335
336}
f9cea31c 337
338//____________________________________________________________________________
339void AliAnaGammaHadron::Exec(Option_t *)
340{
341
342 // Processing of one event
2a1d8a29 343
f9cea31c 344 //Get ESDs
345 Long64_t entry = GetChain()->GetReadEntry() ;
346
347 if (!GetESD()) {
348 AliError("fESD is not connected to the input!") ;
349 return ;
350 }
351
352 if (GetPrintInfo())
353 AliInfo(Form("%s ----> Processing event # %lld", (dynamic_cast<TChain *>(GetChain()))->GetFile()->GetName(), entry)) ;
354
355 //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
356
357 TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
358 TClonesArray * plCTS = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
359 TClonesArray * plNe = new TClonesArray("TParticle",1000); // All particles measured in Jet Calorimeter
360 TClonesArray * plCalo = new TClonesArray("TParticle",1000); // All particles measured in Prompt Gamma calorimeter
361
362
363 TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
364
365
366 Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
367
368 //Fill lists with photons, neutral particles and charged particles
369 //look for the highest energy photon in the event inside fCalorimeter
370 //Fill particle lists
371 AliDebug(2, "Fill particle lists, get prompt gamma");
372
373 //Fill particle lists
374 if(GetCalorimeter() == "PHOS")
375 CreateParticleList(particleList, plCTS,plNe,plCalo);
376 else if(GetCalorimeter() == "EMCAL")
377 CreateParticleList(particleList, plCTS,plCalo,plNe);
378 else
379 AliError("No calorimeter selected");
380
381 //Search highest energy prompt gamma in calorimeter
382 GetPromptGamma(plCalo, plCTS, pGamma, iIsInPHOSorEMCAL) ;
383
384
385 AliDebug(1, Form("Is Gamma in %s? %d",GetCalorimeter().Data(),iIsInPHOSorEMCAL));
386 AliDebug(3,Form("Charged list entries %d, Neutral list entries %d, %s list entries %d",
387 plCTS->GetEntries(),plNe->GetEntries(), GetCalorimeter().Data(),plCalo->GetEntries()));
388
389 //If there is any prompt photon in fCalorimeter,
390 //search jet leading particle
391
392 if(iIsInPHOSorEMCAL){
393
394 if (GetPrintInfo())
395 AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
396
397 AliDebug(2, "Make correlation");
398
399 //Search correlation
400 MakeGammaChargedCorrelation(plCTS, pGamma);
401 MakeGammaNeutralCorrelation(plNe, pGamma);
402
403 }//Gamma in Calo
404
405 AliDebug(2, "End of analysis, delete pointers");
406
407 particleList->Delete() ;
408 plCTS->Delete() ;
409 plNe->Delete() ;
410 plCalo->Delete() ;
411 pGamma->Delete();
412
413 delete plNe ;
414 delete plCalo ;
415 delete plCTS ;
416 delete particleList ;
417 // delete pGamma;
418
419 PostData(0, fOutputContainer);
420}
421
422
423//____________________________________________________________________________
424void AliAnaGammaHadron::MakeGammaChargedCorrelation(TClonesArray * pl, TParticle * pGamma) const
425{
426 //Search for the charged particle with highest with
427 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
428 Double_t ptg = pGamma->Pt();
429 Double_t phig = pGamma->Phi();
430 Double_t pt = -100.;
431 Double_t rat = -100.;
432 Double_t phi = -100. ;
433
434 for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
435
436 TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
437
438 pt = particle->Pt();
439 rat = pt/ptg ;
440 phi = particle->Phi() ;
441
442 AliDebug(3,Form("pt %f, phi %f, phi gamma %f. Cuts: delta phi min %f, max%f, pT min %f",pt,phi,phig,fPhiMinCut,fPhiMaxCut,fMinPtPion));
443
444 fhEtaCharged->Fill(ptg,particle->Eta());
445 fhPhiCharged->Fill(ptg,phi);
446 fhDeltaEtaGammaCharged->Fill(ptg,pGamma->Eta()-particle->Eta());
447 fhDeltaPhiGammaCharged->Fill(ptg,phig-phi);
448 //Selection within angular and energy limits
449 if(((phig-phi)> fPhiMinCut) && ((phig-phi)<fPhiMaxCut) && pt > fMinPtPion){
450 AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
451 fhCorrelationGammaCharged->Fill(ptg,rat);
452 }
453 }//particle loop
454}
455
456//____________________________________________________________________________
457void AliAnaGammaHadron::MakeGammaNeutralCorrelation(TClonesArray * pl, TParticle * pGamma)
458{
459
460 //Search for the neutral pion with highest with
461 //Phi=Phi_gamma-Pi and pT=0.1E_gamma
462 Double_t pt = -100.;
463 Double_t rat = -100.;
464 Double_t phi = -100. ;
465 Double_t ptg = pGamma->Pt();
466 Double_t phig = pGamma->Phi();
467
468 TIter next(pl);
469 TParticle * particle = 0;
470
471 Int_t iPrimary = -1;
472 TLorentzVector gammai,gammaj;
473 Double_t angle = 0., e = 0., invmass = 0.;
474 Int_t ksPdg = 0;
475 Int_t jPrimary=-1;
476
477 while ( (particle = (TParticle*)next()) ) {
478 iPrimary++;
479 ksPdg = particle->GetPdgCode();
480
481 //2 gamma overlapped, found with PID
482 if(ksPdg == 111){
483 pt = particle->Pt();
484 rat = pt/ptg ;
485 phi = particle->Phi() ;
486 fhEtaCharged->Fill(ptg,particle->Eta());
487 fhPhiCharged->Fill(ptg,phi);
488 fhDeltaEtaGammaCharged->Fill(ptg,pGamma->Eta()-particle->Eta());
489 fhDeltaPhiGammaCharged->Fill(ptg,phig-phi);
490 //AliDebug(3,Form("pt %f, phi %f",pt,phi));
491 if (GetPrintInfo())
492 AliInfo(Form("pt %f, phi %f",pt,phi));
493 //Selection within angular and energy limits
494 if((pt> ptg)&& ((phig-phi)>fPhiMinCut)&&((phig-phi)<fPhiMaxCut)){
495 fhCorrelationGammaNeutral ->Fill(ptg,rat);
496 //AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
497 if (GetPrintInfo())
498 AliInfo(Form("Selected: pt %f, phi %f",pt,phi));
499 }// cuts
500 }// pdg = 111
501
502 //Make invariant mass analysis
503 else if(ksPdg == 22){
504 //Search the photon companion in case it comes from a Pi0 decay
505 //Apply several cuts to select the good pair
506 particle->Momentum(gammai);
507 jPrimary=-1;
508 TIter next2(pl);
509 while ( (particle = (TParticle*)next2()) ) {
510 jPrimary++;
511 if(jPrimary>iPrimary){
512 ksPdg = particle->GetPdgCode();
513
514 if(ksPdg == 22){
515 particle->Momentum(gammaj);
516 //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
517 //gammai.Pt(),gammaj.Pt());
518 pt = (gammai+gammaj).Pt();
519 phi = (gammai+gammaj).Phi();
520 if(phi < 0)
521 phi+=TMath::TwoPi();
522 rat = pt/ptg ;
523 invmass = (gammai+gammaj).M();
524 angle = gammaj.Angle(gammai.Vect());
525 e = (gammai+gammaj).E();
526 fhEtaNeutral->Fill(ptg,(gammai+gammaj).Eta());
527 fhPhiNeutral->Fill(ptg,phi);
528 fhDeltaEtaGammaNeutral->Fill(ptg,pGamma->Eta()-(gammai+gammaj).Eta());
529 fhDeltaPhiGammaNeutral->Fill(ptg,phig-phi);
530 // AliDebug(3,Form("pt %f, phi %f",pt,phi));
531 if (GetPrintInfo())
532 AliInfo(Form("pt %f, phi %f",pt,phi));
533
534 //Fill histograms with no cuts applied.
535 fhAnglePairNoCut->Fill(e,angle);
536 fhInvMassPairNoCut->Fill(ptg,invmass);
537
538 //First cut on the energy and azimuth of the pair
539 if((phig-phi) > fPhiMinCut && (phig-phi) < fPhiMaxCut
540 && pt > fMinPtPion){
541 fhAnglePairAzimuthCut ->Fill(e,angle);
542 fhInvMassPairAzimuthCut->Fill(ptg,invmass);
543 AliDebug(3,Form("1st cut: pt %f, phi %f",pt,phi));
544
545 //Second cut on the aperture of the pair
546 if(IsAngleInWindow(angle,e)){
547 fhAnglePairOpeningAngleCut ->Fill(e,angle);
548 fhInvMassPairOpeningAngleCut->Fill(ptg,invmass);
549 AliDebug(3,Form("2nd cut: pt %f, phi %f",pt,phi));
550
551 //Third cut on the invariant mass of the pair
552 if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
553 fhInvMassPairAllCut ->Fill(ptg,invmass);
554 fhAnglePairAllCut ->Fill(e,angle);
555 //Fill correlation histogram
556 fhCorrelationGammaNeutral ->Fill(ptg,rat);
557 //AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
558 if (GetPrintInfo())
559 AliInfo(Form("Selected: pt %f, phi %f",pt,phi));
560 }//(invmass>0.125) && (invmass<0.145)
561 }//Opening angle cut
562 }//Azimuth and pt cut.
563 }//if pdg = 22
564 }//iprimary<jprimary
565 }//while
566 }// if pdg = 22
567 }//while
568}
569
570 //____________________________________________________________________________
2a1d8a29 571void AliAnaGammaHadron::InitParameters()
f9cea31c 572{
573 // Initialisation of branch container
2a1d8a29 574 //AliAnaGammaDirect::InitParameters();
f9cea31c 575
576 //Initialize the parameters of the analysis.
577 //fCalorimeter="PHOS";
578 fAngleMaxParam.Set(4) ;
579 fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
580 fAngleMaxParam.AddAt(-0.25,1) ;
581 fAngleMaxParam.AddAt(0.025,2) ;
582 fAngleMaxParam.AddAt(-2e-4,3) ;
583
584 //fPrintInfo = kTRUE;
585 fInvMassMaxCut = 0.16 ;
586 fInvMassMinCut = 0.11 ;
587 fPhiMaxCut = 4.5;
588 fPhiMinCut = 1.5 ;
589
590 fMinPtPion = 0. ;
591
592 //Fill particle lists when PID is ok
593 // fEMCALPID = kFALSE;
594 // fPHOSPID = kFALSE;
595
f9cea31c 596}
597
598//__________________________________________________________________________-
599Bool_t AliAnaGammaHadron::IsAngleInWindow(const Float_t angle,const Float_t e) {
600 //Check if the opening angle of the candidate pairs is inside
601 //our selection windowd
602
603 Bool_t result = kFALSE;
604 Double_t mpi0 = 0.1349766;
605 Double_t max = fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
606 +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
607 Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
608 Double_t min = 100. ;
609 if(arg>0.)
610 min = TMath::ACos(arg);
611
612 if((angle<max)&&(angle>=min))
613 result = kTRUE;
614
615 return result;
616}
617
f9cea31c 618//____________________________________________________________________________
619void AliAnaGammaHadron::Print(const Option_t * opt) const
620{
621
622 //Print some relevant parameters set for the analysis
623 if(! opt)
624 return;
625
626 Info("Print", "%s %s", GetName(), GetTitle() ) ;
627 printf("pT Pion > %f\n", fMinPtPion) ;
628 printf("Phi Pion < %f\n", fPhiMaxCut) ;
629 printf("Phi Pion > %f\n", fPhiMinCut) ;
630 printf("M_pair < %f\n", fInvMassMaxCut) ;
631 printf("M_pair > %f\n", fInvMassMinCut) ;
632
633}
634
635//__________________________________________
636void AliAnaGammaHadron::Terminate(Option_t *)
637{
638 // The Terminate() function is the last function to be called during
639 // a query. It always runs on the client, it can be used to present
640 // the results graphically or save the results to file.
641
642
643}