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  *
  21  */
  22
  23 //_________________________________________________________________________
  24 // Class for the analysis of gamma correlations (gamma-jet,
  25 // gamma-hadron and isolation cut.
  26 // This class contains 3 main methods: one to fill lists of particles (ESDs) comming
  27 //  from the CTS (ITS+TPC) and the calorimeters;  the other one tags a candidate
  28 //  cluster as isolated;  the last one search in the
  29 //  corresponing calorimeter for the highest energy cluster, identified it as
  30 //  prompt photon;
  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
  41 #include <TFile.h>
  42 #include <TParticle.h>
  43 #include <TH2.h>
  44
  45 #include "AliAnaGammaDirect.h"
  46 #include "AliESD.h"
  47 #include "AliESDtrack.h"
  48 #include "AliESDCaloCluster.h"
  49 #include "Riostream.h"
  50 #include "AliLog.h"
  51
  52 ClassImp(AliAnaGammaDirect)
  53
  54 //____________________________________________________________________________
  55   AliAnaGammaDirect::AliAnaGammaDirect(const char *name) :
  56     AliAnalysisTask(name,""), fChain(0), fESD(0),
  57     fOutputContainer(new TObjArray(100)),
  58     fPrintInfo(0), fMinGammaPt(0.),
  59     fCalorimeter(""), fEMCALPID(0),fPHOSPID(0),
  60     fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0),
  61     fNCones(0),fNPtThres(0),fMakeICMethod(0)
  62 {
  63   //Ctor
  64   TList * list = gDirectory->GetListOfKeys() ;
  65   TIter next(list) ;
  66   TH2F * h = 0 ;
  67   Int_t index ;
  68   for (index = 0 ; index < list->GetSize()-1 ; index++) {
  69     //-1 to avoid GammaJet Task
  70     h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
  71     fOutputContainer->Add(h) ;
  72   }
  73
  74   for(Int_t i = 0; i < 10 ; i++){
  75     fConeSizes[i]=0;
  76     fPtThresholds[i]=0;
  77   }
  78
  79   // Input slot #0 works with an Ntuple
  80   DefineInput(0, TChain::Class());
  81   // Output slot #0 writes into a TH1 container
  82   DefineOutput(0,  TObjArray::Class()) ;
  83
  84 }
  85
  86
  87 //____________________________________________________________________________
  88 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
  89   AliAnalysisTask(g), fChain(g.fChain), fESD(g.fESD),
  90   fOutputContainer(g. fOutputContainer),  fPrintInfo(g.fPrintInfo),
  91   fMinGammaPt(g.fMinGammaPt), fCalorimeter(g.fCalorimeter),
  92   fEMCALPID(g.fEMCALPID),fPHOSPID(g.fPHOSPID),
  93   fConeSize(g.fConeSize),
  94   fPtThreshold(g.fPtThreshold),
  95   fPtSumThreshold(g.fPtSumThreshold),
  96   fNCones(g.fNCones),fNPtThres(g.fNPtThres),
  97   fMakeICMethod(g.fMakeICMethod)
  98 {
  99   // cpy ctor
 100   SetName (g.GetName()) ;
 101   SetTitle(g.GetTitle()) ;
 102
 103   for(Int_t i = 0; i < 10 ; i++){
 104     fConeSizes[i]=  g.fConeSizes[i];
 105     fPtThresholds[i]=   g.fPtThresholds[i];
 106   }
 107 }
 108
 109 //____________________________________________________________________________
 110 AliAnaGammaDirect::~AliAnaGammaDirect()
 111 {
 112   // Remove all pointers
 113   fOutputContainer->Clear() ;
 114   delete fOutputContainer ;
 115   delete fhNGamma    ;
 116   delete fhPhiGamma  ;
 117   delete fhEtaGamma   ;
 118   delete fhPtCandidate ;
 119   delete [] fhPtThresIsolated ;
 120   delete [] fhPtSumIsolated ;
 121
 122 }
 123
 124 //____________________________________________________________________________
 125 void AliAnaGammaDirect::CreateParticleList(TClonesArray * pl,
 126                                         TClonesArray * plCTS,
 127                                         TClonesArray * plEMCAL,
 128                                         TClonesArray * plPHOS){
 129
 130   //Create a list of particles from the ESD. These particles have been measured
 131   //by the Central Tracking system (TPC+ITS), PHOS and EMCAL
 132
 133   Int_t index = pl->GetEntries() ;
 134   Int_t npar  = 0 ;
 135   Float_t *pid = new Float_t[AliPID::kSPECIESN];
 136   AliDebug(3,"Fill particle lists");
 137   if(fPrintInfo)
 138     AliInfo(Form("fCalorimeter %s",fCalorimeter.Data()));
 139
 140   Double_t v[3] ; //vertex ;
 141   fESD->GetVertex()->GetXYZ(v) ;
 142
 143   //########### PHOS ##############
 144   if(fCalorimeter == "PHOS"){
 145
 146     Int_t begphos = fESD->GetFirstPHOSCluster();
 147     Int_t endphos = fESD->GetFirstPHOSCluster() +
 148       fESD->GetNumberOfPHOSClusters() ;
 149     Int_t indexNePHOS = plPHOS->GetEntries() ;
 150     AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));
 151
 152     if(fCalorimeter == "PHOS"){
 153       for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
 154         AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
 155
 156         //Create a TParticle to fill the particle list
 157
 158         Float_t en = clus->GetClusterEnergy() ;
 159         Float_t *p = new Float_t();
 160         clus->GetGlobalPosition(p) ;
 161         TVector3 pos(p[0],p[1],p[2]) ;
 162         Double_t phi  = pos.Phi();
 163         Double_t theta= pos.Theta();
 164         Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
 165         Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
 166         Double_t pz = en*TMath::Cos(theta);
 167
 168         TParticle * particle = new TParticle() ;
 169         particle->SetMomentum(px,py,pz,en) ;
 170         AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
 171
 172         //Select only photons
 173
 174         pid=clus->GetPid();
 175         //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
 176         if( !fPHOSPID)
 177           new((*plPHOS)[indexNePHOS++])   TParticle(*particle) ;
 178         else if( pid[AliPID::kPhoton] > 0.75)
 179           new((*plPHOS)[indexNePHOS++])   TParticle(*particle) ;
 180       }
 181     }
 182   }
 183
 184   //########### CTS (TPC+ITS) #####################
 185   Int_t begtpc   = 0 ;
 186   Int_t endtpc   = fESD->GetNumberOfTracks() ;
 187   Int_t indexCh  = plCTS->GetEntries() ;
 188   AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
 189
 190   for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
 191     AliESDtrack * track = fESD->GetTrack(npar) ; // retrieve track from esd
 192     //We want tracks fitted in the detectors:
 193     ULong_t status=AliESDtrack::kTPCrefit;
 194     status|=AliESDtrack::kITSrefit;
 195
 196     //We want tracks whose PID bit is set:
 197     //     ULong_t status =AliESDtrack::kITSpid;
 198     //     status|=AliESDtrack::kTPCpid;
 199
 200     if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
 201       // Do something with the tracks which were successfully
 202       // re-fitted
 203       Double_t en = 0; //track ->GetTPCsignal() ;
 204       Double_t mom[3];
 205       track->GetPxPyPz(mom) ;
 206       Double_t px = mom[0];
 207       Double_t py = mom[1];
 208       Double_t pz = mom[2]; //Check with TPC people if this is correct.
 209       Int_t pdg = 11; //Give any charged PDG code, in this case electron.
 210       //I just want to tag the particle as charged
 211        TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
 212                                                  px, py, pz, en, v[0], v[1], v[2], 0);
 213
 214       //TParticle * particle = new TParticle() ;
 215       //particle->SetMomentum(px,py,pz,en) ;
 216
 217       new((*plCTS)[indexCh++])       TParticle(*particle) ;
 218       new((*pl)[index++])           TParticle(*particle) ;
 219     }
 220   }
 221
 222   //################ EMCAL ##############
 223
 224   Int_t begem = fESD->GetFirstEMCALCluster();
 225   Int_t endem = fESD->GetFirstEMCALCluster() +
 226     fESD->GetNumberOfEMCALClusters() ;
 227   Int_t indexNe  = plEMCAL->GetEntries() ;
 228
 229   AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
 230
 231     for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
 232       AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
 233       Int_t clustertype= clus->GetClusterType();
 234       if(clustertype == AliESDCaloCluster::kClusterv1){
 235         Float_t en = clus->GetClusterEnergy() ;
 236         Float_t *p = new Float_t();
 237         clus->GetGlobalPosition(p) ;
 238         TVector3 pos(p[0],p[1],p[2]) ;
 239         Double_t phi  = pos.Phi();
 240         Double_t theta= pos.Theta();
 241         Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
 242         Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
 243         Double_t pz = en*TMath::Cos(theta);
 244
 245         pid=clus->GetPid();
 246         if(fCalorimeter == "EMCAL")
 247           {
 248             TParticle * particle = new TParticle() ;
 249             particle->SetMomentum(px,py,pz,en) ;
 250             AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
 251             if(!fEMCALPID) //Only identified particles
 252               new((*plEMCAL)[indexNe++])       TParticle(*particle) ;
 253             else if(pid[AliPID::kPhoton] > 0.75)
 254               new((*plEMCAL)[indexNe++])       TParticle(*particle) ;
 255           }
 256         else
 257           {
 258             Int_t pdg = 0;
 259             if(fEMCALPID)
 260               {
 261                 if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
 262                   pdg = 22;
 263                 else if( pid[AliPID::kPi0] > 0.75)
 264                   pdg = 111;
 265               }
 266             else
 267               pdg = 22; //No PID, assume all photons
 268
 269             TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
 270                                                  px, py, pz, en, v[0], v[1], v[2], 0);
 271             AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
 272
 273             new((*plEMCAL)[indexNe++])       TParticle(*particle) ;
 274             new((*pl)[index++])           TParticle(*particle) ;
 275           }
 276       }
 277     }
 278
 279     AliDebug(3,"Particle lists filled");
 280
 281 }
 282
 283
 284
 285 //____________________________________________________________________________
 286 void AliAnaGammaDirect::Exec(Option_t *)
 287 {
 288
 289   // Processing of one event
 290
 291   //Get ESDs
 292   Long64_t entry = fChain->GetReadEntry() ;
 293
 294   if (!fESD) {
 295     AliError("fESD is not connected to the input!") ;
 296     return ;
 297   }
 298
 299   if (fPrintInfo)
 300      AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(fChain))->GetFile()->GetName(), entry)) ;
 301
 302   //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
 303
 304   TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
 305   TClonesArray * plCTS         = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
 306   TClonesArray * plNe         = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter (EMCAL)
 307   TClonesArray * plPHOS     = new TClonesArray("TParticle",1000);  // All particles measured in PHOS as Gamma calorimeter
 308   TClonesArray * plEMCAL   = new TClonesArray("TParticle",1000);  // All particles measured in EMCAL as Gamma calorimeter
 309
 310   TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
 311   TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
 312
 313   //Fill lists with photons, neutral particles and charged particles
 314   //look for the highest energy photon in the event inside fCalorimeter
 315
 316   AliDebug(2, "Fill particle lists, get prompt gamma");
 317
 318   //Fill particle lists
 319   CreateParticleList(particleList, plCTS,plEMCAL,plPHOS);
 320
 321   if(fCalorimeter == "PHOS")
 322     plNe = plPHOS;
 323   if(fCalorimeter == "EMCAL")
 324     plNe = plEMCAL;
 325
 326
 327   //_______________Analysis 1__________________________
 328   //Look for the prompt photon in the selected calorimeter
 329   if(fMakeICMethod < 3){
 330
 331     Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
 332     //Search highest energy prompt gamma in calorimeter
 333     GetPromptGamma(plNe,  plCTS, pGamma, iIsInPHOSorEMCAL) ;
 334
 335     AliDebug(1, Form("Is Gamma in %s? %d",fCalorimeter.Data(),iIsInPHOSorEMCAL));
 336
 337     //If there is any photon candidate in fCalorimeter
 338     if(iIsInPHOSorEMCAL){
 339       if (fPrintInfo)
 340         AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
 341
 342     }//Gamma in Calo
 343
 344   }//Analysis 1
 345
 346
 347   //_______________Analysis 2__________________________
 348   //Look for the prompt photon in the selected calorimeter
 349   //Isolation Cut Analysis for both methods and different pt cuts and cones
 350   if(fMakeICMethod == 3){
 351
 352     for(Int_t ipr = 0; ipr < plNe->GetEntries() ; ipr ++ ){
 353       TParticle * pCandidate = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
 354
 355       if(pCandidate->Pt() > fMinGammaPt){
 356
 357         Bool_t  icPtThres   = kFALSE;
 358         Bool_t  icPtSum     = kFALSE;
 359
 360         Float_t ptC             = pCandidate->Pt() ;
 361
 362         fhPtCandidate->Fill(ptC);
 363
 364         for(Int_t icone = 0; icone<fNCones; icone++){
 365           fConeSize  = fConeSizes[icone] ;
 366           Float_t coneptsum = 0 ;
 367           for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
 368             fPtThreshold   =  fPtThresholds[ipt] ;
 369             MakeIsolationCut(plCTS,plNe, pCandidate, ipr, icPtThres, icPtSum,coneptsum);
 370             AliDebug(4,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
 371                             pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
 372
 373             fhPtThresIsolated[icone][ipt]->Fill(ptC);
 374           }//pt thresh loop
 375           fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
 376         }//cone size loop
 377       }//min pt candidate
 378     }//candidate loop
 379   }//Analysis 2
 380
 381   AliDebug(2, "End of analysis, delete pointers");
 382
 383   particleList->Delete() ;
 384   plCTS->Delete() ;
 385   plNe->Delete() ;
 386   plPHOS->Delete() ;
 387   pLeading->Delete();
 388   pGamma->Delete();
 389
 390   delete plNe ;
 391   delete plCTS ;
 392   delete particleList ;
 393   //  delete pLeading;
 394   //  delete pGamma;
 395
 396   PostData(0, fOutputContainer);
 397 }
 398
 399
 400 //____________________________________________________________________________
 401 void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &Is) const
 402 {
 403   //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
 404
 405   Double_t pt = 0;
 406   Int_t index = -1;
 407   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
 408     TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
 409
 410     if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
 411       index = ipr ;
 412       pt = particle->Pt();
 413       pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
 414       AliDebug(4,Form("Cluster in calo: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
 415       Is  = kTRUE;
 416     }
 417   }
 418
 419   //Do Isolation?
 420   if(fMakeICMethod && Is)
 421     {
 422       Float_t coneptsum = 0 ;
 423       Bool_t  icPtThres   = kFALSE;
 424       Bool_t  icPtSum     = kFALSE;
 425       MakeIsolationCut(plCTS,pl, pGamma, index,
 426                        icPtThres, icPtSum,coneptsum);
 427       if(fMakeICMethod == 1) //Pt thres method
 428         Is = icPtThres ;
 429       if(fMakeICMethod == 2) //Pt cone sum method
 430         Is = icPtSum ;
 431     }
 432
 433   if(Is){
 434     AliDebug(3,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
 435     AliDebug(3,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
 436     //Fill prompt gamma histograms
 437     fhNGamma  ->Fill(pGamma->Pt());
 438     fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
 439     fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
 440   }
 441   else
 442     AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
 443 }
 444
 445   //____________________________________________________________________________
 446 void AliAnaGammaDirect::Init(const Option_t * )
 447 {
 448   // Initialisation of branch container
 449
 450   AliDebug(2,Form("*** Initialization of %s", GetName())) ;
 451
 452   // Get input data
 453   fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
 454   if (!fChain) {
 455     AliError(Form("Input 0 for %s not found\n", GetName()));
 456     return ;
 457   }
 458
 459   if (!fESD) {
 460     // One should first check if the branch address was taken by some other task
 461     char ** address = (char **)GetBranchAddress(0, "ESD") ;
 462     if (address)
 463       fESD = (AliESD *)(*address) ;
 464     if (!fESD)
 465       fChain->SetBranchAddress("ESD", &fESD) ;
 466   }
 467   // The output objects will be written to
 468   TDirectory * cdir = gDirectory ;
 469   // Open a file for output #0
 470   char outputName[1024] ;
 471   sprintf(outputName, "%s.root", GetName() ) ;
 472   OpenFile(0, outputName , "RECREATE") ;
 473   if (cdir)
 474     cdir->cd() ;
 475
 476  //  //Initialize the parameters of the analysis.
 477   fCalorimeter="PHOS";
 478   fPrintInfo           = kTRUE;
 479   fMinGammaPt  = 10. ;
 480
 481   //Fill particle lists when PID is ok
 482   fEMCALPID = kFALSE;
 483   fPHOSPID = kFALSE;
 484
 485   fConeSize             = 0.2 ;
 486   fPtThreshold         = 2.0;
 487   fPtSumThreshold  = 1.;
 488
 489   fNCones           = 4 ;
 490   fNPtThres         = 4 ;
 491   fConeSizes[0] = 0.1; fConeSizes[0] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
 492   fPtThresholds[0]=1.; fPtThresholds[0]=2.; fPtThresholds[0]=3.; fPtThresholds[0]=4.;
 493
 494   fMakeICMethod = 1; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method, 3 make isolation study
 495
 496   //Initialization of histograms
 497   MakeHistos() ;
 498 }
 499
 500 //__________________________________________________________________
 501 void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
 502                                            TClonesArray * plNe,
 503                                            TParticle * pCandidate,
 504                                            Int_t index,
 505                                            Bool_t  &icmpt,  Bool_t  &icms,
 506                                            Float_t &coneptsum) const
 507 {
 508   //Search in cone around a candidate particle if it is isolated
 509   Float_t phiC  = pCandidate->Phi() ;
 510   Float_t etaC = pCandidate->Eta() ;
 511   Float_t pt     = -100. ;
 512   Float_t eta   = -100.  ;
 513   Float_t phi    = -100.  ;
 514   Float_t rad   = -100 ;
 515   Int_t    n        = 0 ;
 516   TParticle * particle  = new TParticle();
 517
 518   coneptsum = 0;
 519   icmpt = kFALSE;
 520   icms  = kFALSE;
 521
 522   //Check charged particles in cone.
 523   for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
 524     particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
 525     pt    = particle->Pt();
 526     eta  = particle->Eta();
 527     phi  = particle->Phi() ;
 528
 529     //Check if there is any particle inside cone with pt larger than  fPtThreshold
 530     rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
 531                       TMath::Power((phi-phiC),2));
 532     if(rad<fConeSize){
 533       AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
 534       coneptsum+=pt;
 535       if(pt > fPtThreshold ) n++;
 536     }
 537   }// charged particle loop
 538
 539   //Check neutral particles in cone.
 540   for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
 541     if(ipr != index){//Do not count the candidate
 542       particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
 543       pt    = particle->Pt();
 544       eta  = particle->Eta();
 545       phi  = particle->Phi() ;
 546
 547       //Check if there is any particle inside cone with pt larger than  fPtThreshold
 548       rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
 549                         TMath::Power((phi-phiC),2));
 550       if(rad<fConeSize){
 551         AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
 552         coneptsum+=pt;
 553         if(pt > fPtThreshold ) n++;
 554       }
 555     }
 556   }// neutral particle loop
 557
 558   if(n == 0)
 559     icmpt =  kTRUE ;
 560   if(coneptsum < fPtSumThreshold)
 561     icms  =  kTRUE ;
 562
 563 }
 564
 565 //___________________________________________________________________
 566 void AliAnaGammaDirect::MakeHistos()
 567 {
 568   // Create histograms to be saved in output file and
 569   // stores them in fOutputContainer
 570
 571   fOutputContainer = new TObjArray(10000) ;
 572
 573   //Histograms of highest gamma identified in Event
 574   fhNGamma  = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120);
 575   fhNGamma->SetYTitle("N");
 576   fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
 577   fOutputContainer->Add(fhNGamma) ;
 578
 579   fhPhiGamma  = new TH2F
 580     ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
 581   fhPhiGamma->SetYTitle("#phi");
 582   fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
 583   fOutputContainer->Add(fhPhiGamma) ;
 584
 585   fhEtaGamma  = new TH2F
 586     ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
 587   fhEtaGamma->SetYTitle("#eta");
 588   fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
 589   fOutputContainer->Add(fhEtaGamma) ;
 590
 591   if( fMakeICMethod== 3 ){
 592
 593     //Isolation cut histograms
 594     fhPtCandidate  = new TH1F
 595       ("PtCandidate","p_{T} of candidate particles for isolation",240,0,120);
 596     fhPtCandidate->SetXTitle("p_{T} (GeV/c)");
 597     fOutputContainer->Add(fhPtCandidate) ;
 598
 599     char name[128];
 600     char title[128];
 601     for(Int_t icone = 0; icone<fNCones; icone++){
 602       sprintf(name,"PtSumIsolated_Cone_%d",icone);
 603       sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 604       fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
 605       fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 606       fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
 607       fOutputContainer->Add(fhPtSumIsolated[icone]) ;
 608
 609       for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
 610         sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
 611         sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 612         fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
 613         fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 614         fOutputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
 615       }//icone loop
 616   }//ipt loop
 617   }
 618 }
 619
 620 void AliAnaGammaDirect::Print(const Option_t * opt) const
 621 {
 622
 623   //Print some relevant parameters set for the analysis
 624   if(! opt)
 625     return;
 626
 627   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 628   printf("Cone Size               =     %f\n", fConeSize) ;
 629   printf("pT threshold           =     %f\n", fPtThreshold) ;
 630   printf("pT sum threshold   =     %f\n", fPtSumThreshold) ;
 631   printf("Min Gamma pT      =     %f\n", fMinGammaPt) ;
 632   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
 633 }
 634
 635 void AliAnaGammaDirect::Terminate(Option_t *)
 636 {
 637    // The Terminate() function is the last function to be called during
 638    // a query. It always runs on the client, it can be used to present
 639    // the results graphically or save the results to file.
 640
 641
 642 }