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

#include <DK2NuInterface.h>

Inheritance diagram for fluxr::DK2NuInterface:
fluxr::FluxInterface

Public Member Functions

 DK2NuInterface ()
 
 ~DK2NuInterface ()
 
const Long64_t GetEntries ()
 
const int GetRun ()
 
const void SetRun (int run)
 
const float GetPOT ()
 
const TLorentzVector GetNuPosition ()
 
const TLorentzVector GetNuMomentum ()
 
void SetRootFile (TFile *rootFile)
 
bool FillMCFlux (Long64_t ientry, simb::MCFlux &mcflux)
 
bsim::Dk2Nu * GetDk2Nu ()
 
bsim::NuChoice * GetNuChoice ()
 
void Init (fhicl::ParameterSet const &ps)
 
void User2BeamPos (const TLorentzVector &usrxyz, TLorentzVector &beamxyz) const
 
void Beam2UserPos (const TLorentzVector &beamxyz, TLorentzVector &usrxyz) const
 
void Beam2UserP4 (const TLorentzVector &beamp4, TLorentzVector &usrp4) const
 
TVector3 AnglesToAxis (double theta, double phi)
 

Private Attributes

TTree * fDk2NuTree
 
TTree * fDkMetaTree
 
bsim::Dk2Nu * fDk2Nu
 
bsim::DkMeta * fDkMeta
 
bsim::NuChoice * fNuChoice
 
Long64_t fNEntries
 
int fRun
 
float fPOT
 
TLorentzVector fNuPos
 
TLorentzVector fNuMom
 
TRotation fBeamRotXML
 
TRotation fTempRot
 
TLorentzRotation fBeamRot
 
TLorentzRotation fBeamRotInv
 
TVector3 fBeamPosXML
 
TLorentzVector fBeamZero
 
TVector3 fFluxWindowPtUser [3]
 
TLorentzVector fFluxWindowBase
 
TLorentzVector fFluxWindowDir1
 
TLorentzVector fFluxWindowDir2
 
TVector3 fWindowNormal
 
Double_t fFluxWindowLen1
 
Double_t fFluxWindowLen2
 
Double_t fWindowArea
 
TRandom3 fRnd
 

Detailed Description

Definition at line 19 of file DK2NuInterface.h.

Constructor & Destructor Documentation

fluxr::DK2NuInterface::DK2NuInterface ( )

Definition at line 13 of file DK2NuInterface.cxx.

14  {
15  }
fluxr::DK2NuInterface::~DK2NuInterface ( )

Definition at line 17 of file DK2NuInterface.cxx.

18  {
19  }

Member Function Documentation

TVector3 fluxr::DK2NuInterface::AnglesToAxis ( double  theta,
double  phi 
)

Definition at line 190 of file DK2NuInterface.cxx.

191  {
192  //theta,phi in rad
193  double xyz[3];
194  xyz[0] = TMath::Cos(phi)*TMath::Sin(theta);
195  xyz[1] = TMath::Sin(phi)*TMath::Sin(theta);
196  xyz[2] = TMath::Cos(theta);
197  // condition vector to eliminate most floating point errors
198  for (int i=0; i<3; ++i) {
199  const double eps = 1.0e-15;
200  if (TMath::Abs(xyz[i]) < eps ) xyz[i] = 0;
201  if (TMath::Abs(xyz[i]-1) < eps ) xyz[i] = 1;
202  if (TMath::Abs(xyz[i]+1) < eps ) xyz[i] = -1;
203  }
204  return TVector3(xyz[0],xyz[1],xyz[2]);
205  }
void fluxr::DK2NuInterface::Beam2UserP4 ( const TLorentzVector &  beamp4,
TLorentzVector &  usrp4 
) const

Definition at line 184 of file DK2NuInterface.cxx.

186  {
187  usrp4 = fBeamRot*beamp4;
188  }
TLorentzRotation fBeamRot
void fluxr::DK2NuInterface::Beam2UserPos ( const TLorentzVector &  beamxyz,
TLorentzVector &  usrxyz 
) const

Definition at line 179 of file DK2NuInterface.cxx.

