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