20 #include "TLorentzVector.h" 22 #include "TMultiLayerPerceptron.h" 27 #include "DAQChannelMap/DAQChannelMapConstants.h" 34 #include "Utilities/AssociationUtil.h" 161 const double xmass[8] = {0.000511, 0, 0.105658, 0.1349766, 0, 0.938272, 0.939565379, 0.139570};
169 produces<std::vector<jmshower::EID> >();
170 produces< art::Assns<jmshower::EID, rb::Cluster> >();
191 fnameOs +=
"extrahist_fd.root";
194 fnameOs +=
"extrahist_ndos.root";
197 fnameOs +=
"extrahist_ndos.root";
201 <<
" Bad detector id = " << detid;
214 TFile hfileOs(fnameOs.c_str());
215 ht_nue_rw = (TH1D*)hfileOs.Get(
"ht_nue_rw");
222 fnameWeightsEvt +=
"weights_evt_fd.txt";
225 fnameWeightsEvt +=
"weights_evt_ndos.txt";
228 fnameWeightsEvt +=
"weights_evt_ndos.txt";
232 <<
" Bad detector id = " << detid;
249 TTree *trainTree_evt =
new TTree(
"trainTree_evt",
"m_trainTree_evt");
250 trainTree_evt->Branch(
"egLLL", &egLLL,
"egLLL/D");
251 trainTree_evt->Branch(
"egLLT", &egLLT,
"egLLT/D");
252 trainTree_evt->Branch(
"epi0LLL", &epi0LLL,
"epi0LLL/D");
253 trainTree_evt->Branch(
"epi0LLT", &epi0LLT,
"epi0LLT/D");
254 trainTree_evt->Branch(
"epLLL", &epLLL,
"epLLL/D");
255 trainTree_evt->Branch(
"epLLT", &epLLT,
"epLLT/D");
256 trainTree_evt->Branch(
"enLLL", &enLLL,
"enLLL/D");
257 trainTree_evt->Branch(
"enLLT", &enLLT,
"enLLT/D");
258 trainTree_evt->Branch(
"epiLLL", &epiLLL,
"epiLLL/D");
259 trainTree_evt->Branch(
"epiLLT", &epiLLT,
"epiLLT/D");
260 trainTree_evt->Branch(
"pi0mass", &pi0mass,
"pi0mass/D");
261 trainTree_evt->Branch(
"vtxgev", &vtxgev,
"vtxgev/D");
262 trainTree_evt->Branch(
"type", &type,
"type/I");
263 mlp_evt =
new TMultiLayerPerceptron(
"egLLL,egLLT,epi0LLL,epi0LLT,epLLL,epLLT,enLLL,enLLT,epiLLL,epiLLT,pi0mass,vtxgev:18:12:6:type",trainTree_evt);
264 mlp_evt->LoadWeights(fnameWeightsEvt.c_str());
280 fEvent = tfs->
make<TTree>(
"fEvent",
"fEvent");
388 for(
int i=0;
i<3;
i++){
427 if(
fabs(pdg)==12&&ccnc==0&&issig==
true){
430 if(
fabs(pdg)==14&&ccnc==0&&issig==
false) weight = 0.1787396954;
436 double Ee = elecP4.E();
437 double Ehad = etot - Ee;
438 double EhadCor = 0.282525 + 1.0766*Ehad;
439 if(mode!=0&&mode!=1)
return etot;
442 double Enue = (Ee + EhadCor)/0.98343;
471 TLorentzVector evtTrueNuP4(0,0,0,0);
472 TVector3 evtTrueNuVtx(0,0,0);
473 int evtTrueNuCCNC = -1;
474 int evtTrueNuMode = -1;
475 int evtTrueNuPdg = 0;
476 double evtTrueNuEnergy = 0;
511 std::unique_ptr<std::vector<jmshower::EID> > eidcol(
new std::vector<jmshower::EID>);
519 for(
unsigned int sliceIdx = 0; sliceIdx < slicecol->size(); ++sliceIdx){
524 double evtGapTNS = 1e8;
525 for(
unsigned int iic = 0; iic < slicecol->size(); ++iic){
526 if(iic==sliceIdx)
continue;
555 const double wx = (slice.
NXCell() > 0) ? slice.
W(slice.
XCell(0).
get()) : 0;
556 const double wy = (slice.
NYCell() > 0) ? slice.
W(slice.
YCell(0).
get()) : 0;
559 for(
unsigned int i=0;
i<slice.
NCell();
i++){
568 if(rhit.IsCalibrated())gev = rhit.
GeV();
573 std::vector<art::Ptr<jmshower::JMShower> > showercol = jms.at(sliceIdx);
574 int evtNSh = showercol.size();
608 TVector3 v3null(0,0,0);
685 eidcol->push_back(eid);
691 double evtMinX = 99999;
692 double evtMinY = 99999;
693 double evtMinZ = 99999;
694 double evtMaxX = -99999;
695 double evtMaxY = -99999;
696 double evtMaxZ = -99999;
697 for(
unsigned int i=0;
i<showercol.size();
i++){
699 if((showercol[
i]->Start())[0]<evtMinX)evtMinX = (showercol[
i]->Start())[0];
700 if((showercol[
i]->Stop())[0]<evtMinX)evtMinX = (showercol[
i]->Stop())[0];
701 if((showercol[
i]->Start())[1]<evtMinY)evtMinY = (showercol[
i]->Start())[1];
702 if((showercol[
i]->Stop())[1]<evtMinY)evtMinY = (showercol[
i]->Stop())[1];
703 if((showercol[
i]->Start())[2]<evtMinZ)evtMinZ = (showercol[
i]->Start())[2];
704 if((showercol[
i]->Stop())[2]<evtMinZ)evtMinZ = (showercol[
i]->Stop())[2];
706 if((showercol[
i]->Start())[0]>evtMaxX)evtMaxX = (showercol[
i]->Start())[0];
707 if((showercol[
i]->Stop())[0]>evtMaxX)evtMaxX = (showercol[
i]->Stop())[0];
708 if((showercol[
i]->Start())[1]>evtMaxY)evtMaxY = (showercol[
i]->Start())[1];
709 if((showercol[
i]->Stop())[1]>evtMaxY)evtMaxY = (showercol[
i]->Stop())[1];
710 if((showercol[
i]->Start())[2]>evtMaxZ)evtMaxZ = (showercol[
i]->Start())[2];
711 if((showercol[
i]->Stop())[2]>evtMaxZ)evtMaxZ = (showercol[
i]->Stop())[2];
715 std::vector<TLorentzVector> showerpGam;
718 for(
unsigned int i = 0;
i < showercol.size();
i++){
720 double showereraw = showercol[
i]->Energy();;
721 TLorentzVector showerptrk(showercol[
i]->
Dir()[0]*showereraw,showercol[
i]->
Dir()[1]*showereraw,showercol[
i]->
Dir()[2]*showereraw,showereraw);
722 showerpGam.push_back(showerptrk);
730 double sh1eanne = -1;
731 double sh1eannnocos = -1;
732 double sh1eannenocos = -1;
735 double sh1elength = 0;
737 std::vector<TLorentzVector> showerP4Col;
738 std::vector<TLorentzVector> showerP40Col;
739 std::vector<int> showerpid;
741 showerP40Col.clear();
743 double sh1egLLL = -10;
744 double sh1egLLT = -10;
745 double sh1emuLLL = -10;
746 double sh1emuLLT = -10;
747 double sh1epi0LLL = -10;
748 double sh1epi0LLT = -10;
749 double sh1epLLL = -10;
750 double sh1epLLT = -10;
751 double sh1enLLL = -10;
752 double sh1enLLT = -10;
753 double sh1epiLLL = -10;
754 double sh1epiLLT = -10;
755 double sh1vtxgev = 0;
756 double sh1pi0mass = 0;
759 double sh1eLLL = -10;
760 double sh1eLLT = -10;
761 double sh1muLLL = -10;
762 double sh1muLLT = -10;
763 double sh1eelLLL = -10;
764 double sh1eelLLT = -10;
767 int evtNCellToEdge = 10000;
769 for(
unsigned int i = 0;
i < showercol.size();
i++){
770 if(evtNCellToEdge>showercol[
i]->NCellToEdge())evtNCellToEdge=showercol[
i]->NCellToEdge();
771 double showerenergy = showercol[
i]->Energy();
774 double showerann = showercol[
i]->ANN(0);
775 double showeranne = showercol[
i]->ANN(2);
776 double showerannnocos = showercol[
i]->ANN(3);
777 double showerannenocos = showercol[
i]->ANN(4);
779 double ellmax = -9999;
780 double elllmax = -9999;
781 double elltmax = -9999;
784 if(ellmax<showercol[
i]->DedxLongLL(itype)+showercol[
i]->DedxTransLL(itype)){
785 ellmax = showercol[
i]->DedxLongLL(itype)+showercol[
i]->DedxTransLL(itype);
786 elllmax = showercol[
i]->DedxLongLL(itype);
787 elltmax = showercol[
i]->DedxTransLL(itype);
796 double epi0LLL = elllmax - showercol[
i]->DedxLongLL(
int(
jmshower::kPi0));
797 double epi0LLT = elltmax - showercol[
i]->DedxTransLL(
int(
jmshower::kPi0));
808 double vtxgev = showercol[
i]->VtxGeV();
809 double pi0mass =
std::max(showercol[
i]->Pi0Mgg(),0.);
811 if(showercol[
i]->SliceGeV()>1
e-10)shE=showercol[
i]->Energy()/showercol[
i]->SliceGeV();
812 double gap = showercol[
i]->Gap();
823 if(showerann > sh1ann){
826 sh1gev = showerenergy;
830 if(showerenergy > sh1egev){
833 sh1eanne = showeranne;
834 sh1eannnocos = showerannnocos;
835 sh1eannenocos = showerannenocos;
836 sh1egev = showerenergy;
843 sh1epi0LLL = epi0LLL;
844 sh1epi0LLT = epi0LLT;
852 sh1pi0mass = pi0mass;
865 evtetot+=showerenergy;
868 TVector3
p3(showercol[
i]->
Dir()[0]*showerenergy,showercol[
i]->
Dir()[1]*showerenergy,showercol[
i]->
Dir()[2]*showerenergy);
869 double maxll = -9999;
873 if(showercol[
i]->DedxLongLL(itype)+showercol[
i]->DedxTransLL(itype)>maxll){
875 maxll = showercol[
i]->DedxLongLL(itype)+showercol[
i]->DedxTransLL(itype);
880 if(lltype>=0&&lltype<8) m=xmass[lltype];
882 double showerEtot0 =
sqrt(p3.X()*p3.X()+p3.Y()*p3.Y()+p3.Z()*p3.Z());
883 double P3 =
sqrt(showerEtot0*showerEtot0+2*showerEtot0*m);
884 double E = showerEtot0+
m;
885 TLorentzVector showerP40(showercol[
i]->
Dir()[0]*showerEtot0,showercol[
i]->
Dir()[1]*showerEtot0,showercol[
i]->
Dir()[2]*showerEtot0, showerEtot0);
886 TLorentzVector showerP4(showercol[
i]->
Dir()[0]*P3, showercol[
i]->
Dir()[1]*P3, showercol[
i]->
Dir()[2]*P3, E);
889 showerP40Col.push_back(showerP40);
890 showerP4Col.push_back(showerP4);
891 showerpid.push_back(lltype);
900 for(
unsigned int i = 0;
i < showercol.size();
i++){
901 if(
i==(
unsigned int)ish1e)
continue;
902 if(showercol[
i]->
Energy() > sh2gev){
904 sh2gev = showercol[
i]->Energy();
909 TLorentzVector evtSumP4(0,0,0,0);
910 TLorentzVector evtSumP40(0,0,0,0);
912 for(
unsigned int ii=0;ii<showerP4Col.size();ii++){
913 evtSumP4+=showerP4Col[ii];
916 for(
unsigned int ii=0;ii<showerP40Col.size();ii++){
917 evtSumP40+=showerP40Col[ii];
921 double evtSumTheta = evtSumP4.Angle(nuDir);
922 double evtSumPt = evtSumP4.Vect().Mag()*
sin(evtSumTheta);
923 double evtSumP = evtSumP4.Vect().Mag();
924 double evtSumTheta0 = evtSumP40.Angle(nuDir);
925 double evtSumPt0 = evtSumP40.Vect().Mag()*
sin(evtSumTheta0);
926 double evtSumP0 = evtSumP40.Vect().Mag();
929 if(showercol.size()==1){
930 evtetot+=showercol[0]->VtxGeV();
933 bool ismuon = showercol[ish1]->IsMuon();
934 bool ismuone = showercol[ish1e]->IsMuon();
936 if(sh1ann>0.84&&ismuon==
true)
Ncount1++;
939 TVector3 beamDir = evtTrueNuP4.Vect();
942 double sigw =
OscWeight(evtTrueNuPdg,evtTrueNuCCNC, evtTrueNuEnergy,
true);
943 double bkgw =
OscWeight(evtTrueNuPdg,evtTrueNuCCNC, evtTrueNuEnergy,
false);
946 double evtNuEnergy = showercol[ish1e]->NueEnergy();
1014 if(sh1ann>0.84&&sh1gev>0&&ismuon==
false)
Ncount2++;
1022 eid.
SetANNEL(showercol[ish1e]->ANN(5));
1031 eid.
SetVtx(evtTrueNuVtx);
1043 eid.
SetEvtEtot(showercol[ish1e]->SliceEtot());
1072 double evtSh1Sh2Dang = 180*showerP4Col[ish1e].Angle(showerP4Col[ish2].
Vect())/3.14159265;
1097 TVector3 v3null(0,0,0);
1119 eid.
SetDedx0(showercol[ish1e]->PlaneDedx(0));
1120 eid.
SetDedx1(showercol[ish1e]->PlaneDedx(1));
1121 eid.
SetDedx2(showercol[ish1e]->PlaneDedx(2));
1122 eid.
SetDedx3(showercol[ish1e]->PlaneDedx(3));
1123 eid.
SetDedx4(showercol[ish1e]->PlaneDedx(4));
1124 eid.
SetDedx5(showercol[ish1e]->PlaneDedx(5));
1156 eidcol->push_back(eid);
1169 evt.
put(std::move(eidcol));
1170 evt.
put(std::move(assns));
void SetEMuLLT(double emuLLT)
T max(const caf::Proxy< T > &a, T b)
void SetANNEPi0(double ann)
void SetIsMuon(bool ismuon)
SubRunNumber_t subRun() const
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
double m_evtSh1SliceEdge5GeV
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
void SetSh1VtxDoca(double doca)
void SetSh1YNPlane(unsigned int ynplane)
void SetSh1NCellToEdge(int ncell)
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void SetPi0Mass(double pi0mass)
void SetEvtEtot(double etot)
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
fvar< T > fabs(const fvar< T > &x)
void beginRun(art::Run &run)
void SetEPiLLL(double epiLLL)
void SetSh1Gap(double gap)
void SetEPi0LLT(double epi0LLT)
bool LoadHistograms(int detid)
void SetEelLLL(double eelLLL)
unsigned short Plane() const
double DistToEdge(double *point, double detHalfWidth, double detHalfHeight, double detLength)
Find the distance from the given point to the edge of the detector.
void SetSh2TotalLength(double sh1totallength)
void SetSh2Stop(TVector3 sh1stop)
void SetSh1XNPlane(unsigned int xnplane)
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
void SetSh1NCell(unsigned int ncell)
void SetSh1YNCell(unsigned int yncell)
A collection of associated CellHits.
::xsd::cxx::tree::exception< char > exception
void SetNuEnergy(double energy)
void SetSh2NCellToEdge(int ncell)
const PlaneGeo * Plane(unsigned int i) const
void SetVtx(TVector3 vtx)
double m_evtSh1ExclEnergy
DEFINE_ART_MODULE(TestTMapFile)
void SetGLLL(double gLLL)
Horizontal planes which measure Y.
std::string fJMShowerLabel
Label for jmshower::JMShower input.
double m_evtSh1SliceEdge10GeV
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
void SetEPiLLT(double epiLLT)
ProductID put(std::unique_ptr< PROD > &&product)
void SetSh2XNPlane(unsigned int xnplane)
void SetEvtSumP0(double p)
void SetEPi0LLL(double epi0LLL)
void SetEvtGapTNS(double tns)
Far Detector at Ash River, MN.
void SetSh2YNPlane(unsigned int ynplane)
void SetSh2YNCell(unsigned int yncell)
void SetEvtMaxX(double m)
void SetEEdge10Cell(double eedge10cell)
void SetEvtSumP(double p)
void SetSh2NPlane(unsigned int nplane)
NueSel(const fhicl::ParameterSet &pset)
double EnergyEstimator(TVector3 beamDir, TLorentzVector elecP4, double etot, int mode)
void SetEvtMaxZ(double m)
void SetSh2Start(TVector3 sh1start)
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
void SetEvtMinY(double m)
void SetANNEPi0EL(double ann)
void SetSh2Dir(TVector3 sh1dir)
Prototype Near Detector on the surface at FNAL.
void SetSh1NPlane(unsigned int nplane)
T get(std::string const &key) const
void SetELLL(double eLLL)
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
void SetEGLLT(double egLLT)
void SetVtxGeV(double vtxgev)
Near Detector in the NuMI cavern.
double m_evtSh1SliceEdge2GeV
void SetSh1XNCell(unsigned int xncell)
TMultiLayerPerceptron * mlp_evt
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
void SetSh1Sh2Dang(double sh1sh2dang)
double m_evtTrueEnergy[20]
void SetEvtMaxY(double m)
void SetSh2Energy(double sh1energy)
void SetEPLLL(double epLLL)
void SetEEdge2Cell(double eedge2cell)
void SetEGLLL(double egLLL)
unsigned int NYCell() const
Number of cells in the y-view.
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
EventNumber_t event() const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
void SetANNENoCos(double ann)
void SetEvtNCell(int ncell)
T * make(ARGS...args) const
void SetMuLLT(double muLLT)
void SetSh1Dir(TVector3 sh1dir)
void SetANNEL(double ann)
void SetInvGLLL(double invgLLL)
void SetIsInvPhoton(bool isinvphoton)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
void SetBkgOscW(double weight)
void SetDetEnergy(double energy)
void SetSh1ExclEnergy(double sh1exclenergy)
void SetEvtSumCosTheta(double costheta)
void SetEvtSumPt0(double pt)
unsigned int NXCell() const
Number of cells in the x-view.
void SetSh1Start(TVector3 sh1start)
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
void SetEelLLT(double eelLLT)
void SetMuLLL(double muLLL)
void SetEEdge5Cell(double eedge5cell)
void SetSh2VtxDoca(double doca)
void SetEvtMinZ(double m)
void SetELLT(double eLLT)
void SetSh1Energy(double sh1energy)
double OscWeight(int pdg, int ccnc, double energy, bool issig)
void SetSh2XNCell(unsigned int xncell)
void SetENLLL(double enLLL)
void SetEvtNCellToEdge(int ncell)
void SetEvtMinX(double m)
virtual void reconfigure(const fhicl::ParameterSet &pset)
bool IsNoise() const
Is the noise flag set?
void SetANNE(double anne)
virtual void produce(art::Event &evt)
void SetSh2ExclEnergy(double sh1exclenergy)
void SetANNNoCos(double ann)
void SetSh1Stop(TVector3 sh1stop)
void SetEMuLLL(double emuLLL)
double m_evtSh1DistToEdge
void SetEPLLT(double epLLT)
void SetSh1DistToEdge(double dist)
void SetENLLT(double enLLT)
Encapsulate the geometry of one entire detector (near, far, ndos)
const float Theta() const
void SetEvtSumPt(double pt)
void SetSigOscW(double weight)
void SetSh2NCell(unsigned int ncell)
void SetSh2Gap(double gap)
void SetSh1TotalLength(double sh1totallength)