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