181  {
182  usrxyz = fBeamRot*beamxyz + fBeamZero;
183  }
TLorentzVector fBeamZero
TLorentzRotation fBeamRot
bool fluxr::DK2NuInterface::FillMCFlux ( Long64_t  ientry,
simb::MCFlux &  mcflux 
)
virtual

Implements fluxr::FluxInterface.

Definition at line 207 of file DK2NuInterface.cxx.

209  {
210  if (!fDk2NuTree->GetEntry(ientry))
211  return false;
212 
213  TLorentzVector x4beam=fFluxWindowBase+fRnd.Uniform()*fFluxWindowDir1+fRnd.Uniform()*fFluxWindowDir2;
214  double enu,wgt;
215  bsim::calcEnuWgt(fDk2Nu, x4beam.Vect(),enu,wgt);
216 
217  TLorentzVector x4usr;
218  Beam2UserPos(x4beam,x4usr);
219 
220  bsim::NuRay rndnuray=fDk2Nu->nuray[0];
221  fDk2Nu->nuray.clear();
222 
223  TVector3 xyzDk(fDk2Nu->decay.vx,fDk2Nu->decay.vy,fDk2Nu->decay.vz); // origin of decay
224  TVector3 p3beam = enu * (x4beam.Vect()-xyzDk).Unit();
225 
226  //weight due to window being tilted with respect to beam direction
227  double tiltwgt = p3beam.Unit().Dot( fWindowNormal );
228  wgt*=tiltwgt;
229  //weight for the window area and divide by pi (since wgt returned by calcEnuWgt function is flux/(pi*m^2)
230  wgt*=fWindowArea/TMath::Pi();//3.14159;
231 
232  bsim::NuRay anuray(p3beam.x(), p3beam.y(), p3beam.z(), enu, wgt);
233  fDk2Nu->nuray.push_back(rndnuray);
234  fDk2Nu->nuray.push_back(anuray);
235 
236  TLorentzVector p4beam(p3beam,enu);
237  TLorentzVector p4usr;
238  Beam2UserP4(p4beam,p4usr);
239 
240  fNuPos=TLorentzVector(x4usr);
241  fNuMom=TLorentzVector(p4usr);
242 
243  x4usr.SetX(x4usr.X()/100.);
244  x4usr.SetY(x4usr.Y()/100.);
245  x4usr.SetZ(x4usr.Z()/100.);
246 
247  fNuChoice->clear();
248  fNuChoice->pdgNu=fDk2Nu->decay.ntype;
249  fNuChoice->xyWgt=fDk2Nu->nuray[1].wgt;
250  fNuChoice->impWgt=fDk2Nu->decay.nimpwt;
251  fNuChoice->x4NuBeam=x4beam;
252  fNuChoice->p4NuBeam=TLorentzVector(fDk2Nu->nuray[1].px,fDk2Nu->nuray[1].py,fDk2Nu->nuray[1].pz,fDk2Nu->nuray[1].E);
253  //need rotation matrix here to fill these
254  fNuChoice->x4NuUser=x4usr;
255  fNuChoice->p4NuUser=p4usr;
256 
257  flux.fntype = fDk2Nu->decay.ntype;
258  flux.fnimpwt = fDk2Nu->decay.nimpwt;
259  flux.fvx = fDk2Nu->decay.vx;
260  flux.fvy = fDk2Nu->decay.vy;
261  flux.fvz = fDk2Nu->decay.vz;
262  flux.fpdpx = fDk2Nu->decay.pdpx;
263  flux.fpdpy = fDk2Nu->decay.pdpy;
264  flux.fpdpz = fDk2Nu->decay.pdpz;
265  flux.fppdxdz = fDk2Nu->decay.ppdxdz;
266  flux.fppdydz = fDk2Nu->decay.ppdydz;
267  flux.fpppz = fDk2Nu->decay.pppz;
268  flux.fppmedium = fDk2Nu->decay.ppmedium;
269  flux.fptype = fDk2Nu->decay.ptype;
270  flux.fndecay = fDk2Nu->decay.ndecay;
271  flux.fmuparpx = fDk2Nu->decay.muparpx;
272  flux.fmuparpy = fDk2Nu->decay.muparpy;
273  flux.fmuparpz = fDk2Nu->decay.muparpz;
274  flux.fmupare = fDk2Nu->decay.mupare;
275  flux.fnecm = fDk2Nu->decay.necm;
276 
277  flux.ftpx = fDk2Nu->tgtexit.tpx;
278  flux.ftpy = fDk2Nu->tgtexit.tpy;
279  flux.ftpz = fDk2Nu->tgtexit.tpz;
280  flux.ftptype = fDk2Nu->tgtexit.tptype;
281  flux.ftgen = fDk2Nu->tgtexit.tgen;
282 
283  flux.frun = fDk2Nu->job;
284  flux.fevtno = fDk2Nu->potnum;
285  flux.fnenergyn = flux.fnenergyf = enu;
286  flux.fnwtnear = flux.fnwtfar = wgt;
287  flux.fdk2gen = (x4beam.Vect()-xyzDk).Mag();
288  flux.ftgptype = fDk2Nu->ancestor[1].pdg;
289 
290  return true;
291  }
TLorentzVector fNuPos
TLorentzVector fFluxWindowDir2
bsim::NuChoice * fNuChoice
TLorentzVector fNuMom
void Beam2UserPos(const TLorentzVector &beamxyz, TLorentzVector &usrxyz) const
TLorentzVector fFluxWindowBase
void Beam2UserP4(const TLorentzVector &beamp4, TLorentzVector &usrp4) const
TLorentzVector fFluxWindowDir1
bsim::Dk2Nu * fDk2Nu
bsim::Dk2Nu* fluxr::DK2NuInterface::GetDk2Nu ( )
inline

Definition at line 35 of file DK2NuInterface.h.

35 {return fDk2Nu;};
bsim::Dk2Nu * fDk2Nu
const Long64_t fluxr::DK2NuInterface::GetEntries ( )
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 25 of file DK2NuInterface.h.

25 {return fNEntries;};
bsim::NuChoice* fluxr::DK2NuInterface::GetNuChoice ( )
inline

Definition at line 36 of file DK2NuInterface.h.

36 {return fNuChoice;};
bsim::NuChoice * fNuChoice
const TLorentzVector fluxr::DK2NuInterface::GetNuMomentum ( )
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 30 of file DK2NuInterface.h.

30 {return fNuMom;};
TLorentzVector fNuMom
const TLorentzVector fluxr::DK2NuInterface::GetNuPosition ( )
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 29 of file DK2NuInterface.h.

29 {return fNuPos;};
TLorentzVector fNuPos
const float fluxr::DK2NuInterface::GetPOT ( )
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 28 of file DK2NuInterface.h.

28 {return fPOT;};
const int fluxr::DK2NuInterface::GetRun ( )
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 26 of file DK2NuInterface.h.

26 {return fRun;};
void fluxr::DK2NuInterface::Init ( fhicl::ParameterSet const &  ps)

Definition at line 38 of file DK2NuInterface.cxx.

