PWG4/AliAnaGammaDirect.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id$ */
  16
  17 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.9  2007/11/17 16:39:49  gustavo
  21  * removed deleting of not owned data and deleting of histograms which are exported to the output file (MG)
  22  *
  23  * Revision 1.8  2007/10/29 13:48:42  gustavo
  24  * Corrected coding violations
  25  *
  26  * Revision 1.6  2007/08/17 12:40:04  schutz
  27  * New analysis classes by Gustavo Conesa
  28  *
  29  * Revision 1.4.4.4  2007/07/26 10:32:09  schutz
  30  * new analysis classes in the the new analysis framework
  31  *
  32  *
  33  */
  34
  35 //_________________________________________________________________________
  36 // Class for the prompt gamma analysis, isolation cut
  37 //
  38 //  Class created from old AliPHOSGammaJet
  39 //  (see AliRoot versions previous Release 4-09)
  40 //
  41 // -- Author: Gustavo Conesa (LNF-INFN)
  42 //////////////////////////////////////////////////////////////////////////////
  43
  44
  45 // --- ROOT system ---
  46 #include <TParticle.h>
  47 #include <TH2.h>
  48 #include <TList.h>
  49 #include "Riostream.h"
  50
  51 // --- Analysis system ---
  52 #include "AliAnaGammaDirect.h"
  53 #include "AliLog.h"
  54 #include "AliCaloTrackReader.h"
  55 #include "AliIsolationCut.h"
  56 #include "AliNeutralMesonSelection.h"
  57
  58 ClassImp(AliAnaGammaDirect)
  59
  60 //____________________________________________________________________________
  61   AliAnaGammaDirect::AliAnaGammaDirect() :
  62     AliAnaBaseClass(), fDetector(""), fMakeIC(0),  fReMakeIC(0),
  63     fMakeSeveralIC(0), fMakeInvMass(0),
  64     fhPtGamma(0),fhPhiGamma(0),fhEtaGamma(0), fhConeSumPt(0),
  65     //Several IC
  66     fNCones(0),fNPtThresFrac(0), fConeSizes(),  fPtThresholds(),  fPtFractions(),
  67     fhPtThresIsolated(), fhPtFracIsolated(),  fhPtSumIsolated(),
  68     //MC
  69     fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
  70     fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(),  fhPtSumIsolatedPrompt(),
  71     fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
  72     fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(),  fhPtSumIsolatedFragmentation(),
  73     fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
  74     fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(),  fhPtSumIsolatedPi0Decay(),
  75     fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
  76     fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(),  fhPtSumIsolatedOtherDecay(),
  77     fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
  78     fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(),  fhPtSumIsolatedConversion(),
  79     fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0),
  80     fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(),  fhPtSumIsolatedUnknown()
  81 {
  82   //default ctor
  83
  84   //Initialize parameters
  85   InitParameters();
  86
  87   for(Int_t i = 0; i < 5 ; i++){
  88     fConeSizes[i] = 0 ;
  89     fhPtSumIsolated[i] = 0 ;
  90
  91     fhPtSumIsolatedPrompt[i] = 0 ;
  92     fhPtSumIsolatedFragmentation[i] = 0 ;
  93     fhPtSumIsolatedPi0Decay[i] = 0 ;
  94     fhPtSumIsolatedOtherDecay[i] = 0 ;
  95     fhPtSumIsolatedConversion[i] = 0 ;
  96     fhPtSumIsolatedUnknown[i] = 0 ;
  97
  98     for(Int_t j = 0; j < 5 ; j++){
  99       fhPtThresIsolated[i][j] = 0 ;
 100       fhPtFracIsolated[i][j] = 0 ;
 101
 102       fhPtThresIsolatedPrompt[i][j] = 0 ;
 103       fhPtThresIsolatedFragmentation[i][j] = 0 ;
 104       fhPtThresIsolatedPi0Decay[i][j] = 0 ;
 105       fhPtThresIsolatedOtherDecay[i][j] = 0 ;
 106       fhPtThresIsolatedConversion[i][j] = 0 ;
 107       fhPtThresIsolatedUnknown[i][j] = 0 ;
 108
 109       fhPtFracIsolatedPrompt[i][j] = 0 ;
 110       fhPtFracIsolatedFragmentation[i][j] = 0 ;
 111       fhPtFracIsolatedPi0Decay[i][j] = 0 ;
 112       fhPtFracIsolatedOtherDecay[i][j] = 0 ;
 113       fhPtFracIsolatedConversion[i][j] = 0 ;
 114       fhPtFracIsolatedUnknown[i][j] = 0 ;
 115
 116     }
 117   }
 118
 119   for(Int_t i = 0; i < 5 ; i++){
 120     fPtFractions[i]=  0 ;
 121     fPtThresholds[i]= 0 ;
 122   }
 123
 124
 125 }
 126
 127 //____________________________________________________________________________
 128 AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
 129   AliAnaBaseClass(g), fDetector(g.fDetector),
 130   fMakeIC(g.fMakeIC),   fReMakeIC(g.fReMakeIC),
 131   fMakeSeveralIC(g.fMakeSeveralIC),  fMakeInvMass(g.fMakeInvMass),
 132   fhPtGamma(g.fhPtGamma),fhPhiGamma(g.fhPhiGamma),
 133   fhEtaGamma(g.fhEtaGamma), fhConeSumPt(g.fhConeSumPt),
 134   //Several IC
 135   fNCones(g.fNCones),fNPtThresFrac(g.fNPtThresFrac), fConeSizes(), fPtThresholds(),  fPtFractions(),
 136   fhPtThresIsolated(),  fhPtFracIsolated(), fhPtSumIsolated(),
 137   //MC
 138   fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
 139   fhPtThresIsolatedPrompt(), fhPtFracIsolatedPrompt(),  fhPtSumIsolatedPrompt(),
 140   fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
 141   fhPtThresIsolatedFragmentation(), fhPtFracIsolatedFragmentation(),  fhPtSumIsolatedFragmentation(),
 142   fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
 143   fhPtThresIsolatedPi0Decay(), fhPtFracIsolatedPi0Decay(),  fhPtSumIsolatedPi0Decay(),
 144   fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
 145   fhPtThresIsolatedOtherDecay(), fhPtFracIsolatedOtherDecay(),  fhPtSumIsolatedOtherDecay(),
 146   fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
 147   fhPtThresIsolatedConversion(), fhPtFracIsolatedConversion(),  fhPtSumIsolatedConversion(),
 148   fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown),
 149   fhPtThresIsolatedUnknown(), fhPtFracIsolatedUnknown(),  fhPtSumIsolatedUnknown()
 150 {
 151   // cpy ctor
 152
 153   //Several IC
 154   for(Int_t i = 0; i < fNCones ; i++){
 155     fConeSizes[i] =  g.fConeSizes[i];
 156     fhPtSumIsolated[i] = g.fhPtSumIsolated[i];
 157
 158    fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i];
 159    fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i];
 160    fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i];
 161    fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i];
 162    fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i];
 163    fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i];
 164
 165     for(Int_t j = 0; j < fNPtThresFrac ; j++){
 166       fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j];
 167       fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j];
 168
 169       fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j];
 170       fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j];
 171       fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j];
 172       fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j];
 173       fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j];
 174       fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j];
 175
 176       fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j];
 177       fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j];
 178       fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j];
 179       fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j];
 180       fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j];
 181       fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j];
 182
 183     }
 184   }
 185
 186   for(Int_t i = 0; i < fNPtThresFrac ; i++){
 187     fPtFractions[i]=   g.fPtFractions[i];
 188     fPtThresholds[i]=   g.fPtThresholds[i];
 189   }
 190
 191 }
 192
 193 //_________________________________________________________________________
 194 AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & g)
 195 {
 196   // assignment operator
 197
 198   if(&g == this) return *this;
 199
 200   fMakeIC = g.fMakeIC ;
 201   fReMakeIC = g.fReMakeIC ;
 202   fMakeSeveralIC = g.fMakeSeveralIC ;
 203   fMakeInvMass = g.fMakeInvMass ;
 204   fDetector = g.fDetector ;
 205
 206   fhPtGamma = g.fhPtGamma ;
 207   fhPhiGamma = g.fhPhiGamma ;
 208   fhEtaGamma = g.fhEtaGamma ;
 209   fhConeSumPt = g.fhConeSumPt ;
 210
 211   fhPtPrompt = g.fhPtPrompt;
 212   fhPhiPrompt = g.fhPhiPrompt;
 213   fhEtaPrompt = g.fhEtaPrompt;
 214   fhPtFragmentation = g.fhPtFragmentation;
 215   fhPhiFragmentation = g.fhPhiFragmentation;
 216   fhEtaFragmentation = g.fhEtaFragmentation;
 217   fhPtPi0Decay = g.fhPtPi0Decay;
 218   fhPhiPi0Decay = g.fhPhiPi0Decay;
 219   fhEtaPi0Decay = g.fhEtaPi0Decay;
 220   fhPtOtherDecay = g.fhPtOtherDecay;
 221   fhPhiOtherDecay = g.fhPhiOtherDecay;
 222   fhEtaOtherDecay = g.fhEtaOtherDecay;
 223   fhPtConversion = g. fhPtConversion;
 224   fhPhiConversion = g.fhPhiConversion;
 225   fhEtaConversion = g.fhEtaConversion;
 226   fhPtUnknown = g.fhPtUnknown;
 227   fhPhiUnknown = g.fhPhiUnknown;
 228   fhEtaUnknown = g.fhEtaUnknown;
 229
 230   //Several IC
 231   fNCones = g.fNCones ;
 232   fNPtThresFrac = g.fNPtThresFrac ;
 233
 234   for(Int_t i = 0; i < fNCones ; i++){
 235     fConeSizes[i] =  g.fConeSizes[i];
 236     fhPtSumIsolated[i] = g.fhPtSumIsolated[i] ;
 237
 238    fhPtSumIsolatedPrompt[i] = g.fhPtSumIsolatedPrompt[i];
 239    fhPtSumIsolatedFragmentation[i] = g.fhPtSumIsolatedFragmentation[i];
 240    fhPtSumIsolatedPi0Decay[i] = g.fhPtSumIsolatedPi0Decay[i];
 241    fhPtSumIsolatedOtherDecay[i] = g.fhPtSumIsolatedOtherDecay[i];
 242    fhPtSumIsolatedConversion[i] = g.fhPtSumIsolatedConversion[i];
 243    fhPtSumIsolatedUnknown[i] = g.fhPtSumIsolatedUnknown[i];
 244
 245     for(Int_t j = 0; j < fNPtThresFrac ; j++){
 246       fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j] ;
 247       fhPtFracIsolated[i][j] = g.fhPtFracIsolated[i][j] ;
 248
 249       fhPtThresIsolatedPrompt[i][j] = g.fhPtThresIsolatedPrompt[i][j];
 250       fhPtThresIsolatedFragmentation[i][j] = g.fhPtThresIsolatedFragmentation[i][j];
 251       fhPtThresIsolatedPi0Decay[i][j] = g.fhPtThresIsolatedPi0Decay[i][j];
 252       fhPtThresIsolatedOtherDecay[i][j] = g.fhPtThresIsolatedOtherDecay[i][j];
 253       fhPtThresIsolatedConversion[i][j] = g.fhPtThresIsolatedConversion[i][j];
 254       fhPtThresIsolatedUnknown[i][j] = g.fhPtThresIsolatedUnknown[i][j];
 255
 256       fhPtFracIsolatedPrompt[i][j] = g.fhPtFracIsolatedPrompt[i][j];
 257       fhPtFracIsolatedFragmentation[i][j] = g.fhPtFracIsolatedFragmentation[i][j];
 258       fhPtFracIsolatedPi0Decay[i][j] = g.fhPtFracIsolatedPi0Decay[i][j];
 259       fhPtFracIsolatedOtherDecay[i][j] = g.fhPtFracIsolatedOtherDecay[i][j];
 260       fhPtFracIsolatedConversion[i][j] = g.fhPtFracIsolatedConversion[i][j];
 261       fhPtFracIsolatedUnknown[i][j] = g.fhPtFracIsolatedUnknown[i][j];
 262
 263     }
 264   }
 265
 266   for(Int_t i = 0; i < fNPtThresFrac ; i++){
 267     fPtThresholds[i]=   g.fPtThresholds[i];
 268     fPtFractions[i]=   g.fPtFractions[i];
 269   }
 270
 271   return *this;
 272
 273 }
 274
 275 //____________________________________________________________________________
 276 AliAnaGammaDirect::~AliAnaGammaDirect()
 277 {
 278   //dtor
 279   //do not delete histograms
 280
 281  delete [] fConeSizes ;
 282  delete [] fPtThresholds ;
 283  delete [] fPtFractions ;
 284
 285 }
 286 //_________________________________________________________________________
 287 Bool_t AliAnaGammaDirect::CheckInvMass(const Int_t icalo,const TLorentzVector mom, Double_t *vertex, TClonesArray * pl){
 288   //Search if there is a companion decay photon to the candidate
 289   // and discard it in such case
 290   TLorentzVector mom2 ;
 291   for(Int_t jcalo = 0; jcalo < pl->GetEntries(); jcalo++){
 292     if(icalo==jcalo) continue ;
 293     AliAODCaloCluster * calo =  dynamic_cast<AliAODCaloCluster*> (pl->At(jcalo));
 294
 295     //Cluster selection, not charged, with photon id and in fidutial cut, fill TLorentz
 296     if(!SelectCluster(calo, vertex, mom2)) continue ;
 297
 298     //Select good pair (good phit, pt cuts, aperture and invariant mass)
 299     if(GetNeutralMesonSelection()->SelectPair(mom, mom2)){
 300       if(GetDebug()>1)printf("Selected gamma pair: pt %f, phi %f, eta%f",(mom+mom2).Pt(), (mom+mom2).Phi(), (mom+mom2).Eta());
 301       return kTRUE ;
 302     }
 303   }//loop
 304
 305   return kFALSE;
 306
 307 }
 308 //_________________________________________________________________________
 309 Int_t AliAnaGammaDirect::CheckOrigin(const Int_t label){
 310   //Play with the MC stack if available
 311   //Check origin of the candidates, good for PYTHIA
 312
 313   AliStack * stack =  GetMCStack() ;
 314   if(!stack) AliFatal("Stack is not available, check analysis settings in configuration file, STOP!!");
 315
 316   if(label >= 0 && label <  stack->GetNtrack()){
 317     //Mother
 318     TParticle * mom = stack->Particle(label);
 319     Int_t mPdg = TMath::Abs(mom->GetPdgCode());
 320     Int_t mStatus =  mom->GetStatusCode() ;
 321     Int_t iParent =  mom->GetFirstMother() ;
 322     if(label < 8 ) printf("Mother is parton %d\n",iParent);
 323
 324     //GrandParent
 325     TParticle * parent = new TParticle ;
 326     Int_t pPdg = -1;
 327     Int_t pStatus =-1;
 328     if(iParent > 0){
 329       parent = stack->Particle(iParent);
 330       pPdg = TMath::Abs(parent->GetPdgCode());
 331       pStatus = parent->GetStatusCode();
 332     }
 333     else
 334       printf("Parent with label %d\n",iParent);
 335
 336     //return tag
 337     if(mPdg == 22){
 338       if(mStatus == 1){
 339         if(iParent < 8) {
 340           if(pPdg == 22) return kPrompt;
 341           else  return kFragmentation;
 342         }
 343         else if(pStatus == 11){
 344           if(pPdg == 111) return kPi0Decay ;
 345           else if (pPdg == 321)  return kEtaDecay ;
 346           else  return kOtherDecay ;
 347         }
 348       }//Status 1 : Pythia generated
 349       else if(mStatus == 0){
 350         if(pPdg ==22 || pPdg ==11) return kConversion ;
 351         if(pPdg == 111) return kPi0Decay ;
 352         else if (pPdg == 221)  return kEtaDecay ;
 353         else  return kOtherDecay ;
 354       }//status 0 : geant generated
 355     }//Mother gamma
 356     else if(mPdg == 111)  return kPi0 ;
 357     else if(mPdg == 221)  return kEta ;
 358     else if(mPdg ==11){
 359       if(mStatus == 0) return kConversion ;
 360       else return kElectron ;
 361     }
 362     else return kUnknown;
 363   }//Good label value
 364   else{
 365     if(label < 0 ) printf("*** bad label or no stack ***:  label %d \n", label);
 366     if(label >=  stack->GetNtrack()) printf("*** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 367     return kUnknown;
 368   }//Bad label
 369
 370   return kUnknown;
 371
 372 }
 373
 374 //________________________________________________________________________
 375 TList *  AliAnaGammaDirect::GetCreateOutputObjects()
 376 {
 377   // Create histograms to be saved in output file and
 378   // store them in outputContainer
 379   TList * outputContainer = new TList() ;
 380   outputContainer->SetName("DirectGammaHistos") ;
 381
 382   //Histograms of highest gamma identified in Event
 383   fhPtGamma  = new TH1F("hPtGamma","Number of #gamma over calorimeter",240,0,120);
 384   fhPtGamma->SetYTitle("N");
 385   fhPtGamma->SetXTitle("p_{T #gamma}(GeV/c)");
 386   outputContainer->Add(fhPtGamma) ;
 387
 388   fhPhiGamma  = new TH2F
 389     ("hPhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
 390   fhPhiGamma->SetYTitle("#phi");
 391   fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
 392   outputContainer->Add(fhPhiGamma) ;
 393
 394   fhEtaGamma  = new TH2F
 395     ("hEtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
 396   fhEtaGamma->SetYTitle("#eta");
 397   fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
 398   outputContainer->Add(fhEtaGamma) ;
 399
 400   if(!fMakeSeveralIC){
 401     fhConeSumPt  = new TH2F
 402       ("hConePtSum","#Sigma p_{T}  in cone ",200,0,120,100,0,100);
 403     fhConeSumPt->SetYTitle("#Sigma p_{T}");
 404     fhConeSumPt->SetXTitle("p_{T #gamma} (GeV/c)");
 405     outputContainer->Add(fhConeSumPt) ;
 406   }
 407
 408   if(IsDataMC()){
 409
 410     fhPtPrompt  = new TH1F("hPtPrompt","Number of #gamma over calorimeter",240,0,120);
 411     fhPtPrompt->SetYTitle("N");
 412     fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
 413     outputContainer->Add(fhPtPrompt) ;
 414
 415     fhPhiPrompt  = new TH2F
 416       ("hPhiPrompt","#phi_{#gamma}",200,0,120,200,0,7);
 417     fhPhiPrompt->SetYTitle("#phi");
 418     fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 419     outputContainer->Add(fhPhiPrompt) ;
 420
 421     fhEtaPrompt  = new TH2F
 422       ("hEtaPrompt","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
 423     fhEtaPrompt->SetYTitle("#eta");
 424     fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 425     outputContainer->Add(fhEtaPrompt) ;
 426
 427     fhPtFragmentation  = new TH1F("hPtFragmentation","Number of #gamma over calorimeter",240,0,120);
 428     fhPtFragmentation->SetYTitle("N");
 429     fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
 430     outputContainer->Add(fhPtFragmentation) ;
 431
 432     fhPhiFragmentation  = new TH2F
 433       ("hPhiFragmentation","#phi_{#gamma}",200,0,120,200,0,7);
 434     fhPhiFragmentation->SetYTitle("#phi");
 435     fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 436     outputContainer->Add(fhPhiFragmentation) ;
 437
 438     fhEtaFragmentation  = new TH2F
 439       ("hEtaFragmentation","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
 440     fhEtaFragmentation->SetYTitle("#eta");
 441     fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 442     outputContainer->Add(fhEtaFragmentation) ;
 443
 444     fhPtPi0Decay  = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",240,0,120);
 445     fhPtPi0Decay->SetYTitle("N");
 446     fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
 447     outputContainer->Add(fhPtPi0Decay) ;
 448
 449     fhPhiPi0Decay  = new TH2F
 450       ("hPhiPi0Decay","#phi_{#gamma}",200,0,120,200,0,7);
 451     fhPhiPi0Decay->SetYTitle("#phi");
 452     fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 453     outputContainer->Add(fhPhiPi0Decay) ;
 454
 455     fhEtaPi0Decay  = new TH2F
 456       ("hEtaPi0Decay","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
 457     fhEtaPi0Decay->SetYTitle("#eta");
 458     fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 459     outputContainer->Add(fhEtaPi0Decay) ;
 460
 461     fhPtOtherDecay  = new TH1F("hPtOtherDecay","Number of #gamma over calorimeter",240,0,120);
 462     fhPtOtherDecay->SetYTitle("N");
 463     fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
 464     outputContainer->Add(fhPtOtherDecay) ;
 465
 466     fhPhiOtherDecay  = new TH2F
 467       ("hPhiOtherDecay","#phi_{#gamma}",200,0,120,200,0,7);
 468     fhPhiOtherDecay->SetYTitle("#phi");
 469     fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 470     outputContainer->Add(fhPhiOtherDecay) ;
 471
 472     fhEtaOtherDecay  = new TH2F
 473       ("hEtaOtherDecay","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
 474     fhEtaOtherDecay->SetYTitle("#eta");
 475     fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 476     outputContainer->Add(fhEtaOtherDecay) ;
 477
 478     fhPtConversion  = new TH1F("hPtConversion","Number of #gamma over calorimeter",240,0,120);
 479     fhPtConversion->SetYTitle("N");
 480     fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
 481     outputContainer->Add(fhPtConversion) ;
 482
 483     fhPhiConversion  = new TH2F
 484       ("hPhiConversion","#phi_{#gamma}",200,0,120,200,0,7);
 485     fhPhiConversion->SetYTitle("#phi");
 486     fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 487     outputContainer->Add(fhPhiConversion) ;
 488
 489     fhEtaConversion  = new TH2F
 490       ("hEtaConversion","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
 491     fhEtaConversion->SetYTitle("#eta");
 492     fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 493     outputContainer->Add(fhEtaConversion) ;
 494
 495     fhPtUnknown  = new TH1F("hPtUnknown","Number of #gamma over calorimeter",240,0,120);
 496     fhPtUnknown->SetYTitle("N");
 497     fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
 498     outputContainer->Add(fhPtUnknown) ;
 499
 500     fhPhiUnknown  = new TH2F
 501       ("hPhiUnknown","#phi_{#gamma}",200,0,120,200,0,7);
 502     fhPhiUnknown->SetYTitle("#phi");
 503     fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
 504     outputContainer->Add(fhPhiUnknown) ;
 505
 506     fhEtaUnknown  = new TH2F
 507       ("hEtaUnknown","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
 508     fhEtaUnknown->SetYTitle("#eta");
 509     fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
 510     outputContainer->Add(fhEtaUnknown) ;
 511   }//Histos with MC
 512
 513   if(fMakeSeveralIC){
 514     char name[128];
 515     char title[128];
 516     for(Int_t icone = 0; icone<fNCones; icone++){
 517       sprintf(name,"hPtSumIsolated_Cone_%d",icone);
 518       sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 519       fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
 520       fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 521       fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
 522       outputContainer->Add(fhPtSumIsolated[icone]) ;
 523
 524       if(IsDataMC()){
 525         sprintf(name,"hPtSumIsolatedPrompt_Cone_%d",icone);
 526         sprintf(title,"Candidate Prompt cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 527         fhPtSumIsolatedPrompt[icone]  = new TH2F(name, title,240,0,120,120,0,10);
 528         fhPtSumIsolatedPrompt[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 529         fhPtSumIsolatedPrompt[icone]->SetXTitle("p_{T} (GeV/c)");
 530         outputContainer->Add(fhPtSumIsolatedPrompt[icone]) ;
 531
 532         sprintf(name,"hPtSumIsolatedFragmentation_Cone_%d",icone);
 533         sprintf(title,"Candidate Fragmentation cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 534         fhPtSumIsolatedFragmentation[icone]  = new TH2F(name, title,240,0,120,120,0,10);
 535         fhPtSumIsolatedFragmentation[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 536         fhPtSumIsolatedFragmentation[icone]->SetXTitle("p_{T} (GeV/c)");
 537         outputContainer->Add(fhPtSumIsolatedFragmentation[icone]) ;
 538
 539         sprintf(name,"hPtSumIsolatedPi0Decay_Cone_%d",icone);
 540         sprintf(title,"Candidate Pi0Decay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 541         fhPtSumIsolatedPi0Decay[icone]  = new TH2F(name, title,240,0,120,120,0,10);
 542         fhPtSumIsolatedPi0Decay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 543         fhPtSumIsolatedPi0Decay[icone]->SetXTitle("p_{T} (GeV/c)");
 544         outputContainer->Add(fhPtSumIsolatedPi0Decay[icone]) ;
 545
 546         sprintf(name,"hPtSumIsolatedOtherDecay_Cone_%d",icone);
 547         sprintf(title,"Candidate OtherDecay cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 548         fhPtSumIsolatedOtherDecay[icone]  = new TH2F(name, title,240,0,120,120,0,10);
 549         fhPtSumIsolatedOtherDecay[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 550         fhPtSumIsolatedOtherDecay[icone]->SetXTitle("p_{T} (GeV/c)");
 551         outputContainer->Add(fhPtSumIsolatedOtherDecay[icone]) ;
 552
 553         sprintf(name,"hPtSumIsolatedConversion_Cone_%d",icone);
 554         sprintf(title,"Candidate Conversion cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 555         fhPtSumIsolatedConversion[icone]  = new TH2F(name, title,240,0,120,120,0,10);
 556         fhPtSumIsolatedConversion[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 557         fhPtSumIsolatedConversion[icone]->SetXTitle("p_{T} (GeV/c)");
 558         outputContainer->Add(fhPtSumIsolatedConversion[icone]) ;
 559
 560         sprintf(name,"hPtSumIsolatedUnknown_Cone_%d",icone);
 561         sprintf(title,"Candidate Unknown cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
 562         fhPtSumIsolatedUnknown[icone]  = new TH2F(name, title,240,0,120,120,0,10);
 563         fhPtSumIsolatedUnknown[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
 564         fhPtSumIsolatedUnknown[icone]->SetXTitle("p_{T} (GeV/c)");
 565         outputContainer->Add(fhPtSumIsolatedUnknown[icone]) ;
 566
 567       }//Histos with MC
 568
 569       for(Int_t ipt = 0; ipt<fNPtThresFrac;ipt++){
 570         sprintf(name,"hPtThresIsol_Cone_%d_Pt%d",icone,ipt);
 571         sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 572         fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
 573         fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 574         outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
 575
 576         sprintf(name,"hPtFracIsol_Cone_%d_Pt%d",icone,ipt);
 577         sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 578         fhPtFracIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
 579         fhPtFracIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 580         outputContainer->Add(fhPtFracIsolated[icone][ipt]) ;
 581
 582         if(IsDataMC()){
 583           sprintf(name,"hPtThresIsolPrompt_Cone_%d_Pt%d",icone,ipt);
 584           sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 585           fhPtThresIsolatedPrompt[icone][ipt]  = new TH1F(name, title,240,0,120);
 586           fhPtThresIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 587           outputContainer->Add(fhPtThresIsolatedPrompt[icone][ipt]) ;
 588
 589           sprintf(name,"hPtFracIsolPrompt_Cone_%d_Pt%d",icone,ipt);
 590           sprintf(title,"Isolated candidate Prompt p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 591           fhPtFracIsolatedPrompt[icone][ipt]  = new TH1F(name, title,240,0,120);
 592           fhPtFracIsolatedPrompt[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 593           outputContainer->Add(fhPtFracIsolatedPrompt[icone][ipt]) ;
 594
 595           sprintf(name,"hPtThresIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
 596           sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 597           fhPtThresIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,240,0,120);
 598           fhPtThresIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 599           outputContainer->Add(fhPtThresIsolatedFragmentation[icone][ipt]) ;
 600
 601           sprintf(name,"hPtFracIsolFragmentation_Cone_%d_Pt%d",icone,ipt);
 602           sprintf(title,"Isolated candidate Fragmentation p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 603           fhPtFracIsolatedFragmentation[icone][ipt]  = new TH1F(name, title,240,0,120);
 604           fhPtFracIsolatedFragmentation[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 605           outputContainer->Add(fhPtFracIsolatedFragmentation[icone][ipt]) ;
 606
 607           sprintf(name,"hPtThresIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
 608           sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 609           fhPtThresIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,240,0,120);
 610           fhPtThresIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 611           outputContainer->Add(fhPtThresIsolatedPi0Decay[icone][ipt]) ;
 612
 613           sprintf(name,"hPtFracIsolPi0Decay_Cone_%d_Pt%d",icone,ipt);
 614           sprintf(title,"Isolated candidate Pi0Decay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 615           fhPtFracIsolatedPi0Decay[icone][ipt]  = new TH1F(name, title,240,0,120);
 616           fhPtFracIsolatedPi0Decay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 617           outputContainer->Add(fhPtFracIsolatedPi0Decay[icone][ipt]) ;
 618
 619           sprintf(name,"hPtThresIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
 620           sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 621           fhPtThresIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,240,0,120);
 622           fhPtThresIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 623           outputContainer->Add(fhPtThresIsolatedOtherDecay[icone][ipt]) ;
 624
 625           sprintf(name,"hPtFracIsolOtherDecay_Cone_%d_Pt%d",icone,ipt);
 626           sprintf(title,"Isolated candidate OtherDecay p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 627           fhPtFracIsolatedOtherDecay[icone][ipt]  = new TH1F(name, title,240,0,120);
 628           fhPtFracIsolatedOtherDecay[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 629           outputContainer->Add(fhPtFracIsolatedOtherDecay[icone][ipt]) ;
 630
 631           sprintf(name,"hPtThresIsolConversion_Cone_%d_Pt%d",icone,ipt);
 632           sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 633           fhPtThresIsolatedConversion[icone][ipt]  = new TH1F(name, title,240,0,120);
 634           fhPtThresIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 635           outputContainer->Add(fhPtThresIsolatedConversion[icone][ipt]) ;
 636
 637           sprintf(name,"hPtFracIsolConversion_Cone_%d_Pt%d",icone,ipt);
 638           sprintf(title,"Isolated candidate Conversion p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 639           fhPtFracIsolatedConversion[icone][ipt]  = new TH1F(name, title,240,0,120);
 640           fhPtFracIsolatedConversion[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 641           outputContainer->Add(fhPtFracIsolatedConversion[icone][ipt]) ;
 642
 643           sprintf(name,"hPtThresIsolUnknown_Cone_%d_Pt%d",icone,ipt);
 644           sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 645           fhPtThresIsolatedUnknown[icone][ipt]  = new TH1F(name, title,240,0,120);
 646           fhPtThresIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 647           outputContainer->Add(fhPtThresIsolatedUnknown[icone][ipt]) ;
 648
 649           sprintf(name,"hPtFracIsolUnknown_Cone_%d_Pt%d",icone,ipt);
 650           sprintf(title,"Isolated candidate Unknown p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
 651           fhPtFracIsolatedUnknown[icone][ipt]  = new TH1F(name, title,240,0,120);
 652           fhPtFracIsolatedUnknown[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
 653           outputContainer->Add(fhPtFracIsolatedUnknown[icone][ipt]) ;
 654
 655         }//Histos with MC
 656
 657       }//icone loop
 658     }//ipt loop
 659   }
 660
 661   //Keep neutral meson selection histograms if requiered
 662   //Setting done in AliNeutralMesonSelection
 663   if(fMakeInvMass && GetNeutralMesonSelection()){
 664     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
 665     cout<<"NMSHistos "<< nmsHistos<<endl;
 666     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
 667       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
 668   }
 669
 670   return outputContainer ;
 671
 672 }
 673
 674 //__________________________________________________________________
 675 void  AliAnaGammaDirect::MakeAnalysisFillAOD()
 676 {
 677   //Do analysis and fill aods
 678   //Search for the isolated photon in fDetector with pt > GetMinPt()
 679
 680   //Fill CaloClusters if working with ESDs
 681   //if(GetReader()->GetDataType() == AliCaloTrackReader::kESD) ConnectAODCaloClusters();
 682
 683   Int_t n = 0, nfrac = 0;
 684   Bool_t isolated = kFALSE ;
 685   Float_t coneptsum = 0 ;
 686   TClonesArray * pl = new TClonesArray;
 687
 688   //Get vertex for photon momentum calculation
 689   Double_t vertex[]={0,0,0} ; //vertex ;
 690   if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
 691
 692   //Select the detector of the photon
 693   if(fDetector == "PHOS")
 694     pl = GetAODPHOS();
 695   else if (fDetector == "EMCAL")
 696     pl = GetAODEMCAL();
 697   //cout<<"Number of entries "<<pl->GetEntries()<<endl;
 698
 699   //Fill AODCaloClusters and AODParticleCorrelation with PHOS aods
 700   TLorentzVector mom ;
 701   for(Int_t icalo = 0; icalo < pl->GetEntries(); icalo++){
 702     AliAODCaloCluster * calo =  dynamic_cast<AliAODCaloCluster*> (pl->At(icalo));
 703
 704     //Cluster selection, not charged, with photon id and in fidutial cut
 705     if(!SelectCluster(calo,vertex,mom)) continue ;
 706
 707     //If too small pt, skip
 708     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
 709
 710     //Play with the MC stack if available
 711     Int_t tag =-1;
 712     if(IsDataMC()){
 713       //Check origin of the candidates
 714       tag = CheckOrigin(calo->GetLabel(0));
 715       if(GetDebug() > 0) printf("Origin of candidate %d\n",tag);
 716     }//Work with stack also
 717
 718     //Check invariant mass
 719     Bool_t decay = kFALSE ;
 720     if(fMakeInvMass) decay = CheckInvMass(icalo,mom,vertex,pl);
 721     if(decay) continue ;
 722
 723     //Create AOD for analysis
 724     AliAODParticleCorrelation ph = AliAODParticleCorrelation(mom);
 725     ph.SetLabel(calo->GetLabel(0));
 726     ph.SetPdg(AliCaloPID::kPhoton);
 727     ph.SetDetector(fDetector);
 728     ph.SetTag(tag);
 729     if(fMakeIC){
 730       n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
 731       GetIsolationCut()->MakeIsolationCut((TSeqCollection*)GetAODCTS(), (TSeqCollection*)pl,
 732                                           vertex, kTRUE, &ph,icalo,-1,
 733                                           n,nfrac,coneptsum, isolated);
 734       if(isolated){
 735         //cout<<"Isolated : E "<<mom.E()<<" ; Phi"<<mom.Phi()<< " ; Eta "<<mom.Eta()<<endl;
 736         AddAODParticleCorrelation(ph);
 737       }
 738     }
 739     else //Fill all if not isolation requested
 740       AddAODParticleCorrelation(ph);
 741
 742   }//loop
 743
 744   if(GetDebug() > 1) printf("End fill AODs ");
 745
 746 }
 747
 748 //__________________________________________________________________
 749 void  AliAnaGammaDirect::MakeAnalysisFillHistograms()
 750 {
 751   //Do analysis and fill histograms
 752   Int_t n = 0, nfrac = 0;
 753   Bool_t isolated = kFALSE ;
 754   Float_t coneptsum = 0 ;
 755   //cout<<"MakeAnalysisFillHistograms"<<endl;
 756
 757   //Get vertex for photon momentum calculation
 758   Double_t v[]={0,0,0} ; //vertex ;
 759   if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(v);
 760
 761   //Loop on stored AOD photons
 762   Int_t naod = GetAODBranch()->GetEntriesFast();
 763   if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
 764   for(Int_t iaod = 0; iaod < naod ; iaod++){
 765    AliAODParticleCorrelation* ph =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
 766    Int_t pdg = ph->GetPdg();
 767
 768    //Only isolate photons in detector fDetector
 769    if(pdg != AliCaloPID::kPhoton || ph->GetDetector() != fDetector) continue;
 770
 771    if(fMakeSeveralIC) {
 772      //Analysis of multiple IC at same time
 773      MakeSeveralICAnalysis(ph,v);
 774      continue;
 775    }
 776    else if(fReMakeIC){
 777      //In case a more strict IC is needed in the produced AOD
 778      n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
 779      GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(),
 780                                          (TSeqCollection*)ph->GetRefIsolationConeClusters(),
 781                                          v, kFALSE, ph,-1,-1,
 782                                          n,nfrac,coneptsum, isolated);
 783    }
 784
 785    //Fill histograms (normal case)
 786    if(!fReMakeIC || (fReMakeIC && isolated)){
 787
 788      //printf("Isolated Gamma: pt %f, phi %f, eta %f", ph->Pt(),ph->Phi(),ph->Eta()) ;
 789
 790      //Fill prompt gamma histograms
 791      Float_t ptcluster = ph->Pt();
 792      Float_t phicluster = ph->Phi();
 793      Float_t etacluster = ph->Eta();
 794
 795      fhPtGamma   ->Fill(ptcluster);
 796      fhPhiGamma ->Fill(ptcluster,phicluster);
 797      fhEtaGamma ->Fill(ptcluster,etacluster);
 798      fhConeSumPt->Fill(ptcluster,coneptsum);
 799
 800      if(IsDataMC()){
 801        Int_t tag =ph->GetTag();
 802
 803        if(tag == kPrompt){
 804          fhPtPrompt   ->Fill(ptcluster);
 805          fhPhiPrompt ->Fill(ptcluster,phicluster);
 806          fhEtaPrompt ->Fill(ptcluster,etacluster);
 807        }
 808        else if(tag==kFragmentation)
 809          {
 810            fhPtFragmentation   ->Fill(ptcluster);
 811            fhPhiFragmentation ->Fill(ptcluster,phicluster);
 812            fhEtaFragmentation ->Fill(ptcluster,etacluster);
 813          }
 814        else if(tag==kPi0Decay)
 815          {
 816            fhPtPi0Decay   ->Fill(ptcluster);
 817            fhPhiPi0Decay ->Fill(ptcluster,phicluster);
 818            fhEtaPi0Decay ->Fill(ptcluster,etacluster);
 819          }
 820        else if(tag==kEtaDecay || tag==kOtherDecay)
 821          {
 822            fhPtOtherDecay   ->Fill(ptcluster);
 823            fhPhiOtherDecay ->Fill(ptcluster,phicluster);
 824            fhEtaOtherDecay ->Fill(ptcluster,etacluster);
 825          }
 826        else if(tag==kConversion)
 827          {
 828            fhPtConversion   ->Fill(ptcluster);
 829            fhPhiConversion ->Fill(ptcluster,phicluster);
 830            fhEtaConversion ->Fill(ptcluster,etacluster);
 831          }
 832        else{
 833
 834          fhPtUnknown   ->Fill(ptcluster);
 835          fhPhiUnknown ->Fill(ptcluster,phicluster);
 836          fhEtaUnknown ->Fill(ptcluster,etacluster);
 837        }
 838      }//Histograms with MC
 839
 840    }
 841   }// aod loop
 842
 843 }
 844
 845 //____________________________________________________________________________
 846 void AliAnaGammaDirect::InitParameters()
 847 {
 848
 849   //Initialize the parameters of the analysis.
 850
 851   fDetector = "PHOS" ;
 852   fMakeIC = kTRUE;
 853   fReMakeIC = kFALSE ;
 854   fMakeSeveralIC = kFALSE ;
 855   fMakeInvMass = kFALSE ;
 856
 857  //----------- Several IC-----------------
 858   fNCones           = 5 ;
 859   fNPtThresFrac         = 6 ;
 860   fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4; fConeSizes[4] = 0.5;
 861   fPtThresholds[0]=0.; fPtThresholds[1]=1.; fPtThresholds[2]=2.; fPtThresholds[3]=3.; fPtThresholds[4]=4.;fPtThresholds[5]=5.;
 862   fPtFractions[0]=0.05; fPtFractions[1]=0.075; fPtFractions[2]=0.1; fPtFractions[3]=1.25; fPtFractions[4]=1.5;fPtFractions[5]=2.;
 863
 864 }
 865
 866 //__________________________________________________________________
 867 void  AliAnaGammaDirect::MakeSeveralICAnalysis(AliAODParticleCorrelation* ph, Double_t v[3])
 868 {
 869   //Isolation Cut Analysis for both methods and different pt cuts and cones
 870   Float_t ptC = ph->Pt();
 871   Float_t phiC = ph->Phi();
 872   Float_t etaC = ph->Eta();
 873
 874   fhPtGamma->Fill(ptC);
 875   fhPhiGamma->Fill(ptC,ph->Phi());
 876   fhEtaGamma->Fill(ptC,ph->Eta());
 877   Int_t tag =ph->GetTag();
 878
 879   if(IsDataMC()){
 880     if(tag == kPrompt){
 881       fhPtPrompt   ->Fill(ptC);
 882       fhPhiPrompt ->Fill(ptC,phiC);
 883       fhEtaPrompt ->Fill(ptC,etaC);
 884     }
 885     else if(tag==kFragmentation)
 886       {
 887         fhPtFragmentation   ->Fill(ptC);
 888         fhPhiFragmentation ->Fill(ptC,phiC);
 889         fhEtaFragmentation ->Fill(ptC,etaC);
 890       }
 891     else if(tag==kPi0Decay)
 892       {
 893         fhPtPi0Decay   ->Fill(ptC);
 894         fhPhiPi0Decay ->Fill(ptC,phiC);
 895         fhEtaPi0Decay ->Fill(ptC,etaC);
 896       }
 897     else if(tag==kEtaDecay || tag==kOtherDecay)
 898       {
 899         fhPtOtherDecay   ->Fill(ptC);
 900         fhPhiOtherDecay ->Fill(ptC,phiC);
 901         fhEtaOtherDecay ->Fill(ptC,etaC);
 902       }
 903     else if(tag==kConversion)
 904       {
 905         fhPtConversion   ->Fill(ptC);
 906         fhPhiConversion ->Fill(ptC,phiC);
 907         fhEtaConversion ->Fill(ptC,etaC);
 908       }
 909     else{
 910
 911       fhPtUnknown   ->Fill(ptC);
 912       fhPhiUnknown ->Fill(ptC,phiC);
 913       fhEtaUnknown ->Fill(ptC,etaC);
 914     }
 915   }//Histograms with MC
 916   //Keep original setting used when filling AODs, reset at end of analysis
 917   Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
 918   Float_t ptfracorg = GetIsolationCut()->GetPtFraction();
 919   Float_t rorg = GetIsolationCut()->GetConeSize();
 920
 921   Float_t coneptsum = 0 ;
 922   Int_t n[10][10];//[fNCones][fNPtThresFrac];
 923   Int_t nfrac[10][10];//[fNCones][fNPtThresFrac];
 924   Bool_t  isolated   = kFALSE;
 925
 926   for(Int_t icone = 0; icone<fNCones; icone++){
 927     GetIsolationCut()->SetConeSize(fConeSizes[icone]);
 928     coneptsum = 0 ;
 929     for(Int_t ipt = 0; ipt<fNPtThresFrac ;ipt++){
 930       n[icone][ipt]=0;
 931       nfrac[icone][ipt]=0;
 932       GetIsolationCut()->SetPtThreshold(fPtThresholds[ipt]);
 933       GetIsolationCut()->MakeIsolationCut((TSeqCollection*)ph->GetRefIsolationConeTracks(),
 934                                           (TSeqCollection*)ph->GetRefIsolationConeClusters(),
 935                                           v, kFALSE, ph,-1,-1,
 936                                           n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
 937       if(n[icone][ipt] == 0) {
 938         fhPtThresIsolated[icone][ipt]->Fill(ptC);
 939         if(IsDataMC()){
 940           if(tag==kPrompt) fhPtThresIsolatedPrompt[icone][ipt]->Fill(ptC,coneptsum) ;
 941           else if(tag==kConversion) fhPtThresIsolatedConversion[icone][ipt]->Fill(ptC,coneptsum) ;
 942           else if(tag==kFragmentation) fhPtThresIsolatedFragmentation[icone][ipt]->Fill(ptC,coneptsum) ;
 943           else if(tag==kPi0Decay) fhPtThresIsolatedPi0Decay[icone][ipt]->Fill(ptC,coneptsum) ;
 944           else if(tag==kOtherDecay || tag==kEtaDecay) fhPtThresIsolatedOtherDecay[icone][ipt]->Fill(ptC,coneptsum) ;
 945           else  fhPtThresIsolatedUnknown[icone][ipt]->Fill(ptC,coneptsum) ;
 946         }
 947       }
 948       if(nfrac[icone][ipt] == 0) {
 949         fhPtFracIsolated[icone][ipt]->Fill(ptC);
 950         if(IsDataMC()){
 951           if(tag==kPrompt) fhPtFracIsolatedPrompt[icone][ipt]->Fill(ptC,coneptsum) ;
 952           else if(tag==kConversion) fhPtFracIsolatedConversion[icone][ipt]->Fill(ptC,coneptsum) ;
 953           else if(tag==kFragmentation) fhPtFracIsolatedFragmentation[icone][ipt]->Fill(ptC,coneptsum) ;
 954           else if(tag==kPi0Decay) fhPtFracIsolatedPi0Decay[icone][ipt]->Fill(ptC,coneptsum) ;
 955           else if(tag==kOtherDecay || tag==kEtaDecay) fhPtFracIsolatedOtherDecay[icone][ipt]->Fill(ptC,coneptsum) ;
 956           else  fhPtFracIsolatedUnknown[icone][ipt]->Fill(ptC,coneptsum) ;
 957         }
 958       }
 959     }//pt thresh loop
 960     fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
 961     if(IsDataMC()){
 962       if(tag==kPrompt) fhPtSumIsolatedPrompt[icone]->Fill(ptC,coneptsum) ;
 963       else if(tag==kConversion) fhPtSumIsolatedConversion[icone]->Fill(ptC,coneptsum) ;
 964       else if(tag==kFragmentation) fhPtSumIsolatedFragmentation[icone]->Fill(ptC,coneptsum) ;
 965       else if(tag==kPi0Decay) fhPtSumIsolatedPi0Decay[icone]->Fill(ptC,coneptsum) ;
 966       else if(tag==kOtherDecay || tag==kEtaDecay) fhPtSumIsolatedOtherDecay[icone]->Fill(ptC,coneptsum) ;
 967       else  fhPtSumIsolatedUnknown[icone]->Fill(ptC,coneptsum) ;
 968     }
 969
 970   }//cone size loop
 971
 972   //Reset original parameters for AOD analysis
 973   GetIsolationCut()->SetPtThreshold(ptthresorg);
 974   GetIsolationCut()->SetPtFraction(ptfracorg);
 975   GetIsolationCut()->SetConeSize(rorg);
 976
 977 }
 978
 979 //__________________________________________________________________
 980 void AliAnaGammaDirect::Print(const Option_t * opt) const
 981 {
 982
 983   //Print some relevant parameters set for the analysis
 984   if(! opt)
 985     return;
 986
 987   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 988
 989   printf("Min Gamma pT       =     %2.2f\n",  GetMinPt()) ;
 990   printf("Max Gamma pT       =     %3.2f\n",  GetMaxPt()) ;
 991
 992 //   if(IsCaloPIDOn())printf("Check PID \n") ;
 993 //   if(IsCaloPIDRecalculationOn())printf("Recalculate PID \n") ;
 994 //   if(IsFidutialCutOn())printf("Check Fidutial cut \n") ;
 995   printf("Make Isolation     =     %d \n",  fMakeIC) ;
 996   printf("ReMake Isolation   =     %d \n",  fReMakeIC) ;
 997   printf("Make Several Isolation = %d \n",  fMakeSeveralIC) ;
 998
 999   if(fMakeSeveralIC){
1000     printf("N Cone Sizes       =     %d\n", fNCones) ;
1001     printf("Cone Sizes          =    \n") ;
1002     for(Int_t i = 0; i < fNCones; i++)
1003       printf("  %1.2f;",  fConeSizes[i]) ;
1004     printf("    \n") ;
1005
1006     printf("N pT thresholds/fractions = %d\n", fNPtThresFrac) ;
1007     printf(" pT thresholds         =    \n") ;
1008     for(Int_t i = 0; i < fNPtThresFrac; i++)
1009       printf("   %2.2f;",  fPtThresholds[i]) ;
1010
1011     printf("    \n") ;
1012
1013     printf(" pT fractions         =    \n") ;
1014     for(Int_t i = 0; i < fNPtThresFrac; i++)
1015       printf("   %2.2f;",  fPtFractions[i]) ;
1016
1017   }
1018
1019   printf("    \n") ;
1020
1021 }
1022
1023 //____________________________________________________________________________
1024 Bool_t  AliAnaGammaDirect::SelectCluster(AliAODCaloCluster * calo, Double_t vertex[3], TLorentzVector & mom){
1025    //Select cluster depending on its pid and acceptance selections
1026
1027    //Skip matched clusters with tracks
1028    if(calo->GetNTracksMatched() > 0) return kFALSE ;
1029
1030    //Check PID
1031    calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
1032    Int_t pdg = AliCaloPID::kPhoton;
1033    if(IsCaloPIDOn()){
1034      //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
1035      //or redo PID, recommended option for EMCal.
1036      if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
1037        pdg = GetCaloPID()->GetPdg(fDetector,calo->PID(),mom.E());//PID with weights
1038      else
1039        pdg = GetCaloPID()->GetPdg(fDetector,mom,calo->GetM02(),0,0,0,0);//PID recalculated
1040
1041      if(GetDebug() > 1) printf("PDG of identified particle %d\n",pdg);
1042
1043      //If it does not pass pid, skip
1044      if(pdg != AliCaloPID::kPhoton) return kFALSE ;
1045    }
1046
1047    //Check acceptance selection
1048    if(IsFidutialCutOn()){
1049      Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
1050      if(! in ) return kFALSE ;
1051    }
1052
1053    if(GetDebug() > 1) printf("Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
1054
1055    return kTRUE;
1056
1057  }