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.5  2007/08/17 12:40:04  schutz
  21  * New analysis classes by Gustavo Conesa
  22  *
  23  * Revision 1.4.4.2  2007/07/26 10:32:09  schutz
  24  * new analysis classes in the the new analysis framework
  25  *
  26  *
  27  */
  28
  29 //_________________________________________________________________________
  30 // Class for the analysis of gamma - hadron correlations
  31 //*-- Author: Gustavo Conesa (LNF-INFN)
  32 //////////////////////////////////////////////////////////////////////////////
  33
  34
  35 // --- ROOT system ---
  36 #include "Riostream.h"
  37
  38 //---- AliRoot system ----
  39 #include "AliLog.h"
  40 #include "AliNeutralMesonSelection.h"
  41 #include "AliAnaGammaHadron.h"
  42
  43 ClassImp(AliAnaGammaHadron)
  44
  45
  46 //____________________________________________________________________________
  47   AliAnaGammaHadron::AliAnaGammaHadron() :
  48     AliAnaGammaCorrelation(),
  49     fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0),
  50     fhDeltaPhiGammaCharged(0),  fhDeltaPhiGammaNeutral(0),
  51     fhDeltaEtaGammaCharged(0), fhDeltaEtaGammaNeutral(0),
  52     fhDeltaPhiChargedPt(0),
  53     fhCorrelationGammaNeutral(0), fhCorrelationGammaCharged(0)
  54 {
  55   //Default Ctor
  56
  57   SetCorrelationType(kHadron);
  58   //Initialize parameters
  59   InitParameters();
  60 }
  61
  62 //____________________________________________________________________________
  63 AliAnaGammaHadron::AliAnaGammaHadron(const AliAnaGammaHadron & g) :
  64   AliAnaGammaCorrelation(g),
  65   fhPhiCharged(g.fhPhiCharged), fhPhiNeutral(g.fhPhiNeutral),
  66   fhEtaCharged(g.fhEtaCharged), fhEtaNeutral(g.fhEtaNeutral),
  67   fhDeltaPhiGammaCharged(g.fhDeltaPhiGammaCharged),
  68   fhDeltaPhiGammaNeutral(g.fhDeltaPhiGammaNeutral),
  69   fhDeltaEtaGammaCharged(g.fhDeltaEtaGammaCharged),
  70   fhDeltaEtaGammaNeutral(g.fhDeltaEtaGammaNeutral),
  71   fhDeltaPhiChargedPt(g.fhDeltaPhiChargedPt),
  72   fhCorrelationGammaNeutral(g.fhCorrelationGammaNeutral),
  73   fhCorrelationGammaCharged(g.fhCorrelationGammaCharged)
  74 {
  75   // cpy ctor
  76
  77 }
  78
  79 //_________________________________________________________________________
  80 AliAnaGammaHadron & AliAnaGammaHadron::operator = (const AliAnaGammaHadron & source)
  81 {
  82   // assignment operator
  83
  84   if(this == &source)return *this;
  85   ((AliAnaGammaCorrelation *)this)->operator=(source);
  86
  87   fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ;
  88   fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ;
  89   fhDeltaPhiGammaCharged = source.fhDeltaPhiGammaCharged ;
  90   fhDeltaPhiGammaNeutral = source.fhDeltaPhiGammaNeutral ;
  91   fhDeltaEtaGammaCharged = source.fhDeltaEtaGammaCharged ;
  92   fhDeltaEtaGammaNeutral = source.fhDeltaEtaGammaNeutral ;
  93   fhDeltaPhiChargedPt = source.fhDeltaPhiChargedPt ;
  94
  95   fhCorrelationGammaNeutral = source.fhCorrelationGammaNeutral ;
  96   fhCorrelationGammaCharged = source.fhCorrelationGammaCharged ;
  97
  98   return *this;
  99
 100 }
 101
 102 //____________________________________________________________________________
 103 AliAnaGammaHadron::~AliAnaGammaHadron()
 104 {
 105
 106   delete fhPhiCharged  ;
 107   delete fhPhiNeutral   ;
 108   delete fhEtaCharged  ;
 109   delete fhEtaNeutral  ;
 110   delete fhDeltaPhiGammaCharged  ;
 111   delete fhDeltaPhiGammaNeutral   ;
 112   delete fhDeltaEtaGammaCharged  ;
 113   delete fhDeltaEtaGammaNeutral  ;
 114   delete fhDeltaPhiChargedPt  ;
 115
 116   delete fhCorrelationGammaNeutral  ;
 117   delete fhCorrelationGammaCharged  ;
 118
 119 }
 120
 121
 122 //________________________________________________________________________
 123 TList *  AliAnaGammaHadron::GetCreateOutputObjects()
 124 {
 125
 126   // Create histograms to be saved in output file and
 127   // store them in fOutputContainer
 128   TList * outputContainer = new TList() ;
 129   outputContainer->SetName("GammaCorrelationHistos") ;
 130
 131   fhPhiCharged  = new TH2F
 132     ("PhiCharged","#phi_{#pi^{#pm}}  vs p_{T #gamma}",
 133      120,0,120,120,0,7);
 134   fhPhiCharged->SetYTitle("#phi_{#pi^{#pm}} (rad)");
 135   fhPhiCharged->SetXTitle("p_{T #gamma} (GeV/c)");
 136
 137   fhEtaCharged  = new TH2F
 138     ("EtaCharged","#eta_{#pi^{#pm}}  vs p_{T #gamma}",
 139      120,0,120,120,-1,1);
 140   fhEtaCharged->SetYTitle("#eta_{#pi^{#pm}} (rad)");
 141   fhEtaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
 142
 143   fhDeltaPhiGammaCharged  = new TH2F
 144     ("DeltaPhiGammaCharged","#phi_{#gamma} - #phi_{charged #pi} vs p_{T #gamma}",
 145      200,0,120,200,0,6.4);
 146   fhDeltaPhiGammaCharged->SetYTitle("#Delta #phi");
 147   fhDeltaPhiGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
 148
 149   fhDeltaPhiChargedPt  = new TH2F
 150     ("DeltaPhiChargedPt","#phi_{#gamma} - #phi_{charged #pi} vs p_{T #pi}",
 151      200,0,120,200,0,6.4);
 152   fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
 153   fhDeltaPhiChargedPt->SetXTitle("p_{T #pi} (GeV/c)");
 154
 155   fhDeltaEtaGammaCharged  = new TH2F
 156     ("DeltaEtaGammaCharged","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
 157      200,0,120,200,-2,2);
 158   fhDeltaEtaGammaCharged->SetYTitle("#Delta #eta");
 159   fhDeltaEtaGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
 160
 161   fhCorrelationGammaCharged  =
 162     new TH2F("CorrelationGammaCharged","z_{#gamma #pi} = p_{T #pi^{#pm}} / p_{T #gamma}",
 163              240,0.,120.,1000,0.,1.2);
 164   fhCorrelationGammaCharged->SetYTitle("z_{#gamma #pi}");
 165   fhCorrelationGammaCharged->SetXTitle("p_{T #gamma}");
 166
 167   outputContainer->Add(fhPhiCharged) ;
 168   outputContainer->Add(fhEtaCharged) ;
 169   outputContainer->Add(fhDeltaPhiGammaCharged) ;
 170   outputContainer->Add(fhDeltaEtaGammaCharged) ;
 171   outputContainer->Add(fhCorrelationGammaCharged) ;
 172   outputContainer->Add(fhDeltaPhiChargedPt) ;
 173
 174   if(!AreJetsOnlyInCTS()){
 175     //---- kHadron and kJetLeadCone ----
 176     fhPhiNeutral  = new TH2F
 177       ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T #gamma}",
 178        120,0,120,120,0,7);
 179     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
 180     fhPhiNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
 181
 182     fhEtaNeutral  = new TH2F
 183       ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T #gamma}",
 184        120,0,120,120,-1,1);
 185     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
 186     fhEtaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
 187
 188     fhDeltaPhiGammaNeutral  = new TH2F
 189       ("DeltaPhiGammaNeutral","#phi_{#gamma} - #phi_{#pi^{0}} vs p_{T #gamma}",
 190        200,0,120,200,0,6.4);
 191     fhDeltaPhiGammaNeutral->SetYTitle("#Delta #phi");
 192     fhDeltaPhiGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
 193
 194     fhDeltaEtaGammaNeutral  = new TH2F
 195       ("DeltaEtaGammaNeutral","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
 196        200,0,120,200,-2,2);
 197     fhDeltaEtaGammaNeutral->SetYTitle("#Delta #eta");
 198     fhDeltaEtaGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
 199
 200     fhCorrelationGammaNeutral  =
 201       new TH2F("CorrelationGammaNeutral","z_{#gamma #pi} = p_{T #pi^{0}} / p_{T #gamma}",
 202                240,0.,120.,1000,0.,1.2);
 203     fhCorrelationGammaNeutral->SetYTitle("z_{#gamma #pi}");
 204     fhCorrelationGammaNeutral->SetXTitle("p_{T #gamma}");
 205
 206     outputContainer->Add(fhPhiNeutral) ;
 207     outputContainer->Add(fhEtaNeutral) ;
 208     outputContainer->Add(fhDeltaPhiGammaNeutral) ;
 209     outputContainer->Add(fhDeltaEtaGammaNeutral) ;
 210     outputContainer->Add(fhCorrelationGammaNeutral) ;
 211   }
 212
 213   SetOutputContainer(outputContainer);
 214
 215   return outputContainer;
 216 }
 217
 218  //____________________________________________________________________________
 219 void AliAnaGammaHadron::InitParameters()
 220 {
 221
 222   //Initialize the parameters of the analysis.
 223
 224   SetMinPtHadron(2.)   ;
 225   SetDeltaPhiCutRange(1.5,4.5);
 226   SetJetsOnlyInCTS(kFALSE) ;
 227
 228 }
 229
 230 //__________________________________________________________________
 231 void AliAnaGammaHadron::Print(const Option_t * opt) const
 232 {
 233
 234   //Print some relevant parameters set for the analysis
 235   if(! opt)
 236     return;
 237
 238   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 239   printf("Correlation analysis           =     %d\n", kHadron) ;
 240   printf("pT Hadron       >    %f\n", GetMinPtHadron()) ;
 241   printf("Phi gamma-Hadron      <     %f\n", GetDeltaPhiMaxCut()) ;
 242   printf("Phi gamma-Hadron      >     %f\n", GetDeltaPhiMinCut()) ;
 243
 244
 245 }
 246
 247 //____________________________________________________________________________
 248 void  AliAnaGammaHadron::MakeGammaCorrelation(TParticle * pGamma, TClonesArray * plCTS, TClonesArray * plCalo)
 249 {
 250   //Gamma Hadron Correlation Analysis
 251   AliDebug(2, "Make gamma-hadron correlation");
 252
 253   MakeGammaChargedCorrelation(pGamma, plCTS);
 254   if(!AreJetsOnlyInCTS())
 255   MakeGammaNeutralCorrelation(pGamma, plCalo);
 256
 257
 258 }
 259
 260 //____________________________________________________________________________
 261 void  AliAnaGammaHadron::MakeGammaChargedCorrelation(TParticle * pGamma, TClonesArray * pl)
 262 {
 263   //Gamma Charged Hadron Correlation Analysis
 264   AliDebug(2,"Make gamma-charged hadron correlation");
 265
 266   Double_t ptg  = pGamma->Pt();
 267   Double_t phig = pGamma->Phi();
 268   Double_t pt    = -100.;
 269   Double_t rat   = -100.;
 270   Double_t phi   = -100. ;
 271
 272   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
 273
 274     TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
 275
 276     pt    = particle->Pt();
 277     rat   = pt/ptg ;
 278     phi   = particle->Phi() ;
 279
 280     AliDebug(3,Form("pt %f, phi %f, phi gamma %f. Cuts:  delta phi min %f,  max%f, pT min %f",pt,phi,phig,GetDeltaPhiMinCut(),GetDeltaPhiMaxCut(),GetMinPtHadron()));
 281
 282     fhEtaCharged->Fill(ptg,particle->Eta());
 283     fhPhiCharged->Fill(ptg,phi);
 284     fhDeltaEtaGammaCharged->Fill(ptg,pGamma->Eta()-particle->Eta());
 285     fhDeltaPhiGammaCharged->Fill(ptg,phig-phi);
 286     fhDeltaPhiChargedPt->Fill(pt,phig-phi);
 287
 288     //Selection within angular and energy limits
 289     if(((phig-phi)> GetDeltaPhiMinCut()) && ((phig-phi)<GetDeltaPhiMaxCut()) && pt > GetMinPtHadron()){
 290       AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
 291       fhCorrelationGammaCharged->Fill(ptg,rat);
 292     }
 293   }//particle loop
 294 }
 295
 296 //____________________________________________________________________________
 297 void  AliAnaGammaHadron::MakeGammaNeutralCorrelation(TParticle * pGamma, TClonesArray * pl)
 298 {
 299   //Gamma Neutral Hadron Correlation Analysis
 300   AliDebug(2,"Make gamma-neutral hadron correlation");
 301
 302   Double_t pt = -100.;
 303   Double_t rat = -100.;
 304   Double_t phi = -100. ;
 305   Double_t ptg  = pGamma->Pt();
 306   Double_t phig = pGamma->Phi();
 307
 308   TIter next(pl);
 309   TParticle * particlei = 0;
 310   TParticle * particlej = 0;
 311   TLorentzVector gammai;
 312   TLorentzVector gammaj;
 313
 314   Int_t iPrimary = -1;
 315   Int_t ksPdg = 0;
 316   Int_t jPrimary=-1;
 317
 318   while ( (particlei = (TParticle*)next()) ) {
 319     iPrimary++;
 320     ksPdg = particlei->GetPdgCode();
 321     AliDebug(2, Form("neutral particles opposite to gamma: pt %f, pdg %d", particlei->Pt(),ksPdg));
 322     //2 gamma overlapped, found with PID
 323     if(ksPdg == 111){
 324       pt  = particlei->Pt();
 325       rat = pt/ptg ;
 326       phi = particlei->Phi() ;
 327       fhEtaNeutral->Fill(ptg,particlei->Eta());
 328       fhPhiNeutral->Fill(ptg,phi);
 329       fhDeltaEtaGammaNeutral->Fill(ptg,pGamma->Eta()-particlei->Eta());
 330       fhDeltaPhiGammaNeutral->Fill(ptg,phig-phi);
 331
 332       //Selection within angular and energy limits
 333       if( (phig-phi)>GetDeltaPhiMinCut() && (phig-phi)<GetDeltaPhiMaxCut() && pt > GetMinPtHadron()){
 334         fhCorrelationGammaNeutral ->Fill(ptg,rat);
 335         AliDebug(2,Form("Selected pi0: pt %f, phi %f",pt,phi));
 336       }// cuts
 337     }// pdg = 111
 338
 339     //Make invariant mass analysis
 340     else if(ksPdg == 22){//  gamma i
 341
 342       //Search the photon companion in case it comes from  a Pi0 decay
 343       //Apply several cuts to select the good pair;
 344       particlei->Momentum(gammai);
 345        jPrimary=-1;
 346       TIter next2(pl);
 347       while ( (particlej = (TParticle*)next2()) ) {
 348         jPrimary++;
 349         if(jPrimary>iPrimary){
 350           ksPdg = particlej->GetPdgCode();
 351           particlej->Momentum(gammaj);
 352           if(ksPdg == 22 ){
 353
 354             phi = (gammai+gammaj).Phi();
 355             if(phi < 0)
 356               phi+=TMath::TwoPi();
 357             rat          = (gammai+gammaj).Pt()/ptg ;
 358
 359             //Fill histograms
 360             fhEtaNeutral->Fill(ptg,(gammai+gammaj).Eta());
 361             fhPhiNeutral->Fill(ptg,phi);
 362             fhDeltaEtaGammaNeutral->Fill(ptg,pGamma->Eta()-(gammai+gammaj).Eta());
 363             fhDeltaPhiGammaNeutral->Fill(ptg,phig-phi);
 364
 365             //Select good pair (good phit, pt cuts, aperture and invariant mass)
 366             if(GetNeutralMesonSelection()->SelectPair(pGamma, gammai, gammaj)){
 367               AliDebug(2,Form("Selected gamma pair: pt %f, phi %f",(gammai+gammaj).Pt(),phi));
 368               //correlation histogram
 369               fhCorrelationGammaNeutral ->Fill(ptg,rat);
 370             }//pair selection
 371           }//if pair of gammas
 372         }//jPrimary>iPrimary
 373       }//while
 374     }// if pdg = 22
 375   }//while
 376
 377 }