39  {
40  //code to handle flux window stolen from GENIE_R21210/src/FluxDrivers/GNuMIFlux.cxx
41  //rotation matrix
42  try {
43  std::vector<double> rotmat=ps.get<std::vector<double> >("rotmatrix");
44  TVector3 newX,newY,newZ;
45  if (rotmat.size()==9) {
46  std::cout<<"FluxReader: Matrix defined with new axis values"<<std::endl;
47  newX=TVector3(rotmat[0], rotmat[1], rotmat[2]);
48  newY=TVector3(rotmat[3], rotmat[4], rotmat[5]);
49  newZ=TVector3(rotmat[6], rotmat[7], rotmat[8]);
50  } else if (rotmat.size()==6) {
51  std::cout<<"FluxReader: Matrix defined with theta phi rotations"<<std::endl;
52  newX = AnglesToAxis(rotmat[0],rotmat[1]);
53  newY = AnglesToAxis(rotmat[2],rotmat[3]);
54  newZ = AnglesToAxis(rotmat[4],rotmat[5]);
55  fTempRot.RotateAxes(newX,newY,newZ);
56  fBeamRotXML = fTempRot; //.Inverse();
57  }
58  fTempRot.RotateAxes(newX,newY,newZ);
59  fBeamRotXML = fTempRot.Inverse();
60  fBeamRot = TLorentzRotation(fBeamRotXML);
61  fBeamRotInv = fBeamRot.Inverse();
62  } catch (std::exception& e) {
63  std::cout<<"FluxReader: "<<e.what()<<std::endl;
64  std::cout<<"FluxReader: These are not numbers"<<std::endl;
65  abort();
66  }
67 
68  try {
69  fhicl::ParameterSet rotps=ps.get<fhicl::ParameterSet >("rotmatrix");
70  double xrot=rotps.get<double>("x");
71  double yrot=rotps.get<double>("y");
72  double zrot=rotps.get<double>("z");
73  fTempRot.RotateX(xrot);
74  fTempRot.RotateY(yrot);
75  fTempRot.RotateZ(zrot);
76 
77  std::cout<<"FluxReader: Matrix defined with series of rotations"<<std::endl;
78  fBeamRotXML = fTempRot.Inverse();
79  fBeamRot = TLorentzRotation(fBeamRotXML);
80  fBeamRotInv = fBeamRot.Inverse();
81  } catch (std::exception& e) {
82  std::cout<<"FluxReader: "<<e.what()<<std::endl;
83  std::cout<<"FluxReader: Looks like it was numbers"<<std::endl;
84  }
85 
86  int w=10, p=6;
87  std::cout << "FluxReader: fBeamRotXML: " << std::setprecision(p) << std::endl;
88  std::cout << "FluxReader: [ "
89  << std::setw(w) << fBeamRotXML.XX() << " "
90  << std::setw(w) << fBeamRotXML.XY() << " "
91  << std::setw(w) << fBeamRotXML.XZ() << std::endl
92  << " "
93  << std::setw(w) << fBeamRotXML.YX() << " "
94  << std::setw(w) << fBeamRotXML.YY() << " "
95  << std::setw(w) << fBeamRotXML.YZ() << std::endl
96  << " "
97  << std::setw(w) << fBeamRotXML.ZX() << " "
98  << std::setw(w) << fBeamRotXML.ZY() << " "
99  << std::setw(w) << fBeamRotXML.ZZ() << " ] " << std::endl;
100  std::cout << std::endl;
101 
102  std::vector<double> detxyz=ps.get<std::vector<double> >("userbeam");
103  TVector3 userpos,beampos;
104  if (detxyz.size()==3) {
105  userpos=TVector3(0,0,0);
106  beampos=TVector3(detxyz[0], detxyz[1], detxyz[2]);
107  } else if (detxyz.size()==6) {
108  userpos=TVector3(detxyz[0], detxyz[1], detxyz[2]);
109  beampos=TVector3(detxyz[3], detxyz[4], detxyz[5]);
110  } else {
111  std::cout<<"FluxReader: userbeam needs 3 or 6 numbers to be properly defined"<<std::endl;
112  }
113  fBeamPosXML = userpos - fBeamRotXML*beampos;
114  fBeamZero=TLorentzVector(fBeamPosXML,0);
115 
116  w=16; p=10;
117  std::cout << "FluxReader: fBeamPosXML: [ " << std::setprecision(p)
118  << std::setw(w) << fBeamPosXML.X() << " , "
119  << std::setw(w) << fBeamPosXML.Y() << " , "
120  << std::setw(w) << fBeamPosXML.Z() << " ] "
121  << std::endl;
122 
123  std::vector<double> windowBase=ps.get<std::vector<double> >("windowBase");
124  std::vector<double> window1 =ps.get<std::vector<double> >("window1");
125  std::vector<double> window2 =ps.get<std::vector<double> >("window2");
126 
127  fFluxWindowPtUser[0]=TVector3( windowBase[0], windowBase[1], windowBase[2] );
128  fFluxWindowPtUser[1]=TVector3( window1[0] , window1[1] , window1[2] );
129  fFluxWindowPtUser[2]=TVector3( window2[0] , window2[1] , window2[2] );
130 
131  // convert from user to beam coord and from 3 points to base + 2 directions
132  // apply units conversion
133  TLorentzVector ptbm0, ptbm1, ptbm2;
134  User2BeamPos(TLorentzVector(fFluxWindowPtUser[0],0),ptbm0);
135  User2BeamPos(TLorentzVector(fFluxWindowPtUser[1],0),ptbm1);
136  User2BeamPos(TLorentzVector(fFluxWindowPtUser[2],0),ptbm2);
137 
138  fFluxWindowBase = ptbm0;
139  fFluxWindowDir1 = ptbm1 - ptbm0;
140  fFluxWindowDir2 = ptbm2 - ptbm0;
141 
144  fWindowNormal = fFluxWindowDir1.Vect().Cross(fFluxWindowDir2.Vect()).Unit();
145  //in genie flux driver area is divided out when calculating effective POT
146  //here we will keeep the POT of dk2nu file, so boost up weights
147  //convert to m^2 (since window specified in cm)
148  fWindowArea = fFluxWindowDir1.Vect().Cross(fFluxWindowDir2.Vect()).Mag()/10000.;
149 
150  double dot = fFluxWindowDir1.Dot(fFluxWindowDir2);
151  if ( TMath::Abs(dot) > 1.0e-8 )
152  std::cout << "FluxReader: Dot product between window direction vectors was "
153  << dot << "; please check for orthoganality"<<std::endl;
154 
155  //overwrite dkmeta
156  fDkMeta->location.clear();
157  fDkMeta->location.resize(2);
158  fDkMeta->location[0].x=0;
159  fDkMeta->location[0].y=0;
160  fDkMeta->location[0].z=0;
161  TLorentzVector detVec;
162  User2BeamPos(TLorentzVector(0,0,0,0),detVec);
163  fDkMeta->location[1].x=detVec.Vect().X();
164  fDkMeta->location[1].y=detVec.Vect().Y();
165  fDkMeta->location[1].z=detVec.Vect().Z();
166 
167  std::cout<<"FluxReader: Locations "<<std::endl;
168  for (unsigned int i=0;i<fDkMeta->location.size();i++) {
169  std::cout<<i<<"\t"<<fDkMeta->location[i].x<<"\t"
170  <<fDkMeta->location[i].y<<"\t"
171  <<fDkMeta->location[i].z<<std::endl;
172  }
173  }
process_name window2
constexpr auto dot(Vector const &a, Vector const &b)
Return cross product of two vectors.
bsim::DkMeta * fDkMeta
pdgs p
Definition: selectors.fcl:22
TLorentzVector fBeamZero
TLorentzVector fFluxWindowDir2
void User2BeamPos(const TLorentzVector &usrxyz, TLorentzVector &beamxyz) const
TVector3 AnglesToAxis(double theta, double phi)
TVector3 fFluxWindowPtUser[3]
process_name windowBase
TLorentzRotation fBeamRotInv
do i e
process_name window1
TLorentzVector fFluxWindowBase
TLorentzVector fFluxWindowDir1
BEGIN_PROLOG could also be cout
TLorentzRotation fBeamRot
void fluxr::DK2NuInterface::SetRootFile ( TFile *  rootFile)

