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