PWG4/PartCorrDep/AliAnaPi0EbE.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes iGetEntriesFast(s hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 /* $Id: AliAnaPi0EbE.cxx 28688 2008-09-11 15:04:07Z gconesab $ */\r
  16 \r
  17 //_________________________________________________________________________\r
  18 // Class for the analysis of high pT pi0 event by event\r
  19 // Pi0 identified by one of the following:\r
  20 //  -Invariant mass of 2 cluster in calorimeter\r
  21 //  -Shower shape analysis in calorimeter\r
  22 //  -Invariant mass of one cluster in calorimeter and one photon reconstructed in TPC (in near future)\r
  23 //\r
  24 // -- Author: Gustavo Conesa (LNF-INFN) &  Raphaelle Ichou (SUBATECH)\r
  25 //////////////////////////////////////////////////////////////////////////////\r
  26   \r
  27   \r
  28 // --- ROOT system --- \r
  29 #include <TList.h>\r
  30 //#include "Riostream.h"\r
  31 \r
  32 // --- Analysis system --- \r
  33 #include "AliAnaPi0EbE.h" \r
  34 #include "AliLog.h"\r
  35 #include "AliCaloTrackReader.h"\r
  36 #include "AliIsolationCut.h"\r
  37 #include "AliNeutralMesonSelection.h"\r
  38 #include "AliCaloPID.h"\r
  39 #include "AliAODPWG4ParticleCorrelation.h"\r
  40 \r
  41 ClassImp(AliAnaPi0EbE)\r
  42   \r
  43 //____________________________________________________________________________\r
  44 AliAnaPi0EbE::AliAnaPi0EbE() : \r
  45 AliAnaPartCorrBaseClass(),  fAnaType(kIMCalo),fCalorimeter(""),\r
  46 fMinDist(0.),fMinDist2(0.),fMinDist3(0.),       \r
  47 fInputAODGammaConv(0x0),fInputAODGammaConvName(""),\r
  48 fhPtPi0(0),fhPhiPi0(0),fhEtaPi0(0), \r
  49 fhPtMCNoPi0(0),fhPhiMCNoPi0(0),fhEtaMCNoPi0(0), \r
  50 fhPtMCPi0(0),fhPhiMCPi0(0),fhEtaMCPi0(0)\r
  51 {\r
  52         //default ctor\r
  53         \r
  54         //Initialize parameters\r
  55         InitParameters();\r
  56         \r
  57 }\r
  58 \r
  59 //____________________________________________________________________________\r
  60 AliAnaPi0EbE::AliAnaPi0EbE(const AliAnaPi0EbE & p) : \r
  61 AliAnaPartCorrBaseClass(p),  fAnaType(p.fAnaType), fCalorimeter(p.fCalorimeter),\r
  62 fMinDist(p.fMinDist),fMinDist2(p.fMinDist2), fMinDist3(p.fMinDist3),\r
  63 fInputAODGammaConv(new TClonesArray(*p.fInputAODGammaConv)),\r
  64 fInputAODGammaConvName(p.fInputAODGammaConvName),\r
  65 fhPtPi0(p.fhPtPi0),fhPhiPi0(p.fhPhiPi0),fhEtaPi0(p.fhEtaPi0), \r
  66 fhPtMCNoPi0(p.fhPtMCNoPi0),fhPhiMCNoPi0(p.fhPhiMCNoPi0),fhEtaMCNoPi0(p.fhEtaMCNoPi0), \r
  67 fhPtMCPi0(p.fhPtMCPi0),fhPhiMCPi0(p.fhPhiMCPi0),fhEtaMCPi0(p.fhEtaMCPi0) \r
  68 {\r
  69         // cpy ctor\r
  70 }\r
  71 \r
  72 //_________________________________________________________________________\r
  73 AliAnaPi0EbE & AliAnaPi0EbE::operator = (const AliAnaPi0EbE & p)\r
  74 {\r
  75         // assignment operator\r
  76         \r
  77         if(&p == this) return *this;\r
  78         \r
  79         fAnaType     = p.fAnaType ;\r
  80         fCalorimeter = p.fCalorimeter ;\r
  81         fMinDist     = p.fMinDist;\r
  82         fMinDist2    = p.fMinDist2;\r
  83         fMinDist3    = p.fMinDist3;\r
  84         \r
  85         fInputAODGammaConv     = new TClonesArray(*p.fInputAODGammaConv) ;\r
  86         fInputAODGammaConvName = p.fInputAODGammaConvName ;\r
  87         \r
  88         fhPtPi0 = p.fhPtPi0; fhPhiPi0 = p.fhPhiPi0; fhEtaPi0 = p.fhEtaPi0;  \r
  89         fhPtMCNoPi0 = p.fhPtMCNoPi0; fhPhiMCNoPi0 = p.fhPhiMCNoPi0; fhEtaMCNoPi0 = p.fhEtaMCPi0;  \r
  90         fhPtMCPi0 = p.fhPtMCPi0; fhPhiMCPi0 = p.fhPhiMCPi0; fhEtaMCPi0 = p.fhEtaMCPi0;  \r
  91         \r
  92         return *this;\r
  93         \r
  94 }\r
  95 \r
  96 //____________________________________________________________________________\r
  97 AliAnaPi0EbE::~AliAnaPi0EbE() \r
  98 {\r
  99         //dtor\r
 100         if(fInputAODGammaConv){\r
 101                 fInputAODGammaConv->Clear() ; \r
 102                 delete fInputAODGammaConv ;\r
 103         }\r
 104 }\r
 105 \r
 106 //________________________________________________________________________\r
 107 TList *  AliAnaPi0EbE::GetCreateOutputObjects()\r
 108 {  \r
 109         // Create histograms to be saved in output file and \r
 110         // store them in outputContainer\r
 111         TList * outputContainer = new TList() ; \r
 112         outputContainer->SetName("Pi0EbEHistos") ; \r
 113         \r
 114         Int_t nptbins  = GetHistoNPtBins();\r
 115         Int_t nphibins = GetHistoNPhiBins();\r
 116         Int_t netabins = GetHistoNEtaBins();\r
 117         Float_t ptmax  = GetHistoPtMax();\r
 118         Float_t phimax = GetHistoPhiMax();\r
 119         Float_t etamax = GetHistoEtaMax();\r
 120         Float_t ptmin  = GetHistoPtMin();\r
 121         Float_t phimin = GetHistoPhiMin();\r
 122         Float_t etamin = GetHistoEtaMin();      \r
 123         \r
 124         fhPtPi0  = new TH1F("hPtPi0","Number of identified  #pi^{0} decay",nptbins,ptmin,ptmax); \r
 125         fhPtPi0->SetYTitle("N");\r
 126         fhPtPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");\r
 127         outputContainer->Add(fhPtPi0) ; \r
 128         \r
 129         fhPhiPi0  = new TH2F\r
 130     ("hPhiPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
 131         fhPhiPi0->SetYTitle("#phi");\r
 132         fhPhiPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
 133         outputContainer->Add(fhPhiPi0) ; \r
 134         \r
 135         fhEtaPi0  = new TH2F\r
 136     ("hEtaPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
 137         fhEtaPi0->SetYTitle("#eta");\r
 138         fhEtaPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
 139         outputContainer->Add(fhEtaPi0) ;\r
 140         \r
 141         if(IsDataMC()){\r
 142                 \r
 143                 fhPtMCPi0  = new TH1F("hPtMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax); \r
 144                 fhPtMCPi0->SetYTitle("N");\r
 145                 fhPtMCPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");\r
 146                 outputContainer->Add(fhPtMCPi0) ; \r
 147                 \r
 148                 fhPhiMCPi0  = new TH2F\r
 149                 ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
 150                 fhPhiMCPi0->SetYTitle("#phi");\r
 151                 fhPhiMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
 152                 outputContainer->Add(fhPhiMCPi0) ; \r
 153                 \r
 154                 fhEtaMCPi0  = new TH2F\r
 155                 ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
 156                 fhEtaMCPi0->SetYTitle("#eta");\r
 157                 fhEtaMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
 158                 outputContainer->Add(fhEtaMCPi0) ;\r
 159                 \r
 160                 fhPtMCNoPi0  = new TH1F("hPtMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax); \r
 161                 fhPtMCNoPi0->SetYTitle("N");\r
 162                 fhPtMCNoPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");\r
 163                 outputContainer->Add(fhPtMCNoPi0) ; \r
 164                 \r
 165                 fhPhiMCNoPi0  = new TH2F\r
 166                 ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
 167                 fhPhiMCNoPi0->SetYTitle("#phi");\r
 168                 fhPhiMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
 169                 outputContainer->Add(fhPhiMCNoPi0) ; \r
 170                 \r
 171                 fhEtaMCNoPi0  = new TH2F\r
 172                 ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
 173                 fhEtaMCNoPi0->SetYTitle("#eta");\r
 174                 fhEtaMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
 175                 outputContainer->Add(fhEtaMCNoPi0) ;\r
 176                 \r
 177         }//Histos with MC\r
 178     \r
 179         \r
 180         //Keep neutral meson selection histograms if requiered\r
 181         //Setting done in AliNeutralMesonSelection\r
 182         \r
 183         if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){\r
 184                 \r
 185         TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;\r
 186         if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())\r
 187                         for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;\r
 188         }\r
 189         \r
 190         //Save parameters used for analysis\r
 191         TString parList ; //this will be list of parameters used for this analysis.\r
 192         char onePar[255] ;\r
 193         \r
 194         sprintf(onePar,"--- AliAnaPi0EbE ---\n") ;\r
 195         parList+=onePar ;       \r
 196         sprintf(onePar,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;\r
 197         parList+=onePar ;\r
 198         \r
 199         if(fAnaType == kSSCalo){\r
 200                 sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;\r
 201                 parList+=onePar ;\r
 202                 sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;\r
 203                 parList+=onePar ;\r
 204                 sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;\r
 205                 parList+=onePar ;\r
 206                 sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;\r
 207                 parList+=onePar ;\r
 208         }\r
 209         \r
 210         //Get parameters set in base class.\r
 211         parList += GetBaseParametersList() ;\r
 212         \r
 213         //Get parameters set in PID class.\r
 214         if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;\r
 215         \r
 216         TObjString *oString= new TObjString(parList) ;\r
 217         outputContainer->Add(oString);\r
 218         \r
 219         return outputContainer ;\r
 220         \r
 221 }\r
 222 \r
 223 //__________________________________________________________________\r
 224 void  AliAnaPi0EbE::MakeAnalysisFillAOD() \r
 225 {\r
 226         //Do analysis and fill aods\r
 227         \r
 228         switch(fAnaType) \r
 229         {\r
 230                 case kIMCalo:\r
 231                         MakeInvMassInCalorimeter();\r
 232                         break;\r
 233                         \r
 234                 case kSSCalo:\r
 235                         MakeShowerShapeIdentification();\r
 236                         break;\r
 237                         \r
 238                 case kIMCaloTracks:\r
 239                         MakeInvMassInCalorimeterAndCTS();\r
 240                         break;\r
 241                         \r
 242         }\r
 243 }\r
 244 \r
 245 //__________________________________________________________________\r
 246 void  AliAnaPi0EbE::MakeInvMassInCalorimeter() \r
 247 {\r
 248         //Do analysis and fill aods\r
 249         //Search for the photon decay in calorimeters\r
 250         //Read photon list from AOD, produced in class AliAnaPhoton\r
 251         //Check if 2 photons have the mass of the pi0.\r
 252         \r
 253         TLorentzVector mom1;\r
 254         TLorentzVector mom2;\r
 255         TLorentzVector mom ;\r
 256         Int_t tag1 =-1;\r
 257         Int_t tag2 =-1;\r
 258         Int_t tag = -1;\r
 259         \r
 260         if(!GetInputAODBranch())\r
 261                 AliFatal(Form("MakeInvMassInCalo: No input calo photons in AOD with name branch < %s > \n",GetInputAODName().Data()));\r
 262         \r
 263         for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){\r
 264                 AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));\r
 265                 mom1 = *(photon1->Momentum());\r
 266                 \r
 267                 //Play with the MC stack if available\r
 268                 \r
 269                 for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast()-1; jphoton++){\r
 270                         \r
 271                         AliAODPWG4ParticleCorrelation * photon2 =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(jphoton));\r
 272                         mom2 = *(photon2->Momentum());\r
 273                         //Select good pair (good phi, pt cuts, aperture and invariant mass)\r
 274                         if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))\r
 275                         {\r
 276                                 if(GetDebug()>1) printf("Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());\r
 277                                 \r
 278                                 \r
 279                                 if(IsDataMC()){\r
 280                                         //Check origin of the candidates\r
 281                                         tag1 = GetCaloPID()->CheckOrigin(photon1->GetLabel(), GetMCStack());\r
 282                                         tag2 = GetCaloPID()->CheckOrigin(photon2->GetLabel(), GetMCStack());\r
 283                                         if(GetDebug() > 0) printf("Origin of: photon1 %d; photon2 %d \n",tag1, tag2);\r
 284                                         if(tag1 == AliCaloPID::kMCPi0Decay && tag2 == AliCaloPID::kMCPi0Decay) tag = AliCaloPID::kMCPi0;\r
 285                                 }//Work with stack also   \r
 286                                 \r
 287                                 //Create AOD for analysis\r
 288                                 mom = mom1+mom2;\r
 289                                 AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);\r
 290                                 //pi0.SetLabel(calo->GetLabel(0));\r
 291                                 pi0.SetPdg(AliCaloPID::kPi0);\r
 292                                 pi0.SetDetector(photon1->GetDetector());\r
 293                                 pi0.SetTag(tag);  \r
 294                                 //Set the indeces of the original caloclusters  \r
 295                                 pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));\r
 296                                 AddAODParticle(pi0);\r
 297                         }//pi0\r
 298                 }//2n photon loop\r
 299                 \r
 300         }//1st photon loop\r
 301         \r
 302         if(GetDebug() > 1) printf("End fill AODs \n");  \r
 303         \r
 304 }\r
 305 \r
 306 //__________________________________________________________________\r
 307 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() \r
 308 {\r
 309         //Do analysis and fill aods\r
 310         //Search for the photon decay in calorimeters\r
 311         //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion\r
 312         //Check if 2 photons have the mass of the pi0.\r
 313         \r
 314         TLorentzVector mom1;\r
 315         TLorentzVector mom2;\r
 316         TLorentzVector mom ;\r
 317         Int_t tag1 =-1;\r
 318         Int_t tag2 =-1;\r
 319         Int_t tag = -1;\r
 320         \r
 321         if(!GetInputAODBranch())\r
 322                 AliFatal(Form("MakeInvMassInCalo: No input calo photons in AOD branch with name < %s > \n",GetInputAODName().Data()));\r
 323         \r
 324         for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){\r
 325                 AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));\r
 326                 mom1 = *(photon1->Momentum());\r
 327                 \r
 328                 //Play with the MC stack if available\r
 329                 fInputAODGammaConv = (TClonesArray *) GetReader()->GetAOD()->FindListObject(fInputAODGammaConvName);\r
 330                 if(!fInputAODGammaConv) \r
 331                         AliFatal(Form("MakeInvMassInCaloAndCTS: No input gamma conversions in AOD branch with name < %s >",fInputAODGammaConvName.Data()));     \r
 332                 \r
 333                 for(Int_t jphoton = iphoton+1; jphoton < fInputAODGammaConv->GetEntriesFast()-1; jphoton++){\r
 334                         AliAODPWG4ParticleCorrelation * photon2 =  (AliAODPWG4ParticleCorrelation*) (fInputAODGammaConv->At(jphoton));\r
 335                         mom2 = *(photon2->Momentum());\r
 336                         //Select good pair (good phi, pt cuts, aperture and invariant mass)\r
 337                         if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){\r
 338                                 if(GetDebug() > 1) printf("Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());\r
 339                                 \r
 340                                 if(IsDataMC()){\r
 341                                         //Check origin of the candidates\r
 342                                         tag1 = GetCaloPID()->CheckOrigin(photon1->GetLabel(), GetMCStack());\r
 343                                         tag2 = GetCaloPID()->CheckOrigin(photon2->GetLabel(), GetMCStack());\r
 344                                         if(GetDebug() > 0) printf("Origin of: photon1 %d; photon2 %d \n",tag1, tag2);\r
 345                                         if(tag1 == AliCaloPID::kMCPi0Decay && tag2 == AliCaloPID::kMCPi0Decay) tag = AliCaloPID::kMCPi0;\r
 346                                 }//Work with stack also   \r
 347                                 \r
 348                                 //Create AOD for analysis\r
 349                                 mom = mom1+mom2;\r
 350                                 AliAODPWG4ParticleCorrelation pi0 = AliAODPWG4ParticleCorrelation(mom);\r
 351                                 //pi0.SetLabel(calo->GetLabel(0));\r
 352                                 pi0.SetPdg(AliCaloPID::kPi0);\r
 353                                 pi0.SetDetector(photon1->GetDetector());\r
 354                                 pi0.SetTag(tag);\r
 355                                 //Set the indeces of the original tracks or caloclusters  \r
 356                                 pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);\r
 357                                 pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));\r
 358                                 AddAODParticle(pi0);\r
 359                         }//pi0\r
 360                 }//2n photon loop\r
 361                 \r
 362         }//1st photon loop\r
 363         \r
 364         if(GetDebug() > 1) printf("End fill AODs \n");  \r
 365         \r
 366 }\r
 367 \r
 368 \r
 369 //__________________________________________________________________\r
 370 void  AliAnaPi0EbE::MakeShowerShapeIdentification() \r
 371 {\r
 372         //Search for pi0 in fCalorimeter with shower shape analysis \r
 373         \r
 374         TClonesArray * pl = new TClonesArray; \r
 375         \r
 376         //Get vertex for photon momentum calculation\r
 377         Double_t vertex[]={0,0,0} ; //vertex ;\r
 378         if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);\r
 379         \r
 380         //Select the Calorimeter of the photon\r
 381         if(fCalorimeter == "PHOS")\r
 382                 pl = GetAODPHOS();\r
 383         else if (fCalorimeter == "EMCAL")\r
 384                 pl = GetAODEMCAL();\r
 385         //Fill AODCaloClusters and AODParticle with PHOS aods\r
 386         TLorentzVector mom ;\r
 387         for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){\r
 388                 AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));        \r
 389                 \r
 390                 //Cluster selection, not charged, with photon id and in fidutial cut\r
 391                 //Get Momentum vector, \r
 392                 calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line\r
 393                 //If too small or big pt, skip it\r
 394                 if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; \r
 395                 //Check acceptance selection\r
 396                 if(IsFidutialCutOn()){\r
 397                         Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;\r
 398                         if(! in ) continue ;\r
 399                 }\r
 400                 \r
 401                 //Create AOD for analysis\r
 402                 AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);\r
 403                 aodph.SetLabel(calo->GetLabel(0));\r
 404                 aodph.SetDetector(fCalorimeter);\r
 405                 if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Min pt cut and fidutial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());    \r
 406                 \r
 407                 //Check Distance to Bad channel, set bit.\r
 408                 Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel\r
 409                 if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;\r
 410                 if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)\r
 411                         continue ;\r
 412                 \r
 413                 if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Bad channel cut passed %4.2f\n",distBad);\r
 414                 \r
 415                 if(distBad > fMinDist3) aodph.SetDistToBad(2) ;\r
 416                 else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; \r
 417                 else aodph.SetDistToBad(0) ;\r
 418                 \r
 419                 //Check PID\r
 420                 //PID selection or bit setting\r
 421                 if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){\r
 422                         //Get most probable PID, check PID weights (in MC this option is mandatory)\r
 423                         aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights\r
 424                         if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PDG of identified particle %d\n",aodph.GetPdg());\r
 425                         //If primary is not photon, skip it.\r
 426                         if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;\r
 427                 }                                       \r
 428                 else if(IsCaloPIDOn()){\r
 429                         //Skip matched clusters with tracks\r
 430                         if(calo->GetNTracksMatched() > 0) continue ;\r
 431                         \r
 432                         //Get most probable PID, 2 options check PID weights \r
 433                         //or redo PID, recommended option for EMCal.            \r
 434                         if(!IsCaloPIDRecalculationOn())\r
 435                                 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights\r
 436                         else\r
 437                                 aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated\r
 438                         \r
 439                         if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PDG of identified particle %d\n",aodph.GetPdg());\r
 440                         \r
 441                         //If cluster does not pass pid, not photon, skip it.\r
 442                         if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;                    \r
 443                         \r
 444                 }\r
 445                 else{\r
 446                         //Set PID bits for later selection (AliAnaPi0 for example)\r
 447                         //GetPDG already called in SetPIDBits.\r
 448                         GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);\r
 449                         if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PID Bits set \n");            \r
 450                 }\r
 451                 \r
 452                 if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());\r
 453                 \r
 454                 //Play with the MC stack if available\r
 455                 //Check origin of the candidates\r
 456                 if(IsDataMC()){\r
 457                         aodph.SetTag(GetCaloPID()->CheckOrigin(calo->GetLabel(0),GetMCStack()));\r
 458                         if(GetDebug() > 0) printf("AliAnaPi0EbE::FillAOD: Origin of candidate %d\n",aodph.GetTag());\r
 459                 }//Work with stack also   \r
 460                 \r
 461                 //Add AOD with photon object to aod branch\r
 462                 AddAODParticle(aodph);\r
 463                 \r
 464         }//loop\r
 465         \r
 466         if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: End fill AODs \n");  \r
 467         \r
 468 }\r
 469 //__________________________________________________________________\r
 470 void  AliAnaPi0EbE::MakeAnalysisFillHistograms() \r
 471 {\r
 472         //Do analysis and fill histograms\r
 473         \r
 474         if(!GetOutputAODBranch())\r
 475                 AliFatal(Form("MakeInvMassInCalo: No output pi0 in AOD branch with name < %s > \n",GetOutputAODName().Data()));\r
 476         \r
 477         //Loop on stored AOD pi0\r
 478         Int_t naod = GetOutputAODBranch()->GetEntriesFast();\r
 479         if(GetDebug() > 0) printf("pi0 aod branch entries %d\n", naod);\r
 480         \r
 481         for(Int_t iaod = 0; iaod < naod ; iaod++){\r
 482                 \r
 483                 AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));\r
 484                 Int_t pdg = pi0->GetPdg();\r
 485                 \r
 486                 if(pdg != AliCaloPID::kPi0) continue;              \r
 487                 \r
 488                 //Fill pi0 histograms \r
 489                 Float_t pt = pi0->Pt();\r
 490                 Float_t phi = pi0->Phi();\r
 491                 Float_t eta = pi0->Eta();\r
 492                 \r
 493                 \r
 494                 fhPtPi0  ->Fill(pt);\r
 495                 fhPhiPi0 ->Fill(pt,phi);\r
 496                 fhEtaPi0 ->Fill(pt,eta);\r
 497                 \r
 498                 if(IsDataMC()){\r
 499                         if(pi0->GetTag()== AliCaloPID::kMCPi0Decay){\r
 500                                 fhPtMCPi0  ->Fill(pt);\r
 501                                 fhPhiMCPi0 ->Fill(pt,phi);\r
 502                                 fhEtaMCPi0 ->Fill(pt,eta);\r
 503                         }\r
 504                         else{\r
 505                                 fhPtMCNoPi0  ->Fill(pt);\r
 506                                 fhPhiMCNoPi0 ->Fill(pt,phi);\r
 507                                 fhEtaMCNoPi0 ->Fill(pt,eta);\r
 508                         }\r
 509                 }//Histograms with MC\r
 510                 \r
 511         }// aod loop\r
 512         \r
 513 }\r
 514 \r
 515 \r
 516 //____________________________________________________________________________\r
 517 void AliAnaPi0EbE::InitParameters()\r
 518 {\r
 519         \r
 520         //Initialize the parameters of the analysis.\r
 521         SetOutputAODClassName("AliAODPWG4Particle");\r
 522         SetOutputAODName("pi0s");\r
 523         fInputAODGammaConvName = "gammaconv" ;\r
 524         fAnaType = kIMCalo ;\r
 525         fCalorimeter = "PHOS" ;\r
 526         fMinDist  = 2.;\r
 527         fMinDist2 = 4.;\r
 528         fMinDist3 = 5.;\r
 529         \r
 530 }\r
 531 \r
 532 //__________________________________________________________________\r
 533 void AliAnaPi0EbE::Print(const Option_t * opt) const\r
 534 {\r
 535         \r
 536         //Print some relevant parameters set for the analysis\r
 537         if(! opt)\r
 538                 return;\r
 539         \r
 540         printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
 541         AliAnaPartCorrBaseClass::Print("");\r
 542         printf("Analysis Type = %d \n",  fAnaType) ;\r
 543         if(fAnaType == kSSCalo){     \r
 544                 printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;\r
 545                 printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);\r
 546                 printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);\r
 547                 printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3); \r
 548         } \r
 549         printf("    \n") ;\r
 550         \r
 551 } \r