PWG4/PartCorrDep/AliAnaPhoton.cxx

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