PWG4/AliAnaGammaDirect.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.4.4.4  2007/07/26 10:32:09  schutz
  21  * new analysis classes in the the new analysis framework
  22  *
  23  *
  24  */
  25
  26 //_________________________________________________________________________
  27 // Class for the prompt gamma analysis, isolation cut
  28 //
  29 //  Class created from old AliPHOSGammaJet
  30 //  (see AliRoot versions previous Release 4-09)
  31 //
  32 //*-- Author: Gustavo Conesa (LNF-INFN)
  33 //////////////////////////////////////////////////////////////////////////////
  34
  35
  36 // --- ROOT system ---
  37 #include <TParticle.h>
  38 #include <TH2.h>
  39 #include <TList.h>
  40 #include "AliAnaGammaDirect.h"
  41 #include "Riostream.h"
  42 #include "AliLog.h"
  43
  44 ClassImp(AliAnaGammaDirect)
  45
  46 //____________________________________________________________________________
  47   AliAnaGammaDirect::AliAnaGammaDirect() :
  48     TObject(),
  49     fMinGammaPt(0.),
  50     fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0),
  51     fICMethod(0),fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0),
  52     //kSeveralIC
  53     fNCones(0),fNPtThres(0), fConeSizes(),  fPtThresholds(),
  54     fhPtThresIsolated(), fhPtSumIsolated()
  55 {
  56   //default ctor
  57
  58   //Initialize parameters
  59   InitParameters();
  60
  61 }
  62
  63 //____________________________________________________________________________
  64 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
  65   TObject(g),
  66   fMinGammaPt(g.fMinGammaPt),
  67   fConeSize(g.fConeSize),
  68   fPtThreshold(g.fPtThreshold),
  69   fPtSumThreshold(g.fPtSumThreshold),
  70   fICMethod(g.fICMethod),
  71   fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),fhEtaGamma(g.fhEtaGamma),
  72   //kSeveralIC
  73   fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(),
  74   fhPtThresIsolated(), fhPtSumIsolated()
  75 {
  76   // cpy ctor
  77
  78   //kSeveralIC
  79   for(Int_t i = 0; i < fNCones ; i++){
  80     fConeSizes[i] =  g.fConeSizes[i];
  81     fhPtSumIsolated[i] = g.fhPtSumIsolated[i];
  82     for(Int_t j = 0; j < fNPtThres ; j++)
  83       fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j];
  84   }
  85
  86   for(Int_t i = 0; i < fNPtThres ; i++)
  87     fPtThresholds[i]=   g.fPtThresholds[i];
  88 }
  89
  90 //_________________________________________________________________________
  91 AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & source)
  92 {
  93   // assignment operator
  94
  95   if(&source == this) return *this;
  96
  97   fMinGammaPt = source.fMinGammaPt ;
  98   fConeSize = source.fConeSize ;
  99   fPtThreshold = source.fPtThreshold ;
 100   fPtSumThreshold = source.fPtSumThreshold ;
 101   fICMethod = source.fICMethod ;
 102   fhNGamma = source.fhNGamma ;
 103   fhPhiGamma = source.fhPhiGamma ;
 104   fhEtaGamma = source.fhEtaGamma ;
 105
 106   //kSeveralIC
 107   fNCones = source.fNCones ;
 108   fNPtThres = source.fNPtThres ;
 109
 110   for(Int_t i = 0; i < fNCones ; i++){
 111     fConeSizes[i] =  source.fConeSizes[i];
 112     fhPtSumIsolated[i] = source.fhPtSumIsolated[i] ;
 113     for(Int_t j = 0; j < fNPtThres ; j++)
 114       fhPtThresIsolated[i][j] = source.fhPtThresIsolated[i][j] ;
 115   }
 116
 117   for(Int_t i = 0; i < fNPtThres ; i++)
 118     fPtThresholds[i]=   source.fPtThresholds[i];
 119
 120   return *this;
 121
 122 }
 123
 124 //____________________________________________________________________________
 125 AliAnaGammaDirect::~AliAnaGammaDirect()
 126 {
 127   // Remove all pointers
 128
 129   delete fhNGamma    ;
 130   delete fhPhiGamma  ;
 131   delete fhEtaGamma   ;
 132
 133   //kSeveralIC
 134   delete [] fhPtThresIsolated ;
 135   delete [] fhPtSumIsolated ;
 136
 137 }
 138
 139 //________________________________________________________________________
 140 TList *  AliAnaGammaDirect::GetCreateOutputObjects()
 141 {
 142
 143   // Create histograms to be saved in output file and
 144   // store them in outputContainer
 145   TList * outputContainer = new TList() ;
 146   outputContainer->SetName("DirectGammaHistos") ;
 147
 148   //Histograms of highest gamma identified in Event
 149   fhNGamma  = new TH1F("NGamma","Number of #gamma over calorimeter",240,0,120);
 150   fhNGamma->SetYTitle("N");
 151   fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
 152   outputContainer->Add(fhNGamma) ;
 153
 154   fhPhiGamma  = new TH2F
 155     ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
 156   fhPhiGamma->SetYTitle("#phi");
 157   fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
 158   outputContainer->Add(fhPhiGamma) ;
 159
 160   fhEtaGamma  = new TH2F
 161     ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
 162   fhEtaGamma->SetYTitle("#eta");
 163   fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
 164   outputContainer->Add(fhEtaGamma) ;
 165
 166   if(fICMethod == kSeveralIC){
 167     char name[128];
 168     char title[128];
 169     for(Int_t icone = 0; icone<fNCones; icone++){
 170       sprintf(name,"PtSumIsolated_Cone_%d",icone);
 171       sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 172       fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
 173       fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 174       fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
 175       outputContainer->Add(fhPtSumIsolated[icone]) ;
 176
 177       for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
 178         sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
 179         sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 180         fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
 181         fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 182         outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
 183       }//icone loop
 184     }//ipt loop
 185   }
 186
 187   return outputContainer ;
 188
 189 }
 190
 191 //____________________________________________________________________________
 192 void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &found) const
 193 {
 194   //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
 195
 196   Double_t pt = 0;
 197   Int_t index = -1;
 198   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
 199     TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
 200
 201     if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
 202       index = ipr ;
 203       pt = particle->Pt();
 204       pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
 205       found  = kTRUE;
 206     }
 207   }
 208
 209   //Do Isolation?
 210   if( ( fICMethod == kPtIC  ||  fICMethod == kSumPtIC )  && found)
 211     {
 212       Float_t coneptsum = 0 ;
 213       Bool_t  icPtThres   = kFALSE;
 214       Bool_t  icPtSum     = kFALSE;
 215       MakeIsolationCut(plCTS,pl, pGamma, index,
 216                        icPtThres, icPtSum,coneptsum);
 217       if(fICMethod == kPtIC) //Pt thres method
 218         found = icPtThres ;
 219       if(fICMethod == kSumPtIC) //Pt cone sum method
 220         found = icPtSum ;
 221     }
 222
 223   if(found){
 224     AliDebug(1,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
 225     AliDebug(1,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
 226     //Fill prompt gamma histograms
 227     fhNGamma->Fill(pGamma->Pt());
 228     fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
 229     fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
 230   }
 231   else
 232     AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
 233 }
 234
 235   //____________________________________________________________________________
 236 void AliAnaGammaDirect::InitParameters()
 237 {
 238
 239   //Initialize the parameters of the analysis.
 240   fMinGammaPt  = 5. ;
 241
 242   //Fill particle lists when PID is ok
 243   fConeSize             = 0.2 ;
 244   fPtThreshold         = 2.0;
 245   fPtSumThreshold  = 1.;
 246
 247   fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
 248
 249  //-----------kSeveralIC-----------------
 250   fNCones           = 4 ;
 251   fNPtThres         = 4 ;
 252   fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
 253   fPtThresholds[0]=1.; fPtThresholds[1]=2.; fPtThresholds[2]=3.; fPtThresholds[3]=4.;
 254
 255 }
 256
 257 //__________________________________________________________________
 258 void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
 259                                            TClonesArray * plNe,
 260                                            TParticle * pCandidate,
 261                                            Int_t index,
 262                                            Bool_t  &icmpt,  Bool_t  &icms,
 263                                            Float_t &coneptsum) const
 264 {
 265   //Search in cone around a candidate particle if it is isolated
 266   Float_t phiC  = pCandidate->Phi() ;
 267   Float_t etaC = pCandidate->Eta() ;
 268   Float_t pt     = -100. ;
 269   Float_t eta   = -100.  ;
 270   Float_t phi    = -100.  ;
 271   Float_t rad   = -100 ;
 272   Int_t    n        = 0 ;
 273   TParticle * particle  = new TParticle;
 274
 275   coneptsum = 0;
 276   icmpt = kFALSE;
 277   icms  = kFALSE;
 278
 279   //Check charged particles in cone.
 280   for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
 281     particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
 282     pt    = particle->Pt();
 283     eta  = particle->Eta();
 284     phi  = particle->Phi() ;
 285
 286     //Check if there is any particle inside cone with pt larger than  fPtThreshold
 287     rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
 288                       TMath::Power((phi-phiC),2));
 289     if(rad<fConeSize){
 290       AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
 291       coneptsum+=pt;
 292       if(pt > fPtThreshold ) n++;
 293     }
 294   }// charged particle loop
 295
 296   //Check neutral particles in cone.
 297   for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
 298     if(ipr != index){//Do not count the candidate
 299       particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
 300       pt    = particle->Pt();
 301       eta  = particle->Eta();
 302       phi  = particle->Phi() ;
 303
 304       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 305       rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
 306                         TMath::Power((phi-phiC),2));
 307       if(rad<fConeSize){
 308         AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
 309         coneptsum+=pt;
 310         if(pt > fPtThreshold ) n++;
 311       }
 312     }
 313   }// neutral particle loop
 314
 315   if(n == 0)
 316     icmpt =  kTRUE ;
 317   if(coneptsum < fPtSumThreshold)
 318     icms  =  kTRUE ;
 319
 320 }
 321
 322 //__________________________________________________________________
 323 void  AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS)
 324 {
 325   //Isolation Cut Analysis for both methods and different pt cuts and cones
 326
 327   if (fICMethod != kSeveralIC)
 328     AliFatal("Remember to set in config file: directGamma->SetICMethod(kSeveralIC)");
 329
 330   for(Int_t ipr = 0; ipr < plCalo->GetEntries() ; ipr ++ ){
 331     TParticle * pCandidate = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
 332
 333     if(pCandidate->Pt() > fMinGammaPt){
 334
 335       Bool_t  icPtThres   = kFALSE;
 336       Bool_t  icPtSum     = kFALSE;
 337
 338       Float_t ptC      = pCandidate->Pt() ;
 339
 340       fhNGamma->Fill(ptC);
 341       fhPhiGamma->Fill( ptC,pCandidate->Phi());
 342       fhEtaGamma->Fill(ptC,pCandidate->Eta());
 343
 344       for(Int_t icone = 0; icone<fNCones; icone++){
 345         fConeSize=fConeSizes[icone] ;
 346         Float_t coneptsum = 0 ;
 347         for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
 348           fPtThreshold=fPtThresholds[ipt] ;
 349           MakeIsolationCut(plCTS,plCalo, pCandidate, ipr, icPtThres, icPtSum,coneptsum);
 350           AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
 351                           pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
 352
 353           fhPtThresIsolated[icone][ipt]->Fill(ptC);
 354         }//pt thresh loop
 355         fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
 356       }//cone size loop
 357     }//min pt candidate
 358   }//candidate loop
 359 }
 360
 361 void AliAnaGammaDirect::Print(const Option_t * opt) const
 362 {
 363
 364   //Print some relevant parameters set for the analysis
 365   if(! opt)
 366     return;
 367
 368   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 369
 370   printf("Min Gamma pT      =     %f\n",  fMinGammaPt) ;
 371   printf("IC method               =     %d\n", fICMethod) ;
 372   printf("Cone Size               =     %f\n", fConeSize) ;
 373    if(fICMethod == kPtIC) printf("pT threshold           =     %f\n", fPtThreshold) ;
 374    if(fICMethod == kSumPtIC) printf("pT sum threshold   =     %f\n", fPtSumThreshold) ;
 375
 376   if(fICMethod == kSeveralIC){
 377     printf("N Cone Sizes               =     %d\n", fNCones) ;
 378     printf("N pT thresholds           =     %d\n", fNPtThres) ;
 379     printf("Cone Sizes                  =    \n") ;
 380     for(Int_t i = 0; i < fNCones; i++)
 381       printf("   %f;",  fConeSizes[i]) ;
 382     printf("    \n") ;
 383     for(Int_t i = 0; i < fNPtThres; i++)
 384       printf("   %f;",  fPtThresholds[i]) ;
 385   }
 386
 387   printf("    \n") ;
 388
 389 }