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.3  2007/03/08 10:24:32  schutz
  21  * Coding convention
  22  *
  23  * Revision 1.2  2007/02/09 18:40:40  schutz
  24  * New version from Gustavo
  25  *
  26  * Revision 1.1  2007/01/23 17:17:29  schutz
  27  * New Gamma package
  28  *
  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
  61 ClassImp(AliAnaGammaHadron)
  62
  63 //____________________________________________________________________________
  64 AliAnaGammaHadron::AliAnaGammaHadron(const char *name) :
  65   AliAnaGammaDirect(name),
  66   fPhiMaxCut(0.), fPhiMinCut(0.),
  67   fInvMassMaxCut(0.), fInvMassMinCut(0.),
  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
  78 {
  79
  80   // ctor
  81   fAngleMaxParam.Set(4) ;
  82   fAngleMaxParam.Reset(0.);
  83
  84   //Init Parameters
  85   InitParameters();
  86
  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 //____________________________________________________________________________
  96 AliAnaGammaHadron::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)
 110
 111 {
 112   // cpy ctor
 113   SetName (gh.GetName()) ;
 114   SetTitle(gh.GetTitle()) ;
 115
 116 }
 117
 118 //_________________________________________________________________________
 119 AliAnaGammaHadron & 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
 140 //____________________________________________________________________________
 141 AliAnaGammaHadron::~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
 170 //______________________________________________________________________________
 171 void AliAnaGammaHadron::ConnectInputData(const Option_t*)
 172 {
 173   // Initialisation of branch container and histograms
 174   AliAnaGammaDirect::ConnectInputData("");
 175 }
 176
 177 //________________________________________________________________________
 178 void AliAnaGammaHadron::CreateOutputObjects()
 179 {
 180
 181   // Create histograms to be saved in output file and
 182   // stores them in fOutputContainer
 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) ;
 219
 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 }
 337
 338 //____________________________________________________________________________
 339 void AliAnaGammaHadron::Exec(Option_t *)
 340 {
 341
 342   // Processing of one event
 343
 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 //____________________________________________________________________________
 424 void  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 //____________________________________________________________________________
 457 void  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   //____________________________________________________________________________
 571 void AliAnaGammaHadron::InitParameters()
 572 {
 573   // Initialisation of branch container
 574   //AliAnaGammaDirect::InitParameters();
 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
 596 }
 597
 598 //__________________________________________________________________________-
 599 Bool_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
 618 //____________________________________________________________________________
 619 void 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 //__________________________________________
 636 void 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 }