All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Attributes | List of all members
fluxr::BooNEInterface Class Reference

#include <BooNEInterface.h>

Inheritance diagram for fluxr::BooNEInterface:
fluxr::FluxInterface

Public Member Functions

 BooNEInterface ()
 
 ~BooNEInterface ()
 
const Long64_t GetEntries ()
 
const int GetRun ()
 
const void SetRun (int run)
 
const float GetPOT ()
 
const TLorentzVector GetNuPosition ()
 
const TLorentzVector GetNuMomentum ()
 
void SetRootFile (TFile *rootFileName)
 
bool FillMCFlux (Long64_t ientry, simb::MCFlux &mcflux)
 

Private Attributes

TTree * fBNBtree
 
TTree * fWindowtree
 
BooNENtuple fBooneNtp
 
BeamNtuple fBeamNtp
 
double fMaxWeight
 
double fMinWeight
 
double fSumWeights
 
Long64_t fNEntries
 
int fRun
 
float fPOT
 
TLorentzVector fNuPos
 
TLorentzVector fNuMom
 

Detailed Description

Definition at line 51 of file BooNEInterface.h.

Constructor & Destructor Documentation

fluxr::BooNEInterface::BooNEInterface ( )

Definition at line 18 of file BooNEInterface.cxx.

19  {
20  }
fluxr::BooNEInterface::~BooNEInterface ( )

Definition at line 22 of file BooNEInterface.cxx.

23  {
24  }

Member Function Documentation

bool fluxr::BooNEInterface::FillMCFlux ( Long64_t  ientry,
simb::MCFlux &  mcflux 
)
virtual

Implements fluxr::FluxInterface.

Definition at line 88 of file BooNEInterface.cxx.

89  {
90  if (!fBNBtree->GetEntry(ientry))
91  return false;
92 
93  mf::LogDebug("BooNEInterface") << "At entry: " << ientry << std::endl;
94 
95  fWindowtree->GetEntry(0);
96 
97 
98  if ( fBooneNtp.ntp == 1 ) {
99  flux.fntype = 12; //nue
100  }
101  else if ( fBooneNtp.ntp == 2 ) {
102  flux.fntype = -12; //nuebar
103  }
104  else if ( fBooneNtp.ntp == 3 ) {
105  flux.fntype = 14; //numu
106  }
107  else if ( fBooneNtp.ntp == 4 ) {
108  flux.fntype = -14; //numubar
109  }
110  else{
111  mf::LogWarning("BooNEInterface") << "Neutrino type not recognized! ntp = " << fBooneNtp.ntp
112  << std::endl;
113  }
114 
115  flux.fnimpwt = fBooneNtp.beamwgt;
116 
117  // Convert to meters
118  double nu_x = fBooneNtp.ini_pos[0][0] / 100;
119  double nu_y = fBooneNtp.ini_pos[0][1] / 100;
120  double nu_z = fBooneNtp.ini_pos[0][2] / 100;
121 
122  double nu_momx = fBooneNtp.ini_mom[0][0];
123  double nu_momy = fBooneNtp.ini_mom[0][1];
124  double nu_momz = fBooneNtp.ini_mom[0][2];
125 
126  double targ_x = fBeamNtp.targ_pos_beam[0] / 100;
127  double targ_y = fBeamNtp.targ_pos_beam[1] / 100;
128  double targ_z = fBeamNtp.targ_pos_beam[2] / 100;
129 
130  double tank_z = (fBeamNtp.tank_pos_beam[2] + fBeamNtp.windowbase[2]) / 100;
131 
132  // Calculate the neutrino x, y, and z position at the window location
133  double gsimple_vtxx = nu_x + (nu_momx / nu_momz) * (tank_z - nu_z) + targ_x;
134  double gsimple_vtxy = nu_y + (nu_momy / nu_momz) * (tank_z - nu_z) + targ_y;
135  double gsimple_vtxz = tank_z+targ_z;
136 
137  // Calculate the distance from the neutrino production point to the window
138  double dist = std::sqrt(std::pow(nu_x-gsimple_vtxx,2) +
139  std::pow(nu_y-gsimple_vtxy,2) +
140  std::pow(nu_z-tank_z,2));
141  flux.fdk2gen = dist;
142 
143  // Calculate the neutrino time at the window
144  float nu_time = fBooneNtp.ini_t[0]; // ns
145  nu_time += dist / (TMath::C() / 1e9);
146 
147  mf::LogDebug("BooNEInterface") << "weight = " << flux.fnimpwt
148  << ", nu time = " << nu_time
149  << ", nu energy = " << fBooneNtp.ini_eng[0]
150  << ", dist = " << dist
151  << std::endl;
152 
153  fNuPos = TLorentzVector(gsimple_vtxx * 100,
154  gsimple_vtxy * 100,
155  gsimple_vtxz * 100,
156  nu_time);
157  fNuMom = TLorentzVector(fBooneNtp.ini_mom[0][0],
158  fBooneNtp.ini_mom[0][1],
159  fBooneNtp.ini_mom[0][2],
160  fBooneNtp.ini_eng[0]);
161 
162  flux.fnenergyn = flux.fnenergyf = fBooneNtp.ini_eng[0];
163 
164  flux.frun = fRun;
165  flux.fevtno = ientry;
166 
167  int npart = fBooneNtp.npart;
168  flux.ftpx = fBooneNtp.ini_mom[npart-2][0]; // npart-2
169  flux.ftpy = fBooneNtp.ini_mom[npart-2][1]; // npart-2
170  flux.ftpz = fBooneNtp.ini_mom[npart-2][2]; // npart-2
171 
172  flux.fvx = fBooneNtp.ini_pos[0][0]; //0
173  flux.fvy = fBooneNtp.ini_pos[0][1]; //0
174  flux.fvz = fBooneNtp.ini_pos[0][2]; //0
175 
176  flux.fpdpx = fBooneNtp.fin_mom[1][0]; //1 final
177  flux.fpdpy = fBooneNtp.fin_mom[1][1]; //1 final
178  flux.fpdpz = fBooneNtp.fin_mom[1][2]; //1 final
179 
180  flux.fpppz = fBooneNtp.ini_mom[1][2]; //1 init
181  double pppx = fBooneNtp.ini_mom[1][0]; //1 init
182  double pppy = fBooneNtp.ini_mom[1][1]; //1 init
183  double apppz = flux.fpppz;
184  if (TMath::Abs(flux.fpppz) < 1.0e-30) apppz = 1.0e-30;
185  flux.fppdxdz = pppx / apppz;
186  flux.fppdydz = pppy / apppz;
187 
188 
189  flux.fppmedium = 0.; // empty
190 
191 
192  int ptype_input = fBooneNtp.id[1];
193  int tptype_input = fBooneNtp.id[npart-2];
194 
195  if (ptype_input != 0) ptype_input = genie::pdg::GeantToPdg(ptype_input);
196  if (tptype_input != 0) tptype_input= genie::pdg::GeantToPdg(tptype_input);
197 
198  flux.fptype = ptype_input;
199  flux.ftptype = tptype_input;
200 
201 
202  /////
203  // Now need to calculate ndecay
204  /////
205 
206  double Nenergy = fBooneNtp.ini_eng[0];
207  double Ndxdz = fBooneNtp.ini_mom[0][0] / fBooneNtp.ini_mom[0][2];
208  double Ndydz = fBooneNtp.ini_mom[0][1] / fBooneNtp.ini_mom[0][2];
209 
210  double ppenergy = fBooneNtp.ini_eng[1];
211  double pdPx = fBooneNtp.fin_mom[1][0];
212  double pdPy = fBooneNtp.fin_mom[1][1];
213  double pdPz = fBooneNtp.fin_mom[1][2];
214 
215  double ppdxdz = fBooneNtp.ini_mom[1][0] / fBooneNtp.ini_mom[1][2];
216  double ppdydz = fBooneNtp.ini_mom[1][1] / fBooneNtp.ini_mom[1][2];
217  double pppz = fBooneNtp.ini_mom[1][2];
218 
219  // Get the neutrino energy in the parent decay cm
220  double parent_mass = std::sqrt(ppenergy * ppenergy -
221  pppz * pppz * (ppdxdz * ppdxdz +
222  ppdydz * ppdydz +
223  1.));
224 
225  double parent_energy = std::sqrt(pdPx * pdPx +
226  pdPy * pdPy +
227  pdPz * pdPz +
228  parent_mass * parent_mass);
229 
230  double gamma = parent_energy / parent_mass;
231  double beta[3];
232  beta[0] = pdPx / parent_energy;
233  beta[1] = pdPy / parent_energy;
234  beta[2] = pdPz / parent_energy;
235 
236  double partial = fBooneNtp.ini_mom[0][2] * gamma *
237  (beta[0] * Ndxdz +
238  beta[1] * Ndydz +
239  beta[2]);
240 
241  double Necm = gamma * Nenergy - partial;
242 
243  if (fBooneNtp.id[1] == 10 && fBooneNtp.ntp == 1) flux.fndecay = 1;
244  else if (fBooneNtp.id[1] == 10 && fBooneNtp.ntp == 2) flux.fndecay = 2;
245  else if (fBooneNtp.id[1] == 10 && fBooneNtp.ntp == 3) flux.fndecay = 3;
246  else if (fBooneNtp.id[1] == 10 && fBooneNtp.ntp == 4) flux.fndecay = 4;
247  else if (fBooneNtp.id[1] == 11 && fBooneNtp.ntp == 3) {
248  //check if it is a two or three body decay
249  if (fabs((parent_mass*parent_mass-0.105658389*0.105658389)/(2.*parent_mass)-Necm)/Necm <= 0.001)
250  //two body decay (numu + mu+)
251  flux.fndecay = 5;
252  else {
253  //three body decay (numu + pi0 + mu+)
254  flux.fndecay = 7;
255  }
256  }
257  else if (fBooneNtp.id[1] == 11 && fBooneNtp.ntp == 1) flux.fndecay = 6;
258  else if (fBooneNtp.id[1] == 12 && fBooneNtp.ntp == 4) {
259  if (fabs((parent_mass*parent_mass-0.105658389*0.105658389)/(2.*parent_mass)-Necm)/Necm <= 0.001) {
260  //two body decay (numu + mu+)
261  flux.fndecay = 8;
262  }
263  else {
264  //three body decay (numu + pi0 + mu+)
265  flux.fndecay = 10;
266  }
267  }
268  else if (fBooneNtp.id[1] == 12 && fBooneNtp.ntp == 2) flux.fndecay = 9;
269  else if (fBooneNtp.id[1] == 5 ) flux.fndecay = 11;
270  else if (fBooneNtp.id[1] == 6 ) flux.fndecay = 12;
271  else if (fBooneNtp.id[1] == 8 ) flux.fndecay = 13;
272  else if (fBooneNtp.id[1] == 9 ) flux.fndecay = 14;
273 
274  /////
275  // End calculation of ndecay
276  /////
277 
278  double mupare;
279  double muparpx;
280  double muparpy;
281  double muparpz;
282 
283  if ( fBooneNtp.id[1] == 5 || fBooneNtp.id[1] == 6) {
284  mupare = fBooneNtp.ini_eng[2];
285  muparpx = fBooneNtp.fin_mom[2][0];
286  muparpy = fBooneNtp.fin_mom[2][1];
287  muparpz = fBooneNtp.fin_mom[2][2];
288  } else {
289  mupare = -9999.;
290  muparpx = -9999.;
291  muparpy = -9999.;
292  muparpz = -9999.;
293  }
294 
295  flux.fmuparpx = muparpx;
296  flux.fmuparpy = muparpy;
297  flux.fmuparpz = muparpz;
298  flux.fmupare = mupare;
299 
300  flux.fnwtnear = flux.fnwtfar = 1;
301 
302  return true;
303  }
BooNENtuple::ini_t
float ini_t[20]
E of particle in before chain, array length &#39;npart&#39;.
Definition: BooNEInterface.h:27
BooNENtuple::npart
int npart
1,2,3,4: nue, nuebar, numu, numubar
Definition: BooNEInterface.h:20
BooNENtuple::ini_pos
float ini_pos[20][3]
id of each particle in before chain, array length &#39;npart&#39;
Definition: BooNEInterface.h:24
BeamNtuple::windowbase
float windowbase[3]
pot in the file
Definition: BooNEInterface.h:42
BooNENtuple::ini_eng
float ini_eng[20]
3-mom of particle in before chain, array length &#39;npart&#39;
Definition: BooNEInterface.h:26
BooNENtuple::fin_mom
float fin_mom[20][3]
Definition: BooNEInterface.h:29
BooNENtuple::ini_mom
float ini_mom[20][3]
3-pos of particle in before chain, array length &#39;npart&#39;
Definition: BooNEInterface.h:25
BeamNtuple::tank_pos_beam
float tank_pos_beam[3]
Definition: BooNEInterface.h:38
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:563
BeamNtuple::targ_pos_beam
float targ_pos_beam[3]
3-position of the flux window;
Definition: BooNEInterface.h:39
fluxr::BooNEInterface::fNuMom
TLorentzVector fNuMom
Definition: BooNEInterface.h:80
BooNENtuple::ntp
int ntp
Magic Weight: CrossSect * BooBeamNT.
Definition: BooNEInterface.h:19
fluxr::BooNEInterface::fNuPos
TLorentzVector fNuPos
Definition: BooNEInterface.h:79
const Long64_t fluxr::BooNEInterface::GetEntries ( )
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 57 of file BooNEInterface.h.

57 {return fNEntries;};
const TLorentzVector fluxr::BooNEInterface::GetNuMomentum ( )
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 62 of file BooNEInterface.h.

62 {return fNuMom;};
fluxr::BooNEInterface::fNuMom
TLorentzVector fNuMom
Definition: BooNEInterface.h:80
const TLorentzVector fluxr::BooNEInterface::GetNuPosition ( )
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 61 of file BooNEInterface.h.

61 {return fNuPos;};
fluxr::BooNEInterface::fNuPos
TLorentzVector fNuPos
Definition: BooNEInterface.h:79
const float fluxr::BooNEInterface::GetPOT ( )
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 60 of file BooNEInterface.h.

60 {return fPOT;};
const int fluxr::BooNEInterface::GetRun ( )
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 58 of file BooNEInterface.h.

58 {return fRun;};
void fluxr::BooNEInterface::SetRootFile ( TFile *  rootFileName)

Definition at line 26 of file BooNEInterface.cxx.

27  {
28  fBNBtree = dynamic_cast<TTree*>(fluxInputFile->Get("h201"));
29  fBNBtree->SetBranchAddress("beamwgt",&fBooneNtp.beamwgt);
30  fBNBtree->SetBranchAddress("ntp",&fBooneNtp.ntp);
31  fBNBtree->SetBranchAddress("npart",&fBooneNtp.npart);
32  fBNBtree->SetBranchAddress("id",fBooneNtp.id);
33  fBNBtree->SetBranchAddress("ini_pos",&fBooneNtp.ini_pos[0][0]);
34  fBNBtree->SetBranchAddress("ini_mom",&fBooneNtp.ini_mom[0][0]);
35  fBNBtree->SetBranchAddress("ini_eng",fBooneNtp.ini_eng);
36  fBNBtree->SetBranchAddress("ini_t",fBooneNtp.ini_t);
37  fBNBtree->SetBranchAddress("fin_mom",&fBooneNtp.fin_mom[0][0]);
38  fBNBtree->SetBranchAddress("fin_pol",&fBooneNtp.fin_pol[0][0]);
39 
40  fWindowtree = dynamic_cast<TTree*>(fluxInputFile->Get("h220"));
41  fWindowtree->SetBranchAddress("tank_pos_beam",&fBeamNtp.tank_pos_beam[0]);
42  fWindowtree->SetBranchAddress("targ_pos_beam",&fBeamNtp.targ_pos_beam[0]);
43  fWindowtree->SetBranchAddress("pot",&fBeamNtp.pot);
44  fWindowtree->SetBranchAddress("windowbase",&fBeamNtp.windowbase[0]);
45  fWindowtree->SetBranchAddress("windowdir1",&fBeamNtp.windowdir1[0]);
46  fWindowtree->SetBranchAddress("windowdir2",&fBeamNtp.windowdir2[0]);
47 
48 
49  fNEntries = fBNBtree->GetEntries();
50  mf::LogInfo("BooNEInterface") << "Reading " << fNEntries << " events" << std::endl;
51 
52  fMaxWeight = 0;
53  fMinWeight = 100;
54  fSumWeights = 0;
55 
56  std::string run_number = fluxInputFile->GetName();
57  std::regex r("\\d\\d\\d\\d");
58  std::smatch m;
59  std::regex_search(run_number,m,r);
60 
61  std::istringstream buffer(m[0].str().c_str());
62  buffer >> fRun;
63 
64 
65  for ( int iev = 0 ; iev < fNEntries ; iev++ ) {
66  fBNBtree->GetEntry(iev);
67  fSumWeights += fBooneNtp.beamwgt;
68  if( fMaxWeight < fBooneNtp.beamwgt) fMaxWeight = fBooneNtp.beamwgt;
69  if( fMinWeight > fBooneNtp.beamwgt) fMinWeight = fBooneNtp.beamwgt;
70  }
71 
72  mf::LogInfo("BooNEInterface") << "Min Weight: " << fMinWeight
73  << ", Max Weight: " << fMaxWeight
74  << ", Sum Weights: " << fSumWeights
75  << std::endl;
76 
77 
78  fWindowtree->GetEntry(0);
79  fPOT = fBeamNtp.pot;
80 
81  std::cout << "Window:" << std::endl;
82  std::cout << "\tBase: " << fBeamNtp.windowbase[0] << ", " << fBeamNtp.windowbase[1] << ", " << fBeamNtp.windowbase[2] << std::endl;
83  std::cout << "\tDir 1: " << fBeamNtp.windowdir1[0] << ", " << fBeamNtp.windowdir1[1] << ", " << fBeamNtp.windowdir1[2] << std::endl;
84  std::cout << "\tDir 2: " << fBeamNtp.windowdir2[0] << ", " << fBeamNtp.windowdir2[1] << ", " << fBeamNtp.windowdir2[2] << std::endl;
85  }
BooNENtuple::ini_t
float ini_t[20]
E of particle in before chain, array length &#39;npart&#39;.
Definition: BooNEInterface.h:27
BooNENtuple::npart
int npart
1,2,3,4: nue, nuebar, numu, numubar
Definition: BooNEInterface.h:20
BooNENtuple::ini_pos
float ini_pos[20][3]
id of each particle in before chain, array length &#39;npart&#39;
Definition: BooNEInterface.h:24
BeamNtuple::windowbase
float windowbase[3]
pot in the file
Definition: BooNEInterface.h:42
BooNENtuple::ini_eng
float ini_eng[20]
3-mom of particle in before chain, array length &#39;npart&#39;
Definition: BooNEInterface.h:26
BeamNtuple::pot
float pot
Frame conversion from beam to flux frame.
Definition: BooNEInterface.h:40
gen_crt_frags.m
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
Definition: gen_crt_frags.py:1347
BooNENtuple::fin_mom
float fin_mom[20][3]
Definition: BooNEInterface.h:29
BeamNtuple::windowdir2
float windowdir2[3]
Flux window direction vector 1.
Definition: BooNEInterface.h:44
BeamNtuple::windowdir1
float windowdir1[3]
Flux window base vector.
Definition: BooNEInterface.h:43
BooNENtuple::ini_mom
float ini_mom[20][3]
3-pos of particle in before chain, array length &#39;npart&#39;
Definition: BooNEInterface.h:25
BeamNtuple::tank_pos_beam
float tank_pos_beam[3]
Definition: BooNEInterface.h:38
BeamNtuple::targ_pos_beam
float targ_pos_beam[3]
3-position of the flux window;
Definition: BooNEInterface.h:39
BooNENtuple::fin_pol
float fin_pol[20][3]
final 3-mom of particle in before chain, array length &#39;npart&#39;
Definition: BooNEInterface.h:30
r
esac echo uname r
Definition: condor_lar_pubs.sh:698
BooNENtuple::ntp
int ntp
Magic Weight: CrossSect * BooBeamNT.
Definition: BooNEInterface.h:19
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
const void fluxr::BooNEInterface::SetRun ( int  run)
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 59 of file BooNEInterface.h.

59 {fRun = run;};

Member Data Documentation

BeamNtuple fluxr::BooNEInterface::fBeamNtp
private

Definition at line 71 of file BooNEInterface.h.

TTree* fluxr::BooNEInterface::fBNBtree
private

Definition at line 68 of file BooNEInterface.h.

BooNENtuple fluxr::BooNEInterface::fBooneNtp
private

Definition at line 70 of file BooNEInterface.h.

double fluxr::BooNEInterface::fMaxWeight
private

Definition at line 72 of file BooNEInterface.h.

double fluxr::BooNEInterface::fMinWeight
private

Definition at line 73 of file BooNEInterface.h.

Long64_t fluxr::BooNEInterface::fNEntries
private

Definition at line 76 of file BooNEInterface.h.

TLorentzVector fluxr::BooNEInterface::fNuMom
private

Definition at line 80 of file BooNEInterface.h.

TLorentzVector fluxr::BooNEInterface::fNuPos
private

Definition at line 79 of file BooNEInterface.h.

float fluxr::BooNEInterface::fPOT
private

Definition at line 78 of file BooNEInterface.h.

int fluxr::BooNEInterface::fRun
private

Definition at line 77 of file BooNEInterface.h.

double fluxr::BooNEInterface::fSumWeights
private

Definition at line 74 of file BooNEInterface.h.

TTree* fluxr::BooNEInterface::fWindowtree
private

Definition at line 69 of file BooNEInterface.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5