#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),
- fFillWeightHistograms(kFALSE), fFillTMHisto(0), fFillSelectClHisto(0),
+ 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
- fhPtMCNo(0), fhPhiMCNo(0), fhEtaMCNo(0),
- fhPtMC(0), fhPhiMC(0), fhEtaMC(0),
+ fhMCE(), fhMCPt(),
+ fhMCPhi(), fhMCEta(),
+ fhMCEReject(), fhMCPtReject(),
+ fhMCPi0PtGenRecoFraction(0), fhMCEtaPtGenRecoFraction(0),
+ fhMCPi0DecayPt(0), fhMCPi0DecayPtFraction(0),
+ fhMCEtaDecayPt(0), fhMCEtaDecayPtFraction(0),
+ fhMCOtherDecayPt(0),
fhMassPairMCPi0(0), fhMassPairMCEta(0),
fhAnglePairMCPi0(0), fhAnglePairMCEta(0),
// Weight studies
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;
fhEMCLambda0 [i] = 0;
fhEMCLambda0NoTRD [i] = 0;
fhEMCLambda0FracMaxCellCut[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,
Float_t ll0 = 0., ll1 = 0.;
Float_t dispp= 0., dEta = 0., dPhi = 0.;
Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
- if(fCalorimeter == "EMCAL")
+ if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
{
GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
-
+
fhDispEtaE -> Fill(e,dEta);
fhDispPhiE -> Fill(e,dPhi);
fhSumEtaE -> Fill(e,sEta);
if (fAnaType==kSSCalo)
{
// Asymmetry histograms
- fhAsymmetryE ->Fill(e ,asy);
fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
fhELambda0LocMax [indexMax]->Fill(e,l0);
fhELambda1LocMax [indexMax]->Fill(e,l1);
fhEDispersionLocMax[indexMax]->Fill(e,disp);
- if(fCalorimeter=="EMCAL")
+
+ if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
{
fhEDispEtaLocMax [indexMax]-> Fill(e,dEta);
fhEDispPhiLocMax [indexMax]-> Fill(e,dPhi);
if(IsDataMC())
{
- Int_t mcIndex = 0;
-
- if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
- {
- mcIndex = kmcPi0 ;
- }//pi0
- else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
- {
- mcIndex = kmcEta ;
- }//eta
- else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
- GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
- {
- mcIndex = kmcConversion ;
- }//conversion photon
- else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
- {
- mcIndex = kmcPhoton ;
- }//photon no conversion
- else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
- {
- mcIndex = kmcElectron ;
- }//electron
- else
- {
- mcIndex = kmcHadron ;
- }//other particles
+ Int_t mcIndex = GetMCIndex(tag);
fhEMCLambda0[mcIndex] ->Fill(e, l0);
fhEMCLambda1[mcIndex] ->Fill(e, l1);
if(fCalorimeter=="EMCAL" && nSM < 6)
fhEMCLambda0NoTRD[mcIndex]->Fill(e, l0 );
+
if(maxCellFraction < 0.5)
fhEMCLambda0FracMaxCellCut[mcIndex]->Fill(e, l0 );
- if(fCalorimeter == "EMCAL")
+ if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
{
fhMCEDispEta [mcIndex]-> Fill(e,dEta);
fhMCEDispPhi [mcIndex]-> Fill(e,dPhi);
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);
return new TObjString(parList) ;
}
-//__________________________________________________________________
-void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
- AliAODPWG4Particle * photon2,
- Int_t & label, Int_t & tag)
-{
- // Check the labels of pare in case mother was same pi0 or eta
- // Set the new AOD accordingly
-
- Int_t label1 = photon1->GetLabel();
- Int_t label2 = photon2->GetLabel();
-
- if(label1 < 0 || label2 < 0 ) return ;
-
- //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
- //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
- Int_t tag1 = photon1->GetTag();
- Int_t tag2 = photon2->GetTag();
-
- 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) ) ||
- (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
- GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay) )
- )
- {
-
- //Check if pi0/eta 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
- }
- } // STACK
- 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
- }
- }// AOD
-
- //printf("mother1 %d, mother2 %d\n",label1,label2);
- if( label1 == label2 && label1>=0 )
- {
-
- label = label1;
-
- TLorentzVector mom1 = *(photon1->Momentum());
- TLorentzVector mom2 = *(photon2->Momentum());
-
- Double_t angle = mom2.Angle(mom1.Vect());
- Double_t mass = (mom1+mom2).M();
- Double_t epair = (mom1+mom2).E();
-
- if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
- {
- fhMassPairMCPi0 ->Fill(epair,mass);
- fhAnglePairMCPi0->Fill(epair,angle);
- GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
- }
- else
- {
- fhMassPairMCEta ->Fill(epair,mass);
- fhAnglePairMCEta->Fill(epair,angle);
- GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
- }
-
- } // same label
- } // both from eta or pi0 decay
-
-}
-
//_____________________________________________
TList * AliAnaPi0EbE::GetCreateOutputObjects()
{
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);
fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
outputContainer->Add(fhEFracMaxCellNoTRD) ;
-
- fhDispEtaE = new TH2F ("hDispEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
- fhDispEtaE->SetXTitle("E (GeV)");
- fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
- outputContainer->Add(fhDispEtaE);
-
- fhDispPhiE = new TH2F ("hDispPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
- fhDispPhiE->SetXTitle("E (GeV)");
- fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
- outputContainer->Add(fhDispPhiE);
-
- fhSumEtaE = new TH2F ("hSumEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
- fhSumEtaE->SetXTitle("E (GeV)");
- fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
- outputContainer->Add(fhSumEtaE);
-
- fhSumPhiE = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
- nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
- fhSumPhiE->SetXTitle("E (GeV)");
- fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
- outputContainer->Add(fhSumPhiE);
-
- fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
- nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
- fhSumEtaPhiE->SetXTitle("E (GeV)");
- fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
- outputContainer->Add(fhSumEtaPhiE);
-
- fhDispEtaPhiDiffE = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
- nptbins,ptmin,ptmax,200, -10,10);
- fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
- fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
- outputContainer->Add(fhDispEtaPhiDiffE);
-
- fhSphericityE = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
- nptbins,ptmin,ptmax, 200, -1,1);
- fhSphericityE->SetXTitle("E (GeV)");
- fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
- outputContainer->Add(fhSphericityE);
-
- for(Int_t i = 0; i < 7; i++)
+ if(!fFillOnlySimpleSSHisto)
{
- fhDispEtaDispPhi[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
- ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
- fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
- fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
- outputContainer->Add(fhDispEtaDispPhi[i]);
+ fhDispEtaE = new TH2F ("hDispEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
+ fhDispEtaE->SetXTitle("E (GeV)");
+ fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
+ outputContainer->Add(fhDispEtaE);
+
+ fhDispPhiE = new TH2F ("hDispPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
+ fhDispPhiE->SetXTitle("E (GeV)");
+ fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
+ outputContainer->Add(fhDispPhiE);
+
+ fhSumEtaE = new TH2F ("hSumEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E", nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
+ fhSumEtaE->SetXTitle("E (GeV)");
+ fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
+ outputContainer->Add(fhSumEtaE);
- fhLambda0DispEta[i] = new TH2F (Form("hLambda0DispEta_EBin%d",i),Form("#lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
- ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
- fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
- fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
- outputContainer->Add(fhLambda0DispEta[i]);
+ fhSumPhiE = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
+ nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
+ fhSumPhiE->SetXTitle("E (GeV)");
+ fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
+ outputContainer->Add(fhSumPhiE);
- fhLambda0DispPhi[i] = new TH2F (Form("hLambda0DispPhi_EBin%d",i),Form("#lambda^{2}_{0}} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",bin[i],bin[i+1]),
- ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
- fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
- fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
- outputContainer->Add(fhLambda0DispPhi[i]);
+ fhSumEtaPhiE = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
+ nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
+ fhSumEtaPhiE->SetXTitle("E (GeV)");
+ fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
+ outputContainer->Add(fhSumEtaPhiE);
+ fhDispEtaPhiDiffE = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
+ nptbins,ptmin,ptmax,200, -10,10);
+ fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
+ fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
+ outputContainer->Add(fhDispEtaPhiDiffE);
+
+ fhSphericityE = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
+ nptbins,ptmin,ptmax, 200, -1,1);
+ fhSphericityE->SetXTitle("E (GeV)");
+ fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
+ outputContainer->Add(fhSphericityE);
+
+ for(Int_t i = 0; i < 7; i++)
+ {
+ fhDispEtaDispPhi[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
+ ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
+ fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
+ fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+ outputContainer->Add(fhDispEtaDispPhi[i]);
+
+ fhLambda0DispEta[i] = new TH2F (Form("hLambda0DispEta_EBin%d",i),Form("#lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
+ ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
+ fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
+ fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
+ outputContainer->Add(fhLambda0DispEta[i]);
+
+ fhLambda0DispPhi[i] = new TH2F (Form("hLambda0DispPhi_EBin%d",i),Form("#lambda^{2}_{0}} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",bin[i],bin[i+1]),
+ ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
+ fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
+ fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+ outputContainer->Add(fhLambda0DispPhi[i]);
+
+ }
}
}
fhEDispersionLocMax[i]->SetXTitle("E (GeV)");
outputContainer->Add(fhEDispersionLocMax[i]) ;
- if(fCalorimeter == "EMCAL")
+ if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
{
fhEDispEtaLocMax[i] = new TH2F(Form("hEDispEtaLocMax%d",i+1),
Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
fhETime->SetYTitle("t (ns)");
outputContainer->Add(fhETime);
- }// Invariant mass analysis in calorimeters and calorimeter + conversion photons
+ }
if(fAnaType == kIMCalo)
{
outputContainer->Add(fhEOverPNoTRD);
}
- if(IsDataMC())
+ if(IsDataMC() && fFillTMHisto)
{
fhTrackMatchedMCParticle = new TH2F
("hTrackMatchedMCParticle",
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) ");
if(IsDataMC())
{
+ if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
+ {
+ 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) ;
+
+ 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 E primary #gamma / E primary #pi^{0}",
+ nptbins,ptmin,ptmax,100,0,1);
+ fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
+ 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->SetYTitle("N");
+ 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 E primary #gamma / E primary #eta",
+ nptbins,ptmin,ptmax,100,0,1);
+ fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
+ 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);
+ fhMCOtherDecayPt->SetYTitle("N");
+ fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
+ outputContainer->Add(fhMCOtherDecayPt) ;
+
+ }
+
if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
GetReader()->GetDataType() != AliCaloTrackReader::kMC)
{
- fhPtMC = new TH1F("hPtMC","Identified #pi^{0} (#eta) from #pi^{0} (#eta)",nptbins,ptmin,ptmax);
- fhPtMC->SetYTitle("N");
- fhPtMC->SetXTitle("p_{T} (GeV/c)");
- outputContainer->Add(fhPtMC) ;
-
- fhPhiMC = new TH2F
- ("hPhiMC","Identified #pi^{0} (#eta) from #pi^{0} (#eta)",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiMC->SetYTitle("#phi");
- fhPhiMC->SetXTitle("p_{T} (GeV/c)");
- outputContainer->Add(fhPhiMC) ;
-
- fhEtaMC = new TH2F
- ("hEtaMC","Identified #pi^{0} (#eta) from #pi^{0} (#eta)",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaMC->SetYTitle("#eta");
- fhEtaMC->SetXTitle("p_{T} (GeV/c)");
- outputContainer->Add(fhEtaMC) ;
-
- fhPtMCNo = new TH1F("hPtMCNo","Identified #pi^{0} (#eta) not from #pi^{0} (#eta)",nptbins,ptmin,ptmax);
- fhPtMCNo->SetYTitle("N");
- fhPtMCNo->SetXTitle("p_{T} (GeV/c)");
- outputContainer->Add(fhPtMCNo) ;
-
- fhPhiMCNo = new TH2F
- ("hPhiMCNo","Identified #pi^{0} (#eta) not from #pi^{0} (#eta)",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
- fhPhiMCNo->SetYTitle("#phi");
- fhPhiMCNo->SetXTitle("p_{T} (GeV/c)");
- outputContainer->Add(fhPhiMCNo) ;
-
- fhEtaMCNo = new TH2F
- ("hEtaMCNo","Identified #pi^{0} (#eta) not from #pi^{0} (#eta)",nptbins,ptmin,ptmax,netabins,etamin,etamax);
- fhEtaMCNo->SetYTitle("#eta");
- fhEtaMCNo->SetXTitle("p_{T} (GeV/c)");
- outputContainer->Add(fhEtaMCNo) ;
-
fhAnglePairMCPi0 = new TH2F
("AnglePairMCPi0",
"Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
outputContainer->Add(fhAnglePairMCPi0) ;
- fhAnglePairMCEta = new TH2F
- ("AnglePairMCEta",
- "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
- fhAnglePairMCEta->SetYTitle("#alpha (rad)");
- fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
- outputContainer->Add(fhAnglePairMCEta) ;
-
- fhMassPairMCPi0 = new TH2F
- ("MassPairMCPi0",
- "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
- fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
- fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
- outputContainer->Add(fhMassPairMCPi0) ;
-
- fhMassPairMCEta = new TH2F
- ("MassPairMCEta",
- "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
- fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
- fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
- outputContainer->Add(fhMassPairMCEta) ;
-
- if( fFillSelectClHisto )
+ if (fAnaType!= kSSCalo)
{
- for(Int_t i = 0; i < 6; i++)
- {
+ fhAnglePairMCEta = new TH2F
+ ("AnglePairMCEta",
+ "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
+ fhAnglePairMCEta->SetYTitle("#alpha (rad)");
+ fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
+ outputContainer->Add(fhAnglePairMCEta) ;
+
+ fhMassPairMCPi0 = new TH2F
+ ("MassPairMCPi0",
+ "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
+ fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
+ fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
+ outputContainer->Add(fhMassPairMCPi0) ;
+
+ fhMassPairMCEta = new TH2F
+ ("MassPairMCEta",
+ "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
+ fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
+ fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
+ outputContainer->Add(fhMassPairMCEta) ;
+ }
+
+ 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",
+ ptype[i].Data()),
+ nptbins,ptmin,ptmax);
+ fhMCPt[i]->SetYTitle("N");
+ 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()),
+ nptbins,ptmin,ptmax,nphibins,phimin,phimax);
+ fhMCPhi[i]->SetYTitle("#phi");
+ fhMCPhi[i]->SetXTitle("p_{T} (GeV/c)");
+ outputContainer->Add(fhMCPhi[i]) ;
+
+ fhMCEta[i] = new TH2F
+ (Form("hEta_MC%s",pname[i].Data()),
+ Form("Identified as #pi^{0} (#eta), cluster from %s",
+ ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
+ fhMCEta[i]->SetYTitle("#eta");
+ fhMCEta[i]->SetXTitle("p_{T} (GeV/c)");
+ outputContainer->Add(fhMCEta[i]) ;
+
+
+ if( fFillSelectClHisto )
+ {
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);
fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
-
- fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
- Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
- nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
- fhMCEDispEta[i]->SetXTitle("E (GeV)");
- fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
- outputContainer->Add(fhMCEDispEta[i]);
-
- fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
- Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
- nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
- fhMCEDispPhi[i]->SetXTitle("E (GeV)");
- fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
- outputContainer->Add(fhMCEDispPhi[i]);
-
- fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
- Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptype[i].Data()),
- nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
- fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
- fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
- outputContainer->Add(fhMCESumEtaPhi[i]);
-
- fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
- Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
- nptbins,ptmin,ptmax,200,-10,10);
- fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
- fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
- outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
-
- fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
- Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptype[i].Data()),
- nptbins,ptmin,ptmax, 200,-1,1);
- fhMCESphericity[i]->SetXTitle("E (GeV)");
- fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
- outputContainer->Add(fhMCESphericity[i]);
-
- for(Int_t ie = 0; ie < 7; ie++)
+ if(!fFillOnlySimpleSSHisto)
{
- fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
- Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
- ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
- fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
- fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
- outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
+ fhMCEDispEta[i] = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
+ Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
+ nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
+ fhMCEDispEta[i]->SetXTitle("E (GeV)");
+ fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
+ outputContainer->Add(fhMCEDispEta[i]);
+
+ fhMCEDispPhi[i] = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
+ Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
+ nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
+ fhMCEDispPhi[i]->SetXTitle("E (GeV)");
+ fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+ outputContainer->Add(fhMCEDispPhi[i]);
+
+ fhMCESumEtaPhi[i] = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
+ Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptype[i].Data()),
+ nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
+ fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
+ fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
+ outputContainer->Add(fhMCESumEtaPhi[i]);
- fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
- Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
- ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
- fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
- fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
- outputContainer->Add(fhMCLambda0DispEta[ie][i]);
+ fhMCEDispEtaPhiDiff[i] = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
+ Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
+ nptbins,ptmin,ptmax,200,-10,10);
+ fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
+ fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
+ outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
- fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
- Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
- ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
- fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
- fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
- outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
+ fhMCESphericity[i] = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
+ Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptype[i].Data()),
+ nptbins,ptmin,ptmax, 200,-1,1);
+ fhMCESphericity[i]->SetXTitle("E (GeV)");
+ fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
+ outputContainer->Add(fhMCESphericity[i]);
- }
+ for(Int_t ie = 0; ie < 7; ie++)
+ {
+ fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
+ Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
+ ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
+ fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
+ fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+ outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
+
+ fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
+ Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
+ ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
+ fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
+ fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+ outputContainer->Add(fhMCLambda0DispEta[ie][i]);
+
+ fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
+ Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
+ ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
+ fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
+ fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+ outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
+
+ }
+ }
}
fhEMCLambda0FracMaxCellCut[i] = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[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 )
+ 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
}
+//_____________________________________________
+Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
+{
+
+ // Assign mc index depending on MC bit set
+
+ if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) )
+ {
+ return kmcPi0 ;
+ }//pi0
+ else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) )
+ {
+ return kmcEta ;
+ }//eta
+ else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
+ GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
+ {
+ return kmcConversion ;
+ }//conversion photon
+ else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
+ {
+ return kmcPhoton ;
+ }//photon no conversion
+ else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
+ {
+ return kmcElectron ;
+ }//electron
+ else
+ {
+ return kmcHadron ;
+ }//other particles
+
+}
+
+//__________________________________________________________________
+void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
+ AliAODPWG4Particle * photon2,
+ Int_t & label, Int_t & tag)
+{
+ // Check the labels of pare in case mother was same pi0 or eta
+ // Set the new AOD accordingly
+
+ Int_t label1 = photon1->GetLabel();
+ Int_t label2 = photon2->GetLabel();
+
+ if(label1 < 0 || label2 < 0 ) return ;
+
+ //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
+ //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
+ Int_t tag1 = photon1->GetTag();
+ Int_t tag2 = photon2->GetTag();
+
+ 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) ) ||
+ (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
+ GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay) )
+ )
+ {
+
+ //Check if pi0/eta 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
+ }
+ } // STACK
+ 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
+ }
+ }// AOD
+
+ //printf("mother1 %d, mother2 %d\n",label1,label2);
+ if( label1 == label2 && label1>=0 )
+ {
+
+ label = label1;
+
+ TLorentzVector mom1 = *(photon1->Momentum());
+ TLorentzVector mom2 = *(photon2->Momentum());
+
+ Double_t angle = mom2.Angle(mom1.Vect());
+ Double_t mass = (mom1+mom2).M();
+ Double_t epair = (mom1+mom2).E();
+
+ if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
+ {
+ fhMassPairMCPi0 ->Fill(epair,mass);
+ fhAnglePairMCPi0->Fill(epair,angle);
+ GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
+ }
+ else
+ {
+ fhMassPairMCEta ->Fill(epair,mass);
+ fhAnglePairMCEta->Fill(epair,angle);
+ GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
+ }
+
+ } // same label
+ } // both from eta or pi0 decay
+
+}
+
//____________________________________________________________________________
void AliAnaPi0EbE::Init()
{
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 too small or big pt, skip it
- if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
+ if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
//Check acceptance selection
if(IsFiducialCutOn())
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())
{
- if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
- GetReader()->GetDataType() != AliCaloTrackReader::kMC)
- {
- //aodpi0.SetInputFileIndex(input);
- 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())
{
- if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
- GetReader()->GetDataType() != AliCaloTrackReader::kMC){
- if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0))
+ 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 || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
+ {
+ Float_t efracMC = 0;
+ Int_t label = pi0->GetLabel();
+
+ Bool_t ok = kFALSE;
+ TLorentzVector mom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
+ if(!ok) continue;
+
+ if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
+ {
+ TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
+ if(grandmom.E() > 0 && ok)
+ {
+ efracMC = grandmom.E()/ener;
+ fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
+ }
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
{
- fhPtMC ->Fill(pt);
- fhPhiMC ->Fill(pt,phi);
- fhEtaMC ->Fill(pt,eta);
+ fhMCPi0DecayPt->Fill(pt);
+ TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
+ if(grandmom.E() > 0 && ok)
+ {
+ efracMC = mom.E()/grandmom.E();
+ fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
+ }
}
- else
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
+ {
+ TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
+ if(grandmom.E() > 0 && ok)
+ {
+ efracMC = grandmom.E()/ener;
+ fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
+ }
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
{
- fhPtMCNo ->Fill(pt);
- fhPhiMCNo ->Fill(pt,phi);
- fhEtaMCNo ->Fill(pt,eta);
+ fhMCEtaDecayPt->Fill(pt);
+ TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
+ if(grandmom.E() > 0 && ok)
+ {
+ efracMC = mom.E()/grandmom.E();
+ fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
+ }
+ }
+ else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
+ {
+ fhMCOtherDecayPt->Fill(pt);
}
+
}
+
}//Histograms with MC
}// aod loop