Skip to content

Commit 51f0116

Browse files
author
Sabrina Brickner
committed
New MeVPrtl LLP tools
1 parent b62c15d commit 51f0116

3 files changed

Lines changed: 1207 additions & 0 deletions

File tree

‎sbncode/EventGenerator/MeVPrtl/Tools/LLP/CMakeLists.txt‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
cet_build_plugin( Kaon2LLPFlux art::tool
2+
LIBRARIES
3+
art::Framework_Core
4+
art::Framework_Services_Registry
5+
art_root_io::TFileService_service
6+
art::Framework_Principal
7+
art::Persistency_Common
8+
art::Utilities canvas::canvas
9+
cetlib::cetlib cetlib_except::cetlib_except
10+
ROOT::EGPythia6
11+
CLHEP::CLHEP
12+
sbnobj::Common_EventGen_MeVPrtl
13+
nugen::EventGeneratorBase_GENIE
14+
nusimdata::SimulationBase
15+
nurandom::RandomUtils_NuRandomService_service
16+
GSL::gsl
17+
sbncode_EventGenerator_MeVPrtl_Constants
18+
sbnobj::Common_EventGen_MeVPrtl
19+
)
20+
21+
cet_build_plugin( LLPMakeDecay art::tool
22+
LIBRARIES
23+
art::Framework_Core
24+
art::Framework_Services_Registry
25+
art_root_io::TFileService_service
26+
art::Framework_Principal
27+
art::Persistency_Common
28+
art::Utilities canvas::canvas
29+
cetlib::cetlib cetlib_except::cetlib_except
30+
CLHEP::CLHEP
31+
nurandom::RandomUtils_NuRandomService_service
32+
nusimdata::SimulationBase
33+
GSL::gsl
34+
larcorealg::Geometry
35+
sbncode_EventGenerator_MeVPrtl_Dalitz
36+
sbncode_EventGenerator_MeVPrtl_AnThreeBD
37+
sbncode_EventGenerator_MeVPrtl_Constants
38+
sbnobj::Common_EventGen_MeVPrtl
39+
)
40+
41+
42+
install_headers()
43+
install_fhicl()
44+
install_source()

‎sbncode/EventGenerator/MeVPrtl/Tools/LLP/Kaon2LLPFlux_tool.cc‎

