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