#include "DK2NuInterface.h"

Inheritance diagram for fluxr::DK2NuInterface:

Public Member Functions
Long64_t	GetEntries () const

int	GetRun () const

float	GetPOT () const

TLorentzVector	GetNuPosition () const

TLorentzVector	GetNuMomentum () const

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	detAV_rand_user

TLorentzVector	fRandUser

TLorentzVector	fRandBeam

Detailed Description

Definition at line 20 of file DK2NuInterface.h.

Member Function Documentation

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

Definition at line 161 of file DK2NuInterface.cxx.

Referenced by GetNuChoice(), and Init().

   {
     //theta,phi in rad
     double xyz[3];
     xyz[0] = TMath::Cos(phi) * TMath::Sin(theta);
     xyz[1] = TMath::Sin(phi) * TMath::Sin(theta);
     xyz[2] = TMath::Cos(theta);
     // condition vector to eliminate most floating point errors
     for (int i = 0; i < 3; ++i) {
       const double eps = 1.0e-15;
       if (TMath::Abs(xyz[i]) < eps) xyz[i] = 0;
       if (TMath::Abs(xyz[i] - 1) < eps) xyz[i] = 1;
       if (TMath::Abs(xyz[i] + 1) < eps) xyz[i] = -1;
     }
     return TVector3(xyz[0], xyz[1], xyz[2]);
   }

void fluxr::DK2NuInterface::Beam2UserP4	(	const TLorentzVector &	beamp4,
		TLorentzVector &	usrp4
	)		const

Definition at line 156 of file DK2NuInterface.cxx.

References fBeamRot.

Referenced by FillMCFlux(), and GetNuChoice().

   {
     usrp4 = fBeamRot * beamp4;
   }

void fluxr::DK2NuInterface::Beam2UserPos	(	const TLorentzVector &	beamxyz,
		TLorentzVector &	usrxyz
	)		const

Definition at line 152 of file DK2NuInterface.cxx.

References fBeamRot, and fBeamZero.

Referenced by FillMCFlux(), and GetNuChoice().

   {
     usrxyz = fBeamRot * beamxyz + fBeamZero;
   }

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

virtual

Implements fluxr::FluxInterface.

Definition at line 178 of file DK2NuInterface.cxx.

References Beam2UserP4(), Beam2UserPos(), simb::MCFlux::fdk2gen, fDk2Nu, fDk2NuTree, simb::MCFlux::fevtno, simb::MCFlux::fmupare, simb::MCFlux::fmuparpx, simb::MCFlux::fmuparpy, simb::MCFlux::fmuparpz, simb::MCFlux::fndecay, simb::MCFlux::fnecm, simb::MCFlux::fnenergyf, simb::MCFlux::fnenergyn, simb::MCFlux::fnimpwt, simb::MCFlux::fntype, fNuChoice, fNuMom, fNuPos, simb::MCFlux::fnwtfar, simb::MCFlux::fnwtnear, simb::MCFlux::fpdpx, simb::MCFlux::fpdpy, simb::MCFlux::fpdpz, simb::MCFlux::fppdxdz, simb::MCFlux::fppdydz, simb::MCFlux::fppmedium, simb::MCFlux::fpppz, simb::MCFlux::fptype, fRandBeam, simb::MCFlux::frun, simb::MCFlux::ftgen, simb::MCFlux::ftgptype, simb::MCFlux::ftptype, simb::MCFlux::ftpx, simb::MCFlux::ftpy, simb::MCFlux::ftpz, simb::MCFlux::fvx, simb::MCFlux::fvy, and simb::MCFlux::fvz.

