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.8  2007/10/29 13:48:42  gustavo
  21  * Corrected coding violations
  22  *
  23  * Revision 1.6  2007/08/17 12:40:04  schutz
  24  * New analysis classes by Gustavo Conesa
  25  *
  26  * Revision 1.4.4.4  2007/07/26 10:32:09  schutz
  27  * new analysis classes in the the new analysis framework
  28  *
  29  *
  30  */
  31
  32 //_________________________________________________________________________
  33 // Class for the prompt gamma analysis, isolation cut
  34 //
  35 //  Class created from old AliPHOSGammaJet
  36 //  (see AliRoot versions previous Release 4-09)
  37 //
  38 //*-- Author: Gustavo Conesa (LNF-INFN)
  39 //////////////////////////////////////////////////////////////////////////////
  40
  41
  42 // --- ROOT system ---
  43 #include <TParticle.h>
  44 #include <TH2.h>
  45 #include <TList.h>
  46 #include "Riostream.h"
  47 #include "TROOT.h"
  48
  49 // --- AliRoot system ---
  50 #include "AliAnaGammaDirect.h"
  51 #include "AliLog.h"
  52
  53 ClassImp(AliAnaGammaDirect)
  54
  55 //____________________________________________________________________________
  56   AliAnaGammaDirect::AliAnaGammaDirect() :
  57     TObject(),
  58     fMinGammaPt(0.),
  59     fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0),
  60     fICMethod(0), fAnaMC(0), fIsolatePi0(0),
  61     fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0),
  62     fntuplePrompt(0),
  63     //kSeveralIC
  64     fNCones(0),fNPtThres(0), fConeSizes(),  fPtThresholds(),
  65     fhPtThresIsolated(), fhPtSumIsolated(), fntSeveralIC()
  66 {
  67   //default ctor
  68
  69   //Initialize parameters
  70   InitParameters();
  71
  72 }
  73
  74 //____________________________________________________________________________
  75 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
  76   TObject(g),
  77   fMinGammaPt(g.fMinGammaPt),
  78   fConeSize(g.fConeSize),
  79   fPtThreshold(g.fPtThreshold),
  80   fPtSumThreshold(g.fPtSumThreshold),
  81   fICMethod(g.fICMethod),
  82   fAnaMC( g.fAnaMC),
  83   fIsolatePi0(g.fIsolatePi0),
  84   fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),
  85   fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),
  86   fntuplePrompt(g.fntuplePrompt),
  87   //kSeveralIC
  88   fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(),
  89   fhPtThresIsolated(), fhPtSumIsolated()
  90 {
  91   // cpy ctor
  92
  93   //kSeveralIC
  94   for(Int_t i = 0; i < fNCones ; i++){
  95     fConeSizes[i] =  g.fConeSizes[i];
  96     fhPtSumIsolated[i] = g.fhPtSumIsolated[i];
  97     fntSeveralIC[i]= g.fntSeveralIC[i];
  98     for(Int_t j = 0; j < fNPtThres ; j++)
  99       fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j];
 100   }
 101
 102   for(Int_t i = 0; i < fNPtThres ; i++)
 103     fPtThresholds[i]=   g.fPtThresholds[i];
 104
 105
 106 }
 107
 108 //_________________________________________________________________________
 109 AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & source)
 110 {
 111   // assignment operator
 112
 113   if(&source == this) return *this;
 114
 115   fMinGammaPt = source.fMinGammaPt ;
 116   fConeSize = source.fConeSize ;
 117   fPtThreshold = source.fPtThreshold ;
 118   fPtSumThreshold = source.fPtSumThreshold ;
 119   fICMethod = source.fICMethod ;
 120   fAnaMC = source.fAnaMC ;
 121   fIsolatePi0 =  source.fIsolatePi0 ;
 122
 123   fhNGamma = source.fhNGamma ;
 124   fhPhiGamma = source.fhPhiGamma ;
 125   fhEtaGamma = source.fhEtaGamma ;
 126   fhConeSumPt = source.fhConeSumPt ;
 127
 128   fntuplePrompt = source.fntuplePrompt ;
 129
 130   //kSeveralIC
 131   fNCones = source.fNCones ;
 132   fNPtThres = source.fNPtThres ;
 133
 134   for(Int_t i = 0; i < fNCones ; i++){
 135     fConeSizes[i] =  source.fConeSizes[i];
 136     fhPtSumIsolated[i] = source.fhPtSumIsolated[i] ;
 137     fntSeveralIC[i]= source.fntSeveralIC[i];
 138     for(Int_t j = 0; j < fNPtThres ; j++)
 139       fhPtThresIsolated[i][j] = source.fhPtThresIsolated[i][j] ;
 140   }
 141
 142   for(Int_t i = 0; i < fNPtThres ; i++)
 143     fPtThresholds[i]=   source.fPtThresholds[i];
 144
 145   return *this;
 146
 147 }
 148
 149 //____________________________________________________________________________
 150 AliAnaGammaDirect::~AliAnaGammaDirect()
 151 {
 152   // Remove all pointers except analysis output pointers.
 153
 154 }
 155
 156 //________________________________________________________________________
 157 TList *  AliAnaGammaDirect::GetCreateOutputObjects()
 158 {
 159
 160   // Create histograms to be saved in output file and
 161   // store them in outputContainer
 162   TList * outputContainer = new TList() ;
 163   outputContainer->SetName("DirectGammaHistos") ;
 164
 165   //Histograms of highest gamma identified in Event
 166   fhNGamma  = new TH1F("NGamma","Number of #gamma over calorimeter",240,0,120);
 167   fhNGamma->SetYTitle("N");
 168   fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
 169   outputContainer->Add(fhNGamma) ;
 170
 171   fhPhiGamma  = new TH2F
 172     ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
 173   fhPhiGamma->SetYTitle("#phi");
 174   fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
 175   outputContainer->Add(fhPhiGamma) ;
 176
 177   fhEtaGamma  = new TH2F
 178     ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
 179   fhEtaGamma->SetYTitle("#eta");
 180   fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
 181   outputContainer->Add(fhEtaGamma) ;
 182
 183   fhConeSumPt  = new TH2F
 184     ("ConePtSum","#Sigma p_{T}  in cone ",200,0,120,100,0,100);
 185   fhConeSumPt->SetYTitle("#Sigma p_{T}");
 186   fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
 187   outputContainer->Add(fhConeSumPt) ;
 188
 189   //NTUPLE
 190   fntuplePrompt = new TNtuple("ntuplePromptGamma", "Tree of prompt #gamma", "ptcluster:phicluster:etacluster:pdg:status:ptprimary:phiprimary:etaprimary:pdgprimary:statusprimary");
 191   outputContainer->Add(fntuplePrompt) ;
 192
 193   if(fICMethod == kSeveralIC){
 194     char name[128];
 195     char title[128];
 196     for(Int_t icone = 0; icone<fNCones; icone++){
 197       sprintf(name,"PtSumIsolated_Cone_%d",icone);
 198       sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 199       fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
 200       fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 201       fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
 202       outputContainer->Add(fhPtSumIsolated[icone]) ;
 203
 204       sprintf(name,"nt_Cone_%d",icone);
 205       sprintf(title,"ntuple for cone size %d",icone);
 206       fntSeveralIC[icone] = new TNtuple(name, title, "ptcand:phicand:etacand:ptsum:type:ncone0:ncone1:ncone2:ncone3:ncone4:ncone5");
 207       outputContainer->Add(fntSeveralIC[icone]) ;
 208
 209       for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
 210         sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
 211         sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 212         fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
 213         fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 214         outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
 215       }//icone loop
 216     }//ipt loop
 217   }
 218
 219 //  gROOT->cd();
 220
 221   return outputContainer ;
 222
 223 }
 224
 225 //____________________________________________________________________________
 226 void AliAnaGammaDirect::GetPromptGamma(TClonesArray * plCalo, TClonesArray * plCTS, TClonesArray * plPrimCalo, TParticle *pGamma, Bool_t &found) const
 227 {
 228   //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
 229
 230   Double_t pt = 0;
 231   Int_t n = 0;
 232   Int_t index = -1;
 233   Float_t coneptsum = 0 ;
 234
 235   for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
 236     TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
 237
 238     if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt) &&
 239        (particle->GetPdgCode() == 22 || (fIsolatePi0 && particle->GetPdgCode() == 111))){
 240       index = ipr ;
 241       pt = particle->Pt();
 242       pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
 243       pGamma->SetStatusCode(particle->GetStatusCode());
 244       pGamma->SetPdgCode(particle->GetPdgCode());
 245       found  = kTRUE;
 246     }
 247   }
 248
 249   //Do Isolation?
 250   if( ( fICMethod == kPtIC  ||  fICMethod == kSumPtIC )  && found)
 251     {
 252       Bool_t  icPtThres   = kFALSE;
 253       Bool_t  icPtSum     = kFALSE;
 254       MakeIsolationCut(plCTS,plCalo, pGamma, index,n,
 255                        icPtThres, icPtSum,coneptsum);
 256       if(fICMethod == kPtIC) //Pt thres method
 257         found = icPtThres ;
 258       if(fICMethod == kSumPtIC) //Pt cone sum method
 259         found = icPtSum ;
 260     }
 261
 262   if(found){
 263     AliDebug(1,Form("Cluster with p_{T} larger than the pt treshold %f found in calorimeter ", fMinGammaPt)) ;
 264     AliDebug(1,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
 265     //Fill prompt gamma histograms
 266     Float_t ptcluster = pGamma->Pt();
 267     Float_t phicluster = pGamma->Phi();
 268     Float_t etacluster = pGamma->Eta();
 269     Int_t statuscluster = pGamma->GetStatusCode();
 270     Int_t pdgcluster = pGamma->GetPdgCode();
 271
 272     fhNGamma   ->Fill(ptcluster);
 273     fhPhiGamma ->Fill(ptcluster,phicluster);
 274     fhEtaGamma ->Fill(ptcluster,etacluster);
 275     fhConeSumPt->Fill(ptcluster,coneptsum);
 276
 277     Float_t ptprimary = 0 ;
 278     Float_t phiprimary = 0 ;
 279     Float_t etaprimary = 0 ;
 280     Int_t pdgprimary = 0 ;
 281     Int_t statusprimary = 0 ;
 282
 283     if(fAnaMC){
 284       TParticle * primary = dynamic_cast<TParticle *>(plPrimCalo->At(index)) ;
 285       ptprimary = primary->Pt();
 286       phiprimary = primary->Phi();
 287       etaprimary = primary->Eta();
 288       pdgprimary =  TMath::Abs(primary->GetPdgCode()) ;
 289       statusprimary = primary->GetStatusCode(); // = 2 means decay gamma!!!
 290
 291       AliDebug(1, Form("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f",ptcluster,pdgprimary,ptprimary));
 292       //printf("Identified prompt Gamma pT %2.2f; Primary, pdg %d, pT %2.2f \n",ptcluster,pdgprimary,ptprimary);
 293     }
 294
 295     //Fill ntuple with cluster / MC data
 296 //    gROOT->cd();
 297     fntuplePrompt->Fill(ptcluster,phicluster,etacluster,pdgcluster,statuscluster,ptprimary,phiprimary, etaprimary,pdgprimary,statusprimary);
 298   }
 299   else
 300     AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
 301 }
 302
 303   //____________________________________________________________________________
 304 void AliAnaGammaDirect::InitParameters()
 305 {
 306
 307   //Initialize the parameters of the analysis.
 308   fMinGammaPt  = 5. ;
 309
 310   //Fill particle lists when PID is ok
 311   fConeSize             = 0.4 ;
 312   fPtThreshold         = 1.;
 313   fPtSumThreshold  = 1.;
 314
 315   fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
 316   fAnaMC = kFALSE ;
 317   fIsolatePi0 = kFALSE ;
 318  //-----------kSeveralIC-----------------
 319   fNCones           = 5 ;
 320   fNPtThres         = 6 ;
 321   fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4; fConeSizes[4] = 0.5;
 322   fPtThresholds[0]=0.; fPtThresholds[1]=1.; fPtThresholds[2]=2.; fPtThresholds[3]=3.; fPtThresholds[4]=4.;fPtThresholds[5]=5.;
 323
 324 }
 325
 326 //__________________________________________________________________
 327 void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
 328                                           TClonesArray * plNe,
 329                                           TParticle * pCandidate,
 330                                           Int_t index, Int_t & n,
 331                                           Bool_t  &icmpt,  Bool_t  &icms,
 332                                           Float_t &coneptsum) const
 333 {
 334   //Search in cone around a candidate particle if it is isolated
 335   Float_t phiC  = pCandidate->Phi() ;
 336   Float_t etaC = pCandidate->Eta() ;
 337   Float_t pt     = -100. ;
 338   Float_t eta   = -100.  ;
 339   Float_t phi    = -100.  ;
 340   Float_t rad   = -100 ;
 341   TParticle * particle  = new TParticle;
 342
 343   n = 0 ;
 344   coneptsum = 0.;
 345   icmpt = kFALSE;
 346   icms  = kFALSE;
 347
 348   //Check charged particles in cone.
 349   for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
 350     particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
 351     pt    = particle->Pt();
 352     eta  = particle->Eta();
 353     phi  = particle->Phi() ;
 354
 355     //Check if there is any particle inside cone with pt larger than  fPtThreshold
 356     rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
 357     if(rad < fConeSize){
 358       AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
 359       coneptsum+=pt;
 360       if(pt > fPtThreshold ) n++;
 361     }
 362   }// charged particle loop
 363
 364   //Check neutral particles in cone.
 365   for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
 366     if(ipr != index){//Do not count the candidate
 367       particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
 368       pt    = particle->Pt();
 369       eta  = particle->Eta();
 370       phi  = particle->Phi() ;
 371
 372       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 373       rad = TMath::Sqrt((eta-etaC)*(eta-etaC)+ (phi-phiC)*(phi-phiC));
 374       if(rad < fConeSize){
 375         AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
 376         coneptsum+=pt;
 377         if(pt > fPtThreshold ) n++;
 378       }
 379     }
 380   }// neutral particle loop
 381
 382   if(n == 0)
 383     icmpt =  kTRUE ;
 384   if(coneptsum <= fPtSumThreshold)
 385     icms  =  kTRUE ;
 386
 387 }
 388
 389 //__________________________________________________________________
 390 void  AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS)
 391 {
 392   //Isolation Cut Analysis for both methods and different pt cuts and cones
 393   if (fICMethod != kSeveralIC)
 394     AliFatal("Remember to set in config file: directGamma->SetICMethod(kSeveralIC)");
 395   Int_t type = 0;
 396
 397   //Search maximum energy photon in the event
 398   Double_t ptC = 0;
 399   Int_t index = -1;
 400   TParticle * pCandidate = new TParticle();
 401   Bool_t found = kFALSE;
 402
 403   for(Int_t ipr = 0;ipr < plCalo->GetEntries() ; ipr ++ ){
 404     TParticle * particle = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
 405
 406     if((particle->Pt() > fMinGammaPt) && (particle->Pt() > ptC) &&
 407        (particle->GetPdgCode() == 22 ||  (fIsolatePi0 && particle->GetPdgCode() == 111))){
 408       index = ipr ;
 409       ptC = particle->Pt();
 410       pCandidate = particle ;
 411       found  = kTRUE;
 412     }
 413   }
 414
 415   //If there is a large cluster, larger than threshold, study isolation cut
 416   if(found){
 417
 418     fhNGamma->Fill(ptC);
 419     fhPhiGamma->Fill(ptC,pCandidate->Phi());
 420     fhEtaGamma->Fill(ptC,pCandidate->Eta());
 421
 422     Int_t ncone[10][10];//[fNCones][fNPtThres];
 423     Bool_t  icPtThres   = kFALSE;
 424     Bool_t  icPtSum     = kFALSE;
 425
 426     for(Int_t icone = 0; icone<fNCones; icone++){
 427       fConeSize=fConeSizes[icone] ;
 428       Float_t coneptsum = 0 ;
 429       for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
 430         ncone[icone][ipt]=0;
 431         fPtThreshold=fPtThresholds[ipt] ;
 432         MakeIsolationCut(plCTS,plCalo, pCandidate, index,
 433                          ncone[icone][ipt], icPtThres, icPtSum,coneptsum);
 434         AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
 435                         pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
 436 //      if(ptC >15 && ptC < 25 && (icPtThres || icPtSum) && ipt ==0){
 437 //        printf("R %0.1f, ptthres %1.1f, ptsum %1.1f, Candidate pt %2.2f,  eta %2.2f, phi %2.2f, pt in cone %2.2f, Isolated? ICPt %d, ICSum %d\n",
 438 //               fConeSize,  fPtThreshold, fPtSumThreshold, ptC, pCandidate->Eta(), pCandidate->Phi()*TMath::RadToDeg(), coneptsum, icPtThres, icPtSum);
 439 //        //cout<<"mother label "<<pCandidate->GetMother(0)<<endl;
 440 //      }
 441
 442         fhPtThresIsolated[icone][ipt]->Fill(ptC);
 443       }//pt thresh loop
 444       fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
 445 //      gROOT->cd();
 446       fntSeveralIC[icone]->Fill(ptC,pCandidate->Phi(),pCandidate->Eta(), coneptsum,type,ncone[icone][0],ncone[icone][1],ncone[icone][2],ncone[icone][3],ncone[icone][4],ncone[icone][5]);
 447     }//cone size loop
 448   }//found high energy gamma in the event
 449 }
 450
 451 //__________________________________________________________________
 452 void AliAnaGammaDirect::Print(const Option_t * opt) const
 453 {
 454
 455   //Print some relevant parameters set for the analysis
 456   if(! opt)
 457     return;
 458
 459   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 460
 461   printf("Min Gamma pT      =     %f\n",  fMinGammaPt) ;
 462   printf("IC method               =     %d\n", fICMethod) ;
 463   printf("Cone Size               =     %f\n", fConeSize) ;
 464    if(fICMethod == kPtIC) printf("pT threshold           =     %f\n", fPtThreshold) ;
 465    if(fICMethod == kSumPtIC) printf("pT sum threshold   =     %f\n", fPtSumThreshold) ;
 466
 467   if(fICMethod == kSeveralIC){
 468     printf("N Cone Sizes               =     %d\n", fNCones) ;
 469     printf("N pT thresholds           =     %d\n", fNPtThres) ;
 470     printf("Cone Sizes                  =    \n") ;
 471     for(Int_t i = 0; i < fNCones; i++)
 472       printf("   %f;",  fConeSizes[i]) ;
 473     printf("    \n") ;
 474     for(Int_t i = 0; i < fNPtThres; i++)
 475       printf("   %f;",  fPtThresholds[i]) ;
 476   }
 477
 478   printf("    \n") ;
 479
 480 }