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