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