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, providGetMixedEvent()ed 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 "AliFiducialCut.h"
  41 #include "AliVCluster.h"
  42 #include "AliAODMCParticle.h"
  43 #include "AliMixedEvent.h"
  44
  45
  46 ClassImp(AliAnaPhoton)
  47
  48 //____________________________________________________________________________
  49   AliAnaPhoton::AliAnaPhoton() :
  50     AliAnaPartCorrBaseClass(), fCalorimeter(""),
  51     fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
  52     fCheckConversion(kFALSE),fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
  53     fTimeCutMin(-1), fTimeCutMax(9999999), fNCellsCut(0),
  54     fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
  55     //MC
  56     fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
  57     fhPtMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),
  58     fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
  59     fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
  60     fhPtISR(0),fhPhiISR(0),fhEtaISR(0),
  61     fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
  62     fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
  63     fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
  64     fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
  65 {
  66   //default ctor
  67
  68   //Initialize parameters
  69   InitParameters();
  70
  71 }//____________________________________________________________________________
  72 AliAnaPhoton::~AliAnaPhoton()
  73 {
  74   //dtor
  75
  76 }
  77
  78 //________________________________________________________________________
  79 TObjString *  AliAnaPhoton::GetAnalysisCuts()
  80 {
  81   //Save parameters used for analysis
  82   TString parList ; //this will be list of parameters used for this analysis.
  83   const Int_t buffersize = 255;
  84   char onePar[buffersize] ;
  85
  86   snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
  87   parList+=onePar ;
  88   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
  89   parList+=onePar ;
  90   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
  91   parList+=onePar ;
  92   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
  93   parList+=onePar ;
  94   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
  95   parList+=onePar ;
  96   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
  97   parList+=onePar ;
  98
  99   //Get parameters set in base class.
 100   parList += GetBaseParametersList() ;
 101
 102   //Get parameters set in PID class.
 103   parList += GetCaloPID()->GetPIDParametersList() ;
 104
 105   //Get parameters set in FiducialCut class (not available yet)
 106   //parlist += GetFidCut()->GetFidCutParametersList()
 107
 108   return new TObjString(parList) ;
 109 }
 110
 111
 112 //________________________________________________________________________
 113 TList *  AliAnaPhoton::GetCreateOutputObjects()
 114 {
 115   // Create histograms to be saved in output file and
 116   // store them in outputContainer
 117   TList * outputContainer = new TList() ;
 118   outputContainer->SetName("PhotonHistos") ;
 119
 120   Int_t nptbins  = GetHistoPtBins();
 121   Int_t nphibins = GetHistoPhiBins();
 122   Int_t netabins = GetHistoEtaBins();
 123   Float_t ptmax  = GetHistoPtMax();
 124   Float_t phimax = GetHistoPhiMax();
 125   Float_t etamax = GetHistoEtaMax();
 126   Float_t ptmin  = GetHistoPtMin();
 127   Float_t phimin = GetHistoPhiMin();
 128   Float_t etamin = GetHistoEtaMin();
 129
 130   //Histograms of highest Photon identified in Event
 131   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 132   fhPtPhoton->SetYTitle("N");
 133   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
 134   outputContainer->Add(fhPtPhoton) ;
 135
 136   fhPhiPhoton  = new TH2F
 137     ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 138   fhPhiPhoton->SetYTitle("#phi");
 139   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 140   outputContainer->Add(fhPhiPhoton) ;
 141
 142   fhEtaPhoton  = new TH2F
 143     ("hEtaPhoton","#eta_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 144   fhEtaPhoton->SetYTitle("#eta");
 145   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 146   outputContainer->Add(fhEtaPhoton) ;
 147
 148   if(IsDataMC()){
 149     fhDeltaE  = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
 150     fhDeltaE->SetXTitle("#Delta E (GeV)");
 151     outputContainer->Add(fhDeltaE);
 152
 153     fhDeltaPt  = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
 154     fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
 155     outputContainer->Add(fhDeltaPt);
 156
 157     fhRatioE  = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
 158     fhRatioE->SetXTitle("E_{reco}/E_{gen}");
 159     outputContainer->Add(fhRatioE);
 160
 161     fhRatioPt  = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
 162     fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
 163     outputContainer->Add(fhRatioPt);
 164
 165     fh2E  = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 166     fh2E->SetXTitle("E_{rec} (GeV)");
 167     fh2E->SetYTitle("E_{gen} (GeV)");
 168     outputContainer->Add(fh2E);
 169
 170     fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 171     fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
 172     fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
 173     outputContainer->Add(fh2Pt);
 174
 175     fhPtMCPhoton  = new TH1F("hPtMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 176     fhPtMCPhoton->SetYTitle("N");
 177     fhPtMCPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
 178     outputContainer->Add(fhPtMCPhoton) ;
 179
 180     fhPhiMCPhoton  = new TH2F
 181       ("hPhiMCPhoton","#phi_{#gamma}, #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 182     fhPhiMCPhoton->SetYTitle("#phi");
 183     fhPhiMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 184     outputContainer->Add(fhPhiMCPhoton) ;
 185
 186     fhEtaMCPhoton  = new TH2F
 187       ("hEtaMCPhoton","#eta_{#gamma}, #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 188     fhEtaMCPhoton->SetYTitle("#eta");
 189     fhEtaMCPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 190     outputContainer->Add(fhEtaMCPhoton) ;
 191
 192     fhPtPrompt  = new TH1F("hPtMCPrompt","Number of prompt #gamma over calorimeter",nptbins,ptmin,ptmax);
 193     fhPtPrompt->SetYTitle("N");
 194     fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
 195     outputContainer->Add(fhPtPrompt) ;
 196
 197     fhPhiPrompt  = new TH2F
 198       ("hPhiMCPrompt","#phi_{#gamma}, prompt #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 199     fhPhiPrompt->SetYTitle("#phi");
 200     fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 201     outputContainer->Add(fhPhiPrompt) ;
 202
 203     fhEtaPrompt  = new TH2F
 204       ("hEtaMCPrompt","#eta_{#gamma}, prompt #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 205     fhEtaPrompt->SetYTitle("#eta");
 206     fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
 207     outputContainer->Add(fhEtaPrompt) ;
 208
 209     fhPtFragmentation  = new TH1F("hPtMCFragmentation","Number of fragmentation #gamma over calorimeter",nptbins,ptmin,ptmax);
 210     fhPtFragmentation->SetYTitle("N");
 211     fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
 212     outputContainer->Add(fhPtFragmentation) ;
 213
 214     fhPhiFragmentation  = new TH2F
 215       ("hPhiMCFragmentation","#phi_{#gamma}, fragmentation #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 216     fhPhiFragmentation->SetYTitle("#phi");
 217     fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 218     outputContainer->Add(fhPhiFragmentation) ;
 219
 220     fhEtaFragmentation  = new TH2F
 221       ("hEtaMCFragmentation","#eta_{#gamma}, fragmentation #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 222     fhEtaFragmentation->SetYTitle("#eta");
 223     fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
 224     outputContainer->Add(fhEtaFragmentation) ;
 225
 226     fhPtISR  = new TH1F("hPtMCISR","Number of initial state radiation #gamma over calorimeter",nptbins,ptmin,ptmax);
 227     fhPtISR->SetYTitle("N");
 228     fhPtISR->SetXTitle("p_{T #gamma}(GeV/c)");
 229     outputContainer->Add(fhPtISR) ;
 230
 231     fhPhiISR  = new TH2F
 232       ("hPhiMCISR","#phi_{#gamma} initial state radiation",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 233     fhPhiISR->SetYTitle("#phi");
 234     fhPhiISR->SetXTitle("p_{T #gamma} (GeV/c)");
 235     outputContainer->Add(fhPhiISR) ;
 236
 237     fhEtaISR  = new TH2F
 238       ("hEtaMCISR","#eta_{#gamma} initial state radiation",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 239     fhEtaISR->SetYTitle("#eta");
 240     fhEtaISR->SetXTitle("p_{T #gamma} (GeV/c)");
 241     outputContainer->Add(fhEtaISR) ;
 242
 243     fhPtPi0Decay  = new TH1F("hPtMCPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 244     fhPtPi0Decay->SetYTitle("N");
 245     fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
 246     outputContainer->Add(fhPtPi0Decay) ;
 247
 248     fhPhiPi0Decay  = new TH2F
 249       ("hPhiMCPi0Decay","#phi_{#gamma}, #pi^{0} decay #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 250     fhPhiPi0Decay->SetYTitle("#phi");
 251     fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 252     outputContainer->Add(fhPhiPi0Decay) ;
 253
 254     fhEtaPi0Decay  = new TH2F
 255       ("hEtaMCPi0Decay","#eta_{#gamma}, #pi^{0} #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 256     fhEtaPi0Decay->SetYTitle("#eta");
 257     fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
 258     outputContainer->Add(fhEtaPi0Decay) ;
 259
 260     fhPtOtherDecay  = new TH1F("hPtMCOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 261     fhPtOtherDecay->SetYTitle("N");
 262     fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
 263     outputContainer->Add(fhPtOtherDecay) ;
 264
 265     fhPhiOtherDecay  = new TH2F
 266       ("hPhiMCOtherDecay","#phi_{#gamma}, other decay #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 267     fhPhiOtherDecay->SetYTitle("#phi");
 268     fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 269     outputContainer->Add(fhPhiOtherDecay) ;
 270
 271     fhEtaOtherDecay  = new TH2F
 272       ("hEtaMCOtherDecay","#eta_{#gamma}, other decay #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 273     fhEtaOtherDecay->SetYTitle("#eta");
 274     fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
 275     outputContainer->Add(fhEtaOtherDecay) ;
 276
 277     fhPtConversion  = new TH1F("hPtMCConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 278     fhPtConversion->SetYTitle("N");
 279     fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
 280     outputContainer->Add(fhPtConversion) ;
 281
 282     fhPhiConversion  = new TH2F
 283       ("hPhiMCConversion","#phi_{#gamma}, conversion #gamma in MC",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 284     fhPhiConversion->SetYTitle("#phi");
 285     fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 286     outputContainer->Add(fhPhiConversion) ;
 287
 288     fhEtaConversion  = new TH2F
 289       ("hEtaMCConversion","#eta_{#gamma}, conversion #gamma in MC",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 290     fhEtaConversion->SetYTitle("#eta");
 291     fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
 292     outputContainer->Add(fhEtaConversion) ;
 293
 294     fhPtUnknown  = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 295     fhPtUnknown->SetYTitle("N");
 296     fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
 297     outputContainer->Add(fhPtUnknown) ;
 298
 299     fhPhiUnknown  = new TH2F
 300       ("hPhiMCUnknown","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 301     fhPhiUnknown->SetYTitle("#phi");
 302     fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
 303     outputContainer->Add(fhPhiUnknown) ;
 304
 305     fhEtaUnknown  = new TH2F
 306       ("hEtaMCUnknown","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 307     fhEtaUnknown->SetYTitle("#eta");
 308     fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
 309     outputContainer->Add(fhEtaUnknown) ;
 310
 311   }//Histos with MC
 312
 313   return outputContainer ;
 314
 315 }
 316
 317 //____________________________________________________________________________
 318 void AliAnaPhoton::Init()
 319 {
 320
 321   //Init
 322   //Do some checks
 323   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
 324     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 325     abort();
 326   }
 327   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
 328     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 329     abort();
 330   }
 331
 332 }
 333
 334
 335 //____________________________________________________________________________
 336 void AliAnaPhoton::InitParameters()
 337 {
 338
 339   //Initialize the parameters of the analysis.
 340   AddToHistogramsName("AnaPhoton_");
 341
 342   fCalorimeter = "PHOS" ;
 343   fMinDist  = 2.;
 344   fMinDist2 = 4.;
 345   fMinDist3 = 5.;
 346   fMassCut  = 0.03; //30 MeV
 347
 348   fTimeCutMin  = -1;
 349   fTimeCutMax  = 9999999;
 350   fNCellsCut = 0;
 351
 352   fRejectTrackMatch       = kTRUE ;
 353   fCheckConversion        = kFALSE;
 354   fAddConvertedPairsToAOD = kFALSE;
 355
 356 }
 357
 358 //__________________________________________________________________
 359 void  AliAnaPhoton::MakeAnalysisFillAOD()
 360 {
 361   //Do photon analysis and fill aods
 362
 363   //  Double_t vertex2[] = {0,0,0} ; //vertex from second input aod
 364   //Get the vertex and check it is not too large in z, cut for SE
 365   Double_t v[3] = {0,0,0}; //vertex ;
 366   GetReader()->GetVertex(v);
 367   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
 368
 369   //Select the Calorimeter of the photon
 370   TObjArray * pl = 0x0;
 371   if(fCalorimeter == "PHOS")
 372     pl = GetAODPHOS();
 373   else if (fCalorimeter == "EMCAL")
 374     pl = GetAODEMCAL();
 375
 376   if(!pl) {
 377     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
 378     return;
 379   }
 380
 381   //Fill AODParticle with PHOS/EMCAL aods
 382   TLorentzVector mom, mom2 ;
 383   Int_t nCaloClusters = pl->GetEntriesFast();
 384   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
 385   Bool_t * indexConverted = new Bool_t[nCaloClusters];
 386   for (Int_t i = 0; i < nCaloClusters; i++)
 387     indexConverted[i] = kFALSE;
 388
 389   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
 390
 391           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
 392     //printf("calo %d, %f\n",icalo,calo->E());
 393
 394     //Get the index where the cluster comes, to retrieve the corresponding vertex
 395     Int_t evtIndex = 0 ;
 396     if (GetMixedEvent()) {
 397       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 398       //Get the vertex and check it is not too large in z
 399       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue ;
 400     }
 401
 402     //Cluster selection, not charged, with photon id and in fiducial cut
 403
 404     //Input from second AOD?
 405     //Int_t input = 0;
 406     //    if (fCalorimeter == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= icalo)
 407     //      input = 1 ;
 408     //    else if(fCalorimeter == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= icalo)
 409     //      input = 1;
 410
 411     //Get Momentum vector,
 412     //if (input == 0)
 413     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 414       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
 415     else{
 416       Double_t vertex[]={0,0,0};
 417       calo->GetMomentum(mom,vertex) ;
 418     }
 419     //printf("AliAnaPhoton::MakeAnalysisFillAOD(): Vertex : %f,%f,%f\n",GetVertex(evtIndex)[0] ,GetVertex(evtIndex)[1],GetVertex(evtIndex)[2]);
 420
 421     //    else if(input == 1)
 422     //      calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
 423
 424     //If too small or big pt, skip it
 425     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
 426     Double_t tof = calo->GetTOF()*1e9;
 427
 428     if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
 429
 430     if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) continue;
 431
 432     //printf("AliAnaPhoton::Current Event %d; Current File Name : %s, E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",GetReader()->GetEventNumber(),(GetReader()->GetCurrentFileName()).Data(),
 433       //     mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 434
 435     //Check acceptance selection
 436     if(IsFiducialCutOn()){
 437       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 438       if(! in ) continue ;
 439     }
 440     //printf("Fiducial cut passed \n");
 441
 442     //Create AOD for analysis
 443     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
 444     Int_t label = calo->GetLabel();
 445     aodph.SetLabel(label);
 446     //aodph.SetInputFileIndex(input);
 447
 448     //printf("Index %d, Id %d\n",icalo, calo->GetID());
 449     //Set the indeces of the original caloclusters
 450     aodph.SetCaloLabel(calo->GetID(),-1);
 451     aodph.SetDetector(fCalorimeter);
 452     if(GetDebug() > 1)
 453       printf("AliAnaPhoton::MakeAnalysisFillAOD() - Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());
 454
 455     //Check Distance to Bad channel, set bit.
 456     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 457     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 458     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
 459       continue ;
 460
 461     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
 462
 463     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
 464     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
 465     else                         aodph.SetDistToBad(0) ;
 466
 467     //Skip matched clusters with tracks
 468     if(fRejectTrackMatch && IsTrackMatched(calo)) continue ;
 469     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - TrackMatching cut passed \n");
 470
 471     //Check PID
 472     //PID selection or bit setting
 473     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
 474       //Get most probable PID, check PID weights (in MC this option is mandatory)
 475       aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 476       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
 477       //If primary is not photon, skip it.
 478       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
 479     }
 480     else if(IsCaloPIDOn()){
 481
 482       //Get most probable PID, 2 options check PID weights
 483       //or redo PID, recommended option for EMCal.
 484       if(!IsCaloPIDRecalculationOn())
 485         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 486       else
 487         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
 488
 489       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
 490
 491       //If cluster does not pass pid, not photon, skip it.
 492       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
 493
 494     }
 495     else{
 496       //Set PID bits for later selection (AliAnaPi0 for example)
 497       //GetPDG already called in SetPIDBits.
 498       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
 499       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
 500     }
 501
 502     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
 503
 504     //Play with the MC stack if available
 505     //Check origin of the candidates
 506     if(IsDataMC()){
 507
 508       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
 509       if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
 510     }//Work with stack also
 511
 512
 513     // Check if cluster comes from a conversion in the material in front of the calorimeter
 514     // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
 515
 516     if(fCheckConversion && nCaloClusters > 1){
 517       Bool_t bConverted = kFALSE;
 518       Int_t id2 = -1;
 519
 520       //Check if set previously as converted couple, if so skip its use.
 521       if (indexConverted[icalo]) continue;
 522
 523       for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
 524         //Check if set previously as converted couple, if so skip its use.
 525         if (indexConverted[jcalo]) continue;
 526         //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
 527         AliVCluster * calo2 =  (AliVCluster*) (pl->At(jcalo));              //Get cluster kinematics
 528         Int_t evtIndex2 = 0 ;
 529         if (GetMixedEvent()) {
 530           evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
 531
 532         }
 533
 534         if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 535           calo->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
 536         else{
 537           Double_t vertex[]={0,0,0};
 538           calo->GetMomentum(mom2,vertex) ;
 539         }
 540
 541         //Check only certain regions
 542         Bool_t in2 = kTRUE;
 543         if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
 544         if(!in2) continue;
 545
 546         //Get mass of pair, if small, take this pair.
 547         //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut);
 548         if((mom+mom2).M() < fMassCut){
 549           bConverted = kTRUE;
 550           id2 = calo2->GetID();
 551           indexConverted[jcalo]=kTRUE;
 552           break;
 553         }
 554
 555       }//Mass loop
 556
 557       if(bConverted){
 558         if(fAddConvertedPairsToAOD){
 559           //Create AOD of pair analysis
 560           TLorentzVector mpair = mom+mom2;
 561           AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
 562           aodpair.SetLabel(aodph.GetLabel());
 563           //aodpair.SetInputFileIndex(input);
 564
 565           //printf("Index %d, Id %d\n",icalo, calo->GetID());
 566           //Set the indeces of the original caloclusters
 567           aodpair.SetCaloLabel(calo->GetID(),id2);
 568           aodpair.SetDetector(fCalorimeter);
 569           aodpair.SetPdg(aodph.GetPdg());
 570           aodpair.SetTag(aodph.GetTag());
 571
 572           //Add AOD with pair object to aod branch
 573           AddAODParticle(aodpair);
 574           //printf("\t \t both added pair\n");
 575         }
 576
 577         //Do not add the current calocluster
 578         continue;
 579       }//converted pair
 580     }//check conversion
 581     //printf("\t \t added single cluster %d\n",icalo);
 582
 583     //Add AOD with photon object to aod branch
 584     AddAODParticle(aodph);
 585
 586   }//loop
 587
 588   delete [] indexConverted;
 589
 590   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
 591
 592 }
 593
 594 //__________________________________________________________________
 595 void  AliAnaPhoton::MakeAnalysisFillHistograms()
 596 {
 597   //Do analysis and fill histograms
 598
 599         // Access MC information in stack if requested, check that it exists.
 600         AliStack * stack = 0x0;
 601         TParticle * primary = 0x0;
 602         TClonesArray * mcparticles0 = 0x0;
 603         //TClonesArray * mcparticles1 = 0x0;
 604         AliAODMCParticle * aodprimary = 0x0;
 605         if(IsDataMC()){
 606
 607                 if(GetReader()->ReadStack()){
 608                         stack =  GetMCStack() ;
 609                         if(!stack) {
 610                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
 611                                 abort();
 612                         }
 613
 614                 }
 615                 else if(GetReader()->ReadAODMCParticles()){
 616
 617                         //Get the list of MC particles
 618                         mcparticles0 = GetReader()->GetAODMCParticles(0);
 619                         if(!mcparticles0 && GetDebug() > 0)     {
 620                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
 621                         }
 622       //                        if(GetReader()->GetSecondInputAODTree()){
 623       //                                mcparticles1 = GetReader()->GetAODMCParticles(1);
 624       //                                if(!mcparticles1 && GetDebug() > 0)     {
 625       //                                        printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");
 626       //                                }
 627       //                        }
 628
 629                 }
 630         }// is data and MC
 631
 632         //Loop on stored AOD photons
 633         Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 634         if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 635
 636         for(Int_t iaod = 0; iaod < naod ; iaod++){
 637           AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 638           Int_t pdg = ph->GetPdg();
 639
 640           if(GetDebug() > 3)
 641             printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
 642
 643           //If PID used, fill histos with photons in Calorimeter fCalorimeter
 644           if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
 645           if(ph->GetDetector() != fCalorimeter) continue;
 646
 647           if(GetDebug() > 2)
 648             printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
 649
 650           //Fill photon histograms
 651           Float_t ptcluster  = ph->Pt();
 652           Float_t phicluster = ph->Phi();
 653           Float_t etacluster = ph->Eta();
 654           Float_t ecluster   = ph->E();
 655
 656           fhPtPhoton  ->Fill(ptcluster);
 657           fhPhiPhoton ->Fill(ptcluster,phicluster);
 658           fhEtaPhoton ->Fill(ptcluster,etacluster);
 659
 660           //Play with the MC data if available
 661           if(IsDataMC()){
 662
 663             Int_t tag =ph->GetTag();
 664
 665             if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
 666       {
 667         fhPtMCPhoton  ->Fill(ptcluster);
 668         fhPhiMCPhoton ->Fill(ptcluster,phicluster);
 669         fhEtaMCPhoton ->Fill(ptcluster,etacluster);
 670
 671         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
 672         {
 673           fhPtConversion  ->Fill(ptcluster);
 674           fhPhiConversion ->Fill(ptcluster,phicluster);
 675           fhEtaConversion ->Fill(ptcluster,etacluster);
 676         }
 677
 678         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
 679           fhPtPrompt  ->Fill(ptcluster);
 680           fhPhiPrompt ->Fill(ptcluster,phicluster);
 681           fhEtaPrompt ->Fill(ptcluster,etacluster);
 682         }
 683         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
 684         {
 685           fhPtFragmentation  ->Fill(ptcluster);
 686           fhPhiFragmentation ->Fill(ptcluster,phicluster);
 687           fhEtaFragmentation ->Fill(ptcluster,etacluster);
 688         }
 689         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
 690         {
 691           fhPtISR  ->Fill(ptcluster);
 692           fhPhiISR ->Fill(ptcluster,phicluster);
 693           fhEtaISR ->Fill(ptcluster,etacluster);
 694         }
 695         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
 696         {
 697           fhPtPi0Decay  ->Fill(ptcluster);
 698           fhPhiPi0Decay ->Fill(ptcluster,phicluster);
 699           fhEtaPi0Decay ->Fill(ptcluster,etacluster);
 700         }
 701         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
 702         {
 703           fhPtOtherDecay  ->Fill(ptcluster);
 704           fhPhiOtherDecay ->Fill(ptcluster,phicluster);
 705           fhEtaOtherDecay ->Fill(ptcluster,etacluster);
 706         }
 707       }
 708             else{
 709               fhPtUnknown  ->Fill(ptcluster);
 710               fhPhiUnknown ->Fill(ptcluster,phicluster);
 711               fhEtaUnknown ->Fill(ptcluster,etacluster);
 712
 713         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
 714         //                                      ph->GetLabel(),ph->Pt());
 715         //                for(Int_t i = 0; i < 20; i++) {
 716         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
 717         //                }
 718         //                printf("\n");
 719
 720             }
 721
 722
 723             // Access MC information in stack if requested, check that it exists.
 724             Int_t label =ph->GetLabel();
 725             if(label < 0) {
 726               printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
 727               continue;
 728             }
 729
 730             Float_t eprim   = 0;
 731             Float_t ptprim  = 0;
 732             if(GetReader()->ReadStack()){
 733
 734               if(label >=  stack->GetNtrack()) {
 735           if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 736           continue ;
 737               }
 738
 739               primary = stack->Particle(label);
 740               if(!primary){
 741           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
 742           continue;
 743               }
 744               eprim   = primary->Energy();
 745               ptprim  = primary->Pt();
 746
 747             }
 748             else if(GetReader()->ReadAODMCParticles()){
 749               //Check which is the input
 750               if(ph->GetInputFileIndex() == 0){
 751           if(!mcparticles0) continue;
 752           if(label >=  mcparticles0->GetEntriesFast()) {
 753             if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
 754                                        label, mcparticles0->GetEntriesFast());
 755             continue ;
 756           }
 757           //Get the particle
 758           aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
 759
 760               }
 761 //            else {//Second input
 762 //          if(!mcparticles1) continue;
 763 //          if(label >=  mcparticles1->GetEntriesFast()) {
 764 //            if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
 765 //                                       label, mcparticles1->GetEntriesFast());
 766 //            continue ;
 767 //          }
 768 //          //Get the particle
 769 //          aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
 770 //
 771 //            }//second input
 772
 773               if(!aodprimary){
 774           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
 775           continue;
 776               }
 777
 778               eprim   = aodprimary->E();
 779               ptprim  = aodprimary->Pt();
 780
 781             }
 782
 783             fh2E     ->Fill(ecluster, eprim);
 784             fh2Pt    ->Fill(ptcluster, ptprim);
 785             fhDeltaE ->Fill(eprim-ecluster);
 786             fhDeltaPt->Fill(ptprim-ptcluster);
 787             if(eprim > 0)  fhRatioE  ->Fill(ecluster/eprim);
 788             if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
 789
 790           }//Histograms with MC
 791
 792         }// aod loop
 793
 794 }
 795
 796
 797 //__________________________________________________________________
 798 void AliAnaPhoton::Print(const Option_t * opt) const
 799 {
 800   //Print some relevant parameters set for the analysis
 801
 802   if(! opt)
 803     return;
 804
 805   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 806   AliAnaPartCorrBaseClass::Print(" ");
 807
 808   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
 809   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
 810   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
 811   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
 812   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
 813   printf("Check Pair Conversion                = %d\n",fCheckConversion);
 814   printf("Add conversion pair to AOD           = %d\n",fAddConvertedPairsToAOD);
 815   printf("Conversion pair mass cut             = %f\n",fMassCut);
 816   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
 817   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
 818   printf("    \n") ;
 819
 820 }