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