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