Definition at line 21 of file DK2NuInterface.cxx.

22  {
23  fDk2NuTree=dynamic_cast<TTree*>(rootFile->Get("dk2nuTree"));
24  fDkMetaTree=dynamic_cast<TTree*>(rootFile->Get("dkmetaTree"));
25 
26  fDk2Nu = new bsim::Dk2Nu;
27  fDkMeta= new bsim::DkMeta;
28  fNuChoice=new bsim::NuChoice;
29  fDk2NuTree->SetBranchAddress("dk2nu",&fDk2Nu);
30  fDkMetaTree->SetBranchAddress("dkmeta",&fDkMeta);
31 
32  fNEntries=fDk2NuTree->GetEntries();
33  fDkMetaTree->GetEntry(0);
34  fRun=fDkMeta->job;
35  fPOT=fDkMeta->pots;
36  }
bsim::DkMeta * fDkMeta
bsim::NuChoice * fNuChoice
bsim::Dk2Nu * fDk2Nu
const void fluxr::DK2NuInterface::SetRun ( int  run)
inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 27 of file DK2NuInterface.h.

27 {fRun = run;};
void fluxr::DK2NuInterface::User2BeamPos ( const TLorentzVector &  usrxyz,
TLorentzVector &  beamxyz 
) const

Definition at line 174 of file DK2NuInterface.cxx.

