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
 365   //Select the Calorimeter of the photon
 366   TObjArray * pl = 0x0;
 367   if(fCalorimeter == "PHOS")
 368     pl = GetAODPHOS();
 369   else if (fCalorimeter == "EMCAL")
 370     pl = GetAODEMCAL();
 371
 372   if(!pl) {
 373     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
 374     return;
 375   }
 376
 377   //Fill AODParticle with PHOS/EMCAL aods
 378   TLorentzVector mom, mom2 ;
 379   Int_t nCaloClusters = pl->GetEntriesFast();
 380   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
 381   Bool_t * indexConverted = new Bool_t[nCaloClusters];
 382   for (Int_t i = 0; i < nCaloClusters; i++)
 383     indexConverted[i] = kFALSE;
 384
 385   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
 386
 387           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
 388     //printf("calo %d, %f\n",icalo,calo->E());
 389
 390     //Get the index where the cluster comes, to retrieve the corresponding vertex
 391     Int_t evtIndex = 0 ;
 392     if (GetMixedEvent()) {
 393       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 394     }
 395
 396     //Cluster selection, not charged, with photon id and in fiducial cut
 397
 398     //Input from second AOD?
 399     //Int_t input = 0;
 400     //    if (fCalorimeter == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= icalo)
 401     //      input = 1 ;
 402     //    else if(fCalorimeter == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= icalo)
 403     //      input = 1;
 404
 405     //Get Momentum vector,
 406     //if (input == 0)
 407     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 408       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
 409     else{
 410       Double_t vertex[]={0,0,0};
 411       calo->GetMomentum(mom,vertex) ;
 412     }
 413     //printf("AliAnaPhoton::MakeAnalysisFillAOD(): Vertex : %f,%f,%f\n",GetVertex(evtIndex)[0] ,GetVertex(evtIndex)[1],GetVertex(evtIndex)[2]);
 414
 415     //    else if(input == 1)
 416     //      calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line
 417
 418     //If too small or big pt, skip it
 419     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
 420     Double_t tof = calo->GetTOF()*1e9;
 421
 422     if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
 423
 424     if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) continue;
 425
 426     //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(),
 427       //     mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 428
 429     //Check acceptance selection
 430     if(IsFiducialCutOn()){
 431       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 432       if(! in ) continue ;
 433     }
 434     //printf("Fiducial cut passed \n");
 435
 436     //Create AOD for analysis
 437     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
 438     Int_t label = calo->GetLabel();
 439     aodph.SetLabel(label);
 440     //aodph.SetInputFileIndex(input);
 441
 442     //printf("Index %d, Id %d\n",icalo, calo->GetID());
 443     //Set the indeces of the original caloclusters
 444     aodph.SetCaloLabel(calo->GetID(),-1);
 445     aodph.SetDetector(fCalorimeter);
 446     if(GetDebug() > 1)
 447       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());
 448
 449     //Check Distance to Bad channel, set bit.
 450     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 451     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 452     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
 453       continue ;
 454
 455     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
 456
 457     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
 458     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
 459     else                         aodph.SetDistToBad(0) ;
 460
 461     //Skip matched clusters with tracks
 462     if(fRejectTrackMatch && IsTrackMatched(calo)) continue ;
 463     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - TrackMatching cut passed \n");
 464
 465     //Check PID
 466     //PID selection or bit setting
 467     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
 468       //Get most probable PID, check PID weights (in MC this option is mandatory)
 469       aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 470       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
 471       //If primary is not photon, skip it.
 472       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
 473     }
 474     else if(IsCaloPIDOn()){
 475
 476       //Get most probable PID, 2 options check PID weights
 477       //or redo PID, recommended option for EMCal.
 478       if(!IsCaloPIDRecalculationOn())
 479         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
 480       else
 481         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
 482
 483       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
 484
 485       //If cluster does not pass pid, not photon, skip it.
 486       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
 487
 488     }
 489     else{
 490       //Set PID bits for later selection (AliAnaPi0 for example)
 491       //GetPDG already called in SetPIDBits.
 492       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
 493       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
 494     }
 495
 496     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
 497
 498     //Play with the MC stack if available
 499     //Check origin of the candidates
 500     if(IsDataMC()){
 501
 502       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
 503       if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
 504     }//Work with stack also
 505
 506
 507     // Check if cluster comes from a conversion in the material in front of the calorimeter
 508     // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
 509
 510     if(fCheckConversion && nCaloClusters > 1){
 511       Bool_t bConverted = kFALSE;
 512       Int_t id2 = -1;
 513
 514       //Check if set previously as converted couple, if so skip its use.
 515       if (indexConverted[icalo]) continue;
 516
 517       for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
 518         //Check if set previously as converted couple, if so skip its use.
 519         if (indexConverted[jcalo]) continue;
 520         //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
 521         AliVCluster * calo2 =  (AliVCluster*) (pl->At(jcalo));              //Get cluster kinematics
 522         Int_t evtIndex2 = 0 ;
 523         if (GetMixedEvent()) {
 524           evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
 525
 526         }
 527
 528         if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 529           calo->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
 530         else{
 531           Double_t vertex[]={0,0,0};
 532           calo->GetMomentum(mom2,vertex) ;
 533         }
 534
 535         //Check only certain regions
 536         Bool_t in2 = kTRUE;
 537         if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
 538         if(!in2) continue;
 539
 540         //Get mass of pair, if small, take this pair.
 541         //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut);
 542         if((mom+mom2).M() < fMassCut){
 543           bConverted = kTRUE;
 544           id2 = calo2->GetID();
 545           indexConverted[jcalo]=kTRUE;
 546           break;
 547         }
 548
 549       }//Mass loop
 550
 551       if(bConverted){
 552         if(fAddConvertedPairsToAOD){
 553           //Create AOD of pair analysis
 554           TLorentzVector mpair = mom+mom2;
 555           AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
 556           aodpair.SetLabel(aodph.GetLabel());
 557           //aodpair.SetInputFileIndex(input);
 558
 559           //printf("Index %d, Id %d\n",icalo, calo->GetID());
 560           //Set the indeces of the original caloclusters
 561           aodpair.SetCaloLabel(calo->GetID(),id2);
 562           aodpair.SetDetector(fCalorimeter);
 563           aodpair.SetPdg(aodph.GetPdg());
 564           aodpair.SetTag(aodph.GetTag());
 565
 566           //Add AOD with pair object to aod branch
 567           AddAODParticle(aodpair);
 568           //printf("\t \t both added pair\n");
 569         }
 570
 571         //Do not add the current calocluster
 572         continue;
 573       }//converted pair
 574     }//check conversion
 575     //printf("\t \t added single cluster %d\n",icalo);
 576
 577     //Add AOD with photon object to aod branch
 578     AddAODParticle(aodph);
 579
 580   }//loop
 581
 582   delete [] indexConverted;
 583
 584   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
 585
 586 }
 587
 588 //__________________________________________________________________
 589 void  AliAnaPhoton::MakeAnalysisFillHistograms()
 590 {
 591   //Do analysis and fill histograms
 592
 593         // Access MC information in stack if requested, check that it exists.
 594         AliStack * stack = 0x0;
 595         TParticle * primary = 0x0;
 596         TClonesArray * mcparticles0 = 0x0;
 597         //TClonesArray * mcparticles1 = 0x0;
 598         AliAODMCParticle * aodprimary = 0x0;
 599         if(IsDataMC()){
 600
 601                 if(GetReader()->ReadStack()){
 602                         stack =  GetMCStack() ;
 603                         if(!stack) {
 604                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
 605                                 abort();
 606                         }
 607
 608                 }
 609                 else if(GetReader()->ReadAODMCParticles()){
 610
 611                         //Get the list of MC particles
 612                         mcparticles0 = GetReader()->GetAODMCParticles(0);
 613                         if(!mcparticles0 && GetDebug() > 0)     {
 614                                 printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
 615                         }
 616       //                        if(GetReader()->GetSecondInputAODTree()){
 617       //                                mcparticles1 = GetReader()->GetAODMCParticles(1);
 618       //                                if(!mcparticles1 && GetDebug() > 0)     {
 619       //                                        printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");
 620       //                                }
 621       //                        }
 622
 623                 }
 624         }// is data and MC
 625
 626         //Loop on stored AOD photons
 627         Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 628         if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 629
 630         for(Int_t iaod = 0; iaod < naod ; iaod++){
 631           AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 632           Int_t pdg = ph->GetPdg();
 633
 634           if(GetDebug() > 3)
 635             printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
 636
 637           //If PID used, fill histos with photons in Calorimeter fCalorimeter
 638           if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
 639           if(ph->GetDetector() != fCalorimeter) continue;
 640
 641           if(GetDebug() > 2)
 642             printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
 643
 644           //Fill photon histograms
 645           Float_t ptcluster  = ph->Pt();
 646           Float_t phicluster = ph->Phi();
 647           Float_t etacluster = ph->Eta();
 648           Float_t ecluster   = ph->E();
 649
 650           fhPtPhoton  ->Fill(ptcluster);
 651           fhPhiPhoton ->Fill(ptcluster,phicluster);
 652           fhEtaPhoton ->Fill(ptcluster,etacluster);
 653
 654           //Play with the MC data if available
 655           if(IsDataMC()){
 656
 657             Int_t tag =ph->GetTag();
 658
 659             if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
 660       {
 661         fhPtMCPhoton  ->Fill(ptcluster);
 662         fhPhiMCPhoton ->Fill(ptcluster,phicluster);
 663         fhEtaMCPhoton ->Fill(ptcluster,etacluster);
 664
 665         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
 666         {
 667           fhPtConversion  ->Fill(ptcluster);
 668           fhPhiConversion ->Fill(ptcluster,phicluster);
 669           fhEtaConversion ->Fill(ptcluster,etacluster);
 670         }
 671
 672         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
 673           fhPtPrompt  ->Fill(ptcluster);
 674           fhPhiPrompt ->Fill(ptcluster,phicluster);
 675           fhEtaPrompt ->Fill(ptcluster,etacluster);
 676         }
 677         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
 678         {
 679           fhPtFragmentation  ->Fill(ptcluster);
 680           fhPhiFragmentation ->Fill(ptcluster,phicluster);
 681           fhEtaFragmentation ->Fill(ptcluster,etacluster);
 682         }
 683         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
 684         {
 685           fhPtISR  ->Fill(ptcluster);
 686           fhPhiISR ->Fill(ptcluster,phicluster);
 687           fhEtaISR ->Fill(ptcluster,etacluster);
 688         }
 689         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
 690         {
 691           fhPtPi0Decay  ->Fill(ptcluster);
 692           fhPhiPi0Decay ->Fill(ptcluster,phicluster);
 693           fhEtaPi0Decay ->Fill(ptcluster,etacluster);
 694         }
 695         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) || GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
 696         {
 697           fhPtOtherDecay  ->Fill(ptcluster);
 698           fhPhiOtherDecay ->Fill(ptcluster,phicluster);
 699           fhEtaOtherDecay ->Fill(ptcluster,etacluster);
 700         }
 701       }
 702             else{
 703               fhPtUnknown  ->Fill(ptcluster);
 704               fhPhiUnknown ->Fill(ptcluster,phicluster);
 705               fhEtaUnknown ->Fill(ptcluster,etacluster);
 706
 707         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
 708         //                                      ph->GetLabel(),ph->Pt());
 709         //                for(Int_t i = 0; i < 20; i++) {
 710         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
 711         //                }
 712         //                printf("\n");
 713
 714             }
 715
 716
 717             // Access MC information in stack if requested, check that it exists.
 718             Int_t label =ph->GetLabel();
 719             if(label < 0) {
 720               printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
 721               continue;
 722             }
 723
 724             Float_t eprim   = 0;
 725             Float_t ptprim  = 0;
 726             if(GetReader()->ReadStack()){
 727
 728               if(label >=  stack->GetNtrack()) {
 729           if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
 730           continue ;
 731               }
 732
 733               primary = stack->Particle(label);
 734               if(!primary){
 735           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
 736           continue;
 737               }
 738               eprim   = primary->Energy();
 739               ptprim  = primary->Pt();
 740
 741             }
 742             else if(GetReader()->ReadAODMCParticles()){
 743               //Check which is the input
 744               if(ph->GetInputFileIndex() == 0){
 745           if(!mcparticles0) continue;
 746           if(label >=  mcparticles0->GetEntriesFast()) {
 747             if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
 748                                        label, mcparticles0->GetEntriesFast());
 749             continue ;
 750           }
 751           //Get the particle
 752           aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
 753
 754               }
 755 //            else {//Second input
 756 //          if(!mcparticles1) continue;
 757 //          if(label >=  mcparticles1->GetEntriesFast()) {
 758 //            if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
 759 //                                       label, mcparticles1->GetEntriesFast());
 760 //            continue ;
 761 //          }
 762 //          //Get the particle
 763 //          aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
 764 //
 765 //            }//second input
 766
 767               if(!aodprimary){
 768           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
 769           continue;
 770               }
 771
 772               eprim   = aodprimary->E();
 773               ptprim  = aodprimary->Pt();
 774
 775             }
 776
 777             fh2E     ->Fill(ecluster, eprim);
 778             fh2Pt    ->Fill(ptcluster, ptprim);
 779             fhDeltaE ->Fill(eprim-ecluster);
 780             fhDeltaPt->Fill(ptprim-ptcluster);
 781             if(eprim > 0)  fhRatioE  ->Fill(ecluster/eprim);
 782             if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
 783
 784           }//Histograms with MC
 785
 786         }// aod loop
 787
 788 }
 789
 790
 791 //__________________________________________________________________
 792 void AliAnaPhoton::Print(const Option_t * opt) const
 793 {
 794   //Print some relevant parameters set for the analysis
 795
 796   if(! opt)
 797     return;
 798
 799   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 800   AliAnaPartCorrBaseClass::Print(" ");
 801
 802   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
 803   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
 804   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
 805   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
 806   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
 807   printf("Check Pair Conversion                = %d\n",fCheckConversion);
 808   printf("Add conversion pair to AOD           = %d\n",fAddConvertedPairsToAOD);
 809   printf("Conversion pair mass cut             = %f\n",fMassCut);
 810   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
 811   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
 812   printf("    \n") ;
 813
 814 }