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