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