176  {
177  beamxyz = fBeamRotInv*(usrxyz-fBeamZero);
178  }
TLorentzVector fBeamZero
TLorentzRotation fBeamRotInv

Member Data Documentation

TVector3 fluxr::DK2NuInterface::fBeamPosXML
private

Definition at line 59 of file DK2NuInterface.h.

TLorentzRotation fluxr::DK2NuInterface::fBeamRot
private

Definition at line 58 of file DK2NuInterface.h.

TLorentzRotation fluxr::DK2NuInterface::fBeamRotInv
private

Definition at line 58 of file DK2NuInterface.h.

TRotation fluxr::DK2NuInterface::fBeamRotXML
private

Definition at line 57 of file DK2NuInterface.h.

TLorentzVector fluxr::DK2NuInterface::fBeamZero
private

Definition at line 60 of file DK2NuInterface.h.

bsim::Dk2Nu* fluxr::DK2NuInterface::fDk2Nu
private

Definition at line 47 of file DK2NuInterface.h.

TTree* fluxr::DK2NuInterface::fDk2NuTree
private

Definition at line 45 of file DK2NuInterface.h.

bsim::DkMeta* fluxr::DK2NuInterface::fDkMeta
private

Definition at line 48 of file DK2NuInterface.h.

TTree* fluxr::DK2NuInterface::fDkMetaTree
private

Definition at line 46 of file DK2NuInterface.h.

TLorentzVector fluxr::DK2NuInterface::fFluxWindowBase
private

Definition at line 62 of file DK2NuInterface.h.

TLorentzVector fluxr::DK2NuInterface::fFluxWindowDir1
private

Definition at line 62 of file DK2NuInterface.h.

TLorentzVector fluxr::DK2NuInterface::fFluxWindowDir2
private

Definition at line 62 of file DK2NuInterface.h.

Double_t fluxr::DK2NuInterface::fFluxWindowLen1
private

Definition at line 64 of file DK2NuInterface.h.

Double_t fluxr::DK2NuInterface::fFluxWindowLen2
private

Definition at line 64 of file DK2NuInterface.h.

TVector3 fluxr::DK2NuInterface::fFluxWindowPtUser[3]
private

Definition at line 61 of file DK2NuInterface.h.

Long64_t fluxr::DK2NuInterface::fNEntries
private

Definition at line 50 of file DK2NuInterface.h.

bsim::NuChoice* fluxr::DK2NuInterface::fNuChoice
private

Definition at line 49 of file DK2NuInterface.h.

TLorentzVector fluxr::DK2NuInterface::fNuMom
private

Definition at line 55 of file DK2NuInterface.h.

TLorentzVector fluxr::DK2NuInterface::fNuPos
private

Definition at line 54 of file DK2NuInterface.h.

float fluxr::DK2NuInterface::fPOT
private

Definition at line 52 of file DK2NuInterface.h.

TRandom3 fluxr::DK2NuInterface::fRnd
private

Definition at line 66 of file DK2NuInterface.h.

int fluxr::DK2NuInterface::fRun
private

Definition at line 51 of file DK2NuInterface.h.

TRotation fluxr::DK2NuInterface::fTempRot
private

Definition at line 57 of file DK2NuInterface.h.

Double_t fluxr::DK2NuInterface::fWindowArea
private

Definition at line 65 of file DK2NuInterface.h.

TVector3 fluxr::DK2NuInterface::fWindowNormal
private

Definition at line 63 of file DK2NuInterface.h.


The documentation for this class was generated from the following files: