* Contributors are mentioned in the code where appropriate. *
* *
* Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes iGetEntriesFast(s hereby granted *
+ * documentation strictly for non-commercial purposes is hereby granted *
* without fee, provided that the above copyright notice appears in all *
* copies and that both the copyright notice and this permission notice *
* appear in the supporting documentation. The authors make no claims *
// Pi0 identified by one of the following:
// -Invariant mass of 2 cluster in calorimeter
// -Shower shape analysis in calorimeter
-// -Invariant mass of one cluster in calorimeter and one photon reconstructed in TPC (in near future)
+// -Invariant mass of one cluster in calorimeter and one photon reconstructed in CTS
//
// -- Author: Gustavo Conesa (LNF-INFN) & Raphaelle Ichou (SUBATECH)
//////////////////////////////////////////////////////////////////////////////
#include <TList.h>
#include <TClonesArray.h>
#include <TObjString.h>
-#include <TH2F.h>
+#include <TH3F.h>
//#include "Riostream.h"
// --- Analysis system ---
#include "AliNeutralMesonSelection.h"
#include "AliCaloPID.h"
#include "AliMCAnalysisUtils.h"
-#include "AliAODPWG4ParticleCorrelation.h"
#include "AliStack.h"
-#include "AliFidutialCut.h"
+#include "AliFiducialCut.h"
#include "TParticle.h"
-#include "AliAODCaloCluster.h"
+#include "AliVCluster.h"
#include "AliAODEvent.h"
+#include "AliAODMCParticle.h"
ClassImp(AliAnaPi0EbE)
-//____________________________________________________________________________
+//____________________________
AliAnaPi0EbE::AliAnaPi0EbE() :
-AliAnaPartCorrBaseClass(), fAnaType(kIMCalo),fCalorimeter(""),
-fMinDist(0.),fMinDist2(0.),fMinDist3(0.),
-fInputAODGammaConv(0x0),fInputAODGammaConvName(""),
-fhPtPi0(0),fhPhiPi0(0),fhEtaPi0(0),
-fhPtMCNoPi0(0),fhPhiMCNoPi0(0),fhEtaMCNoPi0(0),
-fhPtMCPi0(0),fhPhiMCPi0(0),fhEtaMCPi0(0)
+ AliAnaPartCorrBaseClass(), fAnaType(kIMCalo), fCalorimeter(""),
+ fMinDist(0.),fMinDist2(0.), fMinDist3(0.), fFillWeightHistograms(kFALSE),
+ fInputAODGammaConvName(""),
+ //Histograms
+ fhPtPi0(0), fhEPi0(0), fhEEtaPhiPi0(0),
+ //Shower shape histos
+ fhEDispersion(0), fhELambda0(0), fhELambda1(0),
+ fhELambda0NoTRD(0), fhELambda0FracMaxCellCut(0),
+ fhEFracMaxCell(0), fhEFracMaxCellNoTRD(0),
+ fhENCells(0), fhETime(0), fhEPairDiffTime(0),
+ //MC histos
+ fhPtMCNoPi0(0), fhPhiMCNoPi0(0), fhEtaMCNoPi0(0),
+ fhPtMCPi0(0), fhPhiMCPi0(0), fhEtaMCPi0(0),
+ // Weight studies
+ fhECellClusterRatio(0), fhECellClusterLogRatio(0),
+ fhEMaxCellClusterRatio(0), fhEMaxCellClusterLogRatio(0)
{
//default ctor
+ for(Int_t i = 0; i < 6; i++){
+ fhEMCLambda0[i] = 0;
+ fhEMCLambda0NoTRD[i]= 0;
+ fhEMCLambda0FracMaxCellCut[i]= 0;
+ fhEMCFracMaxCell[i] = 0;
+ fhEMCLambda1[i] = 0;
+ fhEMCDispersion[i] = 0;
+ }
+
+ //Weight studies
+ for(Int_t i =0; i < 14; i++){
+ fhLambda0ForW0[i] = 0;
+ //fhLambda1ForW0[i] = 0;
+ }
+
//Initialize parameters
InitParameters();
}
-//____________________________________________________________________________
-AliAnaPi0EbE::AliAnaPi0EbE(const AliAnaPi0EbE & p) :
-AliAnaPartCorrBaseClass(p), fAnaType(p.fAnaType), fCalorimeter(p.fCalorimeter),
-fMinDist(p.fMinDist),fMinDist2(p.fMinDist2), fMinDist3(p.fMinDist3),
-fInputAODGammaConv(new TClonesArray(*p.fInputAODGammaConv)),
-fInputAODGammaConvName(p.fInputAODGammaConvName),
-fhPtPi0(p.fhPtPi0),fhPhiPi0(p.fhPhiPi0),fhEtaPi0(p.fhEtaPi0),
-fhPtMCNoPi0(p.fhPtMCNoPi0),fhPhiMCNoPi0(p.fhPhiMCNoPi0),fhEtaMCNoPi0(p.fhEtaMCNoPi0),
-fhPtMCPi0(p.fhPtMCPi0),fhPhiMCPi0(p.fhPhiMCPi0),fhEtaMCPi0(p.fhEtaMCPi0)
-{
- // cpy ctor
+//_____________________________________________________________________________________
+void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, const Int_t tag){
+
+ // Fill shower shape, timing and other histograms for selected clusters from decay
+
+ Float_t e = cluster->E();
+ Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
+ Float_t l0 = cluster->GetM02();
+ Float_t l1 = cluster->GetM20();
+ Int_t nSM = GetModuleNumber(cluster);
+
+ AliVCaloCells * cell = 0x0;
+ if(fCalorimeter == "PHOS")
+ cell = GetPHOSCells();
+ else
+ cell = GetEMCALCells();
+
+ Float_t maxCellFraction = 0;
+ GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
+ fhEFracMaxCell->Fill(e,maxCellFraction);
+
+ FillWeightHistograms(cluster);
+
+ fhEDispersion->Fill(e, disp);
+ fhELambda0 ->Fill(e, l0 );
+ fhELambda1 ->Fill(e, l1 );
+
+ if(fCalorimeter=="EMCAL" && nSM < 6) {
+ fhELambda0NoTRD->Fill(e, l0 );
+ fhEFracMaxCellNoTRD->Fill(e,maxCellFraction);
+ }
+
+ if(maxCellFraction < 0.5)
+ fhELambda0FracMaxCellCut->Fill(e, l0 );
+
+ fhETime ->Fill(e, cluster->GetTOF()*1.e9);
+ fhENCells->Fill(e, cluster->GetNCells());
+
+ if(IsDataMC()) {
+ //Photon1
+ if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) ){
+ fhEMCLambda0[kmcPi0] ->Fill(e, l0);
+ fhEMCLambda1[kmcPi0] ->Fill(e, l1);
+ fhEMCDispersion[kmcPi0] ->Fill(e, disp);
+
+ fhEMCFracMaxCell[kmcPi0]->Fill(e,maxCellFraction);
+ if(fCalorimeter=="EMCAL" && nSM < 6)
+ fhEMCLambda0NoTRD[kmcPi0]->Fill(e, l0 );
+ if(maxCellFraction < 0.5)
+ fhEMCLambda0FracMaxCellCut[kmcPi0]->Fill(e, l0 );
+
+ }//pi0
+ else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) ){
+ fhEMCLambda0[kmcEta] ->Fill(e, l0);
+ fhEMCLambda1[kmcEta] ->Fill(e, l1);
+ fhEMCDispersion[kmcEta] ->Fill(e, disp);
+ fhEMCFracMaxCell[kmcEta]->Fill(e,maxCellFraction);
+ if(fCalorimeter=="EMCAL" && nSM < 6)
+ fhEMCLambda0NoTRD[kmcEta]->Fill(e, l0 );
+ if(maxCellFraction < 0.5)
+ fhEMCLambda0FracMaxCellCut[kmcEta]->Fill(e, l0 );
+ }//eta
+ else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
+ GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) ){
+ fhEMCLambda0[kmcConversion] ->Fill(e, l0);
+ fhEMCLambda1[kmcConversion] ->Fill(e, l1);
+ fhEMCDispersion[kmcConversion] ->Fill(e, disp);
+ fhEMCFracMaxCell[kmcConversion]->Fill(e,maxCellFraction);
+ if(fCalorimeter=="EMCAL" && nSM < 6)
+ fhEMCLambda0NoTRD[kmcConversion]->Fill(e, l0 );
+ if(maxCellFraction < 0.5)
+ fhEMCLambda0FracMaxCellCut[kmcConversion]->Fill(e, l0 );
+ }//conversion photon
+ else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) ){
+ fhEMCLambda0[kmcPhoton] ->Fill(e, l0);
+ fhEMCLambda1[kmcPhoton] ->Fill(e, l1);
+ fhEMCDispersion[kmcPhoton] ->Fill(e, disp);
+ fhEMCFracMaxCell[kmcPhoton]->Fill(e,maxCellFraction);
+ if(fCalorimeter=="EMCAL" && nSM < 6)
+ fhEMCLambda0NoTRD[kmcPhoton]->Fill(e, l0 );
+ if(maxCellFraction < 0.5)
+ fhEMCLambda0FracMaxCellCut[kmcPhoton]->Fill(e, l0 );
+ }//photon no conversion
+ else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)){
+ fhEMCLambda0[kmcElectron] ->Fill(e, l0);
+ fhEMCLambda1[kmcElectron] ->Fill(e, l1);
+ fhEMCDispersion[kmcElectron] ->Fill(e, disp);
+ fhEMCFracMaxCell[kmcElectron]->Fill(e,maxCellFraction);
+ if(fCalorimeter=="EMCAL" && nSM < 6)
+ fhEMCLambda0NoTRD[kmcElectron]->Fill(e, l0 );
+ if(maxCellFraction < 0.5)
+ fhEMCLambda0FracMaxCellCut[kmcElectron]->Fill(e, l0 );
+ }//electron
+ else {
+ fhEMCLambda0[kmcHadron] ->Fill(e, l0);
+ fhEMCLambda1[kmcHadron] ->Fill(e, l1);
+ fhEMCDispersion[kmcHadron] ->Fill(e, disp);
+ fhEMCFracMaxCell[kmcHadron]->Fill(e,maxCellFraction);
+ if(fCalorimeter=="EMCAL" && nSM < 6)
+ fhEMCLambda0NoTRD[kmcHadron]->Fill(e, l0 );
+ if(maxCellFraction < 0.5)
+ fhEMCLambda0FracMaxCellCut[kmcHadron]->Fill(e, l0 );
+ }//other particles
+ }//MC
}
-//_________________________________________________________________________
-AliAnaPi0EbE & AliAnaPi0EbE::operator = (const AliAnaPi0EbE & p)
+//________________________________________________________
+void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
{
- // assignment operator
+ // Calculate weights and fill histograms
+
+ if(!fFillWeightHistograms || GetMixedEvent()) return;
- if(&p == this) return *this;
+ AliVCaloCells* cells = 0;
+ if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
+ else cells = GetPHOSCells();
- fAnaType = p.fAnaType ;
- fCalorimeter = p.fCalorimeter ;
- fMinDist = p.fMinDist;
- fMinDist2 = p.fMinDist2;
- fMinDist3 = p.fMinDist3;
+ // First recalculate energy in case non linearity was applied
+ Float_t energy = 0;
+ Float_t ampMax = 0;
+ for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
+
+ Int_t id = clus->GetCellsAbsId()[ipos];
+
+ //Recalibrate cell energy if needed
+ Float_t amp = cells->GetCellAmplitude(id);
+ RecalibrateCellAmplitude(amp,id);
+
+ energy += amp;
+
+ if(amp> ampMax)
+ ampMax = amp;
+
+ } // energy loop
- fInputAODGammaConv = new TClonesArray(*p.fInputAODGammaConv) ;
- fInputAODGammaConvName = p.fInputAODGammaConvName ;
+ if(energy <=0 ) {
+ printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
+ return;
+ }
- fhPtPi0 = p.fhPtPi0; fhPhiPi0 = p.fhPhiPi0; fhEtaPi0 = p.fhEtaPi0;
- fhPtMCNoPi0 = p.fhPtMCNoPi0; fhPhiMCNoPi0 = p.fhPhiMCNoPi0; fhEtaMCNoPi0 = p.fhEtaMCPi0;
- fhPtMCPi0 = p.fhPtMCPi0; fhPhiMCPi0 = p.fhPhiMCPi0; fhEtaMCPi0 = p.fhEtaMCPi0;
+ fhEMaxCellClusterRatio ->Fill(energy,ampMax/energy);
+ fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
- return *this;
+ //Get the ratio and log ratio to all cells in cluster
+ for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
+ Int_t id = clus->GetCellsAbsId()[ipos];
+
+ //Recalibrate cell energy if needed
+ Float_t amp = cells->GetCellAmplitude(id);
+ RecalibrateCellAmplitude(amp,id);
+
+ fhECellClusterRatio ->Fill(energy,amp/energy);
+ fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
+ }
+ //Recalculate shower shape for different W0
+ if(fCalorimeter=="EMCAL"){
+
+ Float_t l0org = clus->GetM02();
+ Float_t l1org = clus->GetM20();
+ Float_t dorg = clus->GetDispersion();
+
+ for(Int_t iw = 0; iw < 14; iw++){
+ GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
+ GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
+
+ fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
+ //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
+
+ } // w0 loop
+
+ // Set the original values back
+ clus->SetM02(l0org);
+ clus->SetM20(l1org);
+ clus->SetDispersion(dorg);
+
+ }// EMCAL
}
-//____________________________________________________________________________
-AliAnaPi0EbE::~AliAnaPi0EbE()
-{
- //dtor
- if(fInputAODGammaConv){
- fInputAODGammaConv->Clear() ;
- delete fInputAODGammaConv ;
+//___________________________________________
+TObjString * AliAnaPi0EbE::GetAnalysisCuts()
+{
+ //Save parameters used for analysis
+ TString parList ; //this will be list of parameters used for this analysis.
+ const Int_t buffersize = 255;
+ char onePar[buffersize] ;
+
+ snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
+ parList+=onePar ;
+
+ if(fAnaType == kSSCalo){
+ snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
+ parList+=onePar ;
+ snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
+ parList+=onePar ;
}
+
+ //Get parameters set in base class.
+ parList += GetBaseParametersList() ;
+
+ //Get parameters set in PID class.
+ if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
+
+ return new TObjString(parList) ;
}
-//________________________________________________________________________
+//_____________________________________________
TList * AliAnaPi0EbE::GetCreateOutputObjects()
{
// Create histograms to be saved in output file and
TList * outputContainer = new TList() ;
outputContainer->SetName("Pi0EbEHistos") ;
- Int_t nptbins = GetHistoNPtBins();
- Int_t nphibins = GetHistoNPhiBins();
- Int_t netabins = GetHistoNEtaBins();
- Float_t ptmax = GetHistoPtMax();
- Float_t phimax = GetHistoPhiMax();
- Float_t etamax = GetHistoEtaMax();
- Float_t ptmin = GetHistoPtMin();
- Float_t phimin = GetHistoPhiMin();
- Float_t etamin = GetHistoEtaMin();
-
+ Int_t nptbins = GetHistoPtBins(); Float_t ptmax = GetHistoPtMax(); Float_t ptmin = GetHistoPtMin();
+ Int_t nphibins = GetHistoPhiBins(); Float_t phimax = GetHistoPhiMax(); Float_t phimin = GetHistoPhiMin();
+ Int_t netabins = GetHistoEtaBins(); Float_t etamax = GetHistoEtaMax(); Float_t etamin = GetHistoEtaMin();
+ Int_t ssbins = GetHistoShowerShapeBins(); Float_t ssmax = GetHistoShowerShapeMax(); Float_t ssmin = GetHistoShowerShapeMin();
+ Int_t tdbins = GetHistoDiffTimeBins() ; Float_t tdmax = GetHistoDiffTimeMax(); Float_t tdmin = GetHistoDiffTimeMin();
+ Int_t tbins = GetHistoTimeBins() ; Float_t tmax = GetHistoTimeMax(); Float_t tmin = GetHistoTimeMin();
+ Int_t nbins = GetHistoNClusterCellBins(); Int_t nmax = GetHistoNClusterCellMax(); Int_t nmin = GetHistoNClusterCellMin();
+
fhPtPi0 = new TH1F("hPtPi0","Number of identified #pi^{0} decay",nptbins,ptmin,ptmax);
fhPtPi0->SetYTitle("N");
fhPtPi0->SetXTitle("p_{T #pi^{0}}(GeV/c)");
outputContainer->Add(fhPtPi0) ;
- fhPhiPi0 = new TH2F
- ("hPhiPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiPi0->SetYTitle("#phi");
- fhPhiPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
- outputContainer->Add(fhPhiPi0) ;
+ fhEPi0 = new TH1F("hEPi0","Number of identified #pi^{0} decay",nptbins,ptmin,ptmax);
+ fhEPi0->SetYTitle("N");
+ fhEPi0->SetXTitle("E #pi^{0}(GeV)");
+ outputContainer->Add(fhEPi0) ;
+
+ fhEEtaPhiPi0 = new TH3F
+ ("hEEtaPhiPi0","Selected #pi^{0} pairs: E vs #eta vs #phi",nptbins,ptmin,ptmax,netabins,etamin,etamax, nphibins,phimin,phimax);
+ fhEEtaPhiPi0->SetZTitle("#phi");
+ fhEEtaPhiPi0->SetYTitle("#eta");
+ fhEEtaPhiPi0->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEEtaPhiPi0) ;
+
+ ////////
+
+ if(fAnaType == kIMCalo || fAnaType == kIMCaloTracks ){
+
+ fhEDispersion = new TH2F
+ ("hEDispersion","Selected #pi^{0} pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEDispersion->SetYTitle("D^{2}");
+ fhEDispersion->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEDispersion) ;
+
+ fhELambda0 = new TH2F
+ ("hELambda0","Selected #pi^{0} pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhELambda0->SetYTitle("#lambda_{0}^{2}");
+ fhELambda0->SetXTitle("E (GeV)");
+ outputContainer->Add(fhELambda0) ;
+
+ fhELambda1 = new TH2F
+ ("hELambda1","Selected #pi^{0} pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhELambda1->SetYTitle("#lambda_{1}^{2}");
+ fhELambda1->SetXTitle("E (GeV)");
+ outputContainer->Add(fhELambda1) ;
+
+ fhELambda0FracMaxCellCut = new TH2F
+ ("hELambda0FracMaxCellCut","Selected #pi^{0} pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
+ fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
+ outputContainer->Add(fhELambda0FracMaxCellCut) ;
+
+ fhEFracMaxCell = new TH2F
+ ("hEFracMaxCell","Selected #pi^{0} pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
+ fhEFracMaxCell->SetYTitle("Fraction");
+ fhEFracMaxCell->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEFracMaxCell) ;
+
+ if(fCalorimeter=="EMCAL"){
+ fhELambda0NoTRD = new TH2F
+ ("hELambda0NoTRD","Selected #pi^{0} pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
+ fhELambda0NoTRD->SetXTitle("E (GeV)");
+ outputContainer->Add(fhELambda0NoTRD) ;
+
+ fhEFracMaxCellNoTRD = new TH2F
+ ("hEFracMaxCellNoTRD","Selected #pi^{0} pairs: E vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
+ fhEFracMaxCellNoTRD->SetYTitle("Fraction");
+ fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEFracMaxCellNoTRD) ;
+ }
+
+ fhENCells = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
+ fhENCells->SetXTitle("E (GeV)");
+ fhENCells->SetYTitle("# of cells in cluster");
+ outputContainer->Add(fhENCells);
+
+ fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
+ fhETime->SetXTitle("E (GeV)");
+ fhETime->SetYTitle(" t (ns)");
+ outputContainer->Add(fhETime);
+
+ fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
+ fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
+ fhEPairDiffTime->SetYTitle("#Delta t (ns)");
+ outputContainer->Add(fhEPairDiffTime);
+
+
+ }// Invariant mass analysis in calorimeters only
- fhEtaPi0 = new TH2F
- ("hEtaPi0","#phi_{#pi^{0}}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaPi0->SetYTitle("#eta");
- fhEtaPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
- outputContainer->Add(fhEtaPi0) ;
+ if(fFillWeightHistograms){
+
+ fhECellClusterRatio = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
+ nptbins,ptmin,ptmax, 100,0,1.);
+ fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
+ fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
+ outputContainer->Add(fhECellClusterRatio);
+
+ fhECellClusterLogRatio = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
+ nptbins,ptmin,ptmax, 100,-10,0);
+ fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
+ fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
+ outputContainer->Add(fhECellClusterLogRatio);
+
+ fhEMaxCellClusterRatio = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
+ nptbins,ptmin,ptmax, 100,0,1.);
+ fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
+ fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
+ outputContainer->Add(fhEMaxCellClusterRatio);
+
+ fhEMaxCellClusterLogRatio = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
+ nptbins,ptmin,ptmax, 100,-10,0);
+ fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
+ fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
+ outputContainer->Add(fhEMaxCellClusterLogRatio);
+
+ for(Int_t iw = 0; iw < 14; iw++){
+ fhLambda0ForW0[iw] = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected decay photons from neutral meson",1+0.5*iw),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
+ fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
+ outputContainer->Add(fhLambda0ForW0[iw]);
+
+// fhLambda1ForW0[iw] = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected decay photons from neutral meson",0.5+0.5*iw),
+// nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+// fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
+// fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
+// outputContainer->Add(fhLambda1ForW0[iw]);
+
+ }
+ }
if(IsDataMC()) {
if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
outputContainer->Add(fhPtMCPi0) ;
fhPhiMCPi0 = new TH2F
- ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
+ ("hPhiMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
fhPhiMCPi0->SetYTitle("#phi");
fhPhiMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
outputContainer->Add(fhPhiMCPi0) ;
fhEtaMCPi0 = new TH2F
- ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
+ ("hEtaMCPi0","Identified pi0 from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
fhEtaMCPi0->SetYTitle("#eta");
fhEtaMCPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
outputContainer->Add(fhEtaMCPi0) ;
outputContainer->Add(fhPtMCNoPi0) ;
fhPhiMCNoPi0 = new TH2F
- ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
+ ("hPhiMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
fhPhiMCNoPi0->SetYTitle("#phi");
fhPhiMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
outputContainer->Add(fhPhiMCNoPi0) ;
fhEtaMCNoPi0 = new TH2F
- ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
+ ("hEtaMCNoPi0","Identified pi0 not from pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax);
fhEtaMCNoPi0->SetYTitle("#eta");
fhEtaMCNoPi0->SetXTitle("p_{T #pi^{0}} (GeV/c)");
outputContainer->Add(fhEtaMCNoPi0) ;
- }
+ if(fAnaType == kIMCalo){
+ TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
+ TString pname[] ={"Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
+ for(Int_t i = 0; i < 6; i++){
+
+ fhEMCLambda0[i] = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
+ Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
+ fhEMCLambda0[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEMCLambda0[i]) ;
+
+ if(fCalorimeter=="EMCAL"){
+ fhEMCLambda0NoTRD[i] = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
+ Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
+ fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
+ }
+
+ fhEMCLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
+ Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
+ fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
+
+ fhEMCFracMaxCell[i] = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
+ Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
+ nptbins,ptmin,ptmax,100,0,1);
+ fhEMCFracMaxCell[i]->SetYTitle("Fraction");
+ fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEMCFracMaxCell[i]) ;
+
+ fhEMCLambda1[i] = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
+ Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
+ fhEMCLambda1[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEMCLambda1[i]) ;
+
+ fhEMCDispersion[i] = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
+ Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
+ nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+ fhEMCDispersion[i]->SetYTitle("D^{2}");
+ fhEMCDispersion[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEMCDispersion[i]) ;
+
+ }//
+
+ }//kIMCalo
+ } //Not MC reader
}//Histos with MC
TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
+ delete nmsHistos;
+
}
- //Save parameters used for analysis
- TString parList ; //this will be list of parameters used for this analysis.
- char onePar[255] ;
-
- sprintf(onePar,"--- AliAnaPi0EbE ---\n") ;
- parList+=onePar ;
- sprintf(onePar,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
- parList+=onePar ;
+ return outputContainer ;
- if(fAnaType == kSSCalo){
- sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
- parList+=onePar ;
- sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
- parList+=onePar ;
- sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
- parList+=onePar ;
- sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
- parList+=onePar ;
+}
+
+//____________________________________________________________________________
+void AliAnaPi0EbE::Init()
+{
+ //Init
+ //Do some checks
+ if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
+ printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
+ abort();
+ }
+ else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
+ printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
+ abort();
}
- //Get parameters set in base class.
- parList += GetBaseParametersList() ;
-
- //Get parameters set in PID class.
- if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
-
- TObjString *oString= new TObjString(parList) ;
- outputContainer->Add(oString);
+}
+
+//____________________________________________________________________________
+void AliAnaPi0EbE::InitParameters()
+{
+ //Initialize the parameters of the analysis.
+ AddToHistogramsName("AnaPi0EbE_");
- return outputContainer ;
+ fInputAODGammaConvName = "PhotonsCTS" ;
+ fAnaType = kIMCalo ;
+ fCalorimeter = "EMCAL" ;
+ fMinDist = 2.;
+ fMinDist2 = 4.;
+ fMinDist3 = 5.;
}
//Do analysis and fill aods
switch(fAnaType)
- {
+ {
case kIMCalo:
MakeInvMassInCalorimeter();
break;
MakeInvMassInCalorimeterAndCTS();
break;
- }
+ }
}
-//__________________________________________________________________
+//____________________________________________
void AliAnaPi0EbE::MakeInvMassInCalorimeter()
{
//Do analysis and fill aods
TLorentzVector mom1;
TLorentzVector mom2;
TLorentzVector mom ;
- Int_t tag1 =-1;
- Int_t tag2 =-1;
- Int_t tag = -1;
+ Int_t tag1 = 0;
+ Int_t tag2 = 0;
+ Int_t tag = 0;
if(!GetInputAODBranch()){
- printf("AliAnaPi0EbE::kIMCalo:: No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
+ printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
abort();
}
- for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
+
+ //Get shower shape information of clusters
+ TObjArray *clusters = 0;
+ if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
+ else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
+
+ for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
+
+ //Vertex cut in case of mixed events
+ Int_t evtIndex1 = 0 ;
+ if(GetMixedEvent())
+ evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
+ if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ; //vertex cut
mom1 = *(photon1->Momentum());
+ //Get original cluster, to recover some information
+ Int_t iclus = -1;
+ AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
- for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast()-1; jphoton++){
+ if(!cluster1){
+ printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
+ return;
+ }
+
+ for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++){
- AliAODPWG4ParticleCorrelation * photon2 = (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(jphoton));
+ AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
+ Int_t evtIndex2 = 0 ;
+ if(GetMixedEvent())
+ evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
+ if(GetMixedEvent() && (evtIndex1 == evtIndex2))
+ continue ;
+ if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ; //vertex cut
mom2 = *(photon2->Momentum());
+ //Get original cluster, to recover some information
+ Int_t iclus2;
+ AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
+
+ if(!cluster2){
+ printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
+ return;
+ }
+
+ Float_t e1 = photon1->E();
+ Float_t e2 = photon2->E();
+
+ //Select clusters with good time window difference
+ Float_t tof1 = cluster1->GetTOF()*1e9;;
+ Float_t tof2 = cluster2->GetTOF()*1e9;;
+ Double_t t12diff = tof1-tof2;
+ fhEPairDiffTime->Fill(e1+e2, t12diff);
+ if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
+
//Select good pair (good phi, pt cuts, aperture and invariant mass)
- if(GetNeutralMesonSelection()->SelectPair(mom1, mom2))
- {
- if(GetDebug()>1)
- 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());
-
- //Play with the MC stack if available
- if(IsDataMC()){
- //Check origin of the candidates
- tag1 = GetMCAnalysisUtils()->CheckOrigin(photon1->GetLabel(), GetMCStack());
- tag2 = GetMCAnalysisUtils()->CheckOrigin(photon2->GetLabel(), GetMCStack());
-
- if(GetDebug() > 0) printf("AliAnaPi0EbE::kIMCalo::Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
- if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){
-
- //Check if pi0 mother is the same
- Int_t label1 = photon1->GetLabel();
- TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
- label1 = mother1->GetFirstMother();
- //mother1 = GetMCStack()->Particle(label1);//pi0
-
- Int_t label2 = photon2->GetLabel();
- TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
- label2 = mother2->GetFirstMother();
- //mother2 = GetMCStack()->Particle(label2);//pi0
-
- //printf("mother1 %d, mother2 %d\n",label1,label2);
- if(label1 == label2)
- tag = AliMCAnalysisUtils::kMCPi0;
- }
- }//Work with stack also
-
- //Create AOD for analysis
- mom = mom1+mom2;
- AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
- //pi0.SetLabel(calo->GetLabel(0));
- pi0.SetPdg(AliCaloPID::kPi0);
- pi0.SetDetector(photon1->GetDetector());
- pi0.SetTag(tag);
- //Set the indeces of the original caloclusters
- pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
- AddAODParticle(pi0);
- }//pi0
+ if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
+ {
+ if(GetDebug()>1)
+ printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
+
+ //Play with the MC stack if available
+ if(IsDataMC()){
+ //Check origin of the candidates
+ Int_t label1 = photon1->GetLabel();
+ Int_t label2 = photon2->GetLabel();
+ if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
+ if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
+
+ if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
+ if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
+ GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
+
+ //Check if pi0 mother is the same
+ if(GetReader()->ReadStack()){
+ if(label1>=0){
+ TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
+ label1 = mother1->GetFirstMother();
+ //mother1 = GetMCStack()->Particle(label1);//pi0
+ }
+ if(label2>=0){
+ TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
+ label2 = mother2->GetFirstMother();
+ //mother2 = GetMCStack()->Particle(label2);//pi0
+ }
+ }
+ else if(GetReader()->ReadAODMCParticles()){//&& (input > -1)){
+ if(label1>=0){
+ AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
+ label1 = mother1->GetMother();
+ //mother1 = GetMCStack()->Particle(label1);//pi0
+ }
+ if(label2>=0){
+ AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
+ label2 = mother2->GetMother();
+ //mother2 = GetMCStack()->Particle(label2);//pi0
+ }
+ }
+
+ //printf("mother1 %d, mother2 %d\n",label1,label2);
+ if(label1 == label2 && label1>=0)
+ GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
+ }
+ }//Work with stack also
+
+
+ //Fill some histograms about shower shape
+ if(clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
+ FillSelectedClusterHistograms(cluster1, tag1);
+ FillSelectedClusterHistograms(cluster2, tag2);
+ }
+
+ // Tag both photons as decay
+ photon1->SetTagged(kTRUE);
+ photon2->SetTagged(kTRUE);
+
+ //Create AOD for analysis
+ mom = mom1+mom2;
+ AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
+ //pi0.SetLabel(calo->GetLabel());
+ pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
+ pi0.SetDetector(photon1->GetDetector());
+ pi0.SetTag(tag);
+ //Set the indeces of the original caloclusters
+ pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
+ //pi0.SetInputFileIndex(input);
+ AddAODParticle(pi0);
+ }//pi0
+
}//2n photon loop
}//1st photon loop
- if(GetDebug() > 1) printf("AliAnaPi0EbE::kIMCalo::End fill AODs \n");
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
}
TLorentzVector mom1;
TLorentzVector mom2;
TLorentzVector mom ;
- Int_t tag1 =-1;
- Int_t tag2 =-1;
- Int_t tag = -1;
+ Int_t tag1 = 0;
+ Int_t tag2 = 0;
+ Int_t tag = 0;
+ Int_t evtIndex = 0;
+ // Check calorimeter input
if(!GetInputAODBranch()){
- printf("AliAnaPi0EbE::kIMCaloCTS: No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
+ printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
abort();
}
+
+ // Get the array with conversion photons
+ TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
+ if(!inputAODGammaConv) {
+
+ inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
+
+ if(!inputAODGammaConv) {
+ printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
+
+ return;
+ }
+ }
+
+ //Get shower shape information of clusters
+ TObjArray *clusters = 0;
+ if (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
+ else if(fCalorimeter=="PHOS") clusters = GetPHOSClusters() ;
+
+ Int_t nCTS = inputAODGammaConv->GetEntriesFast();
+ Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
+ if(nCTS<=0 || nCalo <=0) {
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
+ return;
+ }
+
+ if(GetDebug() > 1)
+ printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
+
+ // Do the loop, first calo, second CTS
for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
AliAODPWG4Particle * photon1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
mom1 = *(photon1->Momentum());
- //Play with the MC stack if available
- fInputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
- if(!fInputAODGammaConv) {
- printf("AliAnaPi0EbE::kIMCaloCTS: No input gamma conversions in AOD branch with name < %s >, STOP",fInputAODGammaConvName.Data());
- abort();
- }
- for(Int_t jphoton = iphoton+1; jphoton < fInputAODGammaConv->GetEntriesFast()-1; jphoton++){
- AliAODPWG4ParticleCorrelation * photon2 = (AliAODPWG4ParticleCorrelation*) (fInputAODGammaConv->At(jphoton));
+ //Get original cluster, to recover some information
+ Int_t iclus = -1;
+ AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
+
+ for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
+ AliAODPWG4Particle * photon2 = (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
+ if(GetMixedEvent())
+ evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
+ if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
+
mom2 = *(photon2->Momentum());
+
//Select good pair (good phi, pt cuts, aperture and invariant mass)
- if(GetNeutralMesonSelection()->SelectPair(mom1, mom2)){
- 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());
-
- if(IsDataMC()){
- //Check origin of the candidates
- tag1 = GetMCAnalysisUtils()->CheckOrigin(photon1->GetLabel(), GetMCStack());
- tag2 = GetMCAnalysisUtils()->CheckOrigin(photon2->GetLabel(), GetMCStack());
- if(GetDebug() > 0) printf("AliAnaPi0EbE::kIMCaloCTS::Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
- if(tag1 == AliMCAnalysisUtils::kMCPi0Decay && tag2 == AliMCAnalysisUtils::kMCPi0Decay){
- //Check if pi0 mother is the same
- Int_t label1 = photon1->GetLabel();
- TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
- label1 = mother1->GetFirstMother();
- //mother1 = GetMCStack()->Particle(label1);//pi0
-
- Int_t label2 = photon2->GetLabel();
- TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
- label2 = mother2->GetFirstMother();
- //mother2 = GetMCStack()->Particle(label2);//pi0
-
- //printf("mother1 %d, mother2 %d\n",label1,label2);
- if(label1 == label2)
- tag = AliMCAnalysisUtils::kMCPi0;
- }
- }//Work with stack also
-
- //Create AOD for analysis
- mom = mom1+mom2;
- AliAODPWG4ParticleCorrelation pi0 = AliAODPWG4ParticleCorrelation(mom);
- //pi0.SetLabel(calo->GetLabel(0));
- pi0.SetPdg(AliCaloPID::kPi0);
- pi0.SetDetector(photon1->GetDetector());
- pi0.SetTag(tag);
- //Set the indeces of the original tracks or caloclusters
- pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
- pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
- AddAODParticle(pi0);
+ if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter)){
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
+
+ if(IsDataMC()){
+ Int_t label1 = photon1->GetLabel();
+ Int_t label2 = photon2->GetLabel();
+ if(label1>=0)tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
+ if(label2>=0)tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
+ if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
+ if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
+ GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
+ //Check if pi0 mother is the same
+
+ if(GetReader()->ReadStack()){
+ if(label1>=0){
+ TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
+ label1 = mother1->GetFirstMother();
+ //mother1 = GetMCStack()->Particle(label1);//pi0
+ }
+ if(label2>=0){
+ TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
+ label2 = mother2->GetFirstMother();
+ //mother2 = GetMCStack()->Particle(label2);//pi0
+ }
+ }
+ else if(GetReader()->ReadAODMCParticles()&& label1>=0 && label2>=0){ //&& (input > -1)){
+ if(label1>=0){
+ AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
+ label1 = mother1->GetMother();
+ //mother1 = GetMCStack()->Particle(label1);//pi0
+ }
+ if(label2>=0){
+ AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
+ label2 = mother2->GetMother();
+ //mother2 = GetMCStack()->Particle(label2);//pi0
+ }
+ }
+
+ //printf("mother1 %d, mother2 %d\n",label1,label2);
+ if(label1 == label2 && label1>=0)
+ GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
+ }
+ }//Work with stack also
+
+ //Fill some histograms about shower shape
+ if(cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
+ FillSelectedClusterHistograms(cluster, tag1);
+ }
+
+ // Tag both photons as decay
+ photon1->SetTagged(kTRUE);
+ photon2->SetTagged(kTRUE);
+
+ //Create AOD for analysis
+ mom = mom1+mom2;
+ AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
+ //pi0.SetLabel(calo->GetLabel());
+ pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
+ pi0.SetDetector(photon1->GetDetector());
+ pi0.SetTag(tag);
+ //Set the indeces of the original tracks or caloclusters
+ pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
+ pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
+ //pi0.SetInputFileIndex(input);
+ AddAODParticle(pi0);
}//pi0
}//2n photon loop
}//1st photon loop
- if(GetDebug() > 1) printf("AliAnaPi0EbE::kIMCaloCTS::End fill AODs \n");
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
}
{
//Search for pi0 in fCalorimeter with shower shape analysis
- TRefArray * pl = new TRefArray;
-
- //Get vertex for photon momentum calculation
- Double_t vertex[]={0,0,0} ; //vertex ;
- if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
-
+ TObjArray * pl = 0x0;
//Select the Calorimeter of the photon
if(fCalorimeter == "PHOS")
- pl = GetAODPHOS();
+ pl = GetPHOSClusters();
else if (fCalorimeter == "EMCAL")
- pl = GetAODEMCAL();
- //Fill AODCaloClusters and AODParticle with PHOS aods
+ pl = GetEMCALClusters();
+
+ if(!pl) {
+ Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
+ return;
+ }
+
TLorentzVector mom ;
for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
- AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));
+ AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
+
+ Int_t evtIndex = 0 ;
+ if (GetMixedEvent()) {
+ evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
+ }
+ if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
- //Cluster selection, not charged, with pi0 id and in fidutial cut
//Get Momentum vector,
- calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
+ if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
+ calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
+ else{
+ Double_t vertex[]={0,0,0};
+ calo->GetMomentum(mom,vertex) ;
+ }
+
//If too small or big pt, skip it
if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
//Check acceptance selection
- if(IsFidutialCutOn()){
- Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;
+ if(IsFiducialCutOn()){
+ Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
if(! in ) continue ;
}
//Create AOD for analysis
AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
- aodpi0.SetLabel(calo->GetLabel(0));
+ aodpi0.SetLabel(calo->GetLabel());
//Set the indeces of the original caloclusters
aodpi0.SetCaloLabel(calo->GetID(),-1);
aodpi0.SetDetector(fCalorimeter);
if(GetDebug() > 1)
- 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());
+ printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodpi0.Pt(),aodpi0.Phi(),aodpi0.Eta());
//Check Distance to Bad channel, set bit.
- Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
+ Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
continue ;
- if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: Bad channel cut passed %4.2f\n",distBad);
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
if(distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
//Check PID
//PID selection or bit setting
- if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
- //Get most probable PID, check PID weights (in MC this option is mandatory)
- aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
- if(GetDebug() > 1)
- printf("AliAnaPi0EbE::kSSCalo::::FillAOD: PDG of identified particle %d\n",aodpi0.GetPdg());
- //If primary is not pi0, skip it.
- if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
- }
- else if(IsCaloPIDOn()){
+ if(IsCaloPIDOn()){
//Skip matched clusters with tracks
- if(calo->GetNTracksMatched() > 0) continue ;
+ if(IsTrackMatched(calo, GetReader()->GetInputEvent())) continue ;
- //Get most probable PID, 2 options check PID weights
- //or redo PID, recommended option for EMCal.
- if(!IsCaloPIDRecalculationOn())
- aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
- else
- aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
+ // Get most probable PID, 2 options check bayesian PID weights or redo PID
+ // By default, redo PID
+
+ aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
- if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: PDG of identified particle %d\n",aodpi0.GetPdg());
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
//If cluster does not pass pid, not pi0, skip it.
- if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
+ if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
}
else{
- //Set PID bits for later selection (AliAnaPi0 for example)
+ //Set PID bits for later selection
//GetPDG already called in SetPIDBits.
- GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0);
- if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: PID Bits set \n");
+ GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodpi0, GetCaloUtils(), GetReader()->GetInputEvent());
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PID Bits set \n");
}
- if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetPdg());
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
//Play with the MC stack if available
//Check origin of the candidates
if(IsDataMC()){
if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
- GetReader()->GetDataType() != AliCaloTrackReader::kMC){
- aodpi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
- if(GetDebug() > 0) printf("AliAnaPi0EbE::kSSCalo::EbE::FillAOD: Origin of candidate %d\n",aodpi0.GetTag());
+ GetReader()->GetDataType() != AliCaloTrackReader::kMC){
+ //aodpi0.SetInputFileIndex(input);
+ Int_t tag =0;
+ tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(), aodpi0.GetInputFileIndex());
+ //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
+ aodpi0.SetTag(tag);
+ if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
}
}//Work with stack also
}//loop
- if(GetDebug() > 1) printf("AliAnaPi0EbE::kSSCalo::::FillAOD: End fill AODs \n");
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
}
//__________________________________________________________________
//Do analysis and fill histograms
if(!GetOutputAODBranch()){
- printf("AliAnaPi0EbE::FillHistos - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
+ printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
abort();
}
//Loop on stored AOD pi0
Int_t naod = GetOutputAODBranch()->GetEntriesFast();
- if(GetDebug() > 0) printf("AliAnaPi0EbE::Histo::pi0 aod branch entries %d\n", naod);
+ if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
for(Int_t iaod = 0; iaod < naod ; iaod++){
AliAODPWG4Particle* pi0 = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
- Int_t pdg = pi0->GetPdg();
-
- if(pdg != AliCaloPID::kPi0) continue;
+ Int_t pdg = pi0->GetIdentifiedParticleType();
+
+ if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
//Fill pi0 histograms
- Float_t pt = pi0->Pt();
- Float_t phi = pi0->Phi();
+ Float_t ener = pi0->E();
+ Float_t pt = pi0->Pt();
+ Float_t phi = pi0->Phi();
+ if(phi < 0) phi+=TMath::TwoPi();
Float_t eta = pi0->Eta();
-
- fhPtPi0 ->Fill(pt);
- fhPhiPi0 ->Fill(pt,phi);
- fhEtaPi0 ->Fill(pt,eta);
+ fhPtPi0 ->Fill(pt);
+ fhEPi0 ->Fill(ener);
+ fhEEtaPhiPi0 ->Fill(ener,eta,phi);
if(IsDataMC()){
if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
- GetReader()->GetDataType() != AliCaloTrackReader::kMC){
- if(pi0->GetTag()== AliMCAnalysisUtils::kMCPi0){
- fhPtMCPi0 ->Fill(pt);
- fhPhiMCPi0 ->Fill(pt,phi);
- fhEtaMCPi0 ->Fill(pt,eta);
- }
- else{
- fhPtMCNoPi0 ->Fill(pt);
- fhPhiMCNoPi0 ->Fill(pt,phi);
- fhEtaMCNoPi0 ->Fill(pt,eta);
- }
+ GetReader()->GetDataType() != AliCaloTrackReader::kMC){
+ if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0)){
+ fhPtMCPi0 ->Fill(pt);
+ fhPhiMCPi0 ->Fill(pt,phi);
+ fhEtaMCPi0 ->Fill(pt,eta);
+ }
+ else{
+ fhPtMCNoPi0 ->Fill(pt);
+ fhPhiMCNoPi0 ->Fill(pt,phi);
+ fhEtaMCNoPi0 ->Fill(pt,eta);
+ }
}
}//Histograms with MC
}
-
-//____________________________________________________________________________
-void AliAnaPi0EbE::Init()
-{
- //Init
- //Do some checks
- if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn()){
- printf("!!ABORT: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
- abort();
- }
- else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn()){
- printf("!!ABORT: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
- abort();
- }
-
-}
-
-//____________________________________________________________________________
-void AliAnaPi0EbE::InitParameters()
-{
- //Initialize the parameters of the analysis.
- SetOutputAODClassName("AliAODPWG4Particle");
- SetOutputAODName("pi0s");
- fInputAODGammaConvName = "gammaconv" ;
- fAnaType = kIMCalo ;
- fCalorimeter = "EMCAL" ;
- fMinDist = 2.;
- fMinDist2 = 4.;
- fMinDist3 = 5.;
-
-}
-
//__________________________________________________________________
void AliAnaPi0EbE::Print(const Option_t * opt) const
{
printf(" \n") ;
}
+
+//___________________________________________________________________________________
+void AliAnaPi0EbE::RecalibrateCellAmplitude(Float_t & amp, const Int_t id)
+{
+ //Recaculate cell energy if recalibration factor
+
+ Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
+ Int_t nModule = GetModuleNumberCellIndexes(id,fCalorimeter, icol, irow, iRCU);
+
+ if (GetCaloUtils()->IsRecalibrationOn()) {
+ if(fCalorimeter == "PHOS") {
+ amp *= GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
+ }
+ else {
+ amp *= GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
+ }
+ }
+}
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff