LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
MCFlux.cxx
Go to the documentation of this file.
1 
11 #include <iostream>
12 #include <iomanip>
13 #include <cmath>
14 
15 namespace simb{
16 
17  //......................................................................
19  : frun(-999)
20  , fevtno(-999)
21  , fndxdz(-999.)
22  , fndydz(-999.)
23  , fnpz(-999.)
24  , fnenergy(-999.)
25  , fndxdznea(-999.)
26  , fndydznea(-999.)
27  , fnenergyn(-999.)
28  , fnwtnear(-999.)
29  , fndxdzfar(-999.)
30  , fndydzfar(-999.)
31  , fnenergyf(-999.)
32  , fnwtfar(-999.)
33  , fnorig(-999)
34  , fndecay(-999)
35  , fntype(-999)
36  , fvx(-999.)
37  , fvy(-999.)
38  , fvz(-999.)
39  , fpdpx(-999.)
40  , fpdpy(-999.)
41  , fpdpz(-999.)
42  , fppdxdz(-999.)
43  , fppdydz(-999.)
44  , fpppz(-999.)
45  , fppenergy(-999.)
46  , fppmedium(-999)
47  , fptype(-999) // converted to PDG
48  , fppvx(-999.)
49  , fppvy(-999.)
50  , fppvz(-999.)
51  , fmuparpx(-999.)
52  , fmuparpy(-999.)
53  , fmuparpz(-999.)
54  , fmupare(-999.)
55  , fnecm(-999.)
56  , fnimpwt(-999.)
57  , fxpoint(-999.)
58  , fypoint(-999.)
59  , fzpoint(-999.)
60  , ftvx(-999.)
61  , ftvy(-999.)
62  , ftvz(-999.)
63  , ftpx(-999.)
64  , ftpy(-999.)
65  , ftpz(-999.)
66  , ftptype(-999) // converted to PDG
67  , ftgen(-999)
68  , ftgptype(-999) // converted to PDG
69  , ftgppx(-999.)
70  , ftgppy(-999.)
71  , ftgppz(-999.)
72  , ftprivx(-999.)
73  , ftprivy(-999.)
74  , ftprivz(-999.)
75  , fbeamx(-999.)
76  , fbeamy(-999.)
77  , fbeamz(-999.)
78  , fbeampx(-999.)
79  , fbeampy(-999.)
80  , fbeampz(-999.)
81  , fFluxType(simb::kGenerator)
82  , fgenx(-999.)
83  , fgeny(-999.)
84  , fgenz(-999.)
85  , fdk2gen(-999.)
86  , fgen2vtx(-999.)
87  {
88 
89  for (int i=0; i<6; ++i) fFluxGen[i]=fFluxPos[i]=fFluxNeg[i]= 0;
90 
91  }
92 
93  //......................................................................
95  {
96  frun = -999;
97  fevtno = -999;
98  fndxdz = -999.;
99  fndydz = -999.;
100  fnpz = -999.;
101  fnenergy = -999.;
102  fndxdznea = -999.;
103  fndydznea = -999.;
104  fnenergyn = -999.;
105  fnwtnear = -999.;
106  fndxdzfar = -999.;
107  fndydzfar = -999.;
108  fnenergyf = -999.;
109  fnwtfar = -999.;
110  fnorig = -999;
111  fndecay = -999;
112  fntype = -999;
113  fvx = -999.;
114  fvy = -999.;
115  fvz = -999.;
116  fpdpx = -999.;
117  fpdpy = -999.;
118  fpdpz = -999.;
119  fppdxdz = -999.;
120  fppdydz = -999.;
121  fpppz = -999.;
122  fppenergy = -999.;
123  fppmedium = -999;
124  fptype = -999; // converted to PDG
125  fppvx = -999.;
126  fppvy = -999.;
127  fppvz = -999.;
128  fmuparpx = -999.;
129  fmuparpy = -999.;
130  fmuparpz = -999.;
131  fmupare = -999.;
132  fnecm = -999.;
133  fnimpwt = -999.;
134  fxpoint = -999.;
135  fypoint = -999.;
136  fzpoint = -999.;
137  ftvx = -999.;
138  ftvy = -999.;
139  ftvz = -999.;
140  ftpx = -999.;
141  ftpy = -999.;
142  ftpz = -999.;
143  ftptype = -999; // converted to PDG
144  ftgen = -999;
145  ftgptype = -999; // converted to PDG
146  ftgppx = -999.;
147  ftgppy = -999.;
148  ftgppz = -999.;
149  ftprivx = -999.;
150  ftprivy = -999.;
151  ftprivz = -999.;
152  fbeamx = -999.;
153  fbeamy = -999.;
154  fbeamz = -999.;
155  fbeampx = -999.;
156  fbeampy = -999.;
157  fbeampz = -999.;
158 
159  fgenx = -999.;
160  fgeny = -999.;
161  fgenz = -999.;
162  fdk2gen = -999.;
163  fgen2vtx = -999.;
164 
165  return;
166  }
167 
168  //......................................................................
176  double MCFlux::Flux(int pdgcode, int which) const
177  {
178  const float* flux = 0;
179  if (which==kGenerator) flux = fFluxGen;
180  if (which==kHistPlusFocus) flux = fFluxPos;
181  if (which==kHistMinusFocus) flux = fFluxNeg;
182  if (pdgcode== 12) return flux[0];
183  if (pdgcode==-12) return flux[1];
184  if (pdgcode== 14) return flux[2];
185  if (pdgcode==-14) return flux[3];
186  if (pdgcode== 16) return flux[4];
187  if (pdgcode==-16) return flux[5];
188  return 0.0;
189  }
190 
191  //......................................................................
192 
193  void MCFlux::SetFluxPos(double nue, double nuebar,
194  double numu, double numubar,
195  double nutau,double nutaubar)
196  {
197  fFluxPos[0] = nue; fFluxPos[1] = nuebar;
198  fFluxPos[2] = numu; fFluxPos[3] = numubar;
199  fFluxPos[4] = nutau; fFluxPos[5] = nutaubar;
200  }
201 
202  //......................................................................
203 
204  void MCFlux::SetFluxNeg(double nue, double nuebar,
205  double numu, double numubar,
206  double nutau,double nutaubar)
207  {
208  fFluxNeg[0] = nue; fFluxNeg[1] = nuebar;
209  fFluxNeg[2] = numu; fFluxNeg[3] = numubar;
210  fFluxNeg[4] = nutau; fFluxNeg[5] = nutaubar;
211  }
212 
213  //......................................................................
214  void MCFlux::SetFluxGen(double nue, double nuebar,
215  double numu, double numubar,
216  double nutau,double nutaubar)
217  {
218  fFluxGen[0] = nue; fFluxGen[1] = nuebar;
219  fFluxGen[2] = numu; fFluxGen[3] = numubar;
220  fFluxGen[4] = nutau; fFluxGen[5] = nutaubar;
221  }
222 
223  //......................................................................
224  void MCFlux::ReDecay(double &newE, double &newW, double x, double y, double z)
225  {
226  //note x,y,z are assumed to be in cm
227  //x,y,z are also assumed to be in the beam reference frame
228  //i.e. 0,0,0 are at the target,
229  //z points along the beamline
230  //right handed coords, x points to the left if you are looking down the beam
231 
232  //these constants should probably be defined elsewhere.
233  //I'll put them here until I figure out where they should properly go
234  const double pimass=.13957;
235  const double mumass=0.105658389;
236  const double kmass=0.49368;
237  const double k0mass=0.49767;
238  double mass=pimass;
239 
240  //determine parent particle mass based on fptype PDG code
241  if(fptype==13||fptype==-13){
242  mass = mumass;
243  }
244  else if(fptype==211||fptype==-211){
245  mass = pimass;
246  }
247  else if(fptype==321||fptype==-321){
248  mass = kmass;
249  }
250  else if(fptype==130){
251  mass = k0mass;
252  }
253  else{
254  std::cerr<<"dont know ptype "<<fptype<<"can not compute new decay weight"<<std::endl;
255  return;
256  }
257 
258  //compute kinematics of parent particle at decay point
259  double p=sqrt(1.*fpdpx*fpdpx+1.*fpdpy*fpdpy+1.*fpdpz*fpdpz);
260  double Eplab=sqrt(1.*fpdpx*fpdpx+1.*fpdpy*fpdpy+1.*fpdpz*fpdpz+1.*mass*mass);
261  double gamma = Eplab/mass;
262  double beta = sqrt((gamma*gamma-1)/(gamma*gamma));
263 
264  //compute components of vector between decay point
265  //and the point you're aiming at
266  double rnx=1.*(x-fvx);
267  double rny=1.*(y-fvy);
268  double rnz=1.*(z-fvz);
269  double rn=sqrt(rnx*rnx+rny*rny+rnz*rnz);
270 
271  //compute angle between parent momentum
272  //and where we want the neutrino to go
273  double rndotp = (rnx*fpdpx+rny*fpdpy+rnz*fpdpz);
274  double costhetan = rndotp/(rn*p);
275 
276  //do some checking of the calculation
277  if(std::abs(costhetan)>1){
278  if(costhetan>0){
279  costhetan = 1;
280  }
281  else{
282  costhetan=-1;
283  }
284  }
285 
286  //now compute the weights
287  double MN=1.;
288  if(p>0){//if it didn't decay at rest
289  //boost
290  MN=1./(gamma*(1-beta*costhetan));
291  }
292 
293 
294  newE = MN*fnecm;
295  // std::cout<<"New E "<<newE<<" old ecm "<<fnecm<<" mn "<<MN<<" gamma "<<gamma<<" beta "<<beta<<" costhetan "<<costhetan<<std::endl;
296 
297  //solid angle
298  // small angle approximation: // double san = 10000./(4*rn*rn);
299  // Alex Radovic's removal of small angle approximation
300  const double kRDET = 100.;
301  double san = (1.0-cos(atan( kRDET / rn )))/2.0;
302 
303  // std::cout<<"san "<<san<<" fvz-z "<<fvz-z<<std::endl;
304  newW = san*MN*MN;
305  // std::cout<<"mn "<<MN<<std::endl;
306 
307 
308  //if its a (polarized) muon decay, we have modify the weight
309  if(fptype==13||fptype==-13){
310  double betav[3]={0.};
311  double partialn=0.;
312  double partial=0.;
313  double P_nun[3]={0.};
314  double P_dcm_nun[4]={0.};
315  double P_pcm_mp[4]={0.};
316 
317  betav[0] = fpdpx/Eplab;
318  betav[1] = fpdpy/Eplab;
319  betav[2] = fpdpz/Eplab;
320 
321  P_nun[0] = rnx*newE/rn;
322  P_nun[1] = rny*newE/rn;
323  P_nun[2] = rnz*newE/rn;
324 
325  partialn =gamma*(betav[0]*P_nun[0]+betav[1]*P_nun[1]+betav[2]*P_nun[2]);
326  partialn = newE - partialn /(gamma+1.);
327 
328  P_dcm_nun[0] = P_nun[0] - betav[0]*gamma*partialn;
329  P_dcm_nun[1] = P_nun[1] - betav[1]*gamma*partialn;
330  P_dcm_nun[2] = P_nun[2] - betav[2]*gamma*partialn;
331  P_dcm_nun[3] = sqrt(pow(P_dcm_nun[0],2)
332  +pow(P_dcm_nun[1],2)
333  +pow(P_dcm_nun[2],2));
334 
335 
336  gamma = fppenergy/mass;
337  betav[0] = fppdxdz*fpppz/fppenergy;
338  betav[1] = fppdydz*fpppz/fppenergy;
339  betav[2] = fpppz/fppenergy;
340 
341  partial = gamma*(betav[0]*fmuparpx +
342  +betav[1]*fmuparpy + betav[2]*fmuparpz);
343  partial = fmupare-partial/(gamma+1.);
344  P_pcm_mp[0] = fmuparpx - betav[0]*gamma*partial;
345  P_pcm_mp[1] = fmuparpy - betav[1]*gamma*partial;
346  P_pcm_mp[2] = fmuparpz - betav[2]*gamma*partial;
347  P_pcm_mp[3] = sqrt(pow(P_pcm_mp[0],2)+
348  pow(P_pcm_mp[1],2)+
349  pow(P_pcm_mp[2],2));
350 
351  //calc new decay angle w.r.t. (anti)spin direction
352  double costhn = 0.;
353  if(P_dcm_nun[3]!=0&&P_pcm_mp[3]!=0){
354  costhn = ( P_dcm_nun[0]*P_pcm_mp[0]+
355  P_dcm_nun[1]*P_pcm_mp[1]+
356  P_dcm_nun[2]*P_pcm_mp[2])/(P_dcm_nun[3]*P_pcm_mp[3]);
357  }
358 
359  if(std::abs(costhn)>1){
360  if(costhn>0){
361  costhn = 1;
362  }
363  else{
364  costhn=-1;
365  }
366  }
367 
368  double wt_ration;
369  if(fntype==14||fntype==-14){
370  double xnu = 2.*fnecm/mass;
371  wt_ration = ( (3.-2.*xnu) - (1.-2.*xnu)*costhn ) / (3.-2.*xnu);
372  }
373  else{
374  wt_ration=1.-costhn;
375  }
376 
377  newW*=wt_ration;
378  }
379 
380  return;
381  }
382 
383  //......................................................................
384  std::ostream& operator<< (std::ostream& output, const simb::MCFlux &mcflux)
385  {
386  output
387  << "MCFlux:" << std::endl
388  // 123456789012
389  << " flux job " << std::setw(11) << mcflux.frun << " "
390  << " pot # " << std::setw(11) << mcflux.fevtno << std::endl
391  << " ntype " << std::setw(11) << mcflux.fntype << " "
392  << " ptype " << std::setw(11) << mcflux.fptype << " "
393  << " impwt " << std::setw(11) << mcflux.fnimpwt << std::endl
394  << " ndecay " << std::setw(11) << mcflux.fndecay << " "
395  << " ppmedium " << std::setw(11) << mcflux.fppmedium << " "
396  << " tptype " << std::setw(11) << mcflux.ftptype << std::endl
397  << " vxyz " << std::setw(11) << mcflux.fvx << " "
398  << std::setw(11) << mcflux.fvy << " "
399  << std::setw(11) << mcflux.fvz << std::endl
400  << " pdpxyz " << std::setw(11) << mcflux.fpdpx << " "
401  << std::setw(11) << mcflux.fpdpy << " "
402  << std::setw(11) << mcflux.fpdpz << std::endl
403  << " tpxyz " << std::setw(11) << mcflux.ftpx << " "
404  << std::setw(11) << mcflux.ftpy << " "
405  << std::setw(11) << mcflux.ftpz << std::endl
406  << " ppd[xy]dz " << std::setw(11) << mcflux.fppdxdz << " "
407  << std::setw(11) << mcflux.fppdydz << std::endl
408  << " pppxyz " << std::setw(11) << mcflux.fppdxdz*mcflux.fpppz << " "
409  << std::setw(11) << mcflux.fppdydz*mcflux.fpppz << " "
410  << std::setw(11) << mcflux.fpppz << std::endl
411  << " muparpxyz " << std::setw(11) << mcflux.fmuparpx << " "
412  << std::setw(11) << mcflux.fmuparpy << " "
413  << std::setw(11) << mcflux.fmuparpz << " "
414  << " mupare " << std::setw(11) << mcflux.fmupare << std::endl
415  << " necm " << std::setw(11) << mcflux.fnecm << " "
416  << " dk2gen " << std::setw(11) << mcflux.fdk2gen << std::endl
417  << " near E " << std::setw(11) << mcflux.fnenergyn << " "
418  << " wgt " << std::setw(11) << mcflux.fnwtnear << std::endl
419  << " far E " << std::setw(11) << mcflux.fnenergyf << " "
420  << " wgt " << std::setw(11) << mcflux.fnwtfar << std::endl;
421  return output;
422  }
423 
424 } // end-of-namespace simb
Float_t x
Definition: compare.C:6
double fgenz
Definition: MCFlux.h:102
int frun
Definition: MCFlux.h:35
double fnenergy
Definition: MCFlux.h:40
double ftgppy
Definition: MCFlux.h:86
double fpdpx
Definition: MCFlux.h:55
double ftpx
Definition: MCFlux.h:79
double Flux(int pdgcode, int which=0) const
Definition: MCFlux.cxx:176
double fgeny
Definition: MCFlux.h:101
double ftprivy
Definition: MCFlux.h:89
int fppmedium
Definition: MCFlux.h:62
double fgen2vtx
distance from ray origin to event vtx
Definition: MCFlux.h:104
Flux for negative horn focus.
Definition: MCFlux.h:19
double fbeamx
Definition: MCFlux.h:91
double fbeampy
Definition: MCFlux.h:95
int ftptype
Definition: MCFlux.h:82
Float_t y
Definition: compare.C:6
double fgenx
origin of ray from flux generator
Definition: MCFlux.h:100
Double_t z
Definition: plot.C:276
double fvx
Definition: MCFlux.h:52
double fnwtfar
Definition: MCFlux.h:48
constexpr auto abs(T v)
Returns the absolute value of the argument.
float fFluxNeg[6]
Fluxes as aboce, for negative horn focus.
Definition: MCFlux.h:109
double fppdxdz
Definition: MCFlux.h:58
A bogus flux assumed by the generator.
Definition: MCFlux.h:20
int ftgen
Definition: MCFlux.h:83
double fdk2gen
distance from decay to ray origin
Definition: MCFlux.h:103
double ftprivx
Definition: MCFlux.h:88
object containing MC flux information
Flux for positive horn focus.
Definition: MCFlux.h:18
double fndydzfar
Definition: MCFlux.h:46
int ftgptype
Definition: MCFlux.h:84
void SetFluxNeg(double nue, double nuebar, double numu, double numubar, double nutau, double nutaubar)
Definition: MCFlux.cxx:204
double fnenergyf
Definition: MCFlux.h:47
void Reset()
Definition: MCFlux.cxx:94
double ftprivz
Definition: MCFlux.h:90
double fbeamz
Definition: MCFlux.h:93
double fnwtnear
Definition: MCFlux.h:44
int fptype
Definition: MCFlux.h:63
double ftvx
Definition: MCFlux.h:76
int fndecay
Definition: MCFlux.h:50
double fndydz
Definition: MCFlux.h:38
int fnorig
Definition: MCFlux.h:49
double fbeamy
Definition: MCFlux.h:92
double ftgppz
Definition: MCFlux.h:87
double fpdpz
Definition: MCFlux.h:57
double fmuparpz
Definition: MCFlux.h:69
double ftgppx
Definition: MCFlux.h:85
int fevtno
Definition: MCFlux.h:36
double fppdydz
Definition: MCFlux.h:59
double ftpz
Definition: MCFlux.h:81
double ftvy
Definition: MCFlux.h:77
double fbeampz
Definition: MCFlux.h:96
ART objects.
double fndxdzfar
Definition: MCFlux.h:45
double fpdpy
Definition: MCFlux.h:56
double fnpz
Definition: MCFlux.h:39
double fzpoint
Definition: MCFlux.h:75
double fvy
Definition: MCFlux.h:53
double fmupare
Definition: MCFlux.h:70
double fxpoint
Definition: MCFlux.h:73
void SetFluxGen(double nue, double nuebar, double numu, double numubar, double nutau, double nutaubar)
Definition: MCFlux.cxx:214
double fndxdznea
Definition: MCFlux.h:41
double fppenergy
Definition: MCFlux.h:61
int fntype
Definition: MCFlux.h:51
double fnecm
Definition: MCFlux.h:71
double fndydznea
Definition: MCFlux.h:42
double fypoint
Definition: MCFlux.h:74
double fmuparpy
Definition: MCFlux.h:68
double fppvx
Definition: MCFlux.h:64
double fpppz
Definition: MCFlux.h:60
double fppvy
Definition: MCFlux.h:65
float fFluxGen[6]
Fluxes as above, assumed by generator.
Definition: MCFlux.h:110
double fbeampx
Definition: MCFlux.h:94
double fppvz
Definition: MCFlux.h:66
double ftpy
Definition: MCFlux.h:80
double fndxdz
Definition: MCFlux.h:37
double fvz
Definition: MCFlux.h:54
void SetFluxPos(double nue, double nuebar, double numu, double numubar, double nutau, double nutaubar)
Definition: MCFlux.cxx:193
friend std::ostream & operator<<(std::ostream &output, const simb::MCFlux &mcflux)
Definition: MCFlux.cxx:384
float fFluxPos[6]
e,ebar,mu,mubar,tau,taubar flux, +horn focus
Definition: MCFlux.h:108
double fmuparpx
Definition: MCFlux.h:67
void ReDecay(double &newE, double &newW, double x, double y, double z)
Definition: MCFlux.cxx:224
double fnenergyn
Definition: MCFlux.h:43
double ftvz
Definition: MCFlux.h:78
double fnimpwt
Definition: MCFlux.h:72