Skip to content
Snippets Groups Projects

New plot to monitor alignment weak modes added to ParticleMassMonitor

Merged New plot to monitor alignment weak modes added to ParticleMassMonitor
miruizdi_monitor_updates into master
All threads resolved!
Merged Miguel Ruiz Diaz requested to merge miruizdi_monitor_updates into master
All threads resolved!
Compare
and
1 file
+ 75
− 65
Compare changes
  • Side-by-side
  • Inline
Phys/ParticleMaker/src/ParticleMassMonitor.cpp
+ 75
− 65
@@ -35,27 +35,17 @@ public:
StatusCode initialize() override {
return Consumer::initialize().andThen( [&] {
/*Assign the names of the daughter particles of the decay. So far only considering D0 -> KPi, Z -> mumu and JPs ->
mumu used for alignment mass constraints*/
if ( m_mother == "D0" ) {
m_daughter1 = "K";
m_daughter2 = "Pi";
} else {
m_daughter1 = "Mu1";
m_daughter2 = "Mu2";
}
// Mass histograms
m_histograms_mass.emplace(
std::piecewise_construct, std::forward_as_tuple( 0 ),
std::forward_as_tuple( this, m_bins, m_mother, m_mass, m_sigma, m_minpt, m_maxpt, m_minp, m_maxp ) );
m_histograms_mass.emplace( std::piecewise_construct, std::forward_as_tuple( 0 ),
+1
std::forward_as_tuple( this, m_bins, m_pdiff_bins, m_particles[0], m_mass, m_sigma, m_minpt,
m_maxpt, m_minp, m_maxp, m_pdiff_range / Gaudi::Units::GeV ) );
// Momentum histograms
std::vector<std::string> particles = {m_mother, m_daughter1, m_daughter2};
for ( int i = 0; i < 3; i++ ) {
for ( unsigned int i = 0; i < m_particles.size(); i++ ) {
m_histograms_momentum.emplace(
std::piecewise_construct, std::forward_as_tuple( i ),
std::forward_as_tuple( this, m_bins, particles[i], m_minpt, m_maxpt, m_minp, m_maxp ) );
std::forward_as_tuple( this, m_bins, m_particles[i], m_minpt, m_maxpt, m_minp, m_maxp ) );
}
} );
}
@@ -64,22 +54,22 @@ public:
private:
// Histogram settings. Default values are the standard ones for D0 -> KPi candidates
Gaudi::Property<unsigned int> m_bins{this, "Bins", 100};
Gaudi::Property<std::string> m_mother{this, "ParticleType", "D0"};
Gaudi::Property<float> m_mass{this, "MassMean", 1864.84 * Gaudi::Units::MeV};
Gaudi::Property<float> m_sigma{this, "MassSigma", 90 * Gaudi::Units::MeV};
Gaudi::Property<float> m_minpt{this, "MinPt", 800 * Gaudi::Units::MeV};
Gaudi::Property<float> m_maxpt{this, "MaxPt", 6000 * Gaudi::Units::MeV};
Gaudi::Property<float> m_minp{this, "MinP", 2000 * Gaudi::Units::MeV};
Gaudi::Property<float> m_maxp{this, "MaxP", 150000 * Gaudi::Units::MeV};
Gaudi::Property<unsigned int> m_bins{this, "Bins", 100};
Gaudi::Property<unsigned int> m_pdiff_bins{this, "PdiffBins", 20};
Gaudi::Property<float> m_mass{this, "MassMean", 1864.84 * Gaudi::Units::MeV};
Gaudi::Property<float> m_sigma{this, "MassSigma", 18 * Gaudi::Units::MeV};
Gaudi::Property<float> m_minpt{this, "MinPt", 800 * Gaudi::Units::MeV};
Gaudi::Property<float> m_maxpt{this, "MaxPt", 5000 * Gaudi::Units::MeV};
Gaudi::Property<float> m_minp{this, "MinP", 2000 * Gaudi::Units::MeV};
Gaudi::Property<float> m_maxp{this, "MaxP", 150000 * Gaudi::Units::MeV};
Gaudi::Property<float> m_pdiff_range{this, "PdiffDelta", 40000 * Gaudi::Units::MeV};
// Vector with the name of the particles in the same order as in the particle combiner
Gaudi::Property<std::vector<std::string>> m_particles{this, "ParticleType", {"D0", "Pi", "K"}};
// Flag for tuple generation
Gaudi::Property<bool> m_generateTuples{this, "GenerateTuples", false};
// Daughter particles
std::string m_daughter1;
std::string m_daughter2;
// Definition of histograms
struct HistoMomentum {
mutable Gaudi::Accumulators::Histogram<1> CandidatePt;
@@ -92,12 +82,14 @@ private:
};
struct HistoMass {
mutable Gaudi::Accumulators::Histogram<1> CandidateMass;
mutable Gaudi::Accumulators::Histogram<2> MassvsPt;
mutable Gaudi::Accumulators::Histogram<2> MassvsP;
HistoMass( const ParticleMassMonitor* owner, unsigned int const& bins, std::string const& name, float const& mass,
float const& sigma, float const& min_pt, float const& max_pt, float const& min_p, float const& max_p )
mutable Gaudi::Accumulators::Histogram<1> CandidateMass;
mutable Gaudi::Accumulators::Histogram<2> MassvsPt;
mutable Gaudi::Accumulators::Histogram<2> MassvsP;
mutable Gaudi::Accumulators::ProfileHistogram<1> MassvsPdiff;
HistoMass( const ParticleMassMonitor* owner, unsigned int const& bins, unsigned int const& bins_pdiff,
std::string const& name, float const& mass, float const& sigma, float const& min_pt, float const& max_pt,
float const& min_p, float const& max_p, float const& range_pdiff )
: CandidateMass{owner,
name + "_mass",
name + " mass",
@@ -109,12 +101,19 @@ private:
, MassvsP{owner,
name + "_MassvsP",
name + " mass vs p",
{{bins, min_p, max_p, "p [MeV/c]"}, {bins, mass - 5 * sigma, mass + 5 * sigma, "m [MeV/c^{2}]"}}} {}
{{bins, min_p, max_p, "p [MeV/c]"}, {bins, mass - 5 * sigma, mass + 5 * sigma, "m [MeV/c^{2}]"}}}
, MassvsPdiff{owner,
name + "_MassvsPdiff",
fmt::format( "m({}) - {:.2f} MeV/c^{{2}} vs daughter p difference", name, mass ),
{bins_pdiff, -range_pdiff, range_pdiff, "Daughter p_{+} - p_{-} [GeV/c]"}} {}
};
// Histogram objects are stored inside of maps
std::map<int, HistoMomentum> m_histograms_momentum;
std::map<int, HistoMass> m_histograms_mass;
// Function to calculate the momentum difference of the daughter particles
double calculate_momentum_difference( LHCb::Particle::ConstVector const& particles ) const;
};
DECLARE_COMPONENT( ParticleMassMonitor )
@@ -122,57 +121,68 @@ DECLARE_COMPONENT( ParticleMassMonitor )
// Algorithm
void ParticleMassMonitor::operator()( LHCb::Particle::Range const& particles ) const {
for ( const auto& p : particles ) {
const double m = p->momentum().M();
const auto daughters = p->daughtersVector();
const double pt = p->pt();
const double mom = p->p();
const double m = p->momentum().M();
const auto daughters = p->daughtersVector();
const double pt = p->pt();
const double mom = p->p();
const double mom_difference = calculate_momentum_difference( daughters );
// Mass plots
auto& mass_histo = m_histograms_mass.at( 0 );
++mass_histo.CandidateMass[m];
++mass_histo.MassvsPt[{pt, m}];
++mass_histo.MassvsP[{mom, m}];
if (m_particles.size() == 3) {
mass_histo.MassvsPdiff[mom_difference / Gaudi::Units::GeV] += m - m_mass; // Momentum in GeV for this plot
}
// Momentum plots
int daughters_counter = 0;
for ( int i = 0; i < 3; i++ ) {
for ( unsigned int i = 0; i < m_particles.size(); i++ ) {
auto& momentum_histo = m_histograms_momentum.at( i );
if ( i == 0 ) {
++momentum_histo.CandidatePt[pt];
++momentum_histo.CandidateMom[mom];
} else {
++momentum_histo.CandidatePt[daughters[daughters_counter]->pt()];
++momentum_histo.CandidateMom[daughters[daughters_counter]->p()];
daughters_counter++;
++momentum_histo.CandidatePt[daughters[i-1]->pt()];
++momentum_histo.CandidateMom[daughters[i-1]->p()];
}
}
// Code for tuple generation
if ( m_generateTuples ) {
// Extra variables to store in the tuples
const auto slopes1 = daughters[0]->slopes();
const auto slopes2 = daughters[1]->slopes();
const auto id = p->particleID().pid();
const auto id1 = daughters[0]->particleID().pid();
const auto id2 = daughters[1]->particleID().pid();
// Creation of the tuple
const auto id = p->particleID().pid();
Tuple particleTuple = nTuple( "Particletuples", "" );
particleTuple->column( m_mother + "_M", m ).ignore();
particleTuple->column( m_mother + "_Pt", pt ).ignore();
particleTuple->column( m_daughter1 + "_Pt", daughters[0]->pt() ).ignore();
particleTuple->column( m_daughter2 + "_Pt", daughters[1]->pt() ).ignore();
particleTuple->column( m_mother + "_ID", id ).ignore();
particleTuple->column( m_daughter1 + "_ID", id1 ).ignore();
particleTuple->column( m_daughter2 + "_ID", id2 ).ignore();
particleTuple->column( m_mother + "_P", mom ).ignore();
particleTuple->column( m_daughter1 + "_P", daughters[0]->p() ).ignore();
particleTuple->column( m_daughter2 + "_P", daughters[1]->p() ).ignore();
particleTuple->column( m_daughter1 + "_PxoverPz", slopes1.X() ).ignore();
particleTuple->column( m_daughter2 + "_PxoverPz", slopes2.X() ).ignore();
particleTuple->column( m_daughter1 + "_PyoverPz", slopes1.Y() ).ignore();
particleTuple->column( m_daughter2 + "_PyoverPz", slopes2.Y() ).ignore();
particleTuple->column( m_particles[0] + "_M", m ).ignore();
particleTuple->column( m_particles[0] + "_Pt", pt ).ignore();
particleTuple->column( m_particles[0] + "_ID", id ).ignore();
particleTuple->column( m_particles[0] + "_P", mom ).ignore();
for ( unsigned int i = 1; i < m_particles.size(); i++ ) {
particleTuple->column( m_particles[i] + "_Pt", daughters[i-1]->pt() ).ignore();
particleTuple->column( m_particles[i] + "_ID", daughters[i-1]->particleID().pid() ).ignore();
particleTuple->column( m_particles[i] + "_P", daughters[i-1]->p() ).ignore();
particleTuple->column( m_particles[i] + "_PxoverPz", daughters[i-1]->slopes().X() ).ignore();
particleTuple->column( m_particles[i] + "_PyoverPz", daughters[i-1]->slopes().Y() ).ignore();
}
if ( m_particles.size() == 3 ){
particleTuple->column( "Daughter_p_plus_p_minus", mom_difference ).ignore();
}
particleTuple->write().ignore();
}
}
}
double ParticleMassMonitor::calculate_momentum_difference( LHCb::Particle::ConstVector const& particles ) const {
if ( particles.size() == 2 ) {
const int charge_p1 = particles[0]->charge();
const int charge_p2 = particles[1]->charge();
const double p1 = particles[0]->p();
const double p2 = particles[1]->p();
const double p_diff = ( charge_p1 > 0 && charge_p2 < 0 ) ? p1 - p2 : p2 - p1;
return p_diff;
} else {
return 0.;
}
}
\ No newline at end of file