diff --git a/sbnanalysis/ana/CMakeLists.txt b/sbnanalysis/ana/CMakeLists.txt index 67eaa213..823de477 100644 --- a/sbnanalysis/ana/CMakeLists.txt +++ b/sbnanalysis/ana/CMakeLists.txt @@ -22,6 +22,7 @@ include_directories( $ENV{MESSAGEFACILITY_INC} ) link_directories( $ENV{MESSAGEFACILITY_LIB} $ENV{FHICLCPP_LIB} + $ENV{BOOST_LIB} $ENV{CANVAS_LIB} $ENV{CANVAS_ROOT_IO_LIB} $ENV{CETLIB_EXCEPT_LIB} diff --git a/sbnanalysis/ana/SBNOsc/Chi2Sensitivity.cxx b/sbnanalysis/ana/SBNOsc/Chi2Sensitivity.cxx index 017d0457..09176d61 100644 --- a/sbnanalysis/ana/SBNOsc/Chi2Sensitivity.cxx +++ b/sbnanalysis/ana/SBNOsc/Chi2Sensitivity.cxx @@ -154,11 +154,6 @@ void Chi2Sensitivity::Initialize(fhicl::ParameterSet* config) { } } - // Uniformly applied weights - if (pconfig.is_key_to_sequence("UniformWeights")) { - fUniformWeights = pconfig.get >("UniformWeights"); - } - // start at 0th event sample fSampleIndex = 0; } @@ -166,10 +161,10 @@ void Chi2Sensitivity::Initialize(fhicl::ParameterSet* config) { Chi2Sensitivity::EventSample::EventSample(const fhicl::ParameterSet& config) { // scaling stuff fName = config.get("name", ""); - fScaleFactor = config.get("scalefactor", 0.); + fScalePOT = config.get("ScalePOT", 0.); + fPOT = 0.; // setup detector stuff - fDistance = config.get("Distance", 0); std::vector xlim = config.get >("DetX"); std::vector ylim = config.get >("DetY"); std::vector zlim = config.get >("DetZ"); @@ -206,13 +201,39 @@ Chi2Sensitivity::EventSample::EventSample(const fhicl::ParameterSet& config) { // Take distance along z-axis as limits int numBinsPerMeter = config.get("NumDistanceBinsPerMeter", 1); double numBinsPerCM = numBinsPerMeter * 0.01; // unit conversion - unsigned n_bins = (unsigned) ((fZlim[1] - fZlim[0]) * numBinsPerCM) + 1; /* round up to be on the safe side */ + + // get beam center and baseline in cm + fBaseline = config.get("Baseline"); + fBeamCenterX = config.get("BeamCenterX", 0.); + fBeamCenterY = config.get("BeamCenterY", 0.); + fBeamFrontZ = config.get("BeamFrontZ", 0.); + + // get full range of distances across detector + double detector_distance = sqrt((fZlim[1] - fZlim[0]) * (fZlim[1] - fZlim[0]) \ + + (fXlim[1] - fXlim[0]) * (fXlim[1] - fXlim[0]) \ + + (fYlim[1] - fYlim[0]) * (fYlim[1] - fYlim[0])); + // get full range of beam + // taken approximately from eyeballing MCFlux information + double beam_width = sqrt( 125. * 125. /* x */ + 125. * 125. /* y */ + 5500. * 5500. /* z */); + unsigned n_bins = (unsigned) ((detector_distance + beam_width) * numBinsPerCM) + 1; unsigned n_limits = n_bins + 1; double dist_binwidth = 1./numBinsPerCM; + + // minimum distance from decayed vertex to detector + // taken approximately from MCFlux info + double min_baseline = fBaseline - 5500.; + for (unsigned i = 0; i < n_limits; i++) { - fDistBins.push_back(fDistance + (i * dist_binwidth) / 100000. /* cm -> km */); + fDistBins.push_back((min_baseline + (i * dist_binwidth)) / 100000. /* cm -> km */); } + // scaling in reco energy bins + fEnergyBinScale = config.get>("energy_bin_scale", {}); + if (fEnergyBinScale.size() == 0) { + fEnergyBinScale = std::vector(fBins.size() - 1, 1.0); + } + else assert(fEnergyBinScale.size() == fBins.size() - 1); + // setup histograms fBkgCounts = new TH1D((fName + " Background").c_str(), (fName + " Background").c_str(), fBins.size() - 1, &fBins[0]); fSignalCounts = new TH3D((fName + " Signal").c_str(), (fName + " Signal").c_str(), @@ -227,6 +248,11 @@ void Chi2Sensitivity::FileCleanup(TTree *eventTree) { // onto the next sample fSampleIndex ++; } + +void Chi2Sensitivity::ProcessSubRun(const SubRun *subrun) { + fCovariance.ProcessSubRun(subrun); + fEventSamples[fSampleIndex].fPOT += subrun->totgoodpot; +} void Chi2Sensitivity::ProcessEvent(const Event *event) { // have the covariance process the event @@ -247,17 +273,15 @@ void Chi2Sensitivity::ProcessEvent(const Event *event) { nuE = event->reco[n].reco_energy; } - // Get distance travelled along z in km - // (To reproduce proposal) - double dist = fEventSamples[fSampleIndex].fDistance + - (event->truth[truth_ind].neutrino.position.Z() - fEventSamples[fSampleIndex].fZlim[0])/100000. /* cm -> km */; - - // Get distance travelled (assuming nu started at (x, y, z) = (0, 0, min_det_zdim - det_dist)) - //double dx = (event->truth[truth_ind].neutrino.position.X() - (fDetDims[sample.fDet][0][1] + fDetDims[sample.fDet][0][0])/2) / 100000 /* cm -> km */, - //dy = (event->truth[truth_ind].neutrino.position.Y() - (fDetDims[sample.fDet][1][1] + fDetDims[sample.fDet][1][0])/2) / 100000 /* cm -> km */, - //dz = (event->truth[truth_ind].neutrino.position.Z() - fDetDims[sample.fDet][2][0]) / 100000 /* cm -> km */; - //double dist = TMath::Sqrt( dx*dx + dy*dy + (fDetDists[sample.fDet] + dz)*(fDetDists[sample.fDet] + dz) ); + // get distance travelled from decay vertex + // parent decay vertex + TVector3 p_decay_vtx = event->reco[n].truth.neutrino.parentDecayVtx; + TVector3 neutrino_vtx = event->reco[n].truth.neutrino.position; + double dz = fEventSamples[fSampleIndex].fBaseline - p_decay_vtx.Z() + neutrino_vtx.Z() - fEventSamples[fSampleIndex].fBeamFrontZ; + double dy = p_decay_vtx.Y() - neutrino_vtx.Y() + fEventSamples[fSampleIndex].fBeamCenterY; + double dx = p_decay_vtx.X() - neutrino_vtx.X() + fEventSamples[fSampleIndex].fBeamCenterX; + double dist = sqrt(dx*dx + dy*dy + dz*dz) / 100000. /* cm -> km */; // check if energy is within bounds if (nuE < fEventSamples[fSampleIndex].fBins[0] || nuE > fEventSamples[fSampleIndex].fBins[fEventSamples[fSampleIndex].fBins.size()-1]) { @@ -278,12 +302,16 @@ void Chi2Sensitivity::ProcessEvent(const Event *event) { // Apply selection (or rejection) efficiencies int isCC = event->truth[truth_ind].neutrino.iscc; - double wgt = isCC*(fSelectionEfficiency) + (1-isCC)*(1 - fBackgroundRejection); + // double wgt = isCC*(fSelectionEfficiency) + (1-isCC)*(1 - fBackgroundRejection); // and scale weight - wgt *= fEventSamples[fSampleIndex].fScaleFactor; + // wgt *= fEventSamples[fSampleIndex].fScaleFactor; // apply uniform weights - for (auto const &key: fUniformWeights) { - wgt *= event->truth[truth_ind].weights.at(key)[0]; + // for (auto const &key: fUniformWeights) { + // wgt *= event->truth[truth_ind].weights.at(key)[0]; + // } + double wgt = event->reco[n].weight; + if (fEventSamples[fSampleIndex].fScalePOT > 0) { + wgt *= fEventSamples[fSampleIndex].fScalePOT / fEventSamples[fSampleIndex].fPOT; } // fill in hitograms @@ -299,6 +327,30 @@ void Chi2Sensitivity::ProcessEvent(const Event *event) { } } +// apply reco energy bin scaling +void Chi2Sensitivity::Scale() { + for (int sample = 0; sample < fEventSamples.size(); sample++) { + for (int i = 1; i < fEventSamples[sample].fSignalCounts->GetNbinsX()+1; i++) { + double scale = fEventSamples[sample].fEnergyBinScale[i-1]; + for (int j = 1; j < fEventSamples[sample].fSignalCounts->GetNbinsY()+1; j++) { + for (int k = 1; k < fEventSamples[sample].fSignalCounts->GetNbinsZ()+1; k++) { + double content = fEventSamples[sample].fSignalCounts->GetBinContent(i, j, k); + //std::cout << "pre sample: " << sample << " bin content: " << fEventSamples[sample].fSignalCounts->GetBinContent(i); + fEventSamples[sample].fSignalCounts->SetBinContent(i, j, k, content * scale); + //std::cout << "post sample: " << sample << " bin content: " << fEventSamples[sample].fSignalCounts->GetBinContent(i); + } + } + } + + for (int i = 1; i < fEventSamples[sample].fBkgCounts->GetNbinsX()+1; i++) { + double scale = fEventSamples[sample].fEnergyBinScale[i-1]; + double content = fEventSamples[sample].fBkgCounts->GetBinContent(i); + fEventSamples[sample].fBkgCounts->SetBinContent(i, content * scale); + } + } +} + + void Chi2Sensitivity::GetChi2() { //// Invert Error Matrix diff --git a/sbnanalysis/ana/SBNOsc/Chi2Sensitivity.h b/sbnanalysis/ana/SBNOsc/Chi2Sensitivity.h index cbea3ae0..4f2d0b33 100644 --- a/sbnanalysis/ana/SBNOsc/Chi2Sensitivity.h +++ b/sbnanalysis/ana/SBNOsc/Chi2Sensitivity.h @@ -27,11 +27,18 @@ class Chi2Sensitivity: public core::PostProcessorBase { // Implement post-processor void Initialize(fhicl::ParameterSet* config); void ProcessEvent(const Event* event); + void ProcessSubRun(const SubRun *subrun); void FileCleanup(TTree* eventTree); - void Finalize() { fCovariance.Finalize(); GetChi2(); GetContours(); Write(); } + void Finalize() { + for (int i = 0; i < fEventSamples.size(); i++) { + std::cout << "POT: " << fEventSamples[i].fPOT << " to: " << fEventSamples[i].fScalePOT << " factor: " << fEventSamples[i].fScalePOT / fEventSamples[i].fPOT << std::endl; + } + + fCovariance.Finalize(); Scale(); GetChi2(); GetContours(); Write(); } // API Functions + void Scale(); void GetChi2(); void GetContours(); void Write(); @@ -58,12 +65,17 @@ class Chi2Sensitivity: public core::PostProcessorBase { // Config std::string fName; - double fDistance; //!< Distance in km + double fBaseline; //!< baseline (along z) in cm + double fBeamCenterX; //!< Center of beam in detector coordinates in x-dir in cm + double fBeamCenterY; //!< Center of beam in detector coordinates in y-dir in cm + double fBeamFrontZ; //!< Front face of detector along z-dir in cm std::array fXlim; //!< Detector size in cm std::array fYlim; //!< Detector size in cm std::array fZlim; //!< Detector size in cm - double fScaleFactor; //!< Factor for POT (etc.) scaling + double fScalePOT; //!< Factor for POT (etc.) scaling + double fPOT; int fOscType; //!< Oscilaltion type: 0 == None, 1 == numu -> nue, 2 == numu -> numu + std::vector fEnergyBinScale; // Storage TH1D *fBkgCounts; //!< Background count Histogram @@ -83,7 +95,6 @@ class Chi2Sensitivity: public core::PostProcessorBase { int fNumDm2; int fNumSin; std::vector fLogDm2Lims, fLogSinLims; - std::vector fUniformWeights; std::string fOutputFile; // whether to save stuff diff --git a/sbnanalysis/ana/SBNOsc/Covariance.cxx b/sbnanalysis/ana/SBNOsc/Covariance.cxx index b01115c8..4be40c48 100644 --- a/sbnanalysis/ana/SBNOsc/Covariance.cxx +++ b/sbnanalysis/ana/SBNOsc/Covariance.cxx @@ -29,7 +29,13 @@ Covariance::EventSample::EventSample(const fhicl::ParameterSet &config, unsigned // get configuration stuff fBins = config.get >("binlims"); fName = config.get("name", ""); - fScaleFactor = config.get("scalefactor", 1.0); + fScalePOT = config.get("ScalePOT", -1); + fPOT = 0.; + fEnergyBinScale = config.get>("energy_bin_scale", {}); + if (fEnergyBinScale.size() == 0) { + fEnergyBinScale = std::vector(fBins.size() - 1, 1.0); + } + else assert(fEnergyBinScale.size() == fBins.size() - 1); // setup histograms std::string cv_title = fName + " Central Value"; @@ -49,18 +55,14 @@ Covariance::EventSample::EventSample(const fhicl::ParameterSet &config, unsigned // Gets scale factors (weights) for different universes -std::vector GetUniWeights(const std::map > &weights, const std::vector &keys, int n_unis) { - - // Tentative format: universe u scale factor is the product of the u-th entries on each vector - // inside the map. For vectors with less than u entries, use the (u - vec_size)-th entry - +std::vector GetUniWeights(const std::map > &weights, const std::vector &keys, int n_unis, int uni_offset) { std::vector uweights; for (int u = 0; u < n_unis; u++) { double weight = 1.; for (auto const &key: keys) { const std::vector& this_weights = weights.at(key); - int wind = u % this_weights.size(); + int wind = u + uni_offset; weight *= this_weights.at(wind); } uweights.push_back(weight); @@ -79,15 +81,36 @@ void Covariance::Initialize(fhicl::ParameterSet* config) { fhicl::ParameterSet pconfig = config->get("Covariance"); fWeightKeys = pconfig.get>>("WeightKey"); + fWeightKeysCC = pconfig.get>>("WeightKeyCC", {}); + fWeightKeysNC = pconfig.get>>("WeightKeyNC", {}); fNVariations = fWeightKeys.size(); - // uniformly applied weights - if (pconfig.is_key_to_sequence("UniformWeights")) { - fUniformWeights = pconfig.get >("UniformWeights"); + // all should have the same size + assert(fWeightKeysCC.size() == 0 || fWeightKeysCC.size() == fWeightKeys.size()); + assert(fWeightKeysNC.size() == 0 || fWeightKeysNC.size() == fWeightKeys.size()); + + // merge weight keys into CC/NC keys + if (fWeightKeysCC.size() == 0) { + fWeightKeysCC = fWeightKeys; + } + else { + for (unsigned i = 0; i < fNVariations; i++) { + fWeightKeysCC[i].insert(fWeightKeysCC[i].end(), fWeightKeys[i].begin(), fWeightKeys[i].end()); + } + } + if (fWeightKeysNC.size() == 0) { + fWeightKeysNC = fWeightKeys; + } + else { + for (unsigned i = 0; i < fNVariations; i++) { + fWeightKeysNC[i].insert(fWeightKeysNC[i].end(), fWeightKeys[i].begin(), fWeightKeys[i].end()); + } } // number of universes to be used fNumAltUnis = pconfig.get("NumAltUnis", 0); + // and offset into the eventweight vector + fAltUniOffset = pconfig.get("AltUniOffset", 0); // Type of energy fEnergyType = pconfig.get("EnergyType", "Reco"); @@ -100,7 +123,7 @@ void Covariance::Initialize(fhicl::ParameterSet* config) { // If an event has any weight larger than this value, it will be // thrown away for the purposes of covariance construction fWeightMax = pconfig.get("WeightMax", -1.); - + // get event samples std::vector configEventSamples = \ config->get >("EventSamples"); @@ -119,6 +142,10 @@ void Covariance::Initialize(fhicl::ParameterSet* config) { fSampleIndex = 0; } +void Covariance::ProcessSubRun(const SubRun *subrun) { + fEventSamples[fSampleIndex].fPOT += subrun->totgoodpot; +} + void Covariance::ProcessEvent(const Event *event) { // iterate over each interaction in the event for (int n = 0; n < event->reco.size(); n++) { @@ -137,20 +164,35 @@ void Covariance::ProcessEvent(const Event *event) { // Apply selection (or rejection) efficiencies int isCC = event->truth[truth_ind].neutrino.iscc; - double wgt = isCC*(fSelectionEfficiency) + (1-isCC)*(1 - fBackgroundRejection); + // double wgt = isCC*(fSelectionEfficiency) + (1-isCC)*(1 - fBackgroundRejection); // apply scaling from fScaleFactor - wgt *= fEventSamples[fSampleIndex].fScaleFactor; + // wgt *= fEventSamples[fSampleIndex].fScaleFactor; // apply uniform weights - for (auto const &key: fUniformWeights) { - wgt *= event->truth[truth_ind].weights.at(key)[0]; + // for (auto const &key: fUniformWeights) { + // wgt *= event->truth[truth_ind].weightmap.at(key)[0]; + // } + double wgt = event->reco[n].weight; + // apply POT scaling if configured + if (fEventSamples[fSampleIndex].fScalePOT > 0) { + wgt *= fEventSamples[fSampleIndex].fScalePOT / fEventSamples[fSampleIndex].fPOT; } - // Get weights for each alternative universe std::vector> uweights; - for (std::vector &weight_keys: fWeightKeys) { - uweights.push_back( - GetUniWeights(event->truth[truth_ind].weights, weight_keys, fNumAltUnis) - ); + // Get weights for each alternative universe + if (isCC) { + for (std::vector &weight_keys: fWeightKeysCC) { + uweights.push_back( + GetUniWeights(event->truth[truth_ind].weightmap, weight_keys, fNumAltUnis, fAltUniOffset) + ); + } + } + else { + for (std::vector &weight_keys: fWeightKeysNC) { + uweights.push_back( + GetUniWeights(event->truth[truth_ind].weightmap, weight_keys, fNumAltUnis, fAltUniOffset) + ); + } + } // see if weight is too big @@ -179,6 +221,25 @@ void Covariance::ProcessEvent(const Event *event) { } } +void Covariance::Scale() { + for (int sample_i = 0; sample_i < fEventSamples.size(); sample_i++) { + for (int i = 1; i < fEventSamples[sample_i].fCentralValue->GetNbinsX()+1; i++) { + //std::cout << "pre sample: " << sample_i << " bin content: " << fEventSamples[sample_i].fCentralValue->GetBinContent(i); + fEventSamples[sample_i].fCentralValue->SetBinContent(i, + fEventSamples[sample_i].fEnergyBinScale[i-1] * fEventSamples[sample_i].fCentralValue->GetBinContent(i)); + //std::cout << "post sample: " << sample_i << " bin content: " << fEventSamples[sample_i].fCentralValue->GetBinContent(i); + } + for (int variation = 0; variation < fNVariations; variation++) { + for (int uni = 0; uni < fNumAltUnis; uni++) { + for (int i = 1; i < fEventSamples[sample_i].fUniverses[variation][uni]->GetNbinsX(); i++) { + double content = fEventSamples[sample_i].fUniverses[variation][uni]->GetBinContent(i); + fEventSamples[sample_i].fUniverses[variation][uni]->SetBinContent(i, content * fEventSamples[sample_i].fEnergyBinScale[i-1]); + } + } + } + } +} + void Covariance::FileCleanup(TTree *eventTree) { // onto the next sample diff --git a/sbnanalysis/ana/SBNOsc/Covariance.h b/sbnanalysis/ana/SBNOsc/Covariance.h index 4d6a718e..c5bd7306 100644 --- a/sbnanalysis/ana/SBNOsc/Covariance.h +++ b/sbnanalysis/ana/SBNOsc/Covariance.h @@ -36,11 +36,13 @@ class Covariance: public core::PostProcessorBase { void FileCleanup(TTree *eventTree); void Initialize(fhicl::ParameterSet *config); void ProcessEvent(const Event *event); - void Finalize() { GetCovs(); Write(); } + void ProcessSubRun(const SubRun *subrun); + void Finalize() { Scale(); GetCovs(); Write(); } // API Functions void GetCovs(); void Write(); + void Scale(); // build the covariance matrix TMatrixDSym CovarianceMatrix(); @@ -57,19 +59,23 @@ class Covariance: public core::PostProcessorBase { /** Constructors. */ EventSample(const fhicl::ParameterSet &config, unsigned nUniverses, unsigned nVariations); - double fScaleFactor; //!< Factor for POT (etc.) scaling + double fPOT; + double fScalePOT; //!< Factor for POT (etc.) scaling std::vector fBins; //!< Energy bin limits TH1D *fCentralValue; //!< central value histogram std::vector> fUniverses; //!< List of histogram per systematic universe. // List has index per covariance matrix to be generated. std::string fName; //!< Name for the sample + std::vector fEnergyBinScale; }; void GetCovPerVariation(unsigned variation); // config std::vector> fWeightKeys; - std::vector fUniformWeights; + std::vector> fWeightKeysCC; + std::vector> fWeightKeysNC; int fNumAltUnis; + int fAltUniOffset; std::string fEnergyType; @@ -89,7 +95,6 @@ class Covariance: public core::PostProcessorBase { // Stored Event Samples std::vector fEventSamples; - }; } // namespace SBNOsc diff --git a/sbnanalysis/ana/SBNOsc/NumuSelection.cxx b/sbnanalysis/ana/SBNOsc/NumuSelection.cxx index bb5df2ab..7bf81b1b 100644 --- a/sbnanalysis/ana/SBNOsc/NumuSelection.cxx +++ b/sbnanalysis/ana/SBNOsc/NumuSelection.cxx @@ -30,7 +30,33 @@ NumuSelection::NumuSelection() : SelectionBase(), _event_counter(0), _nu_count(0), - _interactionInfo(new std::vector) {} + _rand(57 /* seed */), + _interactionInfo(new std::vector) { + // setup the event categories + _eventCategories["CC0pi"] = 0; + _eventCategories["CC1pi"] = 0; + _eventCategories["CC2pi"] = 0; + _eventCategories["CCnpi"] = 0; + + _eventCategories["NC0pi"] = 0; + _eventCategories["NC1pi"] = 0; + _eventCategories["NC2pi"] = 0; + _eventCategories["NCnpi"] = 0; + + // more + _eventCategories["CCnoMU"] = 0; + _eventCategories["NCwiMU"] = 0; + + // more + _eventCategories["CCQE"] = 0; + _eventCategories["CCRES"] = 0; + _eventCategories["CCDIS"] = 0; + _eventCategories["CCCOH"] = 0; + _eventCategories["CCMEC"] = 0; + _eventCategories["CCCOHE"] = 0; + _eventCategories["CCother"] = 0; +} + double aaBoxesMin(const std::vector &boxes, unsigned dim) { return std::min_element(boxes.begin(), boxes.end(), [dim](auto &lhs, auto &rhs) { return lhs.Min()[dim] < rhs.Min()[dim]; })->Min()[dim]; @@ -40,6 +66,20 @@ double aaBoxesMax(const std::vector &boxes, unsigned dim) { return std::max_element(boxes.begin(), boxes.end(), [dim](auto &lhs, auto &rhs) { return lhs.Max()[dim] < rhs.Max()[dim]; })->Max()[dim]; } +geoalgo::AABox shaveVolume(const geoalgo::AABox &select_volume, double delta) { + + double xmin = select_volume.Min()[0] + 5.; + double ymin = select_volume.Min()[1] + 5.; + double zmin = select_volume.Min()[2] + 5.; + + double xmax = select_volume.Max()[0] - 5.; + double ymax = select_volume.Max()[1] - 5.; + double zmax = select_volume.Max()[2] - 5.; + + return geoalgo::AABox(xmin, ymin, zmin, xmax, ymax, zmax); +} + + void NumuSelection::Initialize(fhicl::ParameterSet* config) { if (config) { fhicl::ParameterSet pconfig = config->get("NumuSelection"); @@ -70,7 +110,7 @@ void NumuSelection::Initialize(fhicl::ParameterSet* config) { } _config.doFVCut = pconfig.get("doFVcut", true); - _config.trajPointLength = pconfig.get("trajPointLength", false); + _config.trajPointLength = pconfig.get("trajPointLength", true); _config.vertexDistanceCut = pconfig.get("vertexDistance", -1); _config.minLengthContainedTrack = pconfig.get("minLengthContainedTrack", -1); _config.minLengthExitingTrack = pconfig.get("minLengthExitingTrack", -1); @@ -83,10 +123,13 @@ void NumuSelection::Initialize(fhicl::ParameterSet* config) { _config.acceptShakyTracks = pconfig.get("acceptShakyTracks", false); _config.verbose = pconfig.get("verbose", false); _config.cutKMEC = pconfig.get("cutKMEC", false); + _config.selectMode = pconfig.get("selectMode", -1); + _config.selectCCNC = pconfig.get("selectCCNC", -1); _config.onlyKMEC = pconfig.get("onlyKMEC", false); // setup weight config _config.selectionEfficiency = pconfig.get("selectionEfficiency", 1.0); + _config.backgroundRejection = pconfig.get("backgroundRejection", 0.0); _config.uniformWeights = pconfig.get>("uniformWeights", {}); _config.constantWeight = pconfig.get("constantWeight", 1.0); @@ -113,6 +156,24 @@ void NumuSelection::Initialize(fhicl::ParameterSet* config) { _root_histos[i].h_numu_Vyz = new TH2D(("numu_Vyz_" + cut_names[i]).c_str(), "numu_Vyz", 20, aaBoxesMin(_config.active_volumes, 1), aaBoxesMax(_config.active_volumes, 1), 20, aaBoxesMin(_config.active_volumes, 2), aaBoxesMax(_config.active_volumes, 2)); + + _root_histos[i].h_numu_Vx_sig = new TH1D(("numu_Vx_sig" + cut_names[i]).c_str(), "numu_Vx_sig", 20, + aaBoxesMin(_config.fiducial_volumes, 0), aaBoxesMax(_config.fiducial_volumes, 0)); + _root_histos[i].h_numu_Vy_sig = new TH1D(("numu_Vy_sig" + cut_names[i]).c_str(), "numu_Vy_sig", 20, + aaBoxesMin(_config.fiducial_volumes, 1), aaBoxesMax(_config.fiducial_volumes, 1)); + _root_histos[i].h_numu_Vz_sig = new TH1D(("numu_Vz_sig" + cut_names[i]).c_str(), "numu_Vz_sig", 20, + aaBoxesMin(_config.fiducial_volumes, 2), aaBoxesMax(_config.fiducial_volumes, 2)); + + _root_histos[i].h_numu_Vx_bkg = new TH1D(("numu_Vx_bkg" + cut_names[i]).c_str(), "numu_Vx_bkg", 20, + aaBoxesMin(_config.fiducial_volumes, 0), aaBoxesMax(_config.fiducial_volumes, 0)); + _root_histos[i].h_numu_Vy_bkg = new TH1D(("numu_Vy_bkg" + cut_names[i]).c_str(), "numu_Vy_bkg", 20, + aaBoxesMin(_config.fiducial_volumes, 1), aaBoxesMax(_config.fiducial_volumes, 1)); + _root_histos[i].h_numu_Vz_bkg = new TH1D(("numu_Vz_bkg" + cut_names[i]).c_str(), "numu_Vz_bkg", 20, + aaBoxesMin(_config.fiducial_volumes, 2), aaBoxesMax(_config.fiducial_volumes, 2)); + + _root_histos[i].h_numu_t_is_muon_sig = new TH1D(("t_is_muon_sig_" + cut_names[i]).c_str(), "t_is_muon_sig", 3, -1.5, 1.5); + _root_histos[i].h_numu_t_is_muon_bkg = new TH1D(("t_is_muon_bkg_" + cut_names[i]).c_str(), "t_is_muon_bkg", 3, -1.5, 1.5); + } // set up TGraph keeping track of cut counts @@ -140,8 +201,25 @@ void NumuSelection::Finalize() { _root_histos[i].h_numu_Vxy->Write(); _root_histos[i].h_numu_Vxz->Write(); _root_histos[i].h_numu_Vyz->Write(); + + _root_histos[i].h_numu_Vx_sig->Write(); + _root_histos[i].h_numu_Vy_sig->Write(); + _root_histos[i].h_numu_Vz_sig->Write(); + _root_histos[i].h_numu_Vx_bkg->Write(); + _root_histos[i].h_numu_Vy_bkg->Write(); + _root_histos[i].h_numu_Vz_bkg->Write(); + + _root_histos[i].h_numu_t_is_muon_sig->Write(); + _root_histos[i].h_numu_t_is_muon_bkg->Write(); } _cut_counts->Write(); + + // print out event stats + std::map::iterator it = _eventCategories.begin(); + while (it != _eventCategories.end()) { + std::cout << "Category: " << it->first << " " << it->second << std::endl; + it ++; + } } @@ -151,7 +229,6 @@ bool NumuSelection::ProcessEvent(const gallery::Event& ev, const std::vectorclear(); @@ -168,6 +245,52 @@ bool NumuSelection::ProcessEvent(const gallery::Event& ev, const std::vectorSetPoint(0, 0, _cut_counts->GetY()[0] + mctruths.size()); + // categorize events: + for (unsigned i = 0; i < mctruths.size(); i++) { + auto const &mctruth = mctruths[i]; + Event::Interaction interaction = truth[i]; + double weight = 1.; + // maybe cut MEC + bool pass_kMEC = !(_config.cutKMEC && mctruth.GetNeutrino().Mode() == simb::kMEC) && !(_config.onlyKMEC && mctruth.GetNeutrino().Mode() != simb::kMEC); + if (!pass_kMEC) continue; + for (auto const &key: _config.uniformWeights) { + weight *= interaction.weights.at(key)[0]; + } + if (containedInAV(mctruth.GetNeutrino().Nu().Position().Vect())) { + bool iscc = mctruth.GetNeutrino().CCNC() == 0; + unsigned n_pions = 0; + for (auto const &mctrack: mctracks) { + if (isFromNuVertex(mctruth, mctrack) && abs(mctrack.PdgCode()) == 211 && mctrack.Process() == "primary") { + n_pions += 1; + } + // categorize on number of pions + } + if (iscc) { + if (n_pions == 0) _eventCategories["CC0pi"] += weight; + else if (n_pions == 1) _eventCategories["CC1pi"] += weight; + else if (n_pions == 2) _eventCategories["CC2pi"] += weight; + else _eventCategories["CCnpi"] += weight; + } + else { + if (n_pions == 0) _eventCategories["NC0pi"] += weight; + else if (n_pions == 1) _eventCategories["NC1pi"] += weight; + else if (n_pions == 2) _eventCategories["NC2pi"] += weight; + else _eventCategories["NCnpi"] += weight; + } + // categorize on event mode + if (iscc) { + int mode = mctruth.GetNeutrino().Mode(); + if (mode == simb::kQE) _eventCategories["CCQE"] += weight; + else if (mode == simb::kCoh) _eventCategories["CCCOH"] += weight; + else if (mode == simb::kRes) _eventCategories["CCRES"] += weight; + else if (mode == simb::kMEC) _eventCategories["CCMEC"] += weight; + else if (mode == simb::kDIS) _eventCategories["CCDIS"] += weight; + else if (mode == simb::kCohElastic) _eventCategories["CCCOHE"] += weight; + else _eventCategories["CCother"] += weight; + } + } + } + // Iterate through the neutrinos bool selected = false; for (size_t i=0; iFill(nu.Nu().Vx(), nu.Nu().Vz()); _root_histos[select_i].h_numu_Vyz->Fill(nu.Nu().Vy(), nu.Nu().Vz()); + if (is_signal) { + _root_histos[select_i].h_numu_Vx_sig->Fill(nu.Nu().Vx()); + _root_histos[select_i].h_numu_Vy_sig->Fill(nu.Nu().Vy()); + _root_histos[select_i].h_numu_Vz_sig->Fill(nu.Nu().Vz()); + + _root_histos[select_i].h_numu_t_is_muon_sig->Fill(abs(intInfo.t_pdgid) == 13); + } + else { + _root_histos[select_i].h_numu_Vx_bkg->Fill(nu.Nu().Vx()); + _root_histos[select_i].h_numu_Vy_bkg->Fill(nu.Nu().Vy()); + _root_histos[select_i].h_numu_Vz_bkg->Fill(nu.Nu().Vz()); + + _root_histos[select_i].h_numu_t_is_muon_bkg->Fill(abs(intInfo.t_pdgid) == 13); + } + // also update cut count _cut_counts->SetPoint(select_i+1, select_i+1, _cut_counts->GetY()[select_i+1] + 1); } @@ -296,11 +444,47 @@ NumuSelection::TrackInfo NumuSelection::trackInfo(const sim::MCTrack &track) { pos = track[i].Position(); } } + // length calculation method used in proposal + else if (track.size() != 0 && !_config.trajPointLength) { + TLorentzVector pos = track.Start().Position(); + // get the active volume that the start position is in + int active_volume_index = -1; + // contruct pos Point + geoalgo::Point_t pos_point(pos); + for (int i = 0; i < _config.active_volumes.size(); i++) { + if (_config.active_volumes[i].Contain(pos_point)) { + active_volume_index = i; + } + } + if (active_volume_index >= 0) { + // setup the "shaved" volume used in the proposal + geoalgo::AABox volume = shaveVolume(_config.active_volumes[active_volume_index], 5.); + if (volume.Contain(pos_point)) { + int i = 1; + TLorentzVector this_pos = pos; + for (; i < track.size(); i++) { + this_pos = track[i].Position(); + geoalgo::Point_t p(this_pos.Vect()); + if (!volume.Contain(p)) { + contained_in_AV = false; + break; + } + } + // Length of track is distance from start to last contained position + length = (pos.Vect() - this_pos.Vect()).Mag(); + std::vector volumes {volume}; + contained_length = containedLength(pos.Vect(), this_pos.Vect(), volumes); + } + else contained_in_AV = false; + } + else contained_in_AV = false; + + } // If active volume is misconfigured, then tracks may be generated w/out points. // Optionally, we can accept them. // // Also, use this method if configured - else if (_config.acceptShakyTracks || !_config.trajPointLength) { + else if (_config.acceptShakyTracks) { contained_length = containedLength(track.Start().Position().Vect(), track.End().Position().Vect(), _config.active_volumes); length = (track.Start().Position().Vect() - track.End().Position().Vect()).Mag(); contained_in_AV = containedInAV(track.Start().Position().Vect()) && containedInAV(track.End().Position().Vect()); @@ -315,6 +499,9 @@ NumuSelection::TrackInfo NumuSelection::trackInfo(const sim::MCTrack &track) { //std::cout << "START: " << start.X() << " " << start.Y() << " " << start.Z() << std::endl; //std::cout << "END: " << end.X() << " " << end.Y() << " " << end.Z() << std::endl; } + // else { + // assert(false); //bad + // } return NumuSelection::TrackInfo({contained_in_AV, contained_length, length}); } @@ -392,7 +579,7 @@ NumuSelection::NuMuInteraction NumuSelection::interactionInfo(const gallery::Eve calculator.lepton_contained = false; calculator.lepton_contained_length = -1; calculator.lepton_index = -1; - return NumuSelection::NuMuInteraction({false, -1, -1, -1, -1, -1}); + return NumuSelection::NuMuInteraction({false, -1, -1, -1, -1, -1, TVector3()}); } // otherwise get the track info and energy info else { @@ -405,10 +592,11 @@ NumuSelection::NuMuInteraction NumuSelection::interactionInfo(const gallery::Eve calculator.lepton_index = track_ind; // smear the energy - double smeared_energy = smearLeptonEnergy(mctrack_list[track_ind], calculator); + double smeared_energy = smearLeptonEnergy(_rand, mctrack_list[track_ind], calculator); // truth kinetic energy double truth_energy = (mctrack_list[track_ind].Start().E()) / 1000.; /* MeV -> GeV */ - return NumuSelection::NuMuInteraction({t_info.t_is_contained, t_info.t_contained_length, t_info.t_length, mctrack_list[track_ind].PdgCode(), truth_energy, smeared_energy}); + TVector3 momentum = mctrack_list[track_ind].Start().Momentum().Vect(); + return NumuSelection::NuMuInteraction({t_info.t_is_contained, t_info.t_contained_length, t_info.t_length, mctrack_list[track_ind].PdgCode(), truth_energy, smeared_energy, momentum}); } } @@ -468,6 +656,11 @@ std::array NumuSelection::Select(const gallery::Even // STUDY KMEC: remove MEC events bool pass_kMEC = !(_config.cutKMEC && nu.Mode() == simb::kMEC) && !(_config.onlyKMEC && nu.Mode() != simb::kMEC); + // select another mode if necessary + bool pass_Mode = _config.selectMode < 0 || nu.Mode() == _config.selectMode; + // maybe require cc or nc + bool pass_CCNC = _config.selectCCNC < 0 || nu.CCNC() == _config.selectCCNC; + pass_kMEC = pass_kMEC && pass_Mode && pass_CCNC; // retrun list of cuts return {pass_kMEC, pass_AV && pass_kMEC, pass_valid_track && pass_kMEC && pass_AV, pass_valid_track && pass_kMEC && pass_FV, pass_valid_track && pass_kMEC && pass_FV && pass_min_length}; @@ -496,10 +689,12 @@ bool NumuSelection::passRecoVertex(const TVector3 &truth_v, const TVector3 &reco } bool NumuSelection::passMinLength(double length, bool stop_in_tpc) { + bool pass_stopped = _config.minLengthContainedTrack < 0 || length > _config.minLengthContainedTrack; + bool pass_exiting = _config.minLengthExitingTrack < 0 || length > _config.minLengthExitingTrack; if (!stop_in_tpc) - return _config.minLengthExitingTrack < 0 || length > _config.minLengthExitingTrack; + return pass_stopped && pass_exiting; else - return _config.minLengthContainedTrack < 0 || length > _config.minLengthContainedTrack; + return pass_stopped; } } // namespace SBNOsc diff --git a/sbnanalysis/ana/SBNOsc/NumuSelection.h b/sbnanalysis/ana/SBNOsc/NumuSelection.h index 81d97e04..75428670 100644 --- a/sbnanalysis/ana/SBNOsc/NumuSelection.h +++ b/sbnanalysis/ana/SBNOsc/NumuSelection.h @@ -11,6 +11,7 @@ #include #include +#include #include "canvas/Utilities/InputTag.h" #include "core/SelectionBase.hh" @@ -28,6 +29,8 @@ // take the geobox stuff from uboonecode #include "ubcore/LLBasicTool/GeoAlgo/GeoAABox.h" #include "ubcore/LLBasicTool/GeoAlgo/GeoAlgo.h" +#include "TRandom.h" +#include "TRandomGen.h" class TH2D; @@ -70,6 +73,7 @@ class NumuSelection : public core::SelectionBase { int t_pdgid; //!< PDGID of primary track (muon or pi+) double t_energy_true; //!< True energy of primary track [GeV] double t_energy_smeared; //!< Smeared energy of primary track [GeV] + TVector3 t_momentum; // default constructor -- fills with bogus info /* @@ -103,7 +107,10 @@ class NumuSelection : public core::SelectionBase { // (%) = -A * Log(B * L) where L is the lepton contained length bool cutKMEC; //!< Whether to remove MEC events (useful for studying difference w.r.t. proposal) bool onlyKMEC; //!< Whether to remove all non-MEC events + int selectMode; + int selectCCNC; double selectionEfficiency; //!< Signal efficiency weight applied to signal (charged current) events + double backgroundRejection; //!< Rejection applied to background (NC) events std::vector uniformWeights; //!< Weights taken from "EventWeight" that should be applied to the weight of each event double constantWeight; //!< constant weight to apply uniformly to each event }; @@ -123,6 +130,14 @@ class NumuSelection : public core::SelectionBase { TH2D *h_numu_Vxy; //!< 2D x-y vertex histogram [cm] TH2D *h_numu_Vxz; //!< 2D x-z vertex histogram [cm] TH2D *h_numu_Vyz; //!< 2D y-z vertex histogram [cm] + TH1D *h_numu_Vx_sig; + TH1D *h_numu_Vy_sig; + TH1D *h_numu_Vz_sig; + TH1D *h_numu_Vx_bkg; + TH1D *h_numu_Vy_bkg; + TH1D *h_numu_Vz_bkg; + TH1D *h_numu_t_is_muon_sig; //!< histogram of whether associated track is a muon + TH1D *h_numu_t_is_muon_bkg; //!< histogram of whether associated track is a muon }; // helper struct holding track info -- see NumuInteraction for variable details @@ -207,11 +222,14 @@ class NumuSelection : public core::SelectionBase { return {"MEC", "AV", "Track", "FV", "min_L"}; } + TRandomMT64 _rand; //!< random number generation Config _config; //!< The config std::vector *_interactionInfo; //!< Branch holder RootHistos _root_histos[nCuts]; //!< Histos (one group per cut) + + std::map _eventCategories; }; } // namespace SBNOsc diff --git a/sbnanalysis/ana/SBNOsc/Utilities.cxx b/sbnanalysis/ana/SBNOsc/Utilities.cxx index 8c6c7e44..3243562c 100644 --- a/sbnanalysis/ana/SBNOsc/Utilities.cxx +++ b/sbnanalysis/ana/SBNOsc/Utilities.cxx @@ -1,7 +1,6 @@ #include #include "TDatabasePDG.h" -#include "TRandom.h" #include "nusimdata/SimulationBase/MCTruth.h" #include "nusimdata/SimulationBase/MCNeutrino.h" @@ -14,6 +13,7 @@ #include "ubcore/LLBasicTool/GeoAlgo/GeoLineSegment.h" #include +#include namespace ana { namespace SBNOsc { @@ -160,8 +160,28 @@ double containedLength(const TVector3 &v0, const TVector3 &v1, // none contained -- either must have zero or two intersections if (n_contained == 0) { auto intersections = algo.Intersection(line, box); - assert(intersections.size() == 0 || intersections.size() == 2); - if (intersections.size() == 2) { + if (!(intersections.size() == 0 || intersections.size() == 2)) { + // more floating point error fixes... + // + // figure out which points are near the edge + double tol = 1e-5; + bool p0_edge = algo.SqDist(p0, box) < tol; + bool p1_edge = algo.SqDist(p1, box) < tol; + // and which points are near the intersection + TVector3 vint = intersections.at(0).ToTLorentzVector().Vect(); + bool p0_int = (v0 - vint).Mag() < tol; + bool p1_int = (v1 - vint).Mag() < tol; + // exactly one of them should produce the intersection + assert((p0_int && p0_edge) != (p1_int && p1_edge)); + // both close to edge -- full length is contained + if (p0_edge && p1_edge) { + length += (v0 - v1).Mag(); + } + // otherwise -- one of them is not on an edge, no length is contained + else {} + } + // assert(intersections.size() == 0 || intersections.size() == 2); + else if (intersections.size() == 2) { TVector3 start(intersections.at(0).ToTLorentzVector().Vect()); TVector3 end(intersections.at(1).ToTLorentzVector().Vect()); length += (start - end).Mag(); @@ -172,21 +192,153 @@ double containedLength(const TVector3 &v0, const TVector3 &v1, return length; } +// run the proposal energy smearing w/ MCParticles +double visibleEnergyProposalMCParticles(TRandom &rand, const simb::MCTruth &mctruth, const std::vector mctrack_list, const VisibleEnergyCalculator &calculator) { + double total = 0; + for (int iparticle=0; iparticle GeV */); + double smear_energy = track_energy; + + // In the proposal, different PDG codes got different energies + + // kaons did not have the mass subtracted + if (abs(pdg) == 321) { + track_energy = particle.Momentum().E(); + smear_energy = track_energy; + } + // pions did not have mass subtracted in the smearing + if (abs(pdg) == 211) { + smear_energy = particle.Momentum().E(); + } + + bool skip = false; + + // threshold -- only for protons + if (abs(pdg) == 2212 && track_energy < calculator.track_threshold) { + skip = true; + } + + // add in smeared energy + double this_smeared_energy = rand.Gaus(track_energy, smear_energy * calculator.track_energy_distortion); + // clamp to zero + this_smeared_energy = std::max(this_smeared_energy, 0.); + + // Ignore particles not from nu vertex + if (!isFromNuVertex(mctruth, particle) || particle.Process() != "primary") { + skip = true; + } + // ignore everything except for protons and kaons and pions + if (!(abs(pdg) == 2212 || abs(pdg) == 211 || abs(pdg) == 321)) { + skip = true; + } + // double distance = (mctruth.GetNeutrino().Nu().Trajectory().Position(0).Vect() - mct.Start().Position().Vect()).Mag(); + // double l_distance = (mctruth.GetNeutrino().Nu().Trajectory().Position(0) - mct.Start().Position()).Mag(); + // std::cout << "New particle -- true E: " << track_energy << " smeared E: " << this_smeared_energy << " PDG: " << pdg << " is primary: " << (calculator.lepton_index == ind) << " is skipped: " << skip << " distance: " << distance << " lorentz distance: " << l_distance << " process: " << mct.Process() << std::endl; + if (skip) continue; + + total += this_smeared_energy; + } + + // ...and primary lepton energy (for CC events) + // only add in extra here if identified "lepton" is actually a lepton + if (calculator.lepton_index >= 0 && (abs(mctrack_list[calculator.lepton_index].PdgCode()) == 13 || abs(mctrack_list[calculator.lepton_index].PdgCode()) == 11)) { + double lepton_energy = smearLeptonEnergy(rand, mctrack_list[calculator.lepton_index], calculator); + total += lepton_energy; + } + // std::cout << "Reco Energy: " << total << std::endl; + + return total; + + } + +// copies the exact energy smearing used in the proposal +double visibleEnergyProposal(TRandom &rand, const simb::MCTruth &mctruth, const std::vector &mctrack_list, const VisibleEnergyCalculator &calculator) { + double total = 0.; + + // std::cout << "\n\nNew Interaction" << std::endl; + // std::cout << "True E: " << mctruth.GetNeutrino().Nu().E() << std::endl; + for (unsigned ind = 0; ind < mctrack_list.size(); ind++) { + auto const &mct = mctrack_list[ind]; + int pdg = mct.PdgCode(); + + double track_energy = (mct.Start().E() - PDGMass(pdg)) / 1000. /* MeV -> GeV */; + double smear_energy = track_energy; -double visibleEnergy(const simb::MCTruth &mctruth, const std::vector &mctrack_list, const std::vector &mcshower_list, + // In the proposal, different PDG codes got different energies + + // kaons did not have the mass subtracted + if (abs(pdg) == 321) { + track_energy = mct.Start().E() / 1000.; + smear_energy = track_energy; + } + // pions did not have mass subtracted in the smearing + if (abs(pdg) == 211) { + smear_energy = mct.Start().E() / 1000.; + } + + bool skip = false; + + // threshold -- only for protons + if (abs(pdg) == 2212 && track_energy < calculator.track_threshold) { + skip = true; + } + + // add in smeared energy + double this_smeared_energy = rand.Gaus(track_energy, smear_energy * calculator.track_energy_distortion); + // clamp to zero + this_smeared_energy = std::max(this_smeared_energy, 0.); + + // Ignore particles not from nu vertex + if (!isFromNuVertex(mctruth, mct) || mct.Process() != "primary") { + skip = true; + } + // ignore everything except for protons and kaons and pions + if (!(abs(pdg) == 2212 || abs(pdg) == 211 || abs(pdg) == 321)) { + skip = true; + } + // account for primary track later + if ((abs(pdg) == 13 || abs(pdg) == 11) && calculator.lepton_index == ind) { + skip = true; + } + // double distance = (mctruth.GetNeutrino().Nu().Trajectory().Position(0).Vect() - mct.Start().Position().Vect()).Mag(); + // double l_distance = (mctruth.GetNeutrino().Nu().Trajectory().Position(0) - mct.Start().Position()).Mag(); + // std::cout << "New particle -- true E: " << track_energy << " smeared E: " << this_smeared_energy << " PDG: " << pdg << " is primary: " << (calculator.lepton_index == ind) << " is skipped: " << skip << " distance: " << distance << " lorentz distance: " << l_distance << " process: " << mct.Process() << std::endl; + if (skip) continue; + + total += this_smeared_energy; + } + + // ...and primary lepton energy (for CC events) + // only add in extra here if identified "lepton" is actually a lepton + if (calculator.lepton_index >= 0 && (abs(mctrack_list[calculator.lepton_index].PdgCode()) == 13 || abs(mctrack_list[calculator.lepton_index].PdgCode()) == 11)) { + double lepton_energy = smearLeptonEnergy(rand, mctrack_list[calculator.lepton_index], calculator); + total += lepton_energy; + } + // std::cout << "Reco Energy: " << total << std::endl; + + return total; +} + + +double visibleEnergy(TRandom &rand, const simb::MCTruth &mctruth, const std::vector &mctrack_list, const std::vector &mcshower_list, const VisibleEnergyCalculator &calculator, bool include_showers) { double visible_E = 0; - // set up distortion if need be - TRandom rand; - // primary leptron track const sim::MCTrack *lepton_track = NULL; bool lepton_track_exists = false; // total up visible energy from tracks... - unsigned ind = 0; - for (auto const &mct: mctrack_list) { + double track_visible_energy = 0.; + for (unsigned ind = 0; ind < mctrack_list.size(); ind++) { + auto const &mct = mctrack_list[ind]; // ignore particles not from nu vertex, non primary particles, and uncharged particles if (!isFromNuVertex(mctruth, mct) || abs(PDGCharge(mct.PdgCode())) < 1e-4 || mct.Process() != "primary") continue; @@ -194,21 +346,27 @@ double visibleEnergy(const simb::MCTruth &mctruth, const std::vector 1e-4) { - this_visible_energy = rand.Gaus(this_visible_energy, calculator.track_energy_distortion*this_visible_energy); - // clamp to 0 - this_visible_energy = std::max(this_visible_energy, 0.); - } if (this_visible_energy > calculator.track_threshold) { - visible_E += this_visible_energy; + track_visible_energy += this_visible_energy; } - ind ++; + } + + // do energy smearing + if (calculator.track_energy_distortion > 1e-4) { + track_visible_energy = rand.Gaus(track_visible_energy, track_visible_energy*calculator.track_energy_distortion); + // clamp to 0 + track_visible_energy = std::max(track_visible_energy, 0.); } // ...and showers + double shower_visible_energy = 0.; if (include_showers) { for (auto const &mcs: mcshower_list) { // ignore particles not from nu vertex, non primary particles, and uncharged particles @@ -220,30 +378,38 @@ double visibleEnergy(const simb::MCTruth &mctruth, const std::vector 1e-4) { - this_visible_energy = rand.Gaus(this_visible_energy, calculator.shower_energy_distortion*this_visible_energy); - // clamp to 0 - this_visible_energy = std::max(this_visible_energy, 0.); - } + if (this_visible_energy > calculator.shower_threshold) { - visible_E += this_visible_energy; + shower_visible_energy += this_visible_energy; } } } + // do energy smearing + if (calculator.shower_energy_distortion > 1e-4) { + shower_visible_energy = rand.Gaus(shower_visible_energy, shower_visible_energy*calculator.shower_energy_distortion); + // clamp to 0 + shower_visible_energy = std::max(shower_visible_energy, 0.); + } + + visible_E = track_visible_energy + shower_visible_energy; + // ...and primary lepton energy (for CC events) // only add in extra here if identified "lepton" is actually a lepton if (calculator.lepton_index >= 0 && (abs(mctrack_list[calculator.lepton_index].PdgCode()) == 13 || abs(mctrack_list[calculator.lepton_index].PdgCode()) == 11)) { - visible_E += smearLeptonEnergy(mctrack_list[calculator.lepton_index], calculator); + visible_E += smearLeptonEnergy(rand, mctrack_list[calculator.lepton_index], calculator); } return visible_E; } -double smearLeptonEnergy(const sim::MCTrack &mct, const VisibleEnergyCalculator &calculator) { - // setup distortion - TRandom rand; +double smearLeptonEnergy(TRandom &rand, const sim::MCTrack &mct, const VisibleEnergyCalculator &calculator) { + // if not contained and zero length, return + if (!calculator.lepton_contained && calculator.lepton_contained_length < 1e-4) { + // std::cout << "Out of detector lepton\n"; + return 0; + } double smearing_percentage; if (calculator.lepton_contained) { @@ -259,6 +425,8 @@ double smearLeptonEnergy(const sim::MCTrack &mct, const VisibleEnergyCalculator // clamp to 0 smeared_lepton_visible_energy = std::max(smeared_lepton_visible_energy, 0.); + // std::cout << "Lepton -- is_contained: " << calculator.lepton_contained << " length: " << calculator.lepton_contained_length << " smearing: " << smearing_percentage << " true E: " << lepton_visible_energy << " smeared E: " << smeared_lepton_visible_energy << std::endl; + return smeared_lepton_visible_energy; } @@ -271,6 +439,7 @@ double PDGMass(int pdg) { // regular particle if (pdg < 1000000000) { TParticlePDG* ple = PDGTable->GetParticle(pdg); + if (ple == NULL) return -1; return ple->Mass() * 1000.0; } // ion @@ -296,19 +465,24 @@ double PDGCharge(int pdg) { } } +bool isFromNuVertex(const simb::MCTruth& mc, const simb::MCParticle& mcp, float distance) { + TVector3 nuVtx = mc.GetNeutrino().Nu().Position().Vect(); + TVector3 partStart = mcp.Position().Vect(); + return (partStart - nuVtx).Mag() < distance; +} bool isFromNuVertex(const simb::MCTruth& mc, const sim::MCShower& show, float distance) { - TLorentzVector nuVtx = mc.GetNeutrino().Nu().Trajectory().Position(0); - TLorentzVector showStart = show.Start().Position(); + TVector3 nuVtx = mc.GetNeutrino().Nu().Trajectory().Position(0).Vect(); + TVector3 showStart = show.Start().Position().Vect(); return (showStart - nuVtx).Mag() < distance; } bool isFromNuVertex(const simb::MCTruth& mc, const sim::MCTrack& track, float distance) { - TLorentzVector nuVtx = mc.GetNeutrino().Nu().Trajectory().Position(0); - TLorentzVector trkStart = track.Start().Position(); + TVector3 nuVtx = mc.GetNeutrino().Nu().Trajectory().Position(0).Vect(); + TVector3 trkStart = track.Start().Position().Vect(); return (trkStart - nuVtx).Mag() < distance; } diff --git a/sbnanalysis/ana/SBNOsc/Utilities.h b/sbnanalysis/ana/SBNOsc/Utilities.h index 145d7ea5..dcc98ac2 100644 --- a/sbnanalysis/ana/SBNOsc/Utilities.h +++ b/sbnanalysis/ana/SBNOsc/Utilities.h @@ -22,6 +22,7 @@ #include "core/Event.hh" #include "ubcore/LLBasicTool/GeoAlgo/GeoAABox.h" +#include "TRandom.h" namespace ana { namespace SBNOsc { @@ -107,6 +108,8 @@ bool isFromNuVertex(const simb::MCTruth& mc, const sim::MCShower& show, bool isFromNuVertex(const simb::MCTruth& mc, const sim::MCTrack& track, float distance=5.0); +bool isFromNuVertex(const simb::MCTruth& mc, const simb::MCParticle& mcp, float distance=5.0); + /** * Calculate CCQE energy from associated lepton information (and optional @@ -175,16 +178,18 @@ struct VisibleEnergyCalculator { * * \return Visble energy in GeV. * */ -double visibleEnergy(const simb::MCTruth &mctruth, const std::vector &mctrack_list, const std::vector &mcshower_list, +double visibleEnergy(TRandom &rand, const simb::MCTruth &mctruth, const std::vector &mctrack_list, const std::vector &mcshower_list, const VisibleEnergyCalculator &calculator=VisibleEnergyCalculator(), bool include_showers=true); +double visibleEnergyProposal(TRandom &rand, const simb::MCTruth &mctruth, const std::vector &mctrack_list, const VisibleEnergyCalculator &calculator); +double visibleEnergyProposalMCParticles(TRandom &rand, const simb::MCTruth &mctruth, const std::vector mctrack_list, const VisibleEnergyCalculator &calculator); /** Get the smeared energy from a lepton. * \param mctrack The MCTrack object corresponding to the lepton * \param calculator Struct containing values to be used in energy calculation * * */ -double smearLeptonEnergy(const sim::MCTrack &mct, const VisibleEnergyCalculator &calculator=VisibleEnergyCalculator()); +double smearLeptonEnergy(TRandom &rand, const sim::MCTrack &mct, const VisibleEnergyCalculator &calculator=VisibleEnergyCalculator()); } // namespace SBNOsc } // namespace ana diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfig.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfig.fcl index 34366ac3..8e6fa84a 100644 --- a/sbnanalysis/ana/SBNOsc/config/NumuConfig.fcl +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfig.fcl @@ -1,28 +1,33 @@ BEGIN_PROLOG // standard configurations to be used in proposal sensitivity -eventweight_tags: ["fluxeventweight", "genieeventweight"] +# eventweight_tags: ["fluxeventweight", "genieeventweight"] +# eventweight_tags: ["fluxeventweightA", "genieeventweightA", "genieeventweightB", "fluxeventweightB"] +eventweight_tags: [] SBND_fiducial_volumes: [ { - xmin: -190.9 + xmin: -184.15 ymin: -185 zmin: 15. - xmax: 5.6 + xmax: 184.15 ymax: 185 - zmax: 415. - }, - { - xmin: 10.9 - ymin: -185 - zmin: 15. - xmax: 190.9 - ymax: 185 - zmax: 415. + zmax: 285. } ] SBND_active_volumes: [ + { + xmin: -199.15 + ymin: -200. + zmin: 0.0 + xmax: 199.15 + ymax: 200. + zmax: 365. + } +] + +SBND_Modern_active_volumes: [ { xmin: -199.15 ymin: -200. @@ -34,35 +39,22 @@ SBND_active_volumes: [ ] // fiducial volume of SBND used in proposal -SBND_Proposal_fiducial_volumes: [ +SBND_fiducial_volumes_no_cathode: [ { - xmin: -190.9 - ymin: -185 + xmin: -184.15 + xmax: -1.5 + ymin: -185. zmin: 15. - xmax: 5.6 - ymax: 185 - zmax: 290. + ymax: 185. + zmax: 285. }, { - xmin: 10.9 - ymin: -185 + xmin: 1.5 + xmax: 184.15 + ymin: -185. zmin: 15. - xmax: 190.9 - ymax: 185 - zmax: 290. - } - -] - -// active volume of SBND used in proposal -SBND_Proposal_active_volumes: [ - { - xmin: -199.15 - ymin: -200. - zmin: 0.0 - xmax: 199.15 - ymax: 200. - zmax: 365.0 + ymax: 185. + zmax: 285. } ] @@ -90,60 +82,81 @@ MicroBooNE_active_volumes: [ ICARUS_fiducial_volumes: [ { - xmin: -356.24 + xmin: -349.49 ymin: -158.41 - zmin: -894.950652 - xmax: -224.54 + zmin: -894.951 + xmax: -82.94 ymax: 128.41 - zmax: 799.950652 + zmax: 799.951 }, { - xmin: -207.89 + xmin: 82.94 ymin: -158.41 - zmin: -894.950652 - xmax: -76.19 + zmin: -894.951 + xmax: 349.49 ymax: 128.41 - zmax: 799.950652 + zmax: 799.951 + } +] + +ICARUS_fiducial_volumes_no_cathode: [ + // cathode at: +/- 216.215 + { + xmin: -349.49 + xmax: -217.715 + ymin: -158.41 + zmin: -894.951 + ymax: 128.41 + zmax: 799.951 }, { - xmin: 76.19 + xmin: -217.715 + xmax: -82.94 ymin: -158.41 - zmin: -894.950652 - xmax: 207.89 + zmin: -894.951 ymax: 128.41 - zmax: 799.950652 + zmax: 799.951 }, { - xmin: 224.54 + xmax: 217.715 + xmin: 82.94 + ymin: -158.41 + zmin: -894.951 + ymax: 128.41 + zmax: 799.951 + }, + { + xmax: 349.49 + xmin: 217.715 ymin: -158.41 - zmin: -894.950652 - xmax: 356.24 + zmin: -894.951 ymax: 128.41 - zmax: 799.950652 + zmax: 799.951 } ] ICARUS_active_volumes: [ { xmin: -364.49 - ymin: -173.41 - zmin: -909.950652 xmax: -67.94 + ymin: -173.41 + zmin: -909.951 ymax: 143.41 - zmax: 879.950652 + zmax: 879.951 }, { xmin: 67.94 ymin: -173.41 - zmin: -909.950652 + zmin: -909.951 xmax: 364.49 ymax: 143.41 - zmax: 879.950652 + zmax: 879.951 } ] standard_cuts: { acceptShakyTracks: false + trajPointLength: false doFVCut: true minLengthContainedTrack: 50 // value taken from proposal minLengthExitingTrack: 100 // value taken from proposal @@ -159,6 +172,103 @@ standard_cuts: { standard_cuts_noMEC: { cutKMEC: true acceptShakyTracks: false + trajPointLength: false + doFVCut: true + minLengthContainedTrack: 50 // value taken from proposal + minLengthExitingTrack: 100 // value taken from proposal + trackVisibleEnergyThreshold: 0.021 // [GeV] taken from pospoal + showerEnergyDistortion: 0.05 // not sure about this one + trackEnergyDistortion: 0.05 // taken from proposal + leptonEnergyDistortionContained: 0.02 // needs update + leptonEnergyDistortionLeavingA: 0.102 // from joseph email + leptonEnergyDistortionLeavingB: 0.000612 // from joseph email + selectionEfficiency: 0.8 // from proposal +} + +standard_cuts_noBKG_noMEC: { + cutKMEC: true + acceptShakyTracks: false + trajPointLength: false + doFVCut: true + minLengthExitingTrack: 100 // value taken from proposal + trackVisibleEnergyThreshold: 0.021 // [GeV] taken from pospoal + showerEnergyDistortion: 0.05 // not sure about this one + trackEnergyDistortion: 0.05 // taken from proposal + leptonEnergyDistortionContained: 0.02 // needs update + leptonEnergyDistortionLeavingA: 0.102 // from joseph email + leptonEnergyDistortionLeavingB: 0.000612 // from joseph email + selectionEfficiency: 0.8 // from proposal + backgroundRejection: 1.0 +} + +standard_cuts_noCuts: { + cutKMEC: true + acceptShakyTracks: false + trajPointLength: false + doFVCut: false + minLengthContainedTrack: -1 + minLengthExitingTrack: -1 + trackVisibleEnergyThreshold: 0.021 + showerEnergyDistortion: 0.05 // not sure about this one + trackEnergyDistortion: 0.05 // taken from proposal + leptonEnergyDistortionContained: 0.02 // needs update + leptonEnergyDistortionLeavingA: 0.102 // from joseph email + leptonEnergyDistortionLeavingB: 0.000612 // from joseph email + selectionEfficiency: 1.0 +} + +get_eff_cuts: { + cutKMEC: true + acceptShakyTracks: false + trajPointLength: false + doFVCut: true + minLengthContainedTrack: -1 // value taken from proposal + minLengthExitingTrack: 100 // value taken from proposal + trackVisibleEnergyThreshold: 0.021 // [GeV] taken from pospoal + showerEnergyDistortion: 0.05 // not sure about this one + trackEnergyDistortion: 0.05 // taken from proposal + leptonEnergyDistortionContained: 0.02 // needs update + leptonEnergyDistortionLeavingA: 0.102 // from joseph email + leptonEnergyDistortionLeavingB: 0.000612 // from joseph email + selectionEfficiency: 0.8 // from proposal +} + +standard_cuts_QE: { + # -1: Don't select mode + # 0: QE + # 1: Res + # 2: DIS + # 3: Coh + # 10: MEC + selectMode: 0 + selectCCNC: 0 + cutKMEC: true + acceptShakyTracks: false + trajPointLength: false + doFVCut: true + minLengthContainedTrack: -1 // value taken from proposal + minLengthExitingTrack: -1 // value taken from proposal + trackVisibleEnergyThreshold: 0.021 // [GeV] taken from pospoal + showerEnergyDistortion: 0.05 // not sure about this one + trackEnergyDistortion: 0.05 // taken from proposal + leptonEnergyDistortionContained: 0.02 // needs update + leptonEnergyDistortionLeavingA: 0.102 // from joseph email + leptonEnergyDistortionLeavingB: 0.000612 // from joseph email + selectionEfficiency: 0.8 // from proposal +} + +standard_cuts_RES: { + # -1: Don't select mode + # 0: QE + # 1: Res + # 2: DIS + # 3: Coh + # 10: MEC + selectMode: 1 + selectCCNC: 0 + cutKMEC: true + acceptShakyTracks: false + trajPointLength: false doFVCut: true minLengthContainedTrack: 50 // value taken from proposal minLengthExitingTrack: 100 // value taken from proposal @@ -169,6 +279,7 @@ standard_cuts_noMEC: { leptonEnergyDistortionLeavingA: 0.102 // from joseph email leptonEnergyDistortionLeavingB: 0.000612 // from joseph email selectionEfficiency: 0.8 // from proposal + } standard_uniformWeights: ["bnbcorrection_FluxHist"] diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigGetEffIcarus.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigGetEffIcarus.fcl new file mode 100644 index 00000000..a3fb8153 --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigGetEffIcarus.fcl @@ -0,0 +1,13 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_GetEff_Icarus.root" + +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::ICARUS_fiducial_volumes + active_volumes: @local::ICARUS_active_volumes + @table::get_eff_cuts + uniformWeights: @local::standard_uniformWeights +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigGetEffSBND.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigGetEffSBND.fcl new file mode 100644 index 00000000..46de6850 --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigGetEffSBND.fcl @@ -0,0 +1,15 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_GetEff_SBND.root" +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::SBND_Proposal_fiducial_volumes + active_volumes: @local::SBND_Proposal_active_volumes + @table::get_eff_cuts + # uniformWeights: @local::standard_uniformWeights + uniformWeights: [] + constantWeight: 1.21 // Scale up all event to account for r=110m (now) -> 100m (then) + # verbose: true +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigGetEffUboone.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigGetEffUboone.fcl new file mode 100644 index 00000000..4a341b2d --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigGetEffUboone.fcl @@ -0,0 +1,13 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_GetEff_Uboone.root" + +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::MicroBooNE_fiducial_volumes + active_volumes: @local::MicroBooNE_active_volumes + @table::get_eff_cuts + uniformWeights: @local::standard_uniformWeights +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigModernIcarus.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigModernIcarus.fcl index 2ac9007f..cdf16adb 100644 --- a/sbnanalysis/ana/SBNOsc/config/NumuConfigModernIcarus.fcl +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigModernIcarus.fcl @@ -5,7 +5,7 @@ OutputFile: "output_SBNOsc_NumuSelection_Modern_Icarus.root" MCWeightTags: @local::eventweight_tags NumuSelection: { - fiducial_volumes: @local::ICARUS_fiducial_volumes + fiducial_volumes: @local::ICARUS_fiducial_volumes_no_cathode active_volumes: @local::ICARUS_active_volumes @table::standard_cuts uniformWeights: @local::standard_uniformWeights diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigModernSBND.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigModernSBND.fcl index ac6acb1b..37c8ba2e 100644 --- a/sbnanalysis/ana/SBNOsc/config/NumuConfigModernSBND.fcl +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigModernSBND.fcl @@ -4,9 +4,10 @@ OutputFile: "output_SBNOsc_NumuSelection_Modern_SBND.root" MCWeightTags: @local::eventweight_tags NumuSelection: { - fiducial_volumes: @local::SBND_fiducial_volumes + fiducial_volumes: @local::SBND_fiducial_volumes_no_cathode active_volumes: @local::SBND_active_volumes @table::standard_cuts - uniformWeights: @local::standard_uniformWeights + # uniformWeights: @local::standard_uniformWeights + uniformWeights: [] } diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruIcarus.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruIcarus.fcl new file mode 100644 index 00000000..01eeee12 --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruIcarus.fcl @@ -0,0 +1,13 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_PassThru_Icarus.root" + +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::ICARUS_fiducial_volumes + active_volumes: @local::ICARUS_active_volumes + @table::standard_cuts_noCuts + uniformWeights: @local::standard_uniformWeights +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruModernSBND.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruModernSBND.fcl new file mode 100644 index 00000000..a115676a --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruModernSBND.fcl @@ -0,0 +1,15 @@ +#include "NumuConfig.fcl" + +OutputFile: "SBND_numu_all.root" +MCWeightTags: [] +ExperimentID: 0 + +NumuSelection: { + fiducial_volumes: @local::SBND_Modern_active_volumes + active_volumes: @local::SBND_Modern_active_volumes + @table::standard_cuts_noCuts + # uniformWeights: @local::standard_uniformWeights + uniformWeights: [] + # verbose: true +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruSBND.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruSBND.fcl new file mode 100644 index 00000000..cf95301f --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruSBND.fcl @@ -0,0 +1,15 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_PassThru_SBND.root" +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::SBND_Proposal_fiducial_volumes + active_volumes: @local::SBND_Proposal_active_volumes + @table::standard_cuts_noCuts + # uniformWeights: @local::standard_uniformWeights + uniformWeights: [] + constantWeight: 1.21 // Scale up all event to account for r=110m (now) -> 100m (then) + # verbose: true +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruUboone.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruUboone.fcl new file mode 100644 index 00000000..69696054 --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigPassThruUboone.fcl @@ -0,0 +1,13 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_PassThru_Uboone.root" + +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::MicroBooNE_fiducial_volumes + active_volumes: @local::MicroBooNE_active_volumes + @table::standard_cuts_noCuts + uniformWeights: @local::standard_uniformWeights +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigProposalIcarus.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigProposalIcarus.fcl index df0f611e..7e04c108 100644 --- a/sbnanalysis/ana/SBNOsc/config/NumuConfigProposalIcarus.fcl +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigProposalIcarus.fcl @@ -5,7 +5,7 @@ OutputFile: "output_SBNOsc_NumuSelection_Proposal_Icarus.root" MCWeightTags: @local::eventweight_tags NumuSelection: { - fiducial_volumes: @local::ICARUS_fiducial_volumes + fiducial_volumes: @local::ICARUS_fiducial_volumes active_volumes: @local::ICARUS_active_volumes @table::standard_cuts_noMEC uniformWeights: @local::standard_uniformWeights diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigProposalSBND.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigProposalSBND.fcl index 37decf9c..16b22558 100644 --- a/sbnanalysis/ana/SBNOsc/config/NumuConfigProposalSBND.fcl +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigProposalSBND.fcl @@ -4,10 +4,10 @@ OutputFile: "output_SBNOsc_NumuSelection_Proposal_SBND.root" MCWeightTags: @local::eventweight_tags NumuSelection: { - fiducial_volumes: @local::SBND_Proposal_fiducial_volumes - active_volumes: @local::SBND_Proposal_active_volumes + fiducial_volumes: @local::SBND_fiducial_volumes + active_volumes: @local::SBND_active_volumes @table::standard_cuts_noMEC - uniformWeights: @local::standard_uniformWeights + uniformWeights: [] constantWeight: 1.21 // Scale up all event to account for r=110m (now) -> 100m (then) } diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigProposal_BeamCenteredSBND.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigProposal_BeamCenteredSBND.fcl new file mode 100644 index 00000000..9ef5f11c --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigProposal_BeamCenteredSBND.fcl @@ -0,0 +1,15 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_Proposal_BeamCentered_SBND.root" +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::SBND_Proposal_fiducial_volumes_no_cathode + active_volumes: @local::SBND_Proposal_active_volumes_beam_centered + @table::standard_cuts_noMEC + # uniformWeights: @local::standard_uniformWeights + uniformWeights: [] + constantWeight: 1.21 // Scale up all event to account for r=110m (now) -> 100m (then) + # verbose: true +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigQEIcarus.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigQEIcarus.fcl new file mode 100644 index 00000000..d3a4ca75 --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigQEIcarus.fcl @@ -0,0 +1,13 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_QE_noL_Icarus.root" + +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::ICARUS_fiducial_volumes + active_volumes: @local::ICARUS_active_volumes + @table::standard_cuts_QE + uniformWeights: @local::standard_uniformWeights +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigQESBND.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigQESBND.fcl new file mode 100644 index 00000000..24d3e398 --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigQESBND.fcl @@ -0,0 +1,15 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_QE_noL_SBND.root" +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::SBND_Proposal_fiducial_volumes + active_volumes: @local::SBND_Proposal_active_volumes + @table::standard_cuts_QE + # uniformWeights: @local::standard_uniformWeights + uniformWeights: [] + constantWeight: 1.21 // Scale up all event to account for r=110m (now) -> 100m (then) + # verbose: true +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigQEUboone.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigQEUboone.fcl new file mode 100644 index 00000000..c4fab37d --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigQEUboone.fcl @@ -0,0 +1,13 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_QE_noL_Uboone.root" + +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::MicroBooNE_fiducial_volumes + active_volumes: @local::MicroBooNE_active_volumes + @table::standard_cuts_QE + uniformWeights: @local::standard_uniformWeights +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigRESIcarus.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigRESIcarus.fcl new file mode 100644 index 00000000..e2437a92 --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigRESIcarus.fcl @@ -0,0 +1,13 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_RES_Icarus.root" + +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::ICARUS_fiducial_volumes + active_volumes: @local::ICARUS_active_volumes + @table::standard_cuts_RES + uniformWeights: @local::standard_uniformWeights +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigRESSBND.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigRESSBND.fcl new file mode 100644 index 00000000..4b1b34cd --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigRESSBND.fcl @@ -0,0 +1,15 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_RES_SBND.root" +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::SBND_Proposal_fiducial_volumes + active_volumes: @local::SBND_Proposal_active_volumes + @table::standard_cuts_RES + # uniformWeights: @local::standard_uniformWeights + uniformWeights: [] + constantWeight: 1.21 // Scale up all event to account for r=110m (now) -> 100m (then) + # verbose: true +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfigRESUboone.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfigRESUboone.fcl new file mode 100644 index 00000000..aa4aca2d --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfigRESUboone.fcl @@ -0,0 +1,13 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_RES_Uboone.root" + +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::MicroBooNE_fiducial_volumes + active_volumes: @local::MicroBooNE_active_volumes + @table::standard_cuts_RES + uniformWeights: @local::standard_uniformWeights +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfignoBKGIcarus.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfignoBKGIcarus.fcl new file mode 100644 index 00000000..3ff7c629 --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfignoBKGIcarus.fcl @@ -0,0 +1,13 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_noBKG_Icarus.root" + +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::ICARUS_fiducial_volumes_no_cathode + active_volumes: @local::ICARUS_active_volumes + @table::standard_cuts_noBKG_noMEC + uniformWeights: @local::standard_uniformWeights +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfignoBKGSBND.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfignoBKGSBND.fcl new file mode 100644 index 00000000..400b71ee --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfignoBKGSBND.fcl @@ -0,0 +1,15 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_noBKG_SBND.root" +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::SBND_fiducial_volumes_no_cathode + active_volumes: @local::SBND_active_volumes + @table::standard_cuts_noBKG_noMEC + # uniformWeights: @local::standard_uniformWeights + uniformWeights: [] + constantWeight: 1.21 // Scale up all event to account for r=110m (now) -> 100m (then) + # verbose: true +} + diff --git a/sbnanalysis/ana/SBNOsc/config/NumuConfignoBKGUboone.fcl b/sbnanalysis/ana/SBNOsc/config/NumuConfignoBKGUboone.fcl new file mode 100644 index 00000000..1fa1a26d --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/config/NumuConfignoBKGUboone.fcl @@ -0,0 +1,13 @@ +#include "NumuConfig.fcl" + +OutputFile: "output_SBNOsc_NumuSelection_noBKG_Uboone.root" + +MCWeightTags: @local::eventweight_tags + +NumuSelection: { + fiducial_volumes: @local::MicroBooNE_fiducial_volumes + active_volumes: @local::MicroBooNE_active_volumes + @table::standard_cuts_noBKG_noMEC + uniformWeights: @local::standard_uniformWeights +} + diff --git a/sbnanalysis/ana/SBNOsc/config/SensitivityConfig.fcl b/sbnanalysis/ana/SBNOsc/config/SensitivityConfig.fcl index 0a87e778..710e835e 100644 --- a/sbnanalysis/ana/SBNOsc/config/SensitivityConfig.fcl +++ b/sbnanalysis/ana/SBNOsc/config/SensitivityConfig.fcl @@ -9,10 +9,8 @@ Covariance: { NumAltUnis: 100 UniformWeights: @local::uniform_weights - WeightKey: @local::intra_genie_plus_flux_weights + WeightKey: [@local::inter_genie_plus_flux_weights_proposal, @local::inter_genie_weights_proposal, @local::flux_weights_proposal] EnergyType: @local::proposal_numu_energy_type - SelectionEfficiency: @local::proposal_numu_selection_efficiency - BackgroundRejection: @local::proposal_numu_selection_rejection WeightMax: @local::weight_max } @@ -22,8 +20,6 @@ Sensitivity: { UniformWeights: @local::uniform_weights EnergyType: @local::proposal_numu_energy_type - SelectionEfficiency: @local::proposal_numu_selection_efficiency - BackgroundRejection: @local::proposal_numu_selection_rejection NumDm2: 300 LogDm2Lims: [ -2.0, 2.0 ] diff --git a/sbnanalysis/ana/SBNOsc/config/SensitivityConfigBase.fcl b/sbnanalysis/ana/SBNOsc/config/SensitivityConfigBase.fcl index 39e2b01a..e8865dd1 100644 --- a/sbnanalysis/ana/SBNOsc/config/SensitivityConfigBase.fcl +++ b/sbnanalysis/ana/SBNOsc/config/SensitivityConfigBase.fcl @@ -5,43 +5,46 @@ proposal_binlims: [ 0.2, 0.3, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0 sbnd_event_sample: { name: "SBND #nu_{#mu}" binlims: @local::proposal_binlims - Distance: 0.1 - DetX: [ -190.9, 190.9 ] - DetY: [ -185., 185. ] - DetZ: [ 15., 415. ] + Baseline: 10000. + BeamCenterX: 130. + DetX: [ -200., 200. ] + DetY: [ -200., 200. ] + DetZ: [ 0., 500. ] OscType: 2 TrueELims: [ 0.1, 7.1 ] NumTrueEBins: 25 NumDistanceBinsPerMeter: 10 - scalefactor: 16.67802551725374 // must be updated per input file + ScalePOT: 6.6e20 } uboone_event_sample: { name: "MicroBooNE #nu_{#mu}" binlims: @local::proposal_binlims - Distance: 0.47 - DetX: [ 13.45, 239.8 ] - DetY: [ -100.53, 102.47 ] - DetZ: [ 15.1, 956.9 ] + Baseline: 47000. + BeamCenterX: 130. + DetX: [ -5., 260. ] + DetY: [ -120., 120. ] + DetZ: [ 0., 1050. ] OscType: 2 TrueELims: [ 0.1, 7.1 ] NumTrueEBins: 25 NumDistanceBinsPerMeter: 10 - scalefactor: 1.3426232009357677 + ScalePOT: 1.32e21 } icarus_event_sample: { name: "ICARUS #nu_{#mu}" binlims: @local::proposal_binlims - Distance: 0.6 - DetX: [ -356.24, 356.24 ] - DetY: [ -158.41, 128.41 ] - DetZ: [ -894.950652, 799.950652] + Baseline: 60000. + BeamFrontZ: -909.951 + DetX: [ -365., 365. ] + DetY: [ -175., 150. ] + DetZ: [ -910., 900.] OscType: 2 TrueELims: [ 0.1, 7.1 ] NumTrueEBins: 25 NumDistanceBinsPerMeter: 10 - scalefactor: 5.5516161511654607 + ScalePOT: 6.6e20 } uniform_weights: ["bnbcorrection_FluxHist"] @@ -174,6 +177,117 @@ inter_genie_plus_flux_weights: [ "piplus_PrimaryHadronSWCentralSplineVariation" ] +inter_genie_weights_proposal: [ + "genie_ccresAxial_Genie", + "genie_ncresAxial_Genie", + "genie_qema_Genie", + "genie_NC_Genie", + "genie_NonResRvbarp1pi_Genie", + "genie_NonResRvbarp2pi_Genie", + "genie_NonResRvp1pi_Genie", + "genie_NonResRvp2pi_Genie" +] + +flux_weights_proposal: [ + "horncurrent_FluxUnisim", + "kplus_PrimaryHadronFeynmanScaling", + "kzero_PrimaryHadronSanfordWang", + "nucleoninexsec_FluxUnisim", + "nucleonqexsec_FluxUnisim", + "nucleontotxsec_FluxUnisim", + "piminus_PrimaryHadronSWCentralSplineVariation", + "pioninexsec_FluxUnisim", + "pionqexsec_FluxUnisim", + "piontotxsec_FluxUnisim", + "piplus_PrimaryHadronSWCentralSplineVariation" +] + +inter_genie_plus_flux_weights_proposal: [ + # uncertainties in axial mass + "genie_ccresAxial_Genie", + "genie_ncresAxial_Genie", + "genie_qema_Genie", + + # uncertainties in NC elastic axial mass / eta not considered + # "genie_ncelAxial_Genie", + # "genie_ncelEta_Genie", + + # neutral current normalization + "genie_NC_Genie", + + # non resonant background production + "genie_NonResRvbarp1pi_Genie", + "genie_NonResRvbarp2pi_Genie", + "genie_NonResRvp1pi_Genie", + "genie_NonResRvp2pi_Genie", + + # MISSING: non-resonant neutron interaction uncertainties + # Marco/Furmanski are unsure about this one + + # MISSING: DIS nuclear modificaiton + # Marco/Furmanski say this one is unimportant + # other DIS parameters -- not counted? + # "genie_DISAth_Genie", + # "genie_DISBth_Genie", + # "genie_DISCv1u_Genie", + # "genie_DISCv2u_Genie", + + # FSI reweighting -- not in proposal + # "genie_IntraNukeNabs_Genie", + # "genie_IntraNukeNcex_Genie", + # "genie_IntraNukeNel_Genie", + # "genie_IntraNukeNinel_Genie", + # "genie_IntraNukeNmfp_Genie", + # "genie_IntraNukeNpi_Genie", + # "genie_IntraNukePIabs_Genie", + # "genie_IntraNukePIcex_Genie", + # "genie_IntraNukePIel_Genie", + # "genie_IntraNukePIinel_Genie", + # "genie_IntraNukePImfp_Genie", + # "genie_IntraNukePIpi_Genie", + + # Fermi-gas -- not important? + # not in proposal or Marco/Furmanski + # "genie_FermiGasModelKf_Genie", + # "genie_FermiGasModelSf_Genie", + + # coherent axial mass and R0 + # Furmanski/Marco unsure about it + # not in proposal + # "genie_cohMA_Genie", + # "genie_cohR0_Genie", + + # Resonant decay branching ratios + # not in proposal + # "genie_ResDecayEta_Genie", + # "genie_ResDecayGamma_Genie", + # "genie_ResDecayTheta_Genie", + + # Hadron formation zone + # not in proposal + # "genie_FormZone_Genie", + + # vector mass for CC/NC resonant production + # not in proposal + # "genie_ccresVector_Genie", + # "genie_ncresVector_Genie", + + # all flux weights re in proposal + "expskin_FluxUnisim", + "kminus_PrimaryHadronNormalization", + "horncurrent_FluxUnisim", + "kplus_PrimaryHadronFeynmanScaling", + "kzero_PrimaryHadronSanfordWang", + "nucleoninexsec_FluxUnisim", + "nucleonqexsec_FluxUnisim", + "nucleontotxsec_FluxUnisim", + "piminus_PrimaryHadronSWCentralSplineVariation", + "pioninexsec_FluxUnisim", + "pionqexsec_FluxUnisim", + "piontotxsec_FluxUnisim", + "piplus_PrimaryHadronSWCentralSplineVariation" +] + proposal_numu_selection_efficiency: 0.8 proposal_numu_selection_rejection: 0.0 diff --git a/sbnanalysis/ana/SBNOsc/macro/example_macro.cc b/sbnanalysis/ana/SBNOsc/macro/example_macro.cc new file mode 100644 index 00000000..f077847b --- /dev/null +++ b/sbnanalysis/ana/SBNOsc/macro/example_macro.cc @@ -0,0 +1,42 @@ +#include "TFile.h" +#include "TH1F.h" +#include "TTreeReader.h" +#include "TTreeReaderValue.h" + +void read_numu_file() { + // open the input file + TFile f("/sbnd/data/users/gputnam/NuMu/outputs/sbnd_passthru_v1/SBND_numu_all.root"); + + // histogram which will store neutrino energies + TH1F *hist = new TH1F("", "", 50, 0, 5); + + // maximum number of events to read from file + int max_event = -1; + + // setup a "TTreeReader" to read from the "sbnana" Tree + // Each entry in the branch is a single "Event" in SBND + TTreeReader myReader("sbnana", &f); + + // Access to energy of each neutrino in the event. + // In SBND, there are sometimes multiple neutrinos interacting + // in the detector in each event, so we need to keep track of an + // array of energies. However, most of the time this array will + // just have one entry. + TTreeReaderArray neutrino_energies(myReader, "truth.neutrino.energy"); + + int event_ind = 0; + // Loop over the sbnana Tree + while (myReader.Next()) { + // break if we are done + if (max_event >= 0 && event_ind == max_event) break; + + // fill each neutrino enrgy in each interaction + for (double e: neutrino_energies) { + hist->Fill(e); + } + + event_ind += 1; + } + hist->Draw(); + +} diff --git a/sbnanalysis/core/CMakeLists.txt b/sbnanalysis/core/CMakeLists.txt index f6b5b429..d80bf0bc 100644 --- a/sbnanalysis/core/CMakeLists.txt +++ b/sbnanalysis/core/CMakeLists.txt @@ -20,6 +20,7 @@ include_directories($ENV{LARDATAOBJ_INC}) include_directories($ENV{NUSIMDATA_INC}) include_directories($ENV{FHICLCPP_INC}) include_directories($ENV{LARSIM_INC}) +include_directories($ENV{LARDATA_INC}) include_directories($ENV{SBNDCODE_INC}) include_directories($ENV{UBCORE_INC}) include_directories($ENV{ICARUSCODE_INC}) @@ -48,8 +49,6 @@ ROOT_GENERATE_DICTIONARY(sbn_eventdict Event.hh SubRun.hh LINKDEF linkdef.h) add_library(sbnanalysis_Event SHARED sbn_eventdict.cxx) target_link_libraries( sbnanalysis_Event - larsim_EventWeight_Base - larsim_EventWeight_Base_dict ${ROOT_LIBRARIES} ) diff --git a/sbnanalysis/core/Event.hh b/sbnanalysis/core/Event.hh index 0cfb1d46..32848a1d 100644 --- a/sbnanalysis/core/Event.hh +++ b/sbnanalysis/core/Event.hh @@ -53,15 +53,17 @@ public: public: /** Constructor. */ Neutrino() - : isnc(false), iscc(false), pdg(0), targetPDG(0), genie_intcode(0), - bjorkenX(kUnfilled), inelasticityY(kUnfilled), Q2(kUnfilled), - q0(kUnfilled), modq(kUnfilled), q0_lab(kUnfilled), modq_lab(kUnfilled), + : isnc(false), iscc(false), pdg(0), initpdg(0), targetPDG(0), + genie_intcode(0), bjorkenX(kUnfilled), inelasticityY(kUnfilled), + Q2(kUnfilled), q0(kUnfilled), + modq(kUnfilled), q0_lab(kUnfilled), modq_lab(kUnfilled), w(kUnfilled), t(kUnfilled), energy(kUnfilled), momentum(kUnfilled, kUnfilled, kUnfilled), parentPDG(0), parentDecayMode(0), parentDecayVtx(kUnfilled, kUnfilled, kUnfilled) {} bool isnc; //!< same as LArSoft "ccnc" - 0=CC, 1=NC bool iscc; //!< CC (true) or NC/interference (false) + int initpdg; //!< Initial PDG code of probe neutrino int pdg; //!< PDG code of probe neutrino int targetPDG; //!< PDG code of struck target int genie_intcode; //!< Interaction mode (as for LArSoft MCNeutrino::Mode() ) @@ -99,6 +101,23 @@ public: TVector3 momentum; //!< Three-momentum }; + + /** + * \struct Weight_t; + * \brief Container for an event weight; + */ + class Weight_t { + public: + Weight_t() : param_idx(kUnfilled), value(kUnfilled), weight(kUnfilled), universe(kUnfilled) {} + Weight_t(size_t _param_idx, size_t _universe, float _value, float _weight) + : param_idx(_param_idx), universe(_universe), value(_value), weight(_weight) {} + size_t param_idx; //!< Parameter name + size_t universe; //!< Universe index + float value; //!< Parameter value + float weight; //!< Weight + }; + + /** * \class Event::Interaction * \brief All truth information associated with one neutrino interaction @@ -109,15 +128,27 @@ public: FinalStateParticle lepton; //!< The primary final state lepton /** The other final state particles. */ + size_t nfinalstate; //!< Size of finalstate std::vector finalstate; //!< Final state particles /** * Event weights. * - * This is a map from the weight calculator name to the list of weights - * for all the sampled universes. + * As a map of parameter names to vectors of weights, as in LArSoft. */ - std::map > weights; + std::map > weightmap; + + /** + * Event weights. + * + * To reduce the depth of the event structure, this is a flattened list + * of weight objects that have the weight, the parameter, and the index + * of the universe. + */ + std::vector weights; + size_t nweights; //!< Number of weights + + std::map weight_indices; //!< IDs for parameter names }; /** @@ -161,7 +192,9 @@ public: }; Metadata metadata; //!< Event metadata + size_t ntruth; //!< Size of truth std::vector truth; //!< All truth interactions + size_t nreco; //!< Size of reco std::vector reco; //!< Reconstructed interactions Experiment experiment; //!< Experiment identifier diff --git a/sbnanalysis/core/PostProcessorBase.cxx b/sbnanalysis/core/PostProcessorBase.cxx index a6f43fc9..57a20fba 100644 --- a/sbnanalysis/core/PostProcessorBase.cxx +++ b/sbnanalysis/core/PostProcessorBase.cxx @@ -31,10 +31,19 @@ void PostProcessorBase::Run(std::vector inputFiles) { break; } + FileSetup(fEventTree); + + // process all subruns + fSubRunTree = (TTree *) f.Get("sbnsubrun"); + fSubRunTree->SetBranchAddress("subruns", &fSubRun); + for (int subrun_ind = 0; subrun_ind < fSubRunTree->GetEntries(); subrun_ind++) { + fSubRunTree->GetEntry(subrun_ind); + ProcessSubRun(fSubRun); + } + // set Event fEventTree = (TTree *) f.Get("sbnana"); fEventTree->SetBranchAddress("events", &fEvent); - FileSetup(fEventTree); // process all events for (int event_ind = 0; event_ind < fEventTree->GetEntries(); event_ind++) { diff --git a/sbnanalysis/core/PostProcessorBase.hh b/sbnanalysis/core/PostProcessorBase.hh index fdd150a4..faef6434 100644 --- a/sbnanalysis/core/PostProcessorBase.hh +++ b/sbnanalysis/core/PostProcessorBase.hh @@ -14,6 +14,7 @@ #include "gallery/Event.h" #include "Loader.hh" #include "Event.hh" +#include "SubRun.hh" class TBranch; class TFile; @@ -61,6 +62,8 @@ protected: */ virtual void ProcessEvent(const Event *event) = 0; + virtual void ProcessSubRun(const SubRun *subrun) {} + /** * Setup anything needed per file * @@ -92,6 +95,9 @@ protected: Event* fEvent; TTree* fEventTree; + + SubRun* fSubRun; + TTree* fSubRunTree; }; } // namespace core diff --git a/sbnanalysis/core/ProcessorBase.cxx b/sbnanalysis/core/ProcessorBase.cxx index fd48001a..b4ffd883 100644 --- a/sbnanalysis/core/ProcessorBase.cxx +++ b/sbnanalysis/core/ProcessorBase.cxx @@ -20,6 +20,9 @@ #include "ProcessorBase.hh" #include "ProviderManager.hh" +#include "larsim/MCCheater/BackTrackerService.h" +#include "lardataobj/Simulation/GeneratedParticleInfo.h" + namespace core { ProcessorBase::ProcessorBase() @@ -171,6 +174,19 @@ void ProcessorBase::Teardown() { fOutputFile->Close(); } +void ProcessorBase::SetupServices(gallery::Event& ev) { + if (fProviderManager != NULL) { + // reset the channels of the back tracker + fProviderManager->GetBackTrackerProvider()->ClearEvent(); + fProviderManager->GetBackTrackerProvider()->PrepSimChannels(ev); + + // reset information in particle inventory + fProviderManager->GetParticleInventoryProvider()->ClearEvent(); + fProviderManager->GetParticleInventoryProvider()->PrepParticleList(ev); + fProviderManager->GetParticleInventoryProvider()->PrepMCTruthList(ev); + fProviderManager->GetParticleInventoryProvider()->PrepTrackIdToMCTruthIndex(ev); + } +} void ProcessorBase::BuildEventTree(gallery::Event& ev) { // Add any new subruns to the subrun tree @@ -184,6 +200,11 @@ void ProcessorBase::BuildEventTree(gallery::Event& ev) { ev.getByLabel(fTruthTag, gtruths_handle); bool genie_truth_is_valid = gtruths_handle.isValid(); + // Get MCFlux information + auto const& mcfluxes = \ + *ev.getValidHandle >(fTruthTag); + assert(mctruths.size() == mcfluxes.size()); + // Get MCEventWeight information std::vector > > wghs; @@ -213,6 +234,7 @@ void ProcessorBase::BuildEventTree(gallery::Event& ev) { Event::Interaction interaction; auto const& mctruth = mctruths.at(i); + auto const& mcflux = mcfluxes.at(i); // TODO: What to do with cosmic MC? // For now, ignore them @@ -220,13 +242,29 @@ void ProcessorBase::BuildEventTree(gallery::Event& ev) { // Combine Weights if (!wghs.empty()) { - for (auto const &wgh: wghs) { - // Insert the weights for each individual EventWeight object into the - // Event class "master" weight list - interaction.weights.insert(wgh->at(i).fWeight.begin(), wgh->at(i).fWeight.end()); + size_t wgh_idx = 0; + // Loop through weight generators, which have a list of weights per truth + for (auto const& wgh : wghs) { + interaction.weightmap.insert(wgh->at(i).fWeight.begin(), wgh->at(i).fWeight.end()); + + // Iterate though the weight parameters for this generator, for the + // current MCTruth interaction + for (auto const& it : wgh->at(i).fWeight) { + // Loop through the weights for this parameter, to build the list + for (size_t j=0; jat(i); interaction.neutrino.parentPDG = flux.fptype; @@ -243,6 +281,7 @@ void ProcessorBase::BuildEventTree(gallery::Event& ev) { interaction.neutrino.isnc = nu.CCNC() && (nu.Mode() != simb::kWeakMix); interaction.neutrino.iscc = (!nu.CCNC()) && (nu.Mode() != simb::kWeakMix); interaction.neutrino.pdg = nu.Nu().PdgCode(); + interaction.neutrino.initpdg = mcflux.fntype; interaction.neutrino.targetPDG = nu.Target(); interaction.neutrino.genie_intcode = nu.Mode(); interaction.neutrino.bjorkenX = nu.X(); @@ -284,6 +323,8 @@ void ProcessorBase::BuildEventTree(gallery::Event& ev) { interaction.finalstate.push_back(fsp); } + interaction.nfinalstate = interaction.finalstate.size(); + // GENIE specific if (genie_truth_is_valid) { auto const& gtruth = gtruths_handle->at(i); @@ -299,6 +340,8 @@ void ProcessorBase::BuildEventTree(gallery::Event& ev) { fEvent->truth.push_back(interaction); } + + fEvent->ntruth = fEvent->truth.size(); } } // namespace core diff --git a/sbnanalysis/core/ProcessorBase.hh b/sbnanalysis/core/ProcessorBase.hh index 7ac656f3..e5c7f3da 100644 --- a/sbnanalysis/core/ProcessorBase.hh +++ b/sbnanalysis/core/ProcessorBase.hh @@ -127,6 +127,8 @@ protected: */ void BuildEventTree(gallery::Event& ev); + void SetupServices(gallery::Event& ev); + /** * Update subrun list to include subruns for this event's file. * diff --git a/sbnanalysis/core/ProcessorBlock.cxx b/sbnanalysis/core/ProcessorBlock.cxx index 8a9ab08f..96112d9d 100644 --- a/sbnanalysis/core/ProcessorBlock.cxx +++ b/sbnanalysis/core/ProcessorBlock.cxx @@ -29,6 +29,7 @@ void ProcessorBlock::ProcessFiles(std::vector filenames) { for (gallery::Event ev(filenames); !ev.atEnd(); ev.next()) { for (auto it : fProcessors) { it.first->BuildEventTree(ev); + it.first->SetupServices(ev); bool accept = it.first->ProcessEvent(ev, it.first->fEvent->truth, *it.first->fReco); diff --git a/sbnanalysis/core/ProviderManager.hh b/sbnanalysis/core/ProviderManager.hh index 1a1445a9..10e5fe5d 100644 --- a/sbnanalysis/core/ProviderManager.hh +++ b/sbnanalysis/core/ProviderManager.hh @@ -43,27 +43,27 @@ class ProviderManager { public: ProviderManager(Experiment det, std::string fcl=""); - const geo::GeometryCore* GetGeometryProvider() { + const geo::GeometryCore* GetGeometryProvider() const { return fGeometryProvider.get(); } - const detinfo::LArPropertiesStandard* GetLArPropertiesProvider() { + const detinfo::LArPropertiesStandard* GetLArPropertiesProvider() const { return fLArPropertiesProvider.get(); } - const detinfo::DetectorClocksStandard* GetDetectorClocksProvider() { + const detinfo::DetectorClocksStandard* GetDetectorClocksProvider() const { return fDetectorClocksProvider.get(); } - const detinfo::DetectorPropertiesStandard* GetDetectorPropertiesProvider() { + const detinfo::DetectorPropertiesStandard* GetDetectorPropertiesProvider() const { return fDetectorPropertiesProvider.get(); } - const cheat::BackTracker* GetBackTrackerProvider() { + cheat::BackTracker* GetBackTrackerProvider() const { return fBackTrackerProvider.get(); } - const cheat::ParticleInventory* GetParticleInventoryProvider() { + cheat::ParticleInventory* GetParticleInventoryProvider() const { return fParticleInventoryProvider.get(); } diff --git a/sbnanalysis/core/linkdef.h b/sbnanalysis/core/linkdef.h index 293d2aec..6ab576b3 100644 --- a/sbnanalysis/core/linkdef.h +++ b/sbnanalysis/core/linkdef.h @@ -11,7 +11,9 @@ #pragma link C++ class Event::Interaction+; #pragma link C++ class Event::Neutrino+; #pragma link C++ class Event::FinalStateParticle+; +#pragma link C++ class Event::Weight_t+; #pragma link C++ class std::map >+; +#pragma link C++ class std::map+; #pragma link C++ class std::vector+;