#include "AliFiducialCut.h"
#include "TParticle.h"
#include "AliVCluster.h"
+#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODMCParticle.h"
AliAnaPi0EbE::AliAnaPi0EbE() :
AliAnaCaloTrackCorrBaseClass(),fAnaType(kIMCalo), fCalorimeter(""),
fMinDist(0.),fMinDist2(0.), fMinDist3(0.),
- fTimeCutMin(-10000), fTimeCutMax(10000),
+ fNLMCutMin(-1), fNLMCutMax(10),
+ fUseSplitAsyCut(kFALSE),
+ fTimeCutMin(-10000), fTimeCutMax(10000),
+ fFillPileUpHistograms(0),
fFillWeightHistograms(kFALSE), fFillTMHisto(0),
fFillSelectClHisto(0), fFillOnlySimpleSSHisto(1),
fInputAODGammaConvName(""),
// Histograms
fhPt(0), fhE(0),
fhEEta(0), fhEPhi(0), fhEtaPhi(0),
+ fhPtReject(0), fhEReject(0),
+ fhEEtaReject(0), fhEPhiReject(0), fhEtaPhiReject(0),
+ fhMass(0), fhAsymmetry(0),
+ fhSelectedMass(0), fhSelectedAsymmetry(0),
fhPtDecay(0), fhEDecay(0),
// Shower shape histos
fhEDispersion(0), fhELambda0(0), fhELambda1(0),
fhENCells(0), fhETime(0), fhEPairDiffTime(0),
fhDispEtaE(0), fhDispPhiE(0),
fhSumEtaE(0), fhSumPhiE(0), fhSumEtaPhiE(0),
- fhDispEtaPhiDiffE(0), fhSphericityE(0), fhAsymmetryE(0),
+ fhDispEtaPhiDiffE(0), fhSphericityE(0),
// MC histos
- fhMCPt(), fhMCPhi(), fhMCEta(),
- fhMCPi0PtFraction(0), fhMCEtaPtFraction(0),
+ fhMCE(), fhMCPt(),
+ fhMCPhi(), fhMCEta(),
+ fhMCEReject(), fhMCPtReject(),
+ fhMCPi0PtGenRecoFraction(0), fhMCEtaPtGenRecoFraction(0),
fhMCPi0DecayPt(0), fhMCPi0DecayPtFraction(0),
fhMCEtaDecayPt(0), fhMCEtaDecayPtFraction(0),
fhMCOtherDecayPt(0),
fhTrackMatchedMCParticle(0), fhdEdx(0),
fhEOverP(0), fhEOverPNoTRD(0),
// Number of local maxima in cluster
- fhNLocMax(0)
+ fhNLocMax(0),
+ // PileUp
+ fhTimeENoCut(0), fhTimeESPD(0), fhTimeESPDMulti(0),
+ fhTimeNPileUpVertSPD(0), fhTimeNPileUpVertTrack(0),
+ fhTimeNPileUpVertContributors(0),
+ fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
{
//default ctor
for(Int_t i = 0; i < 6; i++)
{
+ fhMCE [i] = 0;
fhMCPt [i] = 0;
fhMCPhi [i] = 0;
fhMCEta [i] = 0;
fhAsymmetryLambda0 [j] = 0;
fhAsymmetryDispEta [j] = 0;
fhAsymmetryDispPhi [j] = 0;
- }
+
+ fhPtPi0PileUp [j] = 0;
+ }
for(Int_t i = 0; i < 3; i++)
{
}
+//_______________________________________________________________________________
+void AliAnaPi0EbE::FillPileUpHistograms(const Float_t energy, const Float_t time)
+{
+ // Fill some histograms to understand pile-up
+ if(!fFillPileUpHistograms) return;
+
+ //printf("E %f, time %f\n",energy,time);
+ AliVEvent * event = GetReader()->GetInputEvent();
+
+ fhTimeENoCut->Fill(energy,time);
+ if(GetReader()->IsPileUpFromSPD()) fhTimeESPD ->Fill(energy,time);
+ if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
+
+ if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
+
+ AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
+ AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
+
+ // N pile up vertices
+ Int_t nVerticesSPD = -1;
+ Int_t nVerticesTracks = -1;
+
+ if (esdEv)
+ {
+ nVerticesSPD = esdEv->GetNumberOfPileupVerticesSPD();
+ nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
+
+ }//ESD
+ else if (aodEv)
+ {
+ nVerticesSPD = aodEv->GetNumberOfPileupVerticesSPD();
+ nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
+ }//AOD
+
+ fhTimeNPileUpVertSPD ->Fill(time,nVerticesSPD);
+ fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
+
+ //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
+ // GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
+
+ Int_t ncont = -1;
+ Float_t z1 = -1, z2 = -1;
+ Float_t diamZ = -1;
+ for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
+ {
+ if (esdEv)
+ {
+ const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
+ ncont=pv->GetNContributors();
+ z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
+ z2 = pv->GetZ();
+ diamZ = esdEv->GetDiamondZ();
+ }//ESD
+ else if (aodEv)
+ {
+ AliAODVertex *pv=aodEv->GetVertex(iVert);
+ if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
+ ncont=pv->GetNContributors();
+ z1=aodEv->GetPrimaryVertexSPD()->GetZ();
+ z2=pv->GetZ();
+ diamZ = aodEv->GetDiamondZ();
+ }// AOD
+
+ Double_t distZ = TMath::Abs(z2-z1);
+ diamZ = TMath::Abs(z2-diamZ);
+
+ fhTimeNPileUpVertContributors ->Fill(time,ncont);
+ fhTimePileUpMainVertexZDistance->Fill(time,distZ);
+ fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
+
+ }// loop
+}
+
+
+//___________________________________________________________________________________________
+void AliAnaPi0EbE::FillRejectedClusterHistograms(const TLorentzVector mom, const Int_t mctag)
+{
+ // Fill histograms that do not pass the identification (SS case only)
+
+ Float_t ener = mom.E();
+ Float_t pt = mom.Pt();
+ Float_t phi = mom.Phi();
+ if(phi < 0) phi+=TMath::TwoPi();
+ Float_t eta = mom.Eta();
+
+ fhPtReject ->Fill(pt);
+ fhEReject ->Fill(ener);
+
+ fhEEtaReject ->Fill(ener,eta);
+ fhEPhiReject ->Fill(ener,phi);
+ fhEtaPhiReject ->Fill(eta,phi);
+
+ if(IsDataMC())
+ {
+ Int_t mcIndex = GetMCIndex(mctag);
+ fhMCEReject [mcIndex] ->Fill(ener);
+ fhMCPtReject [mcIndex] ->Fill(pt);
+ }
+}
+
//_____________________________________________________________________________________
void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
const Int_t nMaxima,
if (fAnaType==kSSCalo)
{
// Asymmetry histograms
- fhAsymmetryE ->Fill(e ,asy);
fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
if (fAnaType==kSSCalo)
{
- fhMCEAsymmetry [mcIndex]->Fill(e ,asy);
fhMCAsymmetryLambda0[ebin][mcIndex]->Fill(l0 ,asy);
fhMCAsymmetryDispEta[ebin][mcIndex]->Fill(dEta,asy);
fhMCAsymmetryDispPhi[ebin][mcIndex]->Fill(dPhi,asy);
Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();
Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
+ Int_t ntimebins= GetHistogramRanges()->GetHistoTimeBins();
+ Float_t timemax = GetHistogramRanges()->GetHistoTimeMax();
+ Float_t timemin = GetHistogramRanges()->GetHistoTimeMin();
+
TString nlm[] ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
TString pname[] ={"Photon","Conversion", "Pi0", "Eta", "Electron","Hadron"};
fhEtaPhi->SetXTitle("#eta");
outputContainer->Add(fhEtaPhi) ;
+ if(fAnaType == kSSCalo)
+ {
+ fhPtReject = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
+ fhPtReject->SetYTitle("N");
+ fhPtReject->SetXTitle("p_{T} (GeV/c)");
+ outputContainer->Add(fhPtReject) ;
+
+ fhEReject = new TH1F("hEReject","Number of rejected as #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
+ fhEReject->SetYTitle("N");
+ fhEReject->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEReject) ;
+
+ fhEPhiReject = new TH2F
+ ("hEPhiReject","Rejected #pi^{0} (#eta) cluster: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
+ fhEPhiReject->SetYTitle("#phi (rad)");
+ fhEPhiReject->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEPhiReject) ;
+
+ fhEEtaReject = new TH2F
+ ("hEEtaReject","Rejected #pi^{0} (#eta) cluster: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
+ fhEEtaReject->SetYTitle("#eta");
+ fhEEtaReject->SetXTitle("E (GeV)");
+ outputContainer->Add(fhEEtaReject) ;
+
+ fhEtaPhiReject = new TH2F
+ ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
+ fhEtaPhiReject->SetYTitle("#phi (rad)");
+ fhEtaPhiReject->SetXTitle("#eta");
+ outputContainer->Add(fhEtaPhiReject) ;
+ }
+
+ fhMass = new TH2F
+ ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+ fhMass->SetYTitle("mass (GeV/c^{2})");
+ fhMass->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMass) ;
+
+ fhSelectedMass = new TH2F
+ ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+ fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
+ fhSelectedMass->SetXTitle("E (GeV)");
+ outputContainer->Add(fhSelectedMass) ;
+
if(fAnaType != kSSCalo)
{
fhPtDecay = new TH1F("hPtDecay","Number of identified #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
{
if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
{
- fhMCPi0PtFraction = new TH2F("hMCPi0PtFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), pT versus pT / pT mother",
- nptbins,ptmin,ptmax,100,0,1);
- fhMCPi0PtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
- fhMCPi0PtFraction->SetYTitle("E^{gen}_{T} / E^{gen-mother}_{T}");
- outputContainer->Add(fhMCPi0PtFraction) ;
+ fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), pT versus E primary #pi^{0} / E reco",
+ nptbins,ptmin,ptmax,200,0,2);
+ fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
+ fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
+ outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
- fhMCEtaPtFraction = new TH2F("hMCEtaPtFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta), pT versus pT / pT mother",
- nptbins,ptmin,ptmax,100,0,1);
- fhMCEtaPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
- fhMCEtaPtFraction->SetYTitle("E^{gen}_{T} / E^{gen-mother}_{T}");
- outputContainer->Add(fhMCEtaPtFraction) ;
+ fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
+ nptbins,ptmin,ptmax,200,0,2);
+ fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
+ fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
+ outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
fhMCPi0DecayPt->SetYTitle("N");
fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
outputContainer->Add(fhMCPi0DecayPt) ;
- fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta), pT versus pT / pT mother",
+ fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #pi^{0}",
nptbins,ptmin,ptmax,100,0,1);
fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
- fhMCPi0DecayPtFraction->SetYTitle("E^{gen}_{T} / E^{gen-mother}_{T}");
+ fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
outputContainer->Add(fhMCPi0DecayPtFraction) ;
fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
outputContainer->Add(fhMCEtaDecayPt) ;
- fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay identified as #pi^{0} (#eta), pT versus pT / pT mother",
+ fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay identified as #pi^{0} (#eta), pT versus E primary #gamma / E primary #eta",
nptbins,ptmin,ptmax,100,0,1);
fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
- fhMCEtaDecayPtFraction->SetYTitle("E^{gen}_{T} / E^{gen-mother}_{T}");
+ fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
outputContainer->Add(fhMCEtaDecayPtFraction) ;
fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0}) identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
for(Int_t i = 0; i < 6; i++)
{
+ fhMCE[i] = new TH1F
+ (Form("hE_MC%s",pname[i].Data()),
+ Form("Identified as #pi^{0} (#eta), cluster from %s",
+ ptype[i].Data()),
+ nptbins,ptmin,ptmax);
+ fhMCE[i]->SetYTitle("N");
+ fhMCE[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCE[i]) ;
+
fhMCPt[i] = new TH1F
(Form("hPt_MC%s",pname[i].Data()),
Form("Identified as #pi^{0} (#eta), cluster from %s",
fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
outputContainer->Add(fhMCPt[i]) ;
+ if(fAnaType == kSSCalo)
+ {
+ fhMCEReject[i] = new TH1F
+ (Form("hEReject_MC%s",pname[i].Data()),
+ Form("Rejected as #pi^{0} (#eta), cluster from %s",
+ ptype[i].Data()),
+ nptbins,ptmin,ptmax);
+ fhMCEReject[i]->SetYTitle("N");
+ fhMCEReject[i]->SetXTitle("E (GeV)");
+ outputContainer->Add(fhMCEReject[i]) ;
+
+ fhMCPtReject[i] = new TH1F
+ (Form("hPtReject_MC%s",pname[i].Data()),
+ Form("Rejected as #pi^{0} (#eta), cluster from %s",
+ ptype[i].Data()),
+ nptbins,ptmin,ptmax);
+ fhMCPtReject[i]->SetYTitle("N");
+ fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
+ outputContainer->Add(fhMCPtReject[i]) ;
+ }
+
fhMCPhi[i] = new TH2F
(Form("hPhi_MC%s",pname[i].Data()),
Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
} //Not MC reader
}//Histos with MC
+ if(fAnaType==kSSCalo)
+ {
+ fhAsymmetry = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
+ nptbins,ptmin,ptmax, 200, -1,1);
+ fhAsymmetry->SetXTitle("E (GeV)");
+ fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+ outputContainer->Add(fhAsymmetry);
+
+ fhSelectedAsymmetry = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
+ nptbins,ptmin,ptmax, 200, -1,1);
+ fhSelectedAsymmetry->SetXTitle("E (GeV)");
+ fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+ outputContainer->Add(fhSelectedAsymmetry);
+
+ if(IsDataMC())
+ {
+ for(Int_t i = 0; i< 6; i++)
+ {
+ fhMCEAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
+ Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
+ nptbins,ptmin,ptmax, 200,-1,1);
+ fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
+ fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+ outputContainer->Add(fhMCEAsymmetry[i]);
+ }
+ }
+ }
if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
{
- fhAsymmetryE = new TH2F ("hAsymmetryE","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
- nptbins,ptmin,ptmax, 200, -1,1);
- fhAsymmetryE->SetXTitle("E (GeV)");
- fhAsymmetryE->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
- outputContainer->Add(fhAsymmetryE);
for(Int_t i = 0; i< 3; i++)
{
{
for(Int_t i = 0; i< 6; i++)
{
- fhMCEAsymmetry[i] = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
- Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
- nptbins,ptmin,ptmax, 200,-1,1);
- fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
- fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
- outputContainer->Add(fhMCEAsymmetry[i]);
-
for(Int_t ie = 0; ie < 7; ie++)
{
fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
}
}
+ if(fFillPileUpHistograms)
+ {
+
+ TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
+
+ for(Int_t i = 0 ; i < 7 ; i++)
+ {
+ fhPtPi0PileUp[i] = new TH1F(Form("hPtPi0PileUp%s",pileUpName[i].Data()),
+ Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
+ fhPtPi0PileUp[i]->SetXTitle("p_{T} (GeV/c)");
+ outputContainer->Add(fhPtPi0PileUp[i]);
+ }
+
+ fhTimeENoCut = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
+ fhTimeENoCut->SetXTitle("E (GeV)");
+ fhTimeENoCut->SetYTitle("time (ns)");
+ outputContainer->Add(fhTimeENoCut);
+
+ fhTimeESPD = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
+ fhTimeESPD->SetXTitle("E (GeV)");
+ fhTimeESPD->SetYTitle("time (ns)");
+ outputContainer->Add(fhTimeESPD);
+
+ fhTimeESPDMulti = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
+ fhTimeESPDMulti->SetXTitle("E (GeV)");
+ fhTimeESPDMulti->SetYTitle("time (ns)");
+ outputContainer->Add(fhTimeESPDMulti);
+
+ fhTimeNPileUpVertSPD = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
+ fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
+ fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
+ outputContainer->Add(fhTimeNPileUpVertSPD);
+
+ fhTimeNPileUpVertTrack = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
+ fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
+ fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
+ outputContainer->Add(fhTimeNPileUpVertTrack);
+
+ fhTimeNPileUpVertContributors = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
+ fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
+ fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
+ outputContainer->Add(fhTimeNPileUpVertContributors);
+
+ fhTimePileUpMainVertexZDistance = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimebins,timemin,timemax,100,0,50);
+ fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
+ fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
+ outputContainer->Add(fhTimePileUpMainVertexZDistance);
+
+ fhTimePileUpMainVertexZDiamond = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
+ fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
+ fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
+ outputContainer->Add(fhTimePileUpMainVertexZDiamond);
+
+ }
+
//Keep neutral meson selection histograms if requiered
//Setting done in AliNeutralMesonSelection
if(nMaxima1 > 1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
if(nMaxima2 > 1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
+ //Skip events with too few or too many NLM
+ if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
+
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
+
+ //Mass of all pairs
+ fhMass->Fill(epair,(mom1+mom2).M());
+
//Select good pair (good phi, pt cuts, aperture and invariant mass)
if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
{
fhPtDecay->Fill(photon2->Pt());
fhEDecay ->Fill(photon2->E() );
-
+
//Create AOD for analysis
mom = mom1+mom2;
+
+ //Mass of selected pairs
+ fhSelectedMass->Fill(epair,mom.M());
+
+ // Fill histograms to undertand pile-up before other cuts applied
+ // Remember to relax time cuts in the reader
+ FillPileUpHistograms(mom.E(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9) /2);
AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
Int_t nCTS = inputAODGammaConv->GetEntriesFast();
Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
- if(nCTS<=0 || nCalo <=0) {
+ if(nCTS<=0 || nCalo <=0)
+ {
if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
return;
}
else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
else fhMassPairLocMax[2]->Fill(epair,mass);
+ if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
+
//Play with the MC stack if available
if(IsDataMC())
{
HasPairSameMCMother(photon1, photon2, label, tag) ;
}
+ //Mass of selected pairs
+ fhMass->Fill(epair,(mom1+mom2).M());
+
//Select good pair (good phi, pt cuts, aperture and invariant mass)
if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
{
mom = mom1+mom2;
+ //Mass of selected pairs
+ fhSelectedMass->Fill(epair,mom.M());
+
+ // Fill histograms to undertand pile-up before other cuts applied
+ // Remember to relax time cuts in the reader
+ FillPileUpHistograms(mom.E(),cluster->GetTOF()*1e9);
+
AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ; //vertex cut
//Get Momentum vector,
+ Double_t vertex[]={0,0,0};
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(! in ) continue ;
}
- //Create AOD for analysis
- AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
- aodpi0.SetLabel(calo->GetLabel());
-
- //Set the indeces of the original caloclusters
- aodpi0.SetCaloLabel(calo->GetID(),-1);
- aodpi0.SetDetector(fCalorimeter);
if(GetDebug() > 1)
- 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());
-
+ printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",mom.Pt(),mom.Phi(),mom.Eta());
+
//Check Distance to Bad channel, set bit.
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::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
+
//.......................................
// TOF cut, BE CAREFUL WITH THIS CUT
Double_t tof = calo->GetTOF()*1e9;
if(tof < fTimeCutMin || tof > fTimeCutMax) continue ;
+ //Play with the MC stack if available
+ //Check origin of the candidates
+ Int_t tag = 0 ;
+ if(IsDataMC())
+ {
+ tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(),0);
+ //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
+ if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
+ }
- 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) ;
- else aodpi0.SetDistToBad(0) ;
+ //Skip matched clusters with tracks
+ if(IsTrackMatched(calo, GetReader()->GetInputEvent()))
+ {
+ FillRejectedClusterHistograms(mom,tag);
+ continue ;
+ }
+
//Check PID
//PID selection or bit setting
Int_t nMaxima = 0 ;
Double_t mass = 0 , angle = 0;
Double_t e1 = 0 , e2 = 0;
- //Skip matched clusters with tracks
- if(IsTrackMatched(calo, GetReader()->GetInputEvent())) continue ;
-
- // Check if cluster is pi0 via cluster splitting
- aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
- GetVertex(evtIndex),nMaxima,
- mass,angle,e1,e2));
+ Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
+ GetVertex(evtIndex),nMaxima,
+ mass,angle,e1,e2) ;
- if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
+
+
+ //Skip events with too few or too many NLM
+ if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
+ {
+ FillRejectedClusterHistograms(mom,tag);
+ continue ;
+ }
- // If cluster does not pass pid, not pi0, skip it.
- // TO DO, add option for Eta ... or conversions
- if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
+ if(GetDebug() > 1)
+ printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
- if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",
- aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
+ //mass of all clusters
+ fhMass->Fill(mom.E(),mass);
+
+ // Asymmetry of all clusters
+ Float_t asy =-10;
+ if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
+ fhAsymmetry->Fill(mom.E(),asy);
- //Play with the MC stack if available
- //Check origin of the candidates
- Int_t tag = 0 ;
if(IsDataMC())
{
- tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),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
+ Int_t mcIndex = GetMCIndex(tag);
+ fhMCEAsymmetry[mcIndex]->Fill(mom.E(),asy);
+ }
+
+ // If cluster does not pass pid, not pi0/eta, skip it.
+ if (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
+ {
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Pi0\n");
+ FillRejectedClusterHistograms(mom,tag);
+ continue ;
+ }
+
+ else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
+ {
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Eta\n");
+ FillRejectedClusterHistograms(mom,tag);
+ continue ;
+ }
+
+ if(GetDebug() > 1)
+ printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
+ mom.Pt(), idPartType);
+
+ fhSelectedAsymmetry->Fill(mom.E(),asy);
+
+ if( fUseSplitAsyCut && GetCaloPID()->IsInPi0SplitAsymmetryRange(mom.E(),asy,nMaxima) )
+ {
+ if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Too large asymmetry\n");
+ FillRejectedClusterHistograms(mom,tag);
+ continue ;
+ }
+
+ //Mass of selected pairs
+ fhSelectedMass ->Fill(mom.E(),mass);
+
+ //-----------------------
+ //Create AOD for analysis
+
+ AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
+ aodpi0.SetLabel(calo->GetLabel());
+
+ //Set the indeces of the original caloclusters
+ aodpi0.SetCaloLabel(calo->GetID(),-1);
+ aodpi0.SetDetector(fCalorimeter);
+
+ if (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
+ else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
+ else aodpi0.SetDistToBad(0) ;
+
+ // Check if cluster is pi0 via cluster splitting
+ aodpi0.SetIdentifiedParticleType(idPartType);
+
+ // Add number of local maxima to AOD, method name in AOD to be FIXED
+ aodpi0.SetFiducialArea(nMaxima);
+
+ aodpi0.SetTag(tag);
//Fill some histograms about shower shape
if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
{
- Float_t asy =-10;
- if(e1+e2 > 0 ) asy = (e1-e2) / (e1+e2);
FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
- }
+ }
+
+ // Fill histograms to undertand pile-up before other cuts applied
+ // Remember to relax time cuts in the reader
+ FillPileUpHistograms(calo->E(),calo->GetTOF()*1e9);
//Add AOD with pi0 object to aod branch
AddAODParticle(aodpi0);
fhEPhi ->Fill(ener,phi);
fhEtaPhi ->Fill(eta,phi);
+ if(fFillPileUpHistograms)
+ {
+ if(GetReader()->IsPileUpFromSPD()) fhPtPi0PileUp[0]->Fill(pt);
+ if(GetReader()->IsPileUpFromEMCal()) fhPtPi0PileUp[1]->Fill(pt);
+ if(GetReader()->IsPileUpFromSPDOrEMCal()) fhPtPi0PileUp[2]->Fill(pt);
+ if(GetReader()->IsPileUpFromSPDAndEMCal()) fhPtPi0PileUp[3]->Fill(pt);
+ if(GetReader()->IsPileUpFromSPDAndNotEMCal()) fhPtPi0PileUp[4]->Fill(pt);
+ if(GetReader()->IsPileUpFromEMCalAndNotSPD()) fhPtPi0PileUp[5]->Fill(pt);
+ if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPi0PileUp[6]->Fill(pt);
+ }
+
+
if(IsDataMC())
{
Int_t tag = pi0->GetTag();
Int_t mcIndex = GetMCIndex(tag);
+ fhMCE [mcIndex] ->Fill(ener);
fhMCPt [mcIndex] ->Fill(pt);
fhMCPhi[mcIndex] ->Fill(pt,phi);
fhMCEta[mcIndex] ->Fill(pt,eta);
- if(mcIndex==kmcPhoton && fAnaType==kSSCalo)
+ if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
{
Float_t efracMC = 0;
Int_t label = pi0->GetLabel();
TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
if(grandmom.E() > 0 && ok)
{
- efracMC = mom.E()/grandmom.E();
- fhMCPi0PtFraction ->Fill(pt,efracMC);
+ efracMC = grandmom.E()/ener;
+ fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
}
}
else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
{
fhMCPi0DecayPt->Fill(pt);
- TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
+ TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
if(grandmom.E() > 0 && ok)
{
efracMC = mom.E()/grandmom.E();
TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
if(grandmom.E() > 0 && ok)
{
- efracMC = mom.E()/grandmom.E();
- fhMCEtaPtFraction ->Fill(pt,efracMC);
+ efracMC = grandmom.E()/ener;
+ fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
}
}
else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff