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