NuEScatterGen_module.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <fstream>
4 
9 #include "cetlib_except/exception.h"
10 
11 #include "TLorentzVector.h"
12 #include "TRandom3.h"
13 #include "TDatabasePDG.h"
14 #include "TF1.h"
15 
20 
21 namespace evgen {
22  class NuEScatterGen;
23 }
24 
26 public:
27  explicit NuEScatterGen(fhicl::ParameterSet const & p);
28  virtual ~NuEScatterGen();
29 
30  void produce(art::Event & e) override;
31  void beginJob() override;
32  void beginRun(art::Run & run) override;
33  void reconfigure(fhicl::ParameterSet const & p) ;
34 
35  std::vector<simb::MCParticle> GenerateEventKinematics();
36 
37 private:
38 
39  // Energies are in GeV, throughout
40  double fMinEnu;
41  double fMaxEnu;
42 
43  double fGF; // Fermi constant in GeV^-2
44  double fWMA; // sin^2 (Weak mixing angle)
45 
46  // Detector coordinates, in cm
47  double fMinX;
48  double fMaxX;
49  double fMinY;
50  double fMaxY;
51  double fMinZ;
52  double fMaxZ;
53  double fMinT;
54  double fMaxT;
55 
56  double fNueFluxFrac;
57 
58  // Handels the cross section calculation, as a function of outgoing e mom.
59  // Depends on ga, gv, the electron mass, and neutrino energy
60  TF1 *fdsigdT;
61 
62  TRandom3 *fRand;
63 };
64 
65 //------------------------------------------------------------------------------
67 {
68  this->reconfigure(p);
69 
70  produces< std::vector<simb::MCTruth> >();
71  produces< sumdata::RunData, art::InRun >();
72 
73  // There's a prefactor, too, ignoring since I'll just be pulling from dist
74  // [0] = gv [1] = ga [2] = Enu [3] = eMass
75  // ga and gv depend on neutrino flavor
76  fdsigdT = new TF1("xsecform","TMath::Power([0] + [1],2) + TMath::Power([0] - [1],2) * TMath::Power(1-x/[2],2) - ([1]*[1] - [0]*[0] * TMath::Power([3]/[2],2))*x");
77  fRand = new TRandom3(0);
78 }
79 
80 //------------------------------------------------------------------------------
82 {
83 }
84 
85 //------------------------------------------------------------------------------
87 {
88 }
89 
90 //------------------------------------------------------------------------------
92 {
93  // grab the geometry object to see what geometry we are using
95  std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetectorName()));
96 
97  run.put(std::move(runcol));
98 
99  return;
100  }
101 
102 //------------------------------------------------------------------------------
104 {
105  // Set up a vector of ptrs to eventually put into the event
106  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
107  // And we'll fill one
108  simb::MCTruth truth;
109 
110  std::vector<simb::MCParticle> parts = GenerateEventKinematics();
111 
112  assert(parts.size()==2);
113  truth.Add(parts[0]);
114  truth.Add(parts[1]);
115 
116  TVector3 vec0(parts[0].Px(),parts[0].Py(),parts[0].Pz());
117  TVector3 vec1(parts[1].Px(),parts[1].Py(),parts[1].Pz());
118 
119  truthcol->push_back(truth);
120 
121  e.put(std::move(truthcol));
122 
123  return;
124 }
125 
126 //------------------------------------------------------------------------------
128 {
129  fMinEnu = p.get<double>("MinEnu");
130  fMaxEnu = p.get<double>("MaxEnu");
131 
132  fMinX = p.get<double>("MinX");
133  fMaxX = p.get<double>("MaxX");
134  fMinY = p.get<double>("MinY");
135  fMaxY = p.get<double>("MaxY");
136  fMinZ = p.get<double>("MinZ");
137  fMaxZ = p.get<double>("MaxZ");
138  fMinT = p.get<double>("MinT");
139  fMaxT = p.get<double>("MaxT");
140 
141  fNueFluxFrac = p.get<double>("NueFluxFrac");
142 
143  return;
144 }
145 
146 std::vector<simb::MCParticle> evgen::NuEScatterGen::GenerateEventKinematics()
147 {
148  double Enu = fRand->Uniform(fMinEnu,fMaxEnu);
149 
150  int ePDG = 11;
151  static TDatabasePDG pdg;
152  TParticlePDG* pdgp = pdg.GetParticle(ePDG);
153  double eMass = 0;
154  if (pdgp) eMass = pdgp->Mass();
155 
156  // The total cross section for nue-e is A(1-2*ssWMA+4/3*ssWMA*ssWMA)
157  // The total cross section for nue-mu is A(0.25-ssWMA+4/3*ssWMA*ssWMA)
158  // The average nue flux is ~1/3 that of the total flux
159  // Wrap these into a random number pull that decides the type of neutrino
160  assert(fNueFluxFrac >= 0 && fNueFluxFrac <=1);
161  double fracE = fNueFluxFrac * (1-2*fWMA+4./3*fWMA*fWMA) / (fNueFluxFrac * (1-2*fWMA+4./3*fWMA*fWMA) + (1-fNueFluxFrac) * (0.5-fWMA+4./3*fWMA*fWMA));
162  // xsec is same for nutau and numu, for convenience make all neutrinos numu
163  int nuPDG = gRandom->Uniform(0,1) < fracE ? 12 : 14;
164 
165  // ga, gv depend on whether we're scattering nue-e or numu/nutau-e
166  double ga = nuPDG==12 ? 0.5 : -0.5;
167  double gv = nuPDG==12 ? 2*fWMA+0.5 : 2*fWMA-0.5;
168  fdsigdT->SetParameter(0,gv);
169  fdsigdT->SetParameter(1,ga);
170  fdsigdT->SetParameter(2,Enu);
171  fdsigdT->SetParameter(3,eMass);
172  fdsigdT->SetMinimum(0);
173 
174  double tmax = 2*eMass / (TMath::Power(1+eMass/Enu,2) - 1);
175  double eKE = fdsigdT->GetRandom(0,tmax);
176  double eMom = sqrt(eKE*eKE+2*eMass*eKE);
177  std::cout << "Throwing a neutrino of energy " << Enu << " with a maximum electron energy " << tmax << " from cross section " << fdsigdT->GetExpFormula() << ", and got " << eKE << std::endl;
178 
179  // Throw a random neutrino direction
180  double nuCos = fRand->Uniform(-1,1);
181  double nuPhi = gRandom->Uniform(0,2*TMath::Pi());
182  TVector3 nu(sqrt(1-nuCos*nuCos)*sin(nuPhi),
183  sqrt(1-nuCos*nuCos)*cos(nuPhi),
184  nuCos);
185 
186  // It's a bit tricky translating between the electron momentum in the nu
187  // and detector frame. Do it by making a Gram-schmidt basis around the nu
188  // direction, and applying projecting onto it
189  TVector3 zprime = nu.Unit();
190  TVector3 xprime(0,0,0);
191  xprime.SetX(1-nu[0]*nu[0]);
192  xprime.SetY( -nu[0]*nu[1]);
193  xprime.SetZ( -nu[0]*nu[2]);
194  xprime = xprime.Unit();
195  TVector3 yprime(0,0,0);
196  yprime.SetX( -nu[1]*nu[0]-xprime[1]*xprime[0]);
197  yprime.SetY(1-nu[1]*nu[1]-xprime[1]*xprime[1]);
198  yprime.SetZ( -nu[1]*nu[2]-xprime[1]*xprime[2]);
199  yprime = yprime.Unit();
200 
201  double eCos =
202  TMath::Sqrt( (eKE*TMath::Power(1+(eMass/Enu),2) ) / ( 2*eMass + eKE));
203  double ePhi = fRand->Uniform(0,2*TMath::Pi());
204 
205  double eSin = sqrt(1-eCos*eCos);
206  TVector3 e(0,0,0);
207  e.SetX(eCos*zprime[0]+eSin*sin(ePhi)*xprime[0]+eSin*cos(ePhi)*yprime[0]);
208  e.SetY(eCos*zprime[1]+eSin*sin(ePhi)*xprime[1]+eSin*cos(ePhi)*yprime[1]);
209  e.SetZ(eCos*zprime[2]+eSin*sin(ePhi)*xprime[2]+eSin*cos(ePhi)*yprime[2]);
210 
211  e *= eMom;
212  nu *= Enu;
213 
214  TLorentzVector e4d(e,eMass+eKE);
215  TLorentzVector nu4d(nu,Enu);
216 
217  double vtxx = fRand->Uniform(fMinX,fMaxX);
218  double vtxy = fRand->Uniform(fMinY,fMaxY);
219  double vtxz = fRand->Uniform(fMinZ,fMaxZ);
220  double vtxt = fRand->Uniform(fMinT,fMaxT);
221  TLorentzVector vtx(vtxx,vtxy,vtxz,vtxt);
222 
223  // arguments are particle#, pdg code, process, mother, mass, status
224  // status must be 1 to tel GEANT to propogate
225  simb::MCParticle mcNu(0, nuPDG, "primary", 0, 0, 1);
226  mcNu.AddTrajectoryPoint(vtx, e4d);
227  simb::MCParticle mcE(1, 11, "primary", 0, eMass, 1);
228  mcE.AddTrajectoryPoint(vtx, nu4d);
229 
230  std::vector<simb::MCParticle> ret;
231  ret.push_back(mcNu);
232  ret.push_back(mcE);
233 
234  return ret;
235 }
236 
237 
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:257
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
NuEScatterGen(fhicl::ParameterSet const &p)
std::vector< simb::MCParticle > GenerateEventKinematics()
Particle class.
art framework interface to geometry description
Definition: Run.h:17
void reconfigure(fhicl::ParameterSet const &p)
void produce(art::Event &e) override
const double e
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
def move(depos, offset)
Definition: depos.py:107
T get(std::string const &key) const
Definition: ParameterSet.h:271
p
Definition: test.py:223
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
void beginRun(art::Run &run) override
Event generator information.
Definition: MCTruth.h:32
LArSoft geometry interface.
Definition: ChannelGeo.h:16
Event Generation using GENIE, cosmics or single particles.
QTextStream & endl(QTextStream &s)