diff --git a/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx b/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx index e46ae68..5dd7b87 100644 --- a/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx +++ b/src/MiniBooNE/MiniBooNE_NCEL_XSec_Treco_nu.cxx @@ -1,264 +1,264 @@ // 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 "MiniBooNE_NCEL_XSec_Treco_nu.h" #include "TLorentzVector.h" //******************************************************************** MiniBooNE_NCEL_XSec_Treco_nu::MiniBooNE_NCEL_XSec_Treco_nu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MiniBooNE_NCEL_XSec_Treco_nu sample. \n" \ "Target: CH \n" \ "Flux: MiniBooNE Numu Flux \n" \ "Signal: Any event with True NCEL modes \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("T_{reco} (MeV)"); fSettings.SetYTitle("Events/(12 MeV)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/NORM/MASK", "FIX/FULL,DIAG"); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); fSettings.SetEnuRange(0.0, 10.0); double arr_treco[52] = {40.0, 52.0, 63.9, 75.9, 87.8, 99.8, 111.8, 123.7, 135.7, 147.6, 159.6, 171.6, 183.5, 195.5, \ 207.5, 219.4, 231.4, 243.3, 255.3, 267.3, 279.2, 291.2, 303.1, 315.1, 327.1, 339.0, 351.0, 362.9, 374.9, 386.9, \ 398.8, 410.8, 422.7, 434.7, 446.7, 458.6, 470.6, 482.5, 494.5, 506.5, 518.4, 530.4, 542.4, 554.3, 566.3, 578.2, \ 590.2, 602.2, 614.1, 626.1, 638.0, 650.0}; SetDataValues(FitPar::GetDataBase()+"/MiniBooNE/ncqe/input_data.txt", arr_treco); SetCovarMatrix(FitPar::GetDataBase()+"/MiniBooNE/ncqe/ErrorMatrix.tab", 51); SetResponseMatrix(FitPar::GetDataBase()+"/MiniBooNE/ncqe/response_mat.txt", 51, arr_treco); SetFluxHistogram(FitPar::GetDataBase()+"/MiniBooNE/ncqe/flux.txt"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // The scale factor is quite complicated because MB didn't divide by number of targets. // nMolMB is the number of CH_2 molecules in the MB FV (610.6 cm radius sphere) and 0.845 is the published density of the mineral oil. // UPDATE: nMolMB is the number of CH_2 molecules in the MB FV (500 cm radius sphere) and 0.845 is the published density of the mineral oil. // UPDATE UPDATE: They didn't account for fiducial cut, so neither do we. double nMolMB = 6.023E+23*0.845*4.0*M_PI*610.6*610.6*610.6/3.0; double POT = 6.46e20; // Need to update scalefactor to reflect actual flux used fScaleFactor = (this->newFluxHist->Integral("")*POT*(GetEventHistogram()->Integral("width"))*1E-38*14.08/(fNEvents+0.))*nMolMB; fScaleFactor /= GetFluxHistogram()->Integral("width"); // Final setup --------------------------------------------------- FinaliseMeasurement(); // Usually the MCFine histogram is a finer binned version of MC Hist. // In this case we need to use it to save the true distribution before smearing. if (fMCFine) delete fMCFine; fMCFine = new TH1D((this->fName+"_Ttrue").c_str(), (this->fName+this->fPlotTitles).c_str(), 50, 0, 900); }; void MiniBooNE_NCEL_XSec_Treco_nu::Write(std::string arg){ newFluxHist->Write("MB_NCEL_newFlux"); response_mat->Write("MB_NCEL_response matrix"); Measurement1D::Write(arg); return; } void MiniBooNE_NCEL_XSec_Treco_nu::FillEventVariables(FitEvent *event){ double t_raw = 0.0; // Loop and add all Tnucleon for (UInt_t i = 0; i < event->Npart(); i++){ if (event->PartInfo(i)->Status() != kFinalState) continue; int pdg = event->PartInfo(i)->fPID; if (pdg == 2212 || pdg == 2112){ t_raw += FitUtils::T(event->PartInfo(i)->fP)*1.E3; } } fXVar = t_raw; } void MiniBooNE_NCEL_XSec_Treco_nu::ScaleEvents(){ // Now convert Ttrue to Treco... for (int treco = 0; treco < 51; ++treco){ double total = 0.; for (int ttrue = 0; ttrue < 50; ++ttrue) total += fMCFine->GetBinContent(ttrue+1)*response_mat->GetBinContent(ttrue+1, treco+1); fMCHist->SetBinContent(treco+1, total); } // Scale this->fMCHist->Scale(this->fScaleFactor, "width"); this->fMCFine->Scale(this->fScaleFactor, "width"); PlotUtils::ScaleNeutModeArray((TH1**)fMCHist_PDG, fScaleFactor, "width"); // Add in the backgrounds... for (int treco = 0; treco < 51; ++treco){ double total = this->fMCHist->GetBinContent(treco+1) + this->BKGD_other->GetBinContent(treco+1) + this->BKGD_irrid->GetBinContent(treco+1); this->fMCHist->SetBinContent(treco+1, total); } } bool MiniBooNE_NCEL_XSec_Treco_nu::isSignal(FitEvent *event){ // Should put in MB SignalDef eventually if (event->Mode != 51 && event->Mode != 52) return false; // Numu or nue if (event->PDGnu() != 14 && event->PDGnu() != 12) return false; // Enu if (event->Enu() < EnuMin*1000.0 || event->Enu() > EnuMax*1000.0) return false; return true; }; void MiniBooNE_NCEL_XSec_Treco_nu::SetFluxHistogram(std::string dataFile){ this->newFluxHist = PlotUtils::GetTH1DFromFile(dataFile.c_str(), (this->fName+"Real Flux Hist"), "idgaf"); this->newFluxHist->Print(); return; } // Read in the covariance matrix from the file specified in the constructor void MiniBooNE_NCEL_XSec_Treco_nu::SetCovarMatrix(std::string covarFile, int dim){ // Use Utils // // Make a counter to track the line number int row = 0; std::string line; - std::ifstream covar(covarFile.c_str(),ifstream::in); + std::ifstream covar(covarFile.c_str(),std::ifstream::in); this->covar = new TMatrixDSym(dim); if(covar.is_open()) LOG(DEB) << "Reading covariance matrix from file: " << covarFile << std::endl; while(std::getline(covar >> std::ws, line, '\n')){ std::istringstream stream(line); double entry; int column = 0; // Loop over entries and insert them into matrix // Multiply by the errors to get the covariance, rather than the correlation matrix while(stream >> entry){ (*this->covar)(row, column) = entry; if (row == column) this->fDataHist->SetBinError(row+1, sqrt(entry)); column++; } row++; } // // Robust matrix inversion method TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU .Invert().GetMatrixArray(), ""); return; }; // Override the usual function in the base class because this is more complicated for the NCEL sample... void MiniBooNE_NCEL_XSec_Treco_nu::SetDataValues(std::string inputFile, double* arr_treco){ std::string line; - std::ifstream input(inputFile.c_str(),ifstream::in); + std::ifstream input(inputFile.c_str(),std::ifstream::in); if(input.is_open()) LOG(DEB) << "Reading data from file: " << inputFile << std::endl; this->fDataHist = new TH1D((this->fName+"_data").c_str(), (this->fName+this->fPlotTitles).c_str(), 51, arr_treco); this->BKGD_other = new TH1D((this->fName+"_BKGD_other").c_str(), (this->fName+this->fPlotTitles).c_str(), 51, arr_treco); this->BKGD_irrid = new TH1D((this->fName+"_BKGD_irrid").c_str(), (this->fName+this->fPlotTitles).c_str(), 51, arr_treco); // To get the nDOF correct... this->fNDataPointsX= 52; double entry = 0; int xBin = 0; // First line is the MB data std::getline(input >> std::ws, line, '\n'); std::istringstream stream1(line); while(stream1 >> entry){ this->fDataHist->SetBinContent(xBin+1, entry); xBin++; } // Second line is "other" backgrounds std::getline(input >> std::ws, line, '\n'); std::istringstream stream2(line); entry = 0; xBin = 0; while(stream2 >> entry){ this->BKGD_other->SetBinContent(xBin+1, entry); xBin++; } // Third line is the irreducible background std::getline(input >> std::ws, line, '\n'); std::istringstream stream3(line); entry = 0; xBin = 0; while(stream3 >> entry){ this->BKGD_irrid->SetBinContent(xBin+1, entry); xBin++; } }; // Read in the response matrix -- thus far, a response matrix is unique to the NCEL sample void MiniBooNE_NCEL_XSec_Treco_nu::SetResponseMatrix(std::string responseFile, int dim, double* arr_treco){ // Make a counter to track the line number int xBin = 0; std::string line; - std::ifstream response(responseFile.c_str(),ifstream::in); + std::ifstream response(responseFile.c_str(),std::ifstream::in); // Response matrix: x axis is Ttrue, y axis is Treco this->response_mat = new TH2D((this->fName+"_RESPONSE_MATRIX").c_str(), (this->fName+this->fPlotTitles).c_str(), 50, 0, 900, 51, arr_treco); if(response.is_open()) LOG(DEB) << "Reading in the response matrix from file: " << responseFile << std::endl; while(std::getline(response, line, '\n')){ std::istringstream stream(line); double entry; int yBin = 0; // Loop over entries and insert them into matrix // Multiply by the errors to get the covariance, rather than the correlation matrix while(stream >> entry){ this->response_mat->SetBinContent(xBin+1, yBin+1, entry); yBin++; } xBin++; } };