diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx index 31023b2..54593b9 100755 --- a/src/FCN/JointFCN.cxx +++ b/src/FCN/JointFCN.cxx @@ -1,1010 +1,1014 @@ #include "JointFCN.h" #include #include "FitUtils.h" //*************************************************** JointFCN::JointFCN(TFile* outfile) { //*************************************************** fOutputDir = gDirectory; if (outfile) FitPar::Config().out = outfile; std::vector samplekeys = Config::QueryKeys("sample"); LoadSamples(samplekeys); std::vector covarkeys = Config::QueryKeys("covar"); LoadPulls(covarkeys); fCurIter = 0; fMCFilled = false; fIterationTree = false; fDialVals = NULL; fNDials = 0; fUsingEventManager = FitPar::Config().GetParB("EventManager"); fOutputDir->cd(); } //*************************************************** JointFCN::JointFCN(std::vector samplekeys, TFile* outfile) { //*************************************************** fOutputDir = gDirectory; if (outfile) FitPar::Config().out = outfile; LoadSamples(samplekeys); fCurIter = 0; fMCFilled = false; fOutputDir->cd(); fIterationTree = false; fDialVals = NULL; fNDials = 0; fUsingEventManager = FitPar::Config().GetParB("EventManager"); fOutputDir->cd(); } //*************************************************** JointFCN::~JointFCN() { //*************************************************** // Delete Samples for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; delete exp; } for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; delete pull; } // Sort Tree if (fIterationTree) DestroyIterationTree(); if (fDialVals) delete fDialVals; if (fSampleLikes) delete fSampleLikes; }; //*************************************************** void JointFCN::CreateIterationTree(std::string name, FitWeight* rw) { //*************************************************** LOG(FIT) << " Creating new iteration container! " << std::endl; DestroyIterationTree(); fIterationTreeName = name; // Add sample likelihoods and ndof for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; std::string name = exp->GetName(); std::string liketag = name + "_likelihood"; fNameValues.push_back(liketag); fCurrentValues.push_back(0.0); std::string ndoftag = name + "_ndof"; fNameValues.push_back(ndoftag); fCurrentValues.push_back(0.0); } // Add Pull terms for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; std::string name = pull->GetName(); std::string liketag = name + "_likelihood"; fNameValues.push_back(liketag); fCurrentValues.push_back(0.0); std::string ndoftag = name + "_ndof"; fNameValues.push_back(ndoftag); fCurrentValues.push_back(0.0); } // Add Likelihoods fNameValues.push_back("total_likelihood"); fCurrentValues.push_back(0.0); fNameValues.push_back("total_ndof"); fCurrentValues.push_back(0.0); // Setup Containers fSampleN = fSamples.size() + fPulls.size(); fSampleLikes = new double[fSampleN]; fSampleNDOF = new int[fSampleN]; // Add Dials std::vector dials = rw->GetDialNames(); for (size_t i = 0; i < dials.size(); i++){ fNameValues.push_back( dials[i] ); fCurrentValues.push_back( 0.0 ); } fNDials = dials.size(); fDialVals = new double[fNDials]; // Set IterationTree Flag fIterationTree = true; } //*************************************************** void JointFCN::DestroyIterationTree() { //*************************************************** fIterationCount.clear(); fCurrentValues.clear(); fNameValues.clear(); fIterationValues.clear(); } //*************************************************** void JointFCN::WriteIterationTree() { //*************************************************** LOG(FIT) << "Writing iteration tree" << std::endl; // Make a new TTree TTree* itree = new TTree(fIterationTreeName.c_str(), fIterationTreeName.c_str()); double* vals = new double[fNameValues.size()]; int count = 0; itree->Branch("iteration",&count,"Iteration/I"); for (int i = 0; i < fNameValues.size(); i++) { itree->Branch( fNameValues[i].c_str(), &vals[i], (fNameValues[i] + "/D").c_str() ); } // Fill Iterations for (size_t i = 0; i < fIterationValues.size(); i++){ std::vector itervals = fIterationValues[i]; // Fill iteration state count = fIterationCount[i]; for (size_t j = 0; j < itervals.size(); j++){ vals[j] = itervals[j]; } // Save to TTree itree->Fill(); } // Write to file itree->Write(); } //*************************************************** void JointFCN::FillIterationTree(FitWeight* rw) { //*************************************************** // Loop over samples count int count = 0; for (int i = 0; i < fSampleN; i++){ fCurrentValues[count++] = fSampleLikes[i]; fCurrentValues[count++] = double(fSampleNDOF[i]); } // Fill Totals fCurrentValues[count++] = fLikelihood; fCurrentValues[count++] = double(fNDOF); // Loop Over Parameter Counts rw->GetAllDials(fDialVals, fNDials); for (int i = 0; i < fNDials; i++){ fCurrentValues[count++] = double(fDialVals[i]); } // Push Back Into Container fIterationCount.push_back( fCurIter ); fIterationValues.push_back(fCurrentValues); } //*************************************************** double JointFCN::DoEval(const double* x) { //*************************************************** // WEIGHT ENGINE fDialChanged = FitBase::GetRW()->HasRWDialChanged(x); FitBase::GetRW()->UpdateWeightEngine(x); if (fDialChanged) { FitBase::GetRW()->Reconfigure(); FitBase::EvtManager().ResetWeightFlags(); } if (LOG_LEVEL(REC)) { FitBase::GetRW()->Print(); } // SORT SAMPLES ReconfigureSamples(); // GET TEST STAT fLikelihood = GetLikelihood(); fNDOF = GetNDOF(); // PRINT PROGRESS LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood << std::endl; // UPDATE TREE if (fIterationTree) FillIterationTree(FitBase::GetRW()); return fLikelihood; } //*************************************************** int JointFCN::GetNDOF() { //*************************************************** int totaldof = 0; int count = 0; // Total number of Free bins in each MC prediction for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; int dof = exp->GetNDOF(); // Save Seperate DOF if (fIterationTree) { fSampleNDOF[count] = dof; } // Add to total totaldof += dof; count++; } // Loop over pulls for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; double dof = pull->GetLikelihood(); // Save seperate DOF if (fIterationTree) { fSampleNDOF[count] = dof; } // Add to total totaldof += dof; count++; } // Set Data Variable + if (fIterationTree){ fSampleNDOF[count] = totaldof; - + } return totaldof; } //*************************************************** double JointFCN::GetLikelihood() { //*************************************************** LOG(MIN) << std::left << std::setw(43) << "Getting likelihoods..." << " : " << "-2logL" << std::endl; // Loop and add up likelihoods in an uncorrelated way double like = 0.0; int count = 0; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; double newlike = exp->GetLikelihood(); int ndof = exp->GetNDOF(); // Save seperate likelihoods if (fIterationTree) { fSampleLikes[count] = newlike; } LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : " << newlike << "/" << ndof << std::endl; // Add Weight Scaling // like *= FitBase::GetRW()->GetSampleLikelihoodWeight(exp->GetName()); // Add to total like += newlike; count++; } // Loop over pulls for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; double newlike = pull->GetLikelihood(); // Save seperate likelihoods if (fIterationTree) { fSampleLikes[count] = newlike; } // Add to total like += newlike; count++; } // Set Data Variable fLikelihood = like; + if (fIterationTree){ + fSampleLikes[count] = fLikelihood; + } return like; }; void JointFCN::LoadSamples(std::vector samplekeys) { LOG(MIN) << "Loading Samples : " << samplekeys.size() << std::endl; for (size_t i = 0; i < samplekeys.size(); i++) { nuiskey key = samplekeys[i]; // Get Sample Options std::string samplename = key.GetS("name"); std::string samplefile = key.GetS("input"); std::string sampletype = key.GetS("type"); std::string fakeData = ""; LOG(MIN) << "Loading Sample : " << samplename << std::endl; fOutputDir->cd(); MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(key); if (!NewLoadedSample) { ERR(FTL) << "Could not load sample provided: " << samplename << std::endl; ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx" << std::endl; throw; } else { fSamples.push_back(NewLoadedSample); } } } //*************************************************** void JointFCN::LoadPulls(std::vector pullkeys) { //*************************************************** for (size_t i = 0; i < pullkeys.size(); i++) { nuiskey key = pullkeys[i]; std::string pullname = key.GetS("name"); std::string pullfile = key.GetS("input"); std::string pulltype = key.GetS("type"); fOutputDir->cd(); fPulls.push_back(new ParamPull(pullname, pullfile, pulltype)); } } //*************************************************** void JointFCN::ReconfigureSamples(bool fullconfig) { //*************************************************** int starttime = time(NULL); LOG(REC) << "Starting Reconfigure iter. " << this->fCurIter << std::endl; // std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled << // std::endl; // Event Manager Reconf if (fUsingEventManager) { if (!fullconfig and fMCFilled) ReconfigureFastUsingManager(); else ReconfigureUsingManager(); } else { // Loop over all Measurement Classes for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; // If RW Either do signal or full reconfigure. if (fDialChanged or !fMCFilled or fullconfig) { if (!fullconfig and fMCFilled) exp->ReconfigureFast(); else exp->Reconfigure(); // If RW Not needed just do normalisation } else { exp->Renormalise(); } } } // Loop over pulls and update for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; pull->Reconfigure(); } fMCFilled = true; LOG(MIN) << "Finished Reconfigure iter. " << fCurIter << " in " << time(NULL) - starttime << "s" << std::endl; fCurIter++; } //*************************************************** void JointFCN::ReconfigureSignal() { //*************************************************** ReconfigureSamples(false); } //*************************************************** void JointFCN::ReconfigureAllEvents() { //*************************************************** FitBase::GetRW()->Reconfigure(); FitBase::EvtManager().ResetWeightFlags(); ReconfigureSamples(true); } std::vector JointFCN::GetInputList() { std::vector InputList; fIsAllSplines = true; MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); std::vector subsamples = exp->GetSubSamples(); for (size_t i = 0; i < subsamples.size(); i++) { InputHandlerBase* inp = subsamples[i]->GetInput(); if (std::find(InputList.begin(), InputList.end(), inp) == InputList.end()) { if (subsamples[i]->GetInput()->GetType() != kSPLINEPARAMETER) fIsAllSplines = false; InputList.push_back(subsamples[i]->GetInput()); } } } return InputList; } std::vector JointFCN::GetSubSampleList() { std::vector SampleList; MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); std::vector subsamples = exp->GetSubSamples(); for (size_t i = 0; i < subsamples.size(); i++) { SampleList.push_back(subsamples[i]); } } return SampleList; } //*************************************************** void JointFCN::ReconfigureUsingManager() { //*************************************************** // 'Slow' Event Manager Reconfigure LOG(REC) << "Event Manager Reconfigure" << std::endl; int timestart = time(NULL); // Reset all samples MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ResetAll(); } // If we are siving signal, reset all containers. bool savesignal = (FitPar::Config().GetParB("SignalReconfigures")); if (savesignal) { // Reset all of our event signal vectors fSignalEventBoxes.clear(); fSignalEventFlags.clear(); fSampleSignalFlags.clear(); fSignalEventSplines.clear(); } // Make sure we have a list of inputs if (fInputList.empty()) { fInputList = GetInputList(); fSubSampleList = GetSubSampleList(); } // If all inputs are splines make sure the readers are told // they need to be reconfigured. std::vector::iterator inp_iter = fInputList.begin(); if (fIsAllSplines) { for (; inp_iter != fInputList.end(); inp_iter++) { InputHandlerBase* curinput = (*inp_iter); // Tell reader in each BaseEvent it needs a Reconfigure next weight calc. BaseFitEvt* curevent = curinput->FirstBaseEvent(); if (curevent->fSplineRead) { curevent->fSplineRead->SetNeedsReconfigure(true); } } } // MAIN INPUT LOOP ==================== int fillcount = 0; int inputcount = 0; inp_iter = fInputList.begin(); // Loop over each input in manager for (; inp_iter != fInputList.end(); inp_iter++) { InputHandlerBase* curinput = (*inp_iter); // Get event information FitEvent* curevent = curinput->FirstNuisanceEvent(); curinput->CreateCache(); int i = 0; int nevents = curinput->GetNEvents(); int countwidth = nevents / 5; // Start event loop iterating until we get a NULL pointer. while (curevent) { // Get Event Weight curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent); curevent->Weight = curevent->RWWeight * curevent->InputWeight; double rwweight = curevent->Weight; // std::cout << "RWWeight = " << curevent->RWWeight << " " << // curevent->InputWeight << std::endl; // Logging // std::cout << CHECKLOG(1) << std::endl; if (LOGGING(REC)) { if (i % countwidth == 0) { QLOG(REC, curinput->GetName() << " : Processed " << i << " events. [M, W] = [" << curevent->Mode << ", " << rwweight << "]"); } } // Setup flag for if signal found in at least one sample bool foundsignal = false; // Create a new signal bitset for this event std::vector signalbitset(fSubSampleList.size()); // Create a new signal box vector for this event std::vector signalboxes; // Start measurement iterator size_t measitercount = 0; std::vector::iterator meas_iter = fSubSampleList.begin(); // Loop over all subsamples (sub in JointMeas) for (; meas_iter != fSubSampleList.end(); meas_iter++) { MeasurementBase* curmeas = (*meas_iter); // Compare input pointers, to current input, skip if not. // Pointer tells us if it matches without doing ID checks. if (curinput != curmeas->GetInput()) { if (savesignal) { // Set bit to 0 as definitely not signal signalbitset[measitercount] = 0; } // Count up what measurement we are on. measitercount++; // Skip sample as input not signal. continue; } // Fill events for matching inputs. MeasurementVariableBox* box = curmeas->FillVariableBox(curevent); bool signal = curmeas->isSignal(curevent); curmeas->SetSignal(signal); curmeas->FillHistograms(curevent->Weight); // If its Signal tally up fills if (signal) { fillcount++; } // If we are saving signal/splines fill the bitset if (savesignal) { signalbitset[measitercount] = signal; } // If signal save a clone of the event box for use later. if (savesignal and signal) { foundsignal = true; signalboxes.push_back(box->CloneSignalBox()); } // Keep track of Measurement we are on. measitercount++; } // Once we've filled the measurements, if saving signal // push back if any sample flagged this event as signal if (savesignal) { fSignalEventFlags.push_back(foundsignal); } // Save the vector of signal boxes for this event if (savesignal and foundsignal) { fSignalEventBoxes.push_back(signalboxes); fSampleSignalFlags.push_back(signalbitset); } // If all inputs are splines we can save the spline coefficients // for fast in memory reconfigures later. if (fIsAllSplines and savesignal and foundsignal) { // Make temp vector to push back with std::vector coeff; for (size_t l = 0; l < (UInt_t)curevent->fSplineRead->GetNPar(); l++) { coeff.push_back(curevent->fSplineCoeff[l]); } // Push back to signal event splines. Kept in sync with // fSignalEventBoxes size. // int splinecount = fSignalEventSplines.size(); fSignalEventSplines.push_back(coeff); // if (splinecount % 1000 == 0) { // std::cout << "Pushed Back Coeff " << splinecount << " : "; // for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++) // { // std::cout << " " << fSignalEventSplines[splinecount][l]; // } // std::cout << std::endl; // } } // Clean up vectors once done with this event signalboxes.clear(); signalbitset.clear(); // Iterate to the next event. curevent = curinput->NextNuisanceEvent(); i++; } curinput->RemoveCache(); // Keep track of what input we are on. inputcount++; } // End of Event Loop =============================== // Now event loop is finished loop over all Measurements // Converting Binned events to XSec Distributions iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ConvertEventRates(); } // Print out statements on approximate memory usage for profiling. LOG(REC) << "Filled " << fillcount << " signal events." << std::endl; if (savesignal) { int mem = ( // sizeof(fSignalEventBoxes) + // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) + sizeof(MeasurementVariableBox1D) * fillcount) * 1E-6; LOG(REC) << " -> Saved " << fillcount << " signal boxes for faster access. (~" << mem << " MB)" << std::endl; if (fIsAllSplines and !fSignalEventSplines.empty()) { int splmem = sizeof(float) * fSignalEventSplines.size() * fSignalEventSplines[0].size() * 1E-6; LOG(REC) << " -> Saved " << fillcount << " " << fSignalEventSplines.size() << " spline sets into memory. (~" << splmem << " MB)" << std::endl; } } LOG(REC) << "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart << std::endl; // Check SignalReconfigures works for all samples if (savesignal) { double likefull = GetLikelihood(); ReconfigureFastUsingManager(); double likefast = GetLikelihood(); if (fabs(likefull - likefast) > 0.0001) { ERROR(FTL, "Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast); ERROR(FTL, "This means some samples you are using are not setup to use SignalReconfigures=1"); ERROR(FTL, "Please turn OFF signal reconfigures."); throw; } else { LOG(FIT) << "Likelihoods for FULL and FAST match. Will use FAST next time." << std::endl; } } // End of reconfigure return; }; //*************************************************** void JointFCN::ReconfigureFastUsingManager() { //*************************************************** LOG(FIT) << " -> Doing FAST using manager" << std::endl; // Get Start time for profilling int timestart = time(NULL); // Reset all samples MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ResetAll(); } // Check for saved variables if not do a full reconfigure. if (fSignalEventFlags.empty()) { ERR(WRN) << "Signal Flags Empty! Using normal manager." << std::endl; ReconfigureUsingManager(); return; } bool fFillNuisanceEvent = FitPar::Config().GetParB("FullEventOnSignalReconfigure"); // Setup fast vector iterators. std::vector::iterator inpsig_iter = fSignalEventFlags.begin(); std::vector >::iterator box_iter = fSignalEventBoxes.begin(); std::vector >::iterator spline_iter = fSignalEventSplines.begin(); std::vector >::iterator samsig_iter = fSampleSignalFlags.begin(); int splinecount = 0; // Setup stuff for logging int fillcount = 0; int nevents = fSignalEventFlags.size(); int countwidth = nevents / 20; // If All Splines tell splines they need a reconfigure. std::vector::iterator inp_iter = fInputList.begin(); if (fIsAllSplines) { LOG(REC) << "All Spline Inputs so using fast spline loop." << std::endl; for (; inp_iter != fInputList.end(); inp_iter++) { InputHandlerBase* curinput = (*inp_iter); // Tell each fSplineRead in BaseFitEvent to reconf next weight calc BaseFitEvt* curevent = curinput->FirstBaseEvent(); if (curevent->fSplineRead) curevent->fSplineRead->SetNeedsReconfigure(true); } } // Loop over all possible spline inputs double* coreeventweights = new double[fSignalEventBoxes.size()]; splinecount = 0; inp_iter = fInputList.begin(); inpsig_iter = fSignalEventFlags.begin(); spline_iter = fSignalEventSplines.begin(); // Loop over all signal flags // For each valid signal flag add one to splinecount // Get Splines from that count and add to weight // Add splinecount int sigcount = 0; splinecount = 0; // #pragma omp parallel for shared(splinecount,sigcount) for (uint iinput = 0; iinput < fInputList.size(); iinput++) { InputHandlerBase* curinput = fInputList[iinput]; BaseFitEvt* curevent = curinput->FirstBaseEvent(); for (int i = 0; i < curinput->GetNEvents(); i++) { double rwweight = 0.0; if (fSignalEventFlags[sigcount]) { // Get Event Info if (!fIsAllSplines) { if (fFillNuisanceEvent) curinput->GetNuisanceEvent(i); else curevent = curinput->GetBaseEvent(i); } else { curevent->fSplineCoeff = &fSignalEventSplines[splinecount][0]; } curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent); curevent->Weight = curevent->RWWeight * curevent->InputWeight; rwweight = curevent->Weight; coreeventweights[splinecount] = rwweight; if (splinecount % countwidth == 0) { LOG(REC) << "Processed " << splinecount << " event weights. W = " << rwweight << std::endl; } // #pragma omp atomic splinecount++; } // #pragma omp atomic sigcount++; } } LOG(SAM) << "Processed event weights." << std::endl; // #pragma omp barrier // Reset Iterators inpsig_iter = fSignalEventFlags.begin(); spline_iter = fSignalEventSplines.begin(); box_iter = fSignalEventBoxes.begin(); samsig_iter = fSampleSignalFlags.begin(); int nsplineweights = splinecount; splinecount = 0; // Start of Fast Event Loop ============================ // Start input iterators // Loop over number of inputs for (int ispline = 0; ispline < nsplineweights; ispline++) { double rwweight = coreeventweights[ispline]; // Get iterators for this event std::vector::iterator subsamsig_iter = (*samsig_iter).begin(); std::vector::iterator subbox_iter = (*box_iter).begin(); // Loop over all sub measurements. std::vector::iterator meas_iter = fSubSampleList.begin(); for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) { MeasurementBase* curmeas = (*meas_iter); // If event flagged as signal for this sample fill from the box. if (*subsamsig_iter) { curmeas->SetSignal(true); curmeas->FillHistogramsFromBox((*subbox_iter), rwweight); // Move onto next box if there is one. subbox_iter++; fillcount++; } } if (ispline % countwidth == 0) { LOG(REC) << "Filled " << ispline << " sample weights." << std::endl; } // Iterate over the main signal event containers. samsig_iter++; box_iter++; spline_iter++; splinecount++; } // End of Fast Event Loop =================== LOG(SAM) << "Filled sample distributions." << std::endl; // Now loop over all Measurements // Convert Binned events iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ConvertEventRates(); } // Cleanup coreeventweights if (fIsAllSplines) { delete coreeventweights; } // Print some reconfigure profiling. LOG(REC) << "Filled " << fillcount << " signal events." << std::endl; LOG(REC) << "Time taken ReconfigureFastUsingManager() : " << time(NULL) - timestart << std::endl; } //*************************************************** void JointFCN::Write() { //*************************************************** // Save a likelihood/ndof plot LOG(MIN) << "Writing likelihood plot.." << std::endl; std::vector likes; std::vector ndofs; std::vector names; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; double like = exp->GetLikelihood(); double ndof = exp->GetNDOF(); std::string name = exp->GetName(); likes.push_back(like); ndofs.push_back(ndof); names.push_back(name); } TH1D likehist = TH1D("likelihood_hist", "likelihood_hist", likes.size(), 0.0, double(likes.size())); TH1D ndofhist = TH1D("ndof_hist", "ndof_hist", ndofs.size(), 0.0, double(ndofs.size())); TH1D divhist = TH1D("likedivndof_hist", "likedivndof_hist", likes.size(), 0.0, double(likes.size())); for (size_t i = 0; i < likehist.GetNbinsX(); i++) { likehist.SetBinContent(i + 1, likes[i]); ndofhist.SetBinContent(i + 1, ndofs[i]); if (ndofs[i] != 0.0) { divhist.SetBinContent(i + 1, likes[i] / ndofs[i]); } likehist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); ndofhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); divhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); } likehist.Write(); ndofhist.Write(); divhist.Write(); // Loop over individual experiments and call Write LOG(MIN) << "Writing each of the data classes..." << std::endl; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; exp->Write(); } // Save Pull Terms for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; pull->Write(); } if (FitPar::Config().GetParB("EventManager")) { // Get list of inputs std::map fInputs = FitBase::EvtManager().GetInputs(); std::map::const_iterator iterInp; for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) { InputHandlerBase* input = (iterInp->second); input->GetFluxHistogram()->Write(); input->GetXSecHistogram()->Write(); input->GetEventHistogram()->Write(); } } }; //*************************************************** void JointFCN::SetFakeData(std::string fakeinput) { //*************************************************** LOG(MIN) << "Setting fake data from " << fakeinput << std::endl; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; exp->SetFakeDataValues(fakeinput); } return; } //*************************************************** void JointFCN::ThrowDataToy() { //*************************************************** for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; exp->ThrowDataToy(); } return; } diff --git a/src/FitBase/MeasurementBase.cxx b/src/FitBase/MeasurementBase.cxx index a77597e..235b470 100644 --- a/src/FitBase/MeasurementBase.cxx +++ b/src/FitBase/MeasurementBase.cxx @@ -1,561 +1,562 @@ // 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 "MeasurementBase.h" /* Constructor/Destructors */ //******************************************************************** // 2nd Level Constructor (Inherits From MeasurementBase.h) MeasurementBase::MeasurementBase(void) { //******************************************************************** fScaleFactor = 1.0; fMCFilled = false; fNoData = false; fInput = NULL; // Set the default values // After-wards this gets set in SetupMeasurement EnuMin = 0.; EnuMax = 1.E5; fMeasurementSpeciesType = kSingleSpeciesMeasurement; fEventVariables = NULL; fIsJoint = false; }; void MeasurementBase::FinaliseMeasurement() { // Used to setup default data hists, covars, etc. } //******************************************************************** // 2nd Level Destructor (Inherits From MeasurementBase.h) MeasurementBase::~MeasurementBase() { //******************************************************************** }; //******************************************************************** double MeasurementBase::TotalIntegratedFlux(std::string intOpt, double low, double high) { //******************************************************************** // Set Energy Limits if (low == -9999.9) low = this->EnuMin; if (high == -9999.9) high = this->EnuMax; return GetInput()->TotalIntegratedFlux(low, high, intOpt); }; //******************************************************************** double MeasurementBase::PredictedEventRate(std::string intOpt, double low, double high) { //******************************************************************** // Set Energy Limits if (low == -9999.9) low = this->EnuMin; if (high == -9999.9) high = this->EnuMax; return GetInput()->PredictedEventRate(low, high, intOpt) * 1E-38; }; //******************************************************************** void MeasurementBase::SetupInputs(std::string inputfile) { //******************************************************************** // Add this infile to the global manager if (FitPar::Config().GetParB("EventManager")) { fInput = FitBase::AddInput(fName, inputfile); } else { std::vector file_descriptor = GeneralUtils::ParseToStr(inputfile, ":"); if (file_descriptor.size() != 2) { ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfile << "\". expected \"FILETYPE:file.root\"" << std::endl; throw; } InputUtils::InputType inpType = InputUtils::ParseInputType(file_descriptor[0]); fInput = InputUtils::CreateInputHandler(fName, inpType, file_descriptor[1]); } fNEvents = fInput->GetNEvents(); // Expect INPUTTYPE:FileLocation(s) std::vector file_descriptor = GeneralUtils::ParseToStr(inputfile, ":"); if (file_descriptor.size() != 2) { ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfile << "\". expected \"FILETYPE:file.root\"" << std::endl; throw; } fInputType = InputUtils::ParseInputType(file_descriptor[0]); fInputFileName = file_descriptor[1]; if (EnuMin == 0 && EnuMax == 1.E5) { EnuMin = fInput->GetFluxHistogram()->GetBinLowEdge(1); EnuMax = fInput->GetFluxHistogram()->GetBinLowEdge( fInput->GetFluxHistogram()->GetNbinsX() + 1); } fFluxHist = fInput->GetFluxHistogram(); fEventHist = fInput->GetEventHistogram(); } //*********************************************** int MeasurementBase::GetInputID() { //*********************************************** return FitBase::GetInputID(fInputFileName); } //*********************************************** SampleSettings MeasurementBase::LoadSampleSettings(nuiskey samplekey) { //*********************************************** SampleSettings setting = SampleSettings(samplekey); fName = setting.GetS("name"); // Used as an initial setup function incase we need to do anything here. LOG(SAM) << "Loading Sample : " << setting.GetName() << std::endl; if (!fIsJoint) SetupInputs( setting.GetS("input") ); return setting; } //*********************************************** SampleSettings MeasurementBase::LoadSampleSettings(std::string name, std::string input, std::string type) { //*********************************************** nuiskey samplekey = Config::CreateKey("sample"); samplekey.SetS("name",name); samplekey.SetS("input",input); samplekey.SetS("type",type); return LoadSampleSettings(samplekey); } void MeasurementBase::FinaliseSampleSettings() { EnuMin = fSettings.GetD("enu_min"); EnuMax = fSettings.GetD("enu_max"); } //*********************************************** void MeasurementBase::Reconfigure() { //*********************************************** LOG(REC) << " Reconfiguring sample " << fName << std::endl; // Reset Histograms ResetExtraHistograms(); AutoResetExtraTH1(); this->ResetAll(); // FitEvent* cust_event = fInput->GetEventPointer(); int fNEvents = fInput->GetNEvents(); int countwidth = (fNEvents / 5); // MAIN EVENT LOOP FitEvent* cust_event = fInput->FirstNuisanceEvent(); int i = 0; int npassed = 0; while(cust_event){ cust_event->RWWeight = fRW->CalcWeight(cust_event); cust_event->Weight = cust_event->RWWeight * cust_event->InputWeight; Weight = cust_event->Weight; // Initialize fXVar = -999.9; fYVar = -999.9; fZVar = -999.9; Signal = false; Mode = cust_event->Mode; // Extract Measurement Variables this->FillEventVariables(cust_event); Signal = this->isSignal(cust_event); if (Signal) npassed++; GetBox()->SetX(fXVar); GetBox()->SetY(fYVar); GetBox()->SetZ(fZVar); GetBox()->SetMode(Mode); // GetBox()->fSignal = Signal; // Fill Histogram Values GetBox()->FillBoxFromEvent(cust_event); // this->FillExtraHistograms(GetBox(), Weight); this->FillHistogramsFromBox(GetBox(), Weight); // Print Out if (LOG_LEVEL(REC) && countwidth > 0 && !(i % countwidth)) { std::stringstream ss(""); ss.unsetf(std::ios_base::fixed); ss << std::setw(7) << std::right << i << "/" << fNEvents << " events (" << std::setw(2) << double(i) / double(fNEvents) * 100. << std::left << std::setw(5) << "%) " << "[S,X,Y,Z,M,W] = [" << std::fixed << std::setprecision(2) << std::right << Signal << ", " << std::setw(5) << fXVar << ", " << std::setw(5) << fYVar << ", " << std::setw(5) << fYVar << ", " << std::setw(3) << (int)Mode << ", " << std::setw(5) << Weight << "] " << std::endl; LOG(SAM) << ss.str(); } // iterate cust_event = fInput->NextNuisanceEvent(); i++; } LOG(SAM) << npassed << "/" << fNEvents << " passed selection " << std::endl; if (npassed == 0) { LOG(SAM) << "WARNING: NO EVENTS PASSED SELECTION!" << std::endl; } LOG(REC) << std::setw(10) << std::right << NSignal << "/" << fNEvents << " events passed selection + binning after reweight" << std::endl; // Finalise Histograms fMCFilled = true; this->ConvertEventRates(); } void MeasurementBase::FillHistogramsFromBox(MeasurementVariableBox* var, double weight) { fXVar = var->GetX(); fYVar = var->GetY(); fZVar = var->GetZ(); // Signal = var->fSignal; // Mode = var->fMode; Weight = weight; + fEventVariables = var; FillHistograms(); FillExtraHistograms(var, weight); } void MeasurementBase::FillHistograms(double weight){ Weight = weight * GetBox()->GetSampleWeight(); FillHistograms(); FillExtraHistograms(GetBox(), Weight); } MeasurementVariableBox* MeasurementBase::FillVariableBox(FitEvent* event) { GetBox()->Reset(); Mode = event->Mode; Weight = 1.0; //event->Weight; this->FillEventVariables(event); Signal = this->isSignal(event); GetBox()->FillBoxFromEvent(event); GetBox()->SetX(fXVar); GetBox()->SetY(fYVar); GetBox()->SetZ(fZVar); GetBox()->SetMode(event->Mode); GetBox()->SetSampleWeight(Weight); // GetBox()->fSignal = Signal; return GetBox(); } MeasurementVariableBox* MeasurementBase::GetBox() { if (!fEventVariables) fEventVariables = CreateBox(); return fEventVariables; } //*********************************************** void MeasurementBase::ReconfigureFast() { //*********************************************** this->Reconfigure(); } //*********************************************** void MeasurementBase::ConvertEventRates() { //*********************************************** AutoScaleExtraTH1(); ScaleExtraHistograms(GetBox()); this->ScaleEvents(); double normval = fRW->GetSampleNorm(this->fName); if (normval < 0.01 or normval > 10.0){ ERR(WRN) << "Norm Value inside MeasurementBase::ConvertEventRates() looks off!" << std::endl; ERR(WRN) << "It could have become out of sync with the minimizer norm list." << std::endl; ERR(WRN) << "Setting it to 1.0" << std::endl; normval = 1.0; } AutoNormExtraTH1(normval); NormExtraHistograms(GetBox(), normval); this->ApplyNormScale(normval); } //*********************************************** InputHandlerBase* MeasurementBase::GetInput() { //*********************************************** if (!fInput) { ERR(FTL) << "MeasurementBase::fInput not set. Please submit your command " "line options and input cardfile with a bug report to: " "nuisance@projects.hepforge.org" << std::endl; throw; } return fInput; }; //*********************************************** void MeasurementBase::Renormalise() { //*********************************************** // Called when the fitter has changed a measurements normalisation but not any // reweight dials // Means we don't have to call the time consuming reconfigure when this // happens. double norm = fRW->GetDialValue(this->fName + "_norm"); if ((this->fCurrentNorm == 0.0 and norm != 0.0) or not fMCFilled) { this->ReconfigureFast(); return; } if (this->fCurrentNorm == norm) return; this->ApplyNormScale(1.0 / this->fCurrentNorm); this->ApplyNormScale(norm); return; }; //*********************************************** void MeasurementBase::SetSignal(bool sig) { //*********************************************** Signal = sig; } //*********************************************** void MeasurementBase::SetSignal(FitEvent* evt) { //*********************************************** Signal = this->isSignal(evt); } //*********************************************** void MeasurementBase::SetWeight(double wght) { //*********************************************** Weight = wght; } //*********************************************** void MeasurementBase::SetMode(int md) { //*********************************************** Mode = md; } //*********************************************** std::vector MeasurementBase::GetFluxList() { //*********************************************** return GetInput()->GetFluxList(); } //*********************************************** std::vector MeasurementBase::GetEventRateList() { //*********************************************** return GetInput()->GetEventList(); } //*********************************************** std::vector MeasurementBase::GetXSecList() { //*********************************************** return GetInput()->GetXSecList(); } void MeasurementBase::ProcessExtraHistograms(int cmd, MeasurementVariableBox* vars, double weight) { // This should be overriden if we have extra histograms!!! // Add a flag to tell user this... return; } void MeasurementBase::FillExtraHistograms(MeasurementVariableBox* vars, double weight) { ProcessExtraHistograms(kCMD_Fill, vars, weight); } void MeasurementBase::ScaleExtraHistograms(MeasurementVariableBox* vars) { ProcessExtraHistograms(kCMD_Scale, vars, 1.0); } void MeasurementBase::ResetExtraHistograms() { ProcessExtraHistograms(kCMD_Reset, NULL, 1.0); } void MeasurementBase::NormExtraHistograms(MeasurementVariableBox* vars, double norm) { ProcessExtraHistograms(kCMD_Norm, vars, norm); } void MeasurementBase::WriteExtraHistograms() { ProcessExtraHistograms(kCMD_Write, NULL, 1.00); } void MeasurementBase::SetAutoProcessTH1(TH1* hist, int c1, int c2, int c3, int c4, int c5) { FakeStack* fake = new FakeStack(hist); SetAutoProcessTH1(fake, c1, c2, c3, c4, c5); // Need to add a destroy command! } void MeasurementBase::SetAutoProcess(TH1* hist, int c1, int c2, int c3, int c4, int c5) { FakeStack* fake = new FakeStack(hist); SetAutoProcessTH1(fake, c1, c2, c3, c4, c5); // Need to add a destroy command! } void MeasurementBase::SetAutoProcess(TGraph* g, int c1, int c2, int c3, int c4, int c5) { FakeStack* fake = new FakeStack(g); SetAutoProcessTH1(fake, c1, c2, c3, c4, c5); // Need to add a destroy command! } void MeasurementBase::SetAutoProcess(TF1* f, int c1, int c2, int c3, int c4, int c5) { FakeStack* fake = new FakeStack(f); SetAutoProcessTH1(fake, c1, c2, c3, c4, c5); // Need to add a destroy command! } void MeasurementBase::SetAutoProcess(StackBase* hist, int c1, int c2, int c3, int c4, int c5){ SetAutoProcessTH1(hist, c1, c2, c3, c4, c5); } void MeasurementBase::SetAutoProcessTH1(StackBase* hist, int c1, int c2, int c3, int c4, int c5) { // Set Defaults // int ncommands = kCMD_extraplotflags; bool autoflags[5]; autoflags[0] = false; autoflags[1] = false; autoflags[2] = false; autoflags[3] = false; autoflags[4] = false; int givenflags[5]; givenflags[0] = c1; givenflags[1] = c2; givenflags[2] = c3; givenflags[3] = c4; givenflags[4] = c5; fExtraTH1s[hist] = std::vector(5,0); // Setup a default one. if (c1 == -1 && c2 == -1 && c3 == -1 && c4 == -1 && c5 == -1){ fExtraTH1s[hist][kCMD_Reset] = 1; fExtraTH1s[hist][kCMD_Scale] = 1; fExtraTH1s[hist][kCMD_Norm] = 1; fExtraTH1s[hist][kCMD_Write] = 1; } for (int i = 0; i < 5; i++) { switch (givenflags[i]) { // Skip over... case -1: break; case kCMD_Reset: case kCMD_Scale: case kCMD_Norm: case kCMD_Write: fExtraTH1s[hist][givenflags[i]] = 1; break; case kCMD_Fill: ERR(FTL) << "Can't auto fill yet!" << std::endl; autoflags[givenflags[i]] = 1; break; default: break; } } // LOG(SAM) << "AutoProcessing " << hist->GetName() << std::endl; }; void MeasurementBase::AutoFillExtraTH1() { ERR(FTL) << "Can't auto fill yet! it's too inefficent!" << std::endl; return; } void MeasurementBase::AutoResetExtraTH1() { for (std::map >::iterator iter = fExtraTH1s.begin(); iter != fExtraTH1s.end(); iter++) { if (!((*iter).second)[kCMD_Reset]) continue; (*iter).first->Reset(); } }; void MeasurementBase::AutoScaleExtraTH1() { for (std::map >::iterator iter = fExtraTH1s.begin(); iter != fExtraTH1s.end(); iter++) { if (!((*iter).second)[kCMD_Scale]) continue; if (fIsNoWidth){ (*iter).first->Scale(fScaleFactor); } else { (*iter).first->Scale(fScaleFactor, "width"); } } }; void MeasurementBase::AutoNormExtraTH1(double norm) { double sfactor = 0.0; if (norm != 0.0) sfactor = 1.0 / norm; for (std::map >::iterator iter = fExtraTH1s.begin(); iter != fExtraTH1s.end(); iter++) { if (!((*iter).second)[kCMD_Norm]) continue; (*iter).first->Scale(sfactor); } }; void MeasurementBase::AutoWriteExtraTH1() { for (std::map >::iterator iter = fExtraTH1s.begin(); iter != fExtraTH1s.end(); iter++) { if (!(((*iter).second)[kCMD_Write])) continue; (*iter).first->Write(); } }; diff --git a/src/MINERvA/MINERvAVariableBoxes.h b/src/MINERvA/MINERvAVariableBoxes.h index 56288f8..06ed9ed 100644 --- a/src/MINERvA/MINERvAVariableBoxes.h +++ b/src/MINERvA/MINERvAVariableBoxes.h @@ -1,30 +1,57 @@ #ifndef MINERvA_VARIABLES_BOX_H #define MINERvA_VARIABLES_BOX_H #include "MeasurementVariableBox.h" #include "MeasurementVariableBox1D.h" #include "MeasurementVariableBox2D.h" /*! * \addtogroup FitBase * @{ */ /// Custom box used to also save All Pion Tpi for each event. class NTpiVariableBox1D : public MeasurementVariableBox1D { public: - inline NTpiVariableBox1D() { Reset(); }; + inline NTpiVariableBox1D() { }; inline void Reset() { fTpiVect.clear(); } + + inline MeasurementVariableBox* CloneSignalBox(){ + NTpiVariableBox1D* box = new NTpiVariableBox1D(); + box->fX = this->fX; + box->fSampleWeight = this->fSampleWeight; + + box->fTpiVect.clear(); + for (int i = 0; i < this->fTpiVect.size(); i++){ + box->fTpiVect.push_back( this->fTpiVect[i] ); + } + return box; + } + inline void Print(){ + std::cout << "Box Print Size : " << this->fTpiVect.size() << std::endl; + } + std::vector fTpiVect; }; /// Custom box used to also save All Pion Tpi for each event. class NthpiVariableBox1D : public MeasurementVariableBox1D { public: - inline NthpiVariableBox1D() { Reset(); }; + inline NthpiVariableBox1D() { }; inline void Reset() { fthpiVect.clear(); } + inline MeasurementVariableBox* CloneSignalBox(){ + NthpiVariableBox1D* box = new NthpiVariableBox1D(); + box->fX = this->fX; + box->fSampleWeight = this->fSampleWeight; + + box->fthpiVect.clear(); + for (int i = 0; i < this->fthpiVect.size(); i++){ + box->fthpiVect.push_back( this->fthpiVect[i] ); + } + return box; + } std::vector fthpiVect; }; #endif diff --git a/src/MINERvA/MINERvA_CCNpip_XSec_1DTpi_nu.cxx b/src/MINERvA/MINERvA_CCNpip_XSec_1DTpi_nu.cxx index 51fa133..74f9d9a 100644 --- a/src/MINERvA/MINERvA_CCNpip_XSec_1DTpi_nu.cxx +++ b/src/MINERvA/MINERvA_CCNpip_XSec_1DTpi_nu.cxx @@ -1,267 +1,267 @@ // 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 "MINERvA_SignalDef.h" #include "MINERvA_CCNpip_XSec_1DTpi_nu.h" //******************************************************************** MINERvA_CCNpip_XSec_1DTpi_nu::MINERvA_CCNpip_XSec_1DTpi_nu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCNpip_XSec_1DTpi_nu sample. \n" \ "Target: CH \n" \ "Flux: MINERvA Forward Horn Current nue + nuebar \n" \ "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("T_{#pi} (MeV)"); fSettings.SetYTitle("(1/T#Phi) dN_{#pi}/dT_{#pi} (cm^{2}/MeV/nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL"); fSettings.SetEnuRange(1.5, 10.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); fFullPhaseSpace = !fSettings.Found("name", "_20deg"); fFluxCorrection = fSettings.Found("name", "fluxcorr"); fUpdatedData = !fSettings.Found("name", "2015"); fSettings.SetTitle("MINERvA_CCNpip_XSec_1DTpi_nu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- // Full Phase Space if (fFullPhaseSpace) { // 2016 release if (fUpdatedData) { SetDataFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2016/nu-ccNpi+-xsec-pion-kinetic-energy.csv"); // MINERvA has the error quoted as a percentage of the cross-section // Need to make this into an absolute error before we go from correlation // matrix -> covariance matrix since it depends on the error in the ith // bin for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) { fDataHist->SetBinError(i + 1, fDataHist->GetBinContent(i + 1) * (fDataHist->GetBinError(i + 1) / 100.)); } // This is a correlation matrix, not covariance matrix, so needs to be // converted SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2016/nu-ccNpi+-correlation-pion-kinetic-energy.csv"); // 2015 release } else { // If we're doing shape only if (fIsShape) { SetDataFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_shape.txt"); SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_shape_cov.txt"); // If we're doing full cross-section } else { SetDataFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi.txt"); SetCorrelationFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_cov.txt"); } } // Restricted Phase Space } else { // Only 2015 data released restricted muon phase space cross-section // unfortunately if (fUpdatedData) { ERR(FTL) << fName << " has no updated 2016 data for restricted phase space! Using 2015 data." << std::endl; throw; } // If we're using the shape only data if (fIsShape) { SetDataFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_20deg_shape.txt"); SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_20deg_shape_cov.txt"); // Or total cross-section } else { SetDataFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_20deg.txt"); SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_20deg_cov.txt"); } } // Scale the MINERvA data to account for the flux difference // Adjust MINERvA data to flux correction; roughly a 11% normalisation increase in data // Please change when MINERvA releases new data! if (fFluxCorrection) { for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) { fDataHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) * 1.11); } } // Make some auxillary helper plots onePions = (TH1D*)(fDataHist->Clone()); onePions->SetNameTitle((fName + "_1pions").c_str(), (fName + "_1pions" + fPlotTitles).c_str()); SetAutoProcessTH1(onePions, kCMD_Reset, kCMD_Scale, kCMD_Norm); twoPions = (TH1D*)(fDataHist->Clone()); twoPions->SetNameTitle((fName + "_2pions").c_str(), (fName + "_2pions;" + fPlotTitles).c_str()); SetAutoProcessTH1(twoPions, kCMD_Reset, kCMD_Scale, kCMD_Norm); threePions = (TH1D*)(fDataHist->Clone()); threePions->SetNameTitle((fName + "_3pions").c_str(), (fName + "_3pions" + fPlotTitles).c_str()); SetAutoProcessTH1(threePions, kCMD_Reset, kCMD_Scale, kCMD_Norm); morePions = (TH1D*)(fDataHist->Clone()); morePions->SetNameTitle((fName + "_4pions").c_str(), (fName + "_4pions" + fPlotTitles).c_str()); SetAutoProcessTH1(morePions, kCMD_Reset, kCMD_Scale, kCMD_Norm); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** // Here we have to fill for every pion we find in the event void MINERvA_CCNpip_XSec_1DTpi_nu::FillEventVariables(FitEvent *event) { //******************************************************************** if (event->NumFSParticle(211) == 0 && event->NumFSParticle(-211) == 0) return; if (event->NumFSParticle(13) == 0) return; // Need to make this use event boxes // Clear out the vectors GetPionBox()->Reset(); TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; // Loop over the particle stack for (unsigned int j = 2; j < event->Npart(); ++j) { // Only include alive particles if (event->GetParticleState(j) != kFinalState) continue; int PID = (event->PartInfo(j))->fPID; // Pick up the charged pions in the event if (abs(PID) == 211) { double ppi = FitUtils::T(event->PartInfo(j)->fP) * 1000.; GetPionBox()->fTpiVect.push_back(ppi); } } fXVar = 0; return; }; //******************************************************************** // The last bool refers to if we're using restricted phase space or not bool MINERvA_CCNpip_XSec_1DTpi_nu::isSignal(FitEvent *event) { //******************************************************************** // Last false refers to that this is NOT the restricted MINERvA phase space, // in which only forward-going muons are accepted return SignalDef::isCCNpip_MINERvA(event, EnuMin, EnuMax, !fFullPhaseSpace, !fUpdatedData); } //******************************************************************** // Need to override FillHistograms() here because we fill the histogram N_pion // times void MINERvA_CCNpip_XSec_1DTpi_nu::FillHistograms() { //******************************************************************** if (Signal) { - unsigned int nPions = GetPionBox()->fTpiVect.size(); + unsigned int nPions = GetPionBox()->fTpiVect.size(); // Need to loop over all the pions in the sample for (size_t k = 0; k < nPions; ++k) { double tpi = GetPionBox()->fTpiVect[k]; this->fMCHist->Fill(tpi, Weight); this->fMCFine->Fill(tpi, Weight); this->fMCStat->Fill(tpi, 1.0); if (nPions == 1) { onePions->Fill(tpi, Weight); } else if (nPions == 2) { twoPions->Fill(tpi, Weight); } else if (nPions == 3) { threePions->Fill(tpi, Weight); } else if (nPions > 3) { morePions->Fill(tpi, Weight); } if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, tpi, Weight); } } } //******************************************************************** void MINERvA_CCNpip_XSec_1DTpi_nu::ScaleEvents() { //******************************************************************** Measurement1D::ScaleEvents(); onePions->Scale(this->fScaleFactor, "width"); twoPions->Scale(this->fScaleFactor, "width"); threePions->Scale(this->fScaleFactor, "width"); morePions->Scale(this->fScaleFactor, "width"); return; } //******************************************************************** void MINERvA_CCNpip_XSec_1DTpi_nu::Write(std::string drawOpts) { //******************************************************************** Measurement1D::Write(drawOpts); // Make an auto processed pion stack // Draw the npions stack onePions->SetTitle("1#pi"); onePions->SetLineColor(kBlack); // onePions->SetFillStyle(0); onePions->SetFillColor(onePions->GetLineColor()); twoPions->SetTitle("2#pi"); twoPions->SetLineColor(kRed); // twoPions->SetFillStyle(0); twoPions->SetFillColor(twoPions->GetLineColor()); threePions->SetTitle("3#pi"); threePions->SetLineColor(kGreen); // threePions->SetFillStyle(0); threePions->SetFillColor(threePions->GetLineColor()); morePions->SetTitle(">3#pi"); morePions->SetLineColor(kBlue); // morePions->SetFillStyle(0); morePions->SetFillColor(morePions->GetLineColor()); THStack pionStack = THStack((fName + "_pionStack").c_str(), (fName + "_pionStack").c_str()); pionStack.Add(onePions); pionStack.Add(twoPions); pionStack.Add(threePions); pionStack.Add(morePions); pionStack.Write(); return; }