Referenced by GetNuMomentum().

   {
     if (!fDk2NuTree->GetEntry(ientry)) return false;
 
     // error handling
     // Veto the neutrons and unphysical muon decay
     // Returns default mcflux with nu id = -999
     if (fDk2Nu->decay.ptype == 2112) return true;
     if ((fDk2Nu->decay.ptype == 13 || fDk2Nu->decay.ptype == -13) &&
         (fDk2Nu->decay.ntype == 14 || fDk2Nu->decay.ntype == -14)) {
       // Calculate xnu
       double xnu = (2 * fDk2Nu->decay.necm) / 0.105658389; //
       if (xnu > 1) return true;
     }
 
     // bypass flux window method in favor of random point in detector
     TLorentzVector x4beam = fRandBeam;
 
     // enu = resulting energy when the neutrino is redirected
     // wgt = output probability weight
     double enu, wgt;
     bsim::calcEnuWgt(fDk2Nu, x4beam.Vect(), enu, wgt);
 
     TLorentzVector x4usr;
     Beam2UserPos(x4beam, x4usr);
 
     bsim::NuRay rndnuray = fDk2Nu->nuray[0];
     fDk2Nu->nuray.clear();
 
     TVector3 xyzDk(fDk2Nu->decay.vx, fDk2Nu->decay.vy, fDk2Nu->decay.vz); // origin of decay
     TVector3 p3beam = enu * (x4beam.Vect() - xyzDk).Unit();
 
     // divide output wgt by pi so it is in unit area (output wgt is returned as flux/(pi*m^2))
     wgt *= 1 / TMath::Pi();
 
     // with the recalculated energy, compute the momentum components
     bsim::NuRay anuray(p3beam.x(), p3beam.y(), p3beam.z(), enu, wgt);
     fDk2Nu->nuray.push_back(rndnuray);
     fDk2Nu->nuray.push_back(anuray);
 
     TLorentzVector p4beam(p3beam, enu);
     TLorentzVector p4usr;
     Beam2UserP4(p4beam, p4usr);
 
     fNuPos = TLorentzVector(x4usr);
     fNuMom = TLorentzVector(p4usr);
 
     x4usr.SetX(x4usr.X() / 100.); // from cm --> m
     x4usr.SetY(x4usr.Y() / 100.);
     x4usr.SetZ(x4usr.Z() / 100.);
 
     // overwrite nuchoice
     fNuChoice->clear();
     fNuChoice->pdgNu = fDk2Nu->decay.ntype;
     fNuChoice->xyWgt = fDk2Nu->nuray[1].wgt;
     fNuChoice->impWgt = fDk2Nu->decay.nimpwt;
     fNuChoice->x4NuBeam = x4beam;
     fNuChoice->p4NuBeam = TLorentzVector(
       fDk2Nu->nuray[1].px, fDk2Nu->nuray[1].py, fDk2Nu->nuray[1].pz, fDk2Nu->nuray[1].E);
     //need rotation matrix here to fill these
     fNuChoice->x4NuUser = x4usr;
     fNuChoice->p4NuUser = p4usr;
 
     flux.fntype = fDk2Nu->decay.ntype;
     flux.fnimpwt = fDk2Nu->decay.nimpwt;
     flux.fvx = fDk2Nu->decay.vx;
     flux.fvy = fDk2Nu->decay.vy;
     flux.fvz = fDk2Nu->decay.vz;
     flux.fpdpx = fDk2Nu->decay.pdpx;
     flux.fpdpy = fDk2Nu->decay.pdpy;
     flux.fpdpz = fDk2Nu->decay.pdpz;
     flux.fppdxdz = fDk2Nu->decay.ppdxdz;
     flux.fppdydz = fDk2Nu->decay.ppdydz;
     flux.fpppz = fDk2Nu->decay.pppz;
     flux.fppmedium = fDk2Nu->decay.ppmedium;
     flux.fptype = fDk2Nu->decay.ptype;
     flux.fndecay = fDk2Nu->decay.ndecay;
     flux.fmuparpx = fDk2Nu->decay.muparpx;
     flux.fmuparpy = fDk2Nu->decay.muparpy;
     flux.fmuparpz = fDk2Nu->decay.muparpz;
     flux.fmupare = fDk2Nu->decay.mupare;
     flux.fnecm = fDk2Nu->decay.necm;
 
     flux.ftpx = fDk2Nu->tgtexit.tpx;
     flux.ftpy = fDk2Nu->tgtexit.tpy;
     flux.ftpz = fDk2Nu->tgtexit.tpz;
     flux.ftptype = fDk2Nu->tgtexit.tptype;
     flux.ftgen = fDk2Nu->tgtexit.tgen;
 
     flux.frun = fDk2Nu->job;
     flux.fevtno = fDk2Nu->potnum;
     flux.fnenergyn = flux.fnenergyf = enu;
     flux.fnwtnear = flux.fnwtfar = wgt;
     flux.fdk2gen = (x4beam.Vect() - xyzDk).Mag();
     flux.ftgptype = fDk2Nu->ancestor[1].pdg;
 
     return true;
   }

bsim::Dk2Nu* fluxr::DK2NuInterface::GetDk2Nu ( )

