diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx index 3e19373..77f840a 100644 --- a/src/MCStudies/GenericFlux_Vectors.cxx +++ b/src/MCStudies/GenericFlux_Vectors.cxx @@ -1,246 +1,250 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "GenericFlux_Vectors.h" GenericFlux_Vectors::GenericFlux_Vectors(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile) { // Measurement Details fName = name; eventVariables = NULL; // Define our energy range for flux calcs EnuMin = 0.; EnuMax = 100.; // Arbritrarily high energy limit if (Config::HasPar("EnuMin")) { EnuMin = Config::GetParD("EnuMin"); } if (Config::HasPar("EnuMax")) { EnuMax = Config::GetParD("EnuMax"); } // Set default fitter flags fIsDiag = true; fIsShape = false; fIsRawEvents = false; // This function will sort out the input files automatically and parse all the // inputs,flags,etc. // There may be complex cases where you have to do this by hand, but usually // this will do. Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); eventVariables = NULL; // Setup fDataHist as a placeholder this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1); this->SetupDefaultHist(); fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist); covar = StatUtils::GetInvert(fFullCovar); // 1. The generator is organised in SetupMeasurement so it gives the // cross-section in "per nucleon" units. // So some extra scaling for a specific measurement may be required. For // Example to get a "per neutron" measurement on carbon // which we do here, we have to multiple by the number of nucleons 12 and // divide by the number of neutrons 6. // N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is // often included here in other classes that directly integrate the event // histogram. This method is used here as it now respects EnuMin and EnuMax // correctly. this->fScaleFactor = (this->PredictedEventRate("width") / double(fNEvents)) / this->TotalIntegratedFlux(); + this->fScaleFactor_alttest = + (this->PredictedEventRate("width", 0, 1000) / double(fNEvents)) / + this->TotalIntegratedFlux(); LOG(SAM) << " Generic Flux Scaling Factor = " << fScaleFactor << std::endl; if (fScaleFactor <= 0.0) { ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl; } // Setup our TTrees this->AddEventVariablesToTree(); } void GenericFlux_Vectors::AddEventVariablesToTree() { // Setup the TTree to save everything if (!eventVariables) { Config::Get().out->cd(); eventVariables = new TTree((this->fName + "_VARS").c_str(), (this->fName + "_VARS").c_str()); } LOG(SAM) << "Adding Event Variables" << std::endl; eventVariables->Branch("Mode", &Mode, "Mode/I"); eventVariables->Branch("cc", &cc, "cc/B"); eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I"); eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F"); eventVariables->Branch("tgt", &tgt, "tgt/I"); eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I"); eventVariables->Branch("ELep", &ELep, "ELep/F"); eventVariables->Branch("CosLep", &CosLep, "CosLep/F"); // Basic interaction kinematics eventVariables->Branch("Q2", &Q2, "Q2/F"); eventVariables->Branch("q0", &q0, "q0/F"); eventVariables->Branch("q3", &q3, "q3/F"); eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F"); eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F"); eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F"); eventVariables->Branch("W", &W, "W/F"); eventVariables->Branch("x", &x, "x/F"); eventVariables->Branch("y", &y, "y/F"); // Save outgoing particle vectors eventVariables->Branch("nfsp", &nfsp, "nfsp/I"); eventVariables->Branch("px", px, "px[nfsp]/F"); eventVariables->Branch("py", py, "py[nfsp]/F"); eventVariables->Branch("pz", pz, "pz[nfsp]/F"); eventVariables->Branch("E", E, "E[nfsp]/F"); eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I"); // Event Scaling Information eventVariables->Branch("Weight", &Weight, "Weight/D"); eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/D"); eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/D"); eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D"); + eventVariables->Branch("fScaleFactor_alttest", &fScaleFactor_alttest, "fScaleFactor_alttest/D"); return; } void GenericFlux_Vectors::FillEventVariables(FitEvent *event) { // Reset all Function used to extract any variables of interest to the event Mode = PDGnu = tgt = PDGLep = 0; Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y = -999.9; nfsp = 0; for (int i = 0; i < kMAX; ++i) { px[i] = py[i] = pz[i] = E[i] = -999; pdg[i] = 0; } Weight = InputWeight = RWWeight = 0; partList.clear(); // Now fill the information Mode = event->Mode; cc = (abs(event->Mode) < 30); // Get the incoming neutrino and outgoing lepton FitParticle *nu = event->GetNeutrinoIn(); FitParticle *lep = event->GetHMFSAnyLepton(); PDGnu = nu->fPID; Enu_true = nu->fP.E() / 1E3; tgt = event->fTargetPDG; if (lep != NULL) { PDGLep = lep->fPID; ELep = lep->fP.E() / 1E3; CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect())); // Basic interaction kinematics Q2 = -1 * (nu->fP - lep->fP).Mag2() / 1E6; q0 = (nu->fP - lep->fP).E() / 1E3; q3 = (nu->fP - lep->fP).Vect().Mag() / 1E3; // These assume C12 binding from MINERvA... not ideal Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true); Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true); // Get W_true with assumption of initial state nucleon at rest float m_n = (float)PhysConst::mass_proton; W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n); W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n); x = Q2 / (2 * m_n * q0); y = 1 - ELep / Enu_true; } // Loop over the particles and store all the final state particles in a vector for (UInt_t i = 0; i < event->Npart(); ++i) { bool part_alive = event->PartInfo(i)->fIsAlive and event->PartInfo(i)->Status() == kFinalState; if (!part_alive) continue; partList.push_back(event->PartInfo(i)); } // Save outgoing particle vectors nfsp = (int)partList.size(); for (int i = 0; i < nfsp; ++i) { px[i] = partList[i]->fP.X() / 1E3; py[i] = partList[i]->fP.Y() / 1E3; pz[i] = partList[i]->fP.Z() / 1E3; E[i] = partList[i]->fP.E() / 1E3; pdg[i] = partList[i]->fPID; } // Fill event weights Weight = event->RWWeight * event->InputWeight; RWWeight = event->RWWeight; InputWeight = event->InputWeight; // Fill the eventVariables Tree eventVariables->Fill(); return; }; void GenericFlux_Vectors::Write(std::string drawOpt) { // First save the TTree eventVariables->Write(); // Save Flux and Event Histograms too GetInput()->GetFluxHistogram()->Write(); GetInput()->GetEventHistogram()->Write(); return; } // Override functions which aren't really necessary bool GenericFlux_Vectors::isSignal(FitEvent *event) { (void)event; return true; }; void GenericFlux_Vectors::ScaleEvents() { return; } void GenericFlux_Vectors::ApplyNormScale(float norm) { this->fCurrentNorm = norm; return; } void GenericFlux_Vectors::FillHistograms() { return; } void GenericFlux_Vectors::ResetAll() { eventVariables->Reset(); return; } float GenericFlux_Vectors::GetChi2() { return 0.0; } diff --git a/src/MCStudies/GenericFlux_Vectors.h b/src/MCStudies/GenericFlux_Vectors.h index 237e467..613fe63 100644 --- a/src/MCStudies/GenericFlux_Vectors.h +++ b/src/MCStudies/GenericFlux_Vectors.h @@ -1,98 +1,99 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef GenericFlux_Vectors_H_SEEN #define GenericFlux_Vectors_H_SEEN #include "Measurement1D.h" class GenericFlux_Vectors : public Measurement1D { public: GenericFlux_Vectors(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile); virtual ~GenericFlux_Vectors() {}; //! Grab info from event void FillEventVariables(FitEvent *event); //! Fill Custom Histograms void FillHistograms(); //! ResetAll void ResetAll(); //! Scale void ScaleEvents(); //! Norm void ApplyNormScale(float norm); //! Define this samples signal bool isSignal(FitEvent *nvect); //! Write Files void Write(std::string drawOpt); //! Get Chi2 float GetChi2(); void AddEventVariablesToTree(); private: TTree* eventVariables; std::vector partList; int Mode; bool cc; int PDGnu; int tgt; int PDGLep; float ELep; float CosLep; // Basic interaction kinematics float Q2; float q0; float q3; float Enu_QE; float Enu_true; float Q2_QE; float W_nuc_rest; float W; float x; float y; // Save outgoing particle vectors int nfsp; static const int kMAX = 200; float px[kMAX]; float py[kMAX]; float pz[kMAX]; float E[kMAX]; int pdg[kMAX]; // Basic event info double Weight; double InputWeight; double RWWeight; double fScaleFactor; + double fScaleFactor_alttest; }; #endif