Lines changed: 200 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,200 @@
1+
/**
2+
*
3+
*/
4+
5+
// Framework Includes
6+
#include "art/Framework/Core/EDProducer.h"
7+
#include "art/Framework/Principal/Event.h"
8+
#include "art/Framework/Principal/Handle.h"
9+
#include "art/Framework/Services/Registry/ServiceHandle.h"
10+
#include "art/Persistency/Common/PtrMaker.h"
11+
#include "art/Utilities/ToolMacros.h"
12+
#include "cetlib/cpu_timer.h"
13+
#include "fhiclcpp/ParameterSet.h"
14+
#include "messagefacility/MessageLogger/MessageLogger.h"
15+
#include "CLHEP/Random/RandFlat.h"
16+
17+
// local includes
18+
#include "sbncode/EventGenerator/MeVPrtl/Tools/IMeVPrtlFlux.h"
19+
#include "sbncode/EventGenerator/MeVPrtl/Tools/Constants.h"
20+
21+
#include "sbnobj/Common/EventGen/MeVPrtl/MeVPrtlFlux.h"
22+
#include "sbnobj/Common/EventGen/MeVPrtl/MesonParent.h"
23+
24+
// ROOT
25+
#include "TVector3.h"
26+
27+
// std includes
28+
#include <string>
29+
#include <iostream>
30+
#include <memory>
31+
32+
#include "TDatabasePDG.h"
33+
34+
//------------------------------------------------------------------------------------------------------------------------------------------
35+
// implementation follows
36+
37+
namespace evgen {
38+
namespace ldm {
39+
/**
40+
* @brief Kaon2LLPFlux class definiton
41+
*
42+
* Implementation of Kaon->LLP branching ratio taken from:
43+
* arXiv:1912.07622
44+
*/
45+
class Kaon2LLPFlux : public IMeVPrtlFlux
46+
{
47+
public:
48+
/**
49+
* @brief Constructor
50+
*/
51+
Kaon2LLPFlux(fhicl::ParameterSet const &pset);
52+
53+
/**
54+
* @brief Destructor
55+
*/
56+
~Kaon2LLPFlux();
57+
58+
bool MakeFlux(const simb::MCFlux &flux, MeVPrtlFlux &llp, double &weight) override;
59+
void configure(const fhicl::ParameterSet&) override;
60+
61+
double MaxWeight() override;
62+
63+
private:
64+
// config
65+
double fM; //!< Mass of LLP [GeV]
66+
double fTarget2Absorber;
67+
};
68+
69+
Kaon2LLPFlux::Kaon2LLPFlux(fhicl::ParameterSet const &pset):
70+
IMeVPrtlStage("Kaon2LLPFlux"),
71+
IMeVPrtlFlux(pset)
72+
{
73+
this->configure(pset);
74+
75+
}
76+
77+
//------------------------------------------------------------------------------------------------------------------------------------------
78+
79+
Kaon2LLPFlux::~Kaon2LLPFlux()
80+
{
81+
}
82+
83+
// helper functions
84+
//
85+
// kaon -> lepton + LLP
86+
double llp_momentum(double kaon_mass, double lep_mass, double llp_mass) {
87+
if (kaon_mass - lep_mass < llp_mass) return -1.;
88+
89+
return sqrt(kaon_mass * kaon_mass * kaon_mass * kaon_mass
90+
-2 * kaon_mass * kaon_mass * lep_mass * lep_mass
91+
-2 * kaon_mass * kaon_mass * llp_mass * llp_mass
92+
+ lep_mass * lep_mass * lep_mass * lep_mass
93+
+ llp_mass * llp_mass * llp_mass * llp_mass
94+
-2 * lep_mass * lep_mass * llp_mass * llp_mass) / ( 2 * kaon_mass );
95+
}
96+
97+
double SMKaonBR(int kaon_pdg) {
98+
// The Kaons in Dk2nu file only include those that decay to neutrinos.
99+
//
100+
// We want all kaons -- in order to account for this, we divide by the
101+
// branching-ratio of kaons to neutrinos
102+
//
103+
// Taken from:
104+
// /cvmfs/minerva.opensciencegrid.org/minerva/beamsim/x86_64/geant4/source/particles/hadrons/mesons/src/G4KaonPlus.cc
105+
// /cvmfs/minerva.opensciencegrid.org/minerva/beamsim/x86_64/geant4/source/particles/hadrons/mesons/src/G4KaonZerLong.cc
106+
switch (kaon_pdg) {
107+
case 321:
108+
return 0.6339 /* 5 */ + 0.0559 /* 6 */ + 0.0330 /* 7 */;
109+
case -321:
110+
return 0.6339 /* 8 */ + 0.0559 /* 9 */ + 0.0330 /* 10 */;
111+
case 130:
112+
return 0.2020 /* 1 */ + 0.2020 /* 2 */ + 0.1348 /* 3 */ + 0.1348 /* 4 */;
113+
default:
114+
return -1;
115+
}
116+
}
117+
118+
//------------------------------------------------------------------------------------------------------------------------------------------
119+
void Kaon2LLPFlux::configure(fhicl::ParameterSet const &pset)
120+
{
121+
fM = pset.get<double>("M");
122+
123+
fTarget2Absorber = pset.get<double>("Target2Absorber", 5000);
124+
125+
}
126+
127+
double Kaon2LLPFlux::MaxWeight() {
128+
// Weight comes from the NuMi importance weight -- max is 100 (add in an epsilon)
129+
// Scale by the branching ratios here
130+
return 1;
131+
}
132+
133+
bool Kaon2LLPFlux::MakeFlux(const simb::MCFlux &flux, evgen::ldm::MeVPrtlFlux &llp, double &weight) {
134+
// make the kaon parent
135+
evgen::ldm::MesonParent kaon(flux);
136+
if (abs(kaon.meson_pdg) != 321) return false; // Only take charged kaons
137+
138+
TLorentzVector Beam4 = BeamOrigin();
139+
140+
// get position in detector frame
141+
llp.pos_beamcoord = kaon.pos;
142+
llp.pos = kaon.pos;
143+
llp.pos.Transform(fBeam2Det);
144+
llp.pos += Beam4;
145+
146+
// Branch the parent Kaon Decay
147+
double llp_mass = fM;
148+
double br = 1;
149+
bool is_muon = true;
150+
double lep_mass = is_muon ? Constants::Instance().muon_mass : Constants::Instance().elec_mass;
151+
152+
if (fVerbose) std::cout << "BR: " << br << std::endl;
153+
154+
// ignore if we can't make this llp
155+
// Ignore if branching ratio is exactly 0.
156+
if (br == 0.) return false;
157+
158+
// get the momentum direction in the kaon parent rest frame
159+
double p = llp_momentum(Constants::Instance().kplus_mass, lep_mass, llp_mass);
160+
double e = sqrt(p*p + llp_mass * llp_mass);
161+
// Two-body decays are isotropic
162+
llp.mom = TLorentzVector(p*RandomUnitVector(), e);
163+
164+
// boost to lab frame
165+
TLorentzVector mom = llp.mom;
166+
mom.Boost(kaon.mom.BoostVector());
167+
168+
llp.mom_beamcoord = mom;
169+
// rotate to detector frame
170+
llp.mom = mom;
171+
llp.mom.Transform(fBeam2Det);
172+
173+
llp.mmom_beamcoord = kaon.mom;
174+
// also save the kaon momentum in the detector frame
175+
llp.mmom = kaon.mom;
176+
llp.mmom.Transform(fBeam2Det);
177+
178+
// and save the secondary momentum
179+
llp.sec = llp.mmom - llp.mom;
180+
llp.sec_beamcoord = llp.mmom_beamcoord - llp.mom_beamcoord;
181+
182+
// The weight is the importance weight times the branching-ratio weight
183+
weight = kaon.weight * br / SMKaonBR(kaon.meson_pdg);
184+
185+
llp.mass = fM;
186+
187+
llp.meson_pdg = kaon.meson_pdg;
188+
llp.secondary_pdg = (is_muon ? 13 : 11) * (kaon.meson_pdg > 0 ? 1 : -1);
189+
llp.generator = 1; // kLLP
190+
191+
// equivalent neutrino energy
192+
llp.equiv_enu = EnuLab(flux.fnecm, llp.mmom, llp.pos);
193+
194+
return true;
195+
196+
}
197+
DEFINE_ART_CLASS_TOOL(Kaon2LLPFlux)
198+
199+
} // namespace ldm
200+
} // namespace evgen

0 commit comments

Comments
 (0)