PWG4/PartCorrDep/AliAnaPhoton.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 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: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
  16
  17 //_________________________________________________________________________
  18 //
  19 // Class for the photon identification.
  20 // Clusters from calorimeters are identified as photons
  21 // and kept in the AOD. Few histograms produced.
  22 //
  23 // -- Author: Gustavo Conesa (LNF-INFN)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 #include <TH2F.h>
  29 #include <TClonesArray.h>
  30 #include <TObjString.h>
  31 //#include <Riostream.h>
  32 #include "TParticle.h"
  33
  34 // --- Analysis system ---
  35 #include "AliAnaPhoton.h"
  36 #include "AliCaloTrackReader.h"
  37 #include "AliStack.h"
  38 #include "AliCaloPID.h"
  39 #include "AliMCAnalysisUtils.h"
  40 #include "AliFidutialCut.h"
  41 #include "AliAODCaloCluster.h"
  42
  43 ClassImp(AliAnaPhoton)
  44
  45 //____________________________________________________________________________
  46   AliAnaPhoton::AliAnaPhoton() :
  47     AliAnaPartCorrBaseClass(), fCalorimeter(""),
  48     fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
  49         fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
  50     //MC
  51     fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
  52     fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
  53     fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
  54     fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
  55     fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
  56     fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
  57     fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
  58     fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
  59 {
  60   //default ctor
  61
  62   //Initialize parameters
  63   InitParameters();
  64
  65 }
  66
  67 //____________________________________________________________________________
  68 AliAnaPhoton::AliAnaPhoton(const AliAnaPhoton & g) :
  69   AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),
  70    fMinDist(g.fMinDist),fMinDist2(g.fMinDist2), fMinDist3(g.fMinDist3),
  71    fRejectTrackMatch(g.fRejectTrackMatch),
  72    fhPtPhoton(g.fhPtPhoton),fhPhiPhoton(g.fhPhiPhoton),fhEtaPhoton(g.fhEtaPhoton),
  73   //MC
  74   fhDeltaE(g.fhDeltaE), fhDeltaPt(g.fhDeltaPt),
  75   fhRatioE(g.fhRatioE), fhRatioPt(g.fhRatioPt),
  76   fh2E(g.fh2E), fh2Pt(g.fh2Pt),
  77   fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
  78   fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
  79   fhPtISR(g.fhPtISR),fhPhiISR(g.fhPhiISR),fhEtaISR(g.fhEtaISR),
  80   fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
  81   fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
  82   fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
  83   fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown)
  84
  85 {
  86   // cpy ctor
  87
  88 }
  89
  90 //_________________________________________________________________________
  91 AliAnaPhoton & AliAnaPhoton::operator = (const AliAnaPhoton & g)
  92 {
  93   // assignment operator
  94
  95   if(&g == this) return *this;
  96
  97   fCalorimeter = g.fCalorimeter ;
  98   fMinDist  = g.fMinDist;
  99   fMinDist2 = g.fMinDist2;
 100   fMinDist3 = g.fMinDist3;
 101   fRejectTrackMatch = g.fRejectTrackMatch;
 102
 103   fhPtPhoton  = g.fhPtPhoton ;
 104   fhPhiPhoton = g.fhPhiPhoton ;
 105   fhEtaPhoton = g.fhEtaPhoton ;
 106
 107   fhDeltaE   = g.fhDeltaE;
 108   fhDeltaPt  = g.fhDeltaPt;
 109   fhRatioE   = g.fhRatioE;
 110   fhRatioPt  = g.fhRatioPt;
 111   fh2E   = g.fh2E;
 112   fh2Pt  = g.fh2Pt;
 113
 114   fhPtPrompt  = g.fhPtPrompt;
 115   fhPhiPrompt = g.fhPhiPrompt;
 116   fhEtaPrompt = g.fhEtaPrompt;
 117   fhPtFragmentation  = g.fhPtFragmentation;
 118   fhPhiFragmentation = g.fhPhiFragmentation;
 119   fhEtaFragmentation = g.fhEtaFragmentation;
 120   fhPtISR  = g.fhPtISR;
 121   fhPhiISR = g.fhPhiISR;
 122   fhEtaISR = g.fhEtaISR;
 123   fhPtPi0Decay  = g.fhPtPi0Decay;
 124   fhPhiPi0Decay = g.fhPhiPi0Decay;
 125   fhEtaPi0Decay = g.fhEtaPi0Decay;
 126   fhPtOtherDecay  = g.fhPtOtherDecay;
 127   fhPhiOtherDecay = g.fhPhiOtherDecay;
 128   fhEtaOtherDecay = g.fhEtaOtherDecay;
 129   fhPtConversion  = g. fhPtConversion;
 130   fhPhiConversion = g.fhPhiConversion;
 131   fhEtaConversion = g.fhEtaConversion;
 132   fhPtUnknown  = g.fhPtUnknown;
 133   fhPhiUnknown = g.fhPhiUnknown;
 134   fhEtaUnknown = g.fhEtaUnknown;
 135
 136   return *this;
 137
 138 }
 139
 140 //____________________________________________________________________________
 141 AliAnaPhoton::~AliAnaPhoton()
 142 {
 143   //dtor
 144
 145 }
 146
 147
 148 //________________________________________________________________________
 149 TList *  AliAnaPhoton::GetCreateOutputObjects()
 150 {
 151   // Create histograms to be saved in output file and
 152   // store them in outputContainer
 153   TList * outputContainer = new TList() ;
 154   outputContainer->SetName("PhotonHistos") ;
 155
 156   Int_t nptbins  = GetHistoNPtBins();
 157   Int_t nphibins = GetHistoNPhiBins();
 158   Int_t netabins = GetHistoNEtaBins();
 159   Float_t ptmax  = GetHistoPtMax();
 160   Float_t phimax = GetHistoPhiMax();
 161   Float_t etamax = GetHistoEtaMax();
 162   Float_t ptmin  = GetHistoPtMin();
 163   Float_t phimin = GetHistoPhiMin();
 164   Float_t etamin = GetHistoEtaMin();
 165
 166   //Histograms of highest Photon identified in Event
 167   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 168   fhPtPhoton->SetYTitle("N");
 169   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
 170   outputContainer->Add(fhPtPhoton) ;
 171
 172   fhPhiPhoton  = new TH2F
 173     ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 174   fhPhiPhoton->SetYTitle("#phi");
 175   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 176   outputContainer->Add(fhPhiPhoton) ;
 177
 178   fhEtaPhoton  = new TH2F
 179     ("hEtaPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 180   fhEtaPhoton->SetYTitle("#eta");
 181   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 182   outputContainer->Add(fhEtaPhoton) ;
 183
 184   if(IsDataMC()){
 185     fhDeltaE  = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
 186     fhDeltaE->SetXTitle("#Delta E (GeV)");
 187     outputContainer->Add(fhDeltaE);
 188
 189     fhDeltaPt  = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
 190     fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
 191     outputContainer->Add(fhDeltaPt);
 192
 193     fhRatioE  = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
 194     fhRatioE->SetXTitle("E_{reco}/E_{gen}");
 195     outputContainer->Add(fhRatioE);
 196
 197     fhRatioPt  = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
 198     fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
 199     outputContainer->Add(fhRatioPt);
 200
 201     fh2E  = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 202     fh2E->SetYTitle("E_{rec} (GeV)");
 203     fh2E->SetXTitle("E_{gen} (GeV)");
 204     outputContainer->Add(fh2E);
 205
 206     fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 207     fh2Pt->SetYTitle("p_{T,rec} (GeV/c)");
 208     fh2Pt->SetXTitle("p_{T,gen} (GeV/c)");
 209     outputContainer->Add(fh2Pt);
 210
 211     fhPtPrompt  = new TH1F("hPtMCPrompt","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 212     fhPtPrompt->SetYTitle("N");
 213     fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
 214     outputContainer->Add(fhPtPrompt) ;
 215
 216     fhPhiPrompt  = new TH2F
 217       ("hPhiMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 218     fhPhiPrompt->SetYTitle("#phi");
 219     fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 220     outputContainer->Add(fhPhiPrompt) ;
 221
 222     fhEtaPrompt  = new TH2F
 223       ("hEtaMCPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 224     fhEtaPrompt->SetYTitle("#eta");
 225     fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 226     outputContainer->Add(fhEtaPrompt) ;
 227
 228     fhPtFragmentation  = new TH1F("hPtMCFragmentation","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 229     fhPtFragmentation->SetYTitle("N");
 230     fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
 231     outputContainer->Add(fhPtFragmentation) ;
 232
 233     fhPhiFragmentation  = new TH2F
 234       ("hPhiMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 235     fhPhiFragmentation->SetYTitle("#phi");
 236     fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 237     outputContainer->Add(fhPhiFragmentation) ;
 238
 239     fhEtaFragmentation  = new TH2F
 240       ("hEtaMCFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 241     fhEtaFragmentation->SetYTitle("#eta");
 242     fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 243     outputContainer->Add(fhEtaFragmentation) ;
 244
 245     fhPtISR  = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
 246     fhPtISR->SetYTitle("N");
 247     fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
 248     outputContainer->Add(fhPtISR) ;
 249
 250     fhPhiISR  = new TH2F
 251       ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 252     fhPhiISR->SetYTitle("#phi");
 253     fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
 254     outputContainer->Add(fhPhiISR) ;
 255
 256     fhEtaISR  = new TH2F
 257       ("hEtaMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 258     fhEtaISR->SetYTitle("#eta");
 259     fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
 260     outputContainer->Add(fhEtaISR) ;
 261
 262     fhPtPi0Decay  = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 263     fhPtPi0Decay->SetYTitle("N");
 264     fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
 265     outputContainer->Add(fhPtPi0Decay) ;
 266
 267     fhPhiPi0Decay  = new TH2F
 268       ("hPhiMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 269     fhPhiPi0Decay->SetYTitle("#phi");
 270     fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 271     outputContainer->Add(fhPhiPi0Decay) ;
 272
 273     fhEtaPi0Decay  = new TH2F
 274       ("hEtaMCPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 275     fhEtaPi0Decay->SetYTitle("#eta");
 276     fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 277     outputContainer->Add(fhEtaPi0Decay) ;
 278
 279     fhPtOtherDecay  = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 280     fhPtOtherDecay->SetYTitle("N");
 281     fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
 282     outputContainer->Add(fhPtOtherDecay) ;
 283
 284     fhPhiOtherDecay  = new TH2F
 285       ("hPhiMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 286     fhPhiOtherDecay->SetYTitle("#phi");
 287     fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 288     outputContainer->Add(fhPhiOtherDecay) ;
 289
 290     fhEtaOtherDecay  = new TH2F
 291       ("hEtaMCOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 292     fhEtaOtherDecay->SetYTitle("#eta");
 293     fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 294     outputContainer->Add(fhEtaOtherDecay) ;
 295
 296     fhPtConversion  = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 297     fhPtConversion->SetYTitle("N");
 298     fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
 299     outputContainer->Add(fhPtConversion) ;
 300
 301     fhPhiConversion  = new TH2F
 302       ("hPhiMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 303     fhPhiConversion->SetYTitle("#phi");
 304     fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 305     outputContainer->Add(fhPhiConversion) ;
 306
 307     fhEtaConversion  = new TH2F
 308       ("hEtaMCConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 309     fhEtaConversion->SetYTitle("#eta");
 310     fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 311     outputContainer->Add(fhEtaConversion) ;
 312
 313     fhPtUnknown  = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 314     fhPtUnknown->SetYTitle("N");
 315     fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
 316     outputContainer->Add(fhPtUnknown) ;
 317
 318     fhPhiUnknown  = new TH2F
 319       ("hPhiMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 320     fhPhiUnknown->SetYTitle("#phi");
 321     fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
 322     outputContainer->Add(fhPhiUnknown) ;
 323
 324     fhEtaUnknown  = new TH2F
 325       ("hEtaMCUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 326     fhEtaUnknown->SetYTitle("#eta");
 327     fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
 328     outputContainer->Add(fhEtaUnknown) ;
 329
 330   }//Histos with MC
 331
 332   //Save parameters used for analysis
 333   TString parList ; //this will be list of parameters used for this analysis.
 334   char onePar[255] ;
 335
 336   sprintf(onePar,"--- AliAnaPhoton ---\n") ;
 337   parList+=onePar ;
 338   sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 339   parList+=onePar ;
 340   sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 341   parList+=onePar ;
 342   sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 343   parList+=onePar ;
 344   sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 345   parList+=onePar ;
 346   sprintf(onePar,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
 347   parList+=onePar ;
 348
 349   //Get parameters set in base class.
 350   parList += GetBaseParametersList() ;
 351
 352   //Get parameters set in PID class.
 353   parList += GetCaloPID()->GetPIDParametersList() ;
 354
 355   //Get parameters set in FidutialCut class (not available yet)
 356   //parlist += GetFidCut()->GetFidCutParametersList()
 357
 358   TObjString *oString= new TObjString(parList) ;
 359   outputContainer->Add(oString);
 360
 361   return outputContainer ;
 362
 363 }
 364
 365 //____________________________________________________________________________
 366 void AliAnaPhoton::Init()
 367 {
 368
 369   //Init
 370   //Do some checks
 371   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn()){
 372     printf("AliAnaPhoton::Init() - !!ABORT: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 373     abort();
 374   }
 375   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
 376     printf("AliAnaPhoton::Init() - !!ABORT: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 377     abort();
 378   }
 379
 380 }
 381
 382
 383 //____________________________________________________________________________
 384 void AliAnaPhoton::InitParameters()
 385 {
 386
 387   //Initialize the parameters of the analysis.
 388   SetOutputAODClassName("AliAODPWG4Particle");
 389   SetOutputAODName("PWG4Particle");
 390
 391   AddToHistogramsName("AnaPhoton_");
 392
 393   fCalorimeter = "PHOS" ;
 394   fMinDist  = 2.;
 395   fMinDist2 = 4.;
 396   fMinDist3 = 5.;
 397   fRejectTrackMatch = kTRUE ;
 398 }
 399
 400 //__________________________________________________________________
 401 void  AliAnaPhoton::MakeAnalysisFillAOD()
 402 {
 403   //Do analysis and fill aods
 404   //Search for photons in fCalorimeter
 405   TRefArray * pl = new TRefArray();
 406
 407   //Get vertex for photon momentum calculation
 408   Double_t vertex[]={0,0,0} ; //vertex ;
 409   if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
 410   //Select the Calorimeter of the photon
 411   if(fCalorimeter == "PHOS")
 412     pl = GetAODPHOS();
 413   else if (fCalorimeter == "EMCAL")
 414     pl = GetAODEMCAL();
 415
 416   //Fill AODCaloClusters and AODParticle with PHOS aods
 417   TLorentzVector mom ;
 418
 419   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
 420     AliAODCaloCluster * calo =  (AliAODCaloCluster*) (pl->At(icalo));
 421     //Cluster selection, not charged, with photon id and in fidutial cut
 422     //Get Momentum vector,
 423     calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
 424     //If too small or big pt, skip it
 425     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
 426     //Check acceptance selection
 427     if(IsFidutialCutOn()){
 428       Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;
 429       if(! in ) continue ;
 430     }
 431
 432     //Create AOD for analysis
 433     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
 434     aodph.SetLabel(calo->GetLabel(0));
 435     //printf("Index %d, Id %d\n",icalo, calo->GetID());
 436     //Set the indeces of the original caloclusters
 437     aodph.SetCaloLabel(calo->GetID(),-1);
 438     aodph.SetDetector(fCalorimeter);
 439     if(GetDebug() > 1)
 440       printf("AliAnaPhoton::MakeAnalysisFillAOD() - Min pt cut and fidutial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());
 441
 442     //Check Distance to Bad channel, set bit.
 443     Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
 444     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 445     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
 446       continue ;
 447
 448     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
 449
 450     if(distBad > fMinDist3) aodph.SetDistToBad(2) ;
 451     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
 452     else aodph.SetDistToBad(0) ;
 453
 454     //Check PID
 455     //PID selection or bit setting
 456     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
 457       //Get most probable PID, check PID weights (in MC this option is mandatory)
 458       aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
 459       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
 460       //If primary is not photon, skip it.
 461       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
 462     }
 463     else if(IsCaloPIDOn()){
 464       //Skip matched clusters with tracks
 465       if(fRejectTrackMatch && calo->GetNTracksMatched() > 0) continue ;
 466
 467       //Get most probable PID, 2 options check PID weights
 468       //or redo PID, recommended option for EMCal.
 469       if(!IsCaloPIDRecalculationOn())
 470         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
 471       else
 472         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
 473
 474       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
 475
 476       //If cluster does not pass pid, not photon, skip it.
 477       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
 478
 479     }
 480     else{
 481       //Set PID bits for later selection (AliAnaPi0 for example)
 482       //GetPDG already called in SetPIDBits.
 483       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
 484       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
 485     }
 486
 487     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
 488
 489     //Play with the MC stack if available
 490     //Check origin of the candidates
 491     if(IsDataMC()){
 492       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
 493       if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate %d\n",aodph.GetTag());
 494     }//Work with stack also
 495
 496     //Add AOD with photon object to aod branch
 497     AddAODParticle(aodph);
 498
 499   }//loop
 500
 501   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs \n");
 502
 503 }
 504
 505 //__________________________________________________________________
 506 void  AliAnaPhoton::MakeAnalysisFillHistograms()
 507 {
 508     //Do analysis and fill histograms
 509
 510     // Access MC information in stack if requested, check that it exists.
 511     AliStack * stack = 0x0;
 512     TParticle * primary = 0x0;
 513     if(IsDataMC()){
 514           stack =  GetMCStack() ;
 515           if(!stack) {
 516                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
 517                 abort();
 518           }
 519     }
 520
 521         //Loop on stored AOD photons
 522         Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 523         if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 524
 525         for(Int_t iaod = 0; iaod < naod ; iaod++){
 526     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 527     Int_t pdg = ph->GetPdg();
 528
 529     if(GetDebug() > 3)
 530       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
 531
 532     //If PID used, fill histos with photons in Calorimeter fCalorimeter
 533     if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
 534     if(ph->GetDetector() != fCalorimeter) continue;
 535
 536     if(GetDebug() > 1)
 537       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
 538
 539     //Fill photon histograms
 540     Float_t ptcluster  = ph->Pt();
 541     Float_t phicluster = ph->Phi();
 542     Float_t etacluster = ph->Eta();
 543     Float_t ecluster   = ph->E();
 544
 545     fhPtPhoton  ->Fill(ptcluster);
 546     fhPhiPhoton ->Fill(ptcluster,phicluster);
 547     fhEtaPhoton ->Fill(ptcluster,etacluster);
 548
 549   //Play with the MC stack if available
 550     if(IsDataMC()){
 551
 552       Int_t tag =ph->GetTag();
 553       if(tag == AliMCAnalysisUtils::kMCPrompt){
 554           fhPtPrompt  ->Fill(ptcluster);
 555           fhPhiPrompt ->Fill(ptcluster,phicluster);
 556           fhEtaPrompt ->Fill(ptcluster,etacluster);
 557       }
 558       else if(tag==AliMCAnalysisUtils::kMCFragmentation)
 559         {
 560           fhPtFragmentation  ->Fill(ptcluster);
 561           fhPhiFragmentation ->Fill(ptcluster,phicluster);
 562           fhEtaFragmentation ->Fill(ptcluster,etacluster);
 563         }
 564       else if(tag==AliMCAnalysisUtils::kMCISR)
 565         {
 566           fhPtISR  ->Fill(ptcluster);
 567           fhPhiISR ->Fill(ptcluster,phicluster);
 568           fhEtaISR ->Fill(ptcluster,etacluster);
 569         }
 570       else if(tag==AliMCAnalysisUtils::kMCPi0Decay)
 571         {
 572           fhPtPi0Decay  ->Fill(ptcluster);
 573           fhPhiPi0Decay ->Fill(ptcluster,phicluster);
 574           fhEtaPi0Decay ->Fill(ptcluster,etacluster);
 575         }
 576       else if(tag==AliMCAnalysisUtils::kMCEtaDecay || tag==AliMCAnalysisUtils::kMCOtherDecay)
 577         {
 578           fhPtOtherDecay  ->Fill(ptcluster);
 579           fhPhiOtherDecay ->Fill(ptcluster,phicluster);
 580           fhEtaOtherDecay ->Fill(ptcluster,etacluster);
 581         }
 582       else if(tag==AliMCAnalysisUtils::kMCConversion)
 583         {
 584           fhPtConversion  ->Fill(ptcluster);
 585           fhPhiConversion ->Fill(ptcluster,phicluster);
 586           fhEtaConversion ->Fill(ptcluster,etacluster);
 587
 588         }
 589       else{
 590         fhPtUnknown  ->Fill(ptcluster);
 591         fhPhiUnknown ->Fill(ptcluster,phicluster);
 592         fhEtaUnknown ->Fill(ptcluster,etacluster);
 593       }
 594
 595         Int_t label =ph->GetLabel();
 596         if(label < 0) {
 597                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
 598                 continue;
 599         }
 600
 601         if(label >=  stack->GetNtrack()) {
 602                 if(GetDebug() > 2)  printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 603                 continue ;
 604         }
 605
 606         primary = stack->Particle(label);
 607         if(!primary){
 608                 printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
 609                 continue;
 610         }
 611
 612         fh2E->Fill(primary->Energy(),ecluster);
 613         fh2Pt->Fill(primary->Pt(), ptcluster);
 614         fhDeltaE->Fill(primary->Energy()-ecluster);
 615         fhDeltaPt->Fill(primary->Pt()-ptcluster);
 616         if(primary->Energy()!=0) fhRatioE  ->Fill(ecluster/primary->Energy());
 617         if(primary->Pt()!=0) fhRatioPt ->Fill(ptcluster/primary->Pt());
 618
 619     }//Histograms with MC
 620
 621   }// aod loop
 622
 623 }
 624
 625
 626 //__________________________________________________________________
 627 void AliAnaPhoton::Print(const Option_t * opt) const
 628 {
 629   //Print some relevant parameters set for the analysis
 630
 631   if(! opt)
 632     return;
 633
 634   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 635   AliAnaPartCorrBaseClass::Print(" ");
 636
 637   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
 638   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
 639   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
 640   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
 641   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
 642
 643   printf("    \n") ;
 644
 645 }