inline

Definition at line 31 of file DK2NuInterface.h.

References fDk2Nu.

31 { return fDk2Nu; };

fluxr::DK2NuInterface::fDk2Nu

bsim::Dk2Nu * fDk2Nu

Definition: DK2NuInterface.h:43

Long64_t fluxr::DK2NuInterface::GetEntries ( ) const

inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 22 of file DK2NuInterface.h.

References fNEntries.

22 { return fNEntries; };

fluxr::DK2NuInterface::fNEntries

Long64_t fNEntries

Definition: DK2NuInterface.h:46

bsim::NuChoice* fluxr::DK2NuInterface::GetNuChoice ( )

inline

Definition at line 32 of file DK2NuInterface.h.

References AnglesToAxis(), Beam2UserP4(), Beam2UserPos(), fNuChoice, Init(), and User2BeamPos().

32 { return fNuChoice; };

fluxr::DK2NuInterface::fNuChoice

bsim::NuChoice * fNuChoice

Definition: DK2NuInterface.h:45

TLorentzVector fluxr::DK2NuInterface::GetNuMomentum ( ) const

inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 26 of file DK2NuInterface.h.

References FillMCFlux(), fNuMom, and SetRootFile().

26 { return fNuMom; };

fluxr::DK2NuInterface::fNuMom

TLorentzVector fNuMom

Definition: DK2NuInterface.h:51

TLorentzVector fluxr::DK2NuInterface::GetNuPosition ( ) const

inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 25 of file DK2NuInterface.h.

References fNuPos.

25 { return fNuPos; };

fluxr::DK2NuInterface::fNuPos

TLorentzVector fNuPos

Definition: DK2NuInterface.h:50

float fluxr::DK2NuInterface::GetPOT ( ) const

inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 24 of file DK2NuInterface.h.

References fPOT.

24 { return fPOT; };

fluxr::DK2NuInterface::fPOT

float fPOT

Definition: DK2NuInterface.h:48

int fluxr::DK2NuInterface::GetRun ( ) const

inlinevirtual

Implements fluxr::FluxInterface.

Definition at line 23 of file DK2NuInterface.h.

References fRun.

23 { return fRun; };

fluxr::DK2NuInterface::fRun

int fRun

Definition: DK2NuInterface.h:47

void fluxr::DK2NuInterface::Init ( fhicl::ParameterSet const & ps )

Definition at line 33 of file DK2NuInterface.cxx.

References AnglesToAxis(), detAV_rand_user, fBeamPosXML, fBeamRot, fBeamRotInv, fBeamRotXML, fBeamZero, fDkMeta, fRandBeam, fRandUser, fTempRot, fhicl::ParameterSet::get(), fhicl::ParameterSet::is_key_to_table(), User2BeamPos(), and w.

Referenced by GetNuChoice().

   {
     // Code to handle flux window adapted from
     // GENIE_R21210/src/FluxDrivers/GNuMIFlux.cxx
     if (ps.is_key_to_table("rotmatrix")) {
       // Matrix defined with series of rotations
       fhicl::ParameterSet rotps = ps.get<fhicl::ParameterSet>("rotmatrix");
       fTempRot.RotateX(rotps.get<double>("x"));
       fTempRot.RotateY(rotps.get<double>("y"));
       fTempRot.RotateZ(rotps.get<double>("z"));
     }
     else {
       std::vector<double> rotmat = ps.get<std::vector<double>>("rotmatrix");
       TVector3 newX, newY, newZ;
 
       if (rotmat.size() == 9) {
         // Matrix defined with new axis values
         newX = TVector3(rotmat[0], rotmat[1], rotmat[2]);
         newY = TVector3(rotmat[3], rotmat[4], rotmat[5]);
         newZ = TVector3(rotmat[6], rotmat[7], rotmat[8]);
       }
       else if (rotmat.size() == 6) {
         // Matrix defined with theta phi rotations
         newX = AnglesToAxis(rotmat[0], rotmat[1]);
         newY = AnglesToAxis(rotmat[2], rotmat[3]);
         newZ = AnglesToAxis(rotmat[4], rotmat[5]);
         fTempRot.RotateAxes(newX, newY, newZ);
         fBeamRotXML = fTempRot;
       }
 
       fTempRot.RotateAxes(newX, newY, newZ);
     }
 
     fBeamRotXML = fTempRot.Inverse();
     fBeamRot = TLorentzRotation(fBeamRotXML);
     fBeamRotInv = fBeamRot.Inverse();
 
     int w = 10, p = 6;
     std::cout << " fBeamRotXML: " << std::setprecision(p) << std::endl;
     std::cout << " [ " << std::setw(w) << fBeamRotXML.XX() << " " << std::setw(w)
               << fBeamRotXML.XY() << " " << std::setw(w) << fBeamRotXML.XZ() << std::endl
               << "   " << std::setw(w) << fBeamRotXML.YX() << " " << std::setw(w)
               << fBeamRotXML.YY() << " " << std::setw(w) << fBeamRotXML.YZ() << std::endl
               << "   " << std::setw(w) << fBeamRotXML.ZX() << " " << std::setw(w)
               << fBeamRotXML.ZY() << " " << std::setw(w) << fBeamRotXML.ZZ() << " ] " << std::endl;
     std::cout << std::endl;
 
     std::vector<double> detxyz = ps.get<std::vector<double>>("userbeam");
     if (detxyz.size() == 3) { fBeamPosXML = TVector3(detxyz[0], detxyz[1], detxyz[2]); }
     else if (detxyz.size() == 6) {
       TVector3 userpos = TVector3(detxyz[0], detxyz[1], detxyz[2]);
       TVector3 beampos = TVector3(detxyz[3], detxyz[4], detxyz[5]);
       fBeamPosXML = userpos - fBeamRotXML * beampos;
     }
     else {
       std::cout << "userbeam needs 3 or 6 numbers to be properly defined" << std::endl;
     }
 
     fBeamZero = TLorentzVector(fBeamPosXML, 0);
 
     w = 16;
     p = 10;
     std::cout << " fBeamPosXML: [ " << std::setprecision(p) << std::setw(w) << fBeamPosXML.X()
               << " , " << std::setw(w) << fBeamPosXML.Y() << " , " << std::setw(w)
               << fBeamPosXML.Z() << " ] " << std::endl;
 
     //
     // calculation from flux window
     // discontinued in favor of random point in active volume
 
     //
 
     // calculation from random point in detector
 
     std::vector<double> x_detAV = ps.get<std::vector<double>>("x_detAV");
     std::vector<double> y_detAV = ps.get<std::vector<double>>("y_detAV");
     std::vector<double> z_detAV = ps.get<std::vector<double>>("z_detAV");
 
     // assign random point
 
     detAV_rand_user[0] = gRandom->Uniform(x_detAV[0], x_detAV[1]);
     detAV_rand_user[1] = gRandom->Uniform(y_detAV[0], y_detAV[1]);
     detAV_rand_user[2] = gRandom->Uniform(z_detAV[0], z_detAV[1]);
 
     // convert to beam coordinates
     fRandUser = TLorentzVector(detAV_rand_user, 0);
     User2BeamPos(fRandUser, fRandBeam);
 
     // this will serve as the input point for calcEnuWgt function
 
 
     //overwrite dkmeta
     fDkMeta->location.clear();
     fDkMeta->location.resize(2);
     fDkMeta->location[0].x = 0;
     fDkMeta->location[0].y = 0;
     fDkMeta->location[0].z = 0;
     TLorentzVector detVec;
     User2BeamPos(TLorentzVector(0, 0, 0, 0), detVec);
     fDkMeta->location[1].x = detVec.Vect().X();
     fDkMeta->location[1].y = detVec.Vect().Y();
     fDkMeta->location[1].z = detVec.Vect().Z();
 
     std::cout << "Locations " << std::endl;
     for (unsigned int i = 0; i < fDkMeta->location.size(); i++) {
       std::cout << i << "\t" << fDkMeta->location[i].x << "\t" << fDkMeta->location[i].y << "\t"
                 << fDkMeta->location[i].z << std::endl;
     }
   }

void fluxr::DK2NuInterface::SetRootFile ( TFile * rootFile )

Definition at line 16 of file DK2NuInterface.cxx.

References fDk2Nu, fDk2NuTree, fDkMeta, fDkMetaTree, fNEntries, fNuChoice, fPOT, and fRun.

Referenced by GetNuMomentum().

   {
     fDk2NuTree = dynamic_cast<TTree*>(rootFile->Get("dk2nuTree"));
     fDkMetaTree = dynamic_cast<TTree*>(rootFile->Get("dkmetaTree"));
 
     fDk2Nu = new bsim::Dk2Nu;
     fDkMeta = new bsim::DkMeta;
     fNuChoice = new bsim::NuChoice;
     fDk2NuTree->SetBranchAddress("dk2nu", &fDk2Nu);
     fDkMetaTree->SetBranchAddress("dkmeta", &fDkMeta);
 
     fNEntries = fDk2NuTree->GetEntries();
     fDkMetaTree->GetEntry(0);
     fRun = fDkMeta->job;
     fPOT = fDkMeta->pots;
   }

void fluxr::DK2NuInterface::User2BeamPos	(	const TLorentzVector &	usrxyz,
		TLorentzVector &	beamxyz
	)		const

Definition at line 148 of file DK2NuInterface.cxx.

References fBeamRotInv, and fBeamZero.

Referenced by GetNuChoice(), and Init().

   {
     beamxyz = fBeamRotInv * (usrxyz - fBeamZero);
   }

Member Data Documentation

TVector3 fluxr::DK2NuInterface::detAV_rand_user

private

Definition at line 57 of file DK2NuInterface.h.

Referenced by Init().

TVector3 fluxr::DK2NuInterface::fBeamPosXML

private

Definition at line 55 of file DK2NuInterface.h.

Referenced by Init().

TLorentzRotation fluxr::DK2NuInterface::fBeamRot

private

Definition at line 54 of file DK2NuInterface.h.

Referenced by Beam2UserP4(), Beam2UserPos(), and Init().

TLorentzRotation fluxr::DK2NuInterface::fBeamRotInv

private

Definition at line 54 of file DK2NuInterface.h.

Referenced by Init(), and User2BeamPos().

TRotation fluxr::DK2NuInterface::fBeamRotXML

private

Definition at line 53 of file DK2NuInterface.h.

Referenced by Init().

TLorentzVector fluxr::DK2NuInterface::fBeamZero

private

Definition at line 56 of file DK2NuInterface.h.

Referenced by Beam2UserPos(), Init(), and User2BeamPos().

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

private

Definition at line 43 of file DK2NuInterface.h.

Referenced by FillMCFlux(), GetDk2Nu(), and SetRootFile().

TTree* fluxr::DK2NuInterface::fDk2NuTree

private

Definition at line 41 of file DK2NuInterface.h.

Referenced by FillMCFlux(), and SetRootFile().

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

private

Definition at line 44 of file DK2NuInterface.h.

Referenced by Init(), and SetRootFile().

TTree* fluxr::DK2NuInterface::fDkMetaTree

private

Definition at line 42 of file DK2NuInterface.h.

Referenced by SetRootFile().

Long64_t fluxr::DK2NuInterface::fNEntries

private

Definition at line 46 of file DK2NuInterface.h.

Referenced by GetEntries(), and SetRootFile().

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

private

Definition at line 45 of file DK2NuInterface.h.

Referenced by FillMCFlux(), GetNuChoice(), and SetRootFile().

TLorentzVector fluxr::DK2NuInterface::fNuMom

private

Definition at line 51 of file DK2NuInterface.h.

Referenced by FillMCFlux(), and GetNuMomentum().

TLorentzVector fluxr::DK2NuInterface::fNuPos

private

Definition at line 50 of file DK2NuInterface.h.

Referenced by FillMCFlux(), and GetNuPosition().

float fluxr::DK2NuInterface::fPOT

private

Definition at line 48 of file DK2NuInterface.h.

Referenced by GetPOT(), and SetRootFile().

TLorentzVector fluxr::DK2NuInterface::fRandBeam

private

Definition at line 58 of file DK2NuInterface.h.

Referenced by FillMCFlux(), and Init().

TLorentzVector fluxr::DK2NuInterface::fRandUser

private

Definition at line 58 of file DK2NuInterface.h.

Referenced by Init().

int fluxr::DK2NuInterface::fRun

private

Definition at line 47 of file DK2NuInterface.h.

Referenced by GetRun(), and SetRootFile().

TRotation fluxr::DK2NuInterface::fTempRot

private

Definition at line 53 of file DK2NuInterface.h.

Referenced by Init().

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

larsim/v09_43_00/source/larsim/PPFXFluxReader/DK2NuInterface.h
larsim/v09_43_00/source/larsim/PPFXFluxReader/DK2NuInterface.cxx

Public Member Functions

Private Attributes

Detailed Description

Member Function Documentation

Member Data Documentation