diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx index 6ce228b..85ed5a1 100755 --- a/src/FCN/JointFCN.cxx +++ b/src/FCN/JointFCN.cxx @@ -1,1123 +1,1154 @@ #include "JointFCN.h" #include "FitUtils.h" #include //*************************************************** JointFCN::JointFCN(TFile *outfile) { //*************************************************** fOutputDir = gDirectory; if (outfile) Config::Get().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) Config::Get().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) { //*************************************************** NUIS_LOG(FIT, " Creating new iteration container! "); 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() { //*************************************************** NUIS_LOG(FIT, "Writing iteration tree"); // 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 (size_t 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) { //*************************************************** + double *par_vals = new double[fNPars]; + + for (int i = 0; i < fNPars; ++i) { + if (fMirroredParams.count(i)) { + if (!fMirroredParams[i].mirror_above && + (x[i] < fMirroredParams[i].mirror_value)) { + double xabove = fMirroredParams[i].mirror_value - x[i]; + par_vals[i] = fMirroredParams[i].mirror_value + xabove; + std::cout << "\t--Parameter " << i << " mirrored from " << x[i] + << " -> " << par_vals[i] << std::endl; + } else if (fMirroredParams[i].mirror_above && + (x[i] >= fMirroredParams[i].mirror_value)) { + double xabove = x[i] - fMirroredParams[i].mirror_value; + par_vals[i] = fMirroredParams[i].mirror_value - xabove; + std::cout << "\t--Parameter " << i << " mirrored from " << x[i] + << " -> " << par_vals[i] << std::endl; + } else { + par_vals[i] = x[i]; + } + } else { + par_vals[i] = x[i]; + } + } + // WEIGHT ENGINE - fDialChanged = FitBase::GetRW()->HasRWDialChanged(x); - FitBase::GetRW()->UpdateWeightEngine(x); + fDialChanged = FitBase::GetRW()->HasRWDialChanged(par_vals); + FitBase::GetRW()->UpdateWeightEngine(par_vals); 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 - NUIS_LOG(FIT, "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood); + NUIS_LOG(FIT, + "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood); // UPDATE TREE if (fIterationTree) FillIterationTree(FitBase::GetRW()); + delete[] par_vals; + 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() { //*************************************************** NUIS_LOG(MIN, std::left << std::setw(43) << "Getting likelihoods..." - << " : " - << "-2logL"); + << " : " + << "-2logL"); // 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; } NUIS_LOG(MIN, "-> " << std::left << std::setw(40) << exp->GetName() << " : " - << newlike << "/" << ndof); + << newlike << "/" << ndof); // 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) { NUIS_LOG(MIN, "Loading Samples : " << samplekeys.size()); 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 = ""; NUIS_LOG(MIN, "Loading Sample : " << samplename); fOutputDir->cd(); MeasurementBase *NewLoadedSample = SampleUtils::CreateSample(key); if (!NewLoadedSample) { NUIS_ERR(FTL, "Could not load sample provided: " << samplename); NUIS_ERR(FTL, "Check spelling with that in src/FCN/SampleList.cxx"); 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); NUIS_LOG(REC, "Starting Reconfigure iter. " << this->fCurIter); - // std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled << - // std::endl; - // Event Manager Reconf + // std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled + // << std::endl; Event Manager Reconf if (fUsingEventManager) { if (!fullconfig && 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; NUIS_LOG(MIN, "Finished Reconfigure iter. " << fCurIter << " in " - << time(NULL) - starttime << "s"); + << time(NULL) - starttime << "s"); 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 NUIS_LOG(REC, "Event Manager Reconfigure"); 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. + // 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 / 10; // Start event loop iterating until we get a NULL pointer. while (curevent) { // Get Event Weight // The reweighting weight curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent); // The Custom weight and reweight curevent->Weight = curevent->RWWeight * curevent->InputWeight * curevent->CustomWeight; if (LOGGING(REC)) { if (countwidth && (i % countwidth == 0)) { NUIS_LOG(REC, curinput->GetName() - << " : Processed " << i << " events. [M, W] = [" - << curevent->Mode << ", " << curevent->Weight << "]"); + << " : Processed " << i << " events. [M, W] = [" + << curevent->Mode << ", " << curevent->Weight + << "]"); } } // 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 && 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 && savesignal && 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++) + // 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. NUIS_LOG(REC, "Filled " << fillcount << " signal events."); if (savesignal) { int mem = ( // sizeof(fSignalEventBoxes) + // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) + sizeof(MeasurementVariableBox1D) * fillcount) * 1E-6; - NUIS_LOG(REC, " -> Saved " << fillcount << " signal boxes for faster access. (~" - << mem << " MB)"); + NUIS_LOG(REC, " -> Saved " << fillcount + << " signal boxes for faster access. (~" << mem + << " MB)"); if (fIsAllSplines and !fSignalEventSplines.empty()) { int splmem = sizeof(float) * fSignalEventSplines.size() * fSignalEventSplines[0].size() * 1E-6; - NUIS_LOG(REC, " -> Saved " << fillcount << " " << fSignalEventSplines.size() - << " spline sets into memory. (~" << splmem - << " MB)"); + NUIS_LOG(REC, " -> Saved " + << fillcount << " " << fSignalEventSplines.size() + << " spline sets into memory. (~" << splmem << " MB)"); } } NUIS_LOG(REC, - "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart); + "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart); // Check SignalReconfigures works for all samples if (savesignal) { double likefull = GetLikelihood(); ReconfigureFastUsingManager(); double likefast = GetLikelihood(); if (fabs(likefull - likefast) > 0.0001) { NUIS_ERR(FTL, "Fast and Full Likelihoods DIFFER! : " << likefull << " : " - << likefast); - NUIS_ERR(FTL, "This means some samples you are using are not setup to use " - "SignalReconfigures=1"); + << likefast); + NUIS_ERR(FTL, + "This means some samples you are using are not setup to use " + "SignalReconfigures=1"); NUIS_ERR(FTL, "Please turn OFF signal reconfigures."); throw; } else { NUIS_LOG(FIT, - "Likelihoods for FULL and FAST match. Will use FAST next time."); + "Likelihoods for FULL and FAST match. Will use FAST next time."); } } }; //*************************************************** void JointFCN::ReconfigureFastUsingManager() { //*************************************************** NUIS_LOG(FIT, " -> Doing FAST using manager"); // 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()) { NUIS_ERR(WRN, "Signal Flags Empty! Using normal manager."); 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; // This is just the total number of events // int nevents = fSignalEventFlags.size(); // This is the number of events that are signal int nevents = fSignalEventBoxes.size(); int countwidth = nevents / 10; // If All Splines tell splines they need a reconfigure. std::vector::iterator inp_iter = fInputList.begin(); if (fIsAllSplines) { NUIS_LOG(REC, "All Spline Inputs so using fast spline loop."); 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; // #pragma omp parallel for shared(splinecount,sigcount) for (uint iinput = 0; iinput < fInputList.size(); iinput++) { InputHandlerBase *curinput = fInputList[iinput]; BaseFitEvt *curevent = curinput->FirstBaseEvent(); // Loop over the events in each input for (int i = 0; i < curinput->GetNEvents(); i++) { double rwweight = 0.0; // If the event is a signal event if (fSignalEventFlags[sigcount]) { // Get Event Info if (!fIsAllSplines) { if (fFillNuisanceEvent) { curevent = 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 * curevent->CustomWeight; rwweight = curevent->Weight; coreeventweights[splinecount] = rwweight; if (countwidth && ((splinecount % countwidth) == 0)) { NUIS_LOG(REC, curinput->GetName() - << " : Processed " << i - << " events. W = " << curevent->Weight << std::endl); + << " : Processed " << i << " events. W = " + << curevent->Weight << std::endl); } // #pragma omp atomic splinecount++; } // #pragma omp atomic sigcount++; } } NUIS_LOG(SAM, "Processed event weights."); // #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) { NUIS_LOG(REC, "Filled " << ispline << " sample weights."); } // Iterate over the main signal event containers. samsig_iter++; box_iter++; spline_iter++; splinecount++; } // End of Fast Event Loop =================== NUIS_LOG(SAM, "Filled sample distributions."); // Now loop over all Measurements // Convert Binned events iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase *exp = (*iterSam); exp->ConvertEventRates(); } // Cleanup coreeventweights delete coreeventweights; // Print some reconfigure profiling. NUIS_LOG(REC, "Filled " << fillcount << " signal events."); - NUIS_LOG(REC, - "Time taken ReconfigureFastUsingManager() : " << time(NULL) - timestart); + NUIS_LOG(REC, "Time taken ReconfigureFastUsingManager() : " << time(NULL) - + timestart); } //*************************************************** void JointFCN::Write() { //*************************************************** // Save a likelihood/ndof plot NUIS_LOG(MIN, "Writing likelihood plot..."); 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); } if (likes.size()) { 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 (int 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 NUIS_LOG(MIN, "Writing each of the data classes..."); 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) { //*************************************************** NUIS_LOG(MIN, "Setting fake data from " << fakeinput); 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; } //*************************************************** std::vector JointFCN::GetAllNames() { //*************************************************** // Vect of all likelihoods and total std::vector namevect; // Loop over samples first for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase *exp = *iter; // Get Likelihoods and push to vector namevect.push_back(exp->GetName()); } // Loop over pulls second for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull *pull = *iter; // Push back to vector namevect.push_back(pull->GetName()); } // Finally add the total namevect.push_back("total"); return namevect; } //*************************************************** std::vector JointFCN::GetAllLikelihoods() { //*************************************************** // Vect of all likelihoods and total std::vector likevect; double total_likelihood = 0.0; NUIS_LOG(MIN, "Likelihoods : "); // Loop over samples first for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase *exp = *iter; // Get Likelihoods and push to vector double singlelike = exp->GetLikelihood(); likevect.push_back(singlelike); total_likelihood += singlelike; // Print Out NUIS_LOG(MIN, "-> " << std::left << std::setw(40) << exp->GetName() << " : " - << singlelike); + << singlelike); } // Loop over pulls second for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull *pull = *iter; // Push back to vector double singlelike = pull->GetLikelihood(); likevect.push_back(singlelike); total_likelihood += singlelike; // Print Out - NUIS_LOG(MIN, "-> " << std::left << std::setw(40) << pull->GetName() << " : " - << singlelike); + NUIS_LOG(MIN, "-> " << std::left << std::setw(40) << pull->GetName() + << " : " << singlelike); } // Finally add the total likelihood likevect.push_back(total_likelihood); return likevect; } //*************************************************** std::vector JointFCN::GetAllNDOF() { //*************************************************** // Vect of all ndof and total std::vector ndofvect; int total_ndof = 0; // Loop over samples first for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase *exp = *iter; // Get Likelihoods and push to vector int singlendof = exp->GetNDOF(); ndofvect.push_back(singlendof); total_ndof += singlendof; } // Loop over pulls second for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull *pull = *iter; // Push back to vector int singlendof = pull->GetNDOF(); ndofvect.push_back(singlendof); total_ndof += singlendof; } // Finally add the total ndof ndofvect.push_back(total_ndof); return ndofvect; } diff --git a/src/FCN/JointFCN.h b/src/FCN/JointFCN.h index 1f47f28..895bb9c 100755 --- a/src/FCN/JointFCN.h +++ b/src/FCN/JointFCN.h @@ -1,180 +1,196 @@ #ifndef _JOINT_FCN_H_ #define _JOINT_FCN_H_ /*! * \addtogroup FCN * @{ */ #include #include #include #include // ROOT headers #include "TTree.h" #include "TH1D.h" #include "TH2D.h" #include #include "TGraphErrors.h" #include #include #include "FitUtils.h" // External fitter headers #include "MeasurementBase.h" #include "SampleList.h" #define GetCurrentDir getcwd #include "EventManager.h" #include "Math/Functor.h" #include "ParamPull.h" #include "NuisConfig.h" #include "NuisKey.h" #include "MeasurementVariableBox.h" #include "MeasurementVariableBox1D.h" using namespace FitUtils; using namespace FitBase; //! Main FCN Class which ROOT's joint function needs to evaulate the chi2 at each stage of the fit. class JointFCN { public: //! Constructor //! cardfile = Path to input card file listing samples JointFCN(std::vector samplekeys, TFile* outfile = NULL); JointFCN(TFile* outfile = NULL); // Loads from global config //! Destructor ~JointFCN(); //! Create sample list from cardfile void LoadSamples(std::vector samplekeys); void LoadPulls(std::vector pullkeys); //! Main Likelihood evaluation FCN double DoEval(const double *x); //! Func Wrapper for ROOT inline double operator() (const std::vector & x) { double* x_array = new double[x.size()]; return this->DoEval(x_array); }; //! Func Wrapper for ROOT inline double operator() (const double *x) { return this->DoEval(x); }; //! Create a TTree to save all dial value iterations for this FCN void CreateIterationTree(std::string name, FitWeight* rw); //! Fills dial values and sample likelihoods into tree void FillIterationTree(FitWeight* rw); //! Writes TTree to fOutput directory void WriteIterationTree(); //! Deletes TTree void DestroyIterationTree(); //! Get Degrees of Freedom for samples (NBins) int GetNDOF(); //! Return NDOF wrapper inline unsigned int NDim() {return this->GetNDOF();}; //! Reconfigure samples where we force all events to be looped over. void ReconfigureAllEvents() ; //! Call Reconfigure on samples. //! Choice of reconfigure type depends on whether dials have changed //! and the MC has been filled. void ReconfigureSamples(bool fullconfig = false); //! Call reconfigure on only signal events (defaults to all events if CurIter=0) void ReconfigureSignal(); //! Gets likelihood for all samples in FCN (assuming uncorrelated) double GetLikelihood(); //! Returns list of pointers to the samples inline std::list GetSampleList() { return fSamples; } //! Return list of pointers to all the pulls inline std::list GetPullList() { return fPulls; }; //! Write all samples to output DIR void Write(); //! Set Fake data from file/MC void SetFakeData(std::string fakeinput); //! Reconfigure looping over duplicate inputs void ReconfigureUsingManager(); //! Reconfigure Fast looping over duplicate inputs void ReconfigureFastUsingManager(); /// Throws data according to current stats void ThrowDataToy(); std::vector GetSubSampleList(); std::vector GetInputList(); std::vector GetAllNames(); std::vector GetAllLikelihoods(); std::vector GetAllNDOF(); + + void SetVariableMirrored(int ipar, double mirror_value, bool mirror_above){ + mirror_param mir_par; + mir_par.mirror_value = mirror_value; + mir_par.mirror_above = mirror_above; + fMirroredParams[ipar] = mir_par; + } + void SetNParams(int npar){ + fNPars = npar; + } + private: //! Append the experiments to include in the fit to this list std::list fSamples; //! Append parameter pull terms to include penalties in the fit to this list std::list fPulls; TDirectory *fOutputDir; //!< Directory to save contents std::string fCardFile; //!< Input Card text file bool fDialChanged; //!< Flag for if RW dials changed UInt_t fCurIter; //!< Counter for how many times reconfigure called bool fMCFilled; //!< Check MC has at least been filled once bool fIterationTree; //!< Tree to save RW values on each function call int fNDials; //!< Number of RW Dials in FitWeight double* fDialVals; //!< Current Values of RW Dials double fLikelihood; //!< Current likelihood for joint sample likelihood double fNDOF; //!< Total NDOF double* fSampleLikes; //!< Likelihoods for each individual measurement in list int * fSampleNDOF; //!< NDOF for each individual measurement in list bool fUsingEventManager; //!< Flag for doing joint comparisons std::vector< std::vector > fSignalEventSplines; std::vector< std::vector > fSignalEventBoxes; std::vector< bool > fSignalEventFlags; std::vector< std::vector > fSampleSignalFlags; std::vector fInputList; std::vector fSubSampleList; bool fIsAllSplines; std::vector< int > fIterationCount; std::vector< double > fCurrentValues; std::vector< std::string > fNameValues; std::vector< std::vector > fIterationValues; int fSampleN; std::string fIterationTreeName; - - + struct mirror_param { + double mirror_value; + bool mirror_above; + }; + std::map fMirroredParams; + //the number of pars added to the minimizer, should be the same as fNDials + int fNPars; }; /*! @} */ #endif // _JOINT_FCN_H_ diff --git a/src/Routines/ComparisonRoutines.cxx b/src/Routines/ComparisonRoutines.cxx index 78807a2..a651d24 100755 --- a/src/Routines/ComparisonRoutines.cxx +++ b/src/Routines/ComparisonRoutines.cxx @@ -1,526 +1,637 @@ // 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 "ComparisonRoutines.h" /* Constructor/Destructor */ //************************ void ComparisonRoutines::Init() { //************************ fOutputFile = ""; fOutputRootFile = NULL; fStrategy = "Compare"; fRoutines.clear(); fCardFile = ""; fFakeDataInput = ""; fSampleFCN = NULL; fAllowedRoutines = ("Compare"); }; //************************************* ComparisonRoutines::~ComparisonRoutines() { //************************************* delete fOutputRootFile; }; /* Input Functions */ //************************************* ComparisonRoutines::ComparisonRoutines(int argc, char *argv[]) { //************************************* // Initialise Defaults Init(); nuisconfig configuration = Config::Get(); // Default containers std::string cardfile = ""; std::string maxevents = "-1"; std::string skipevents = "0"; int errorcount = 0; int verbocount = 0; std::vector xmlcmds; std::vector configargs; // Make easier to handle arguments. std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv); ParserUtils::ParseArgument(args, "-c", fCardFile, true); ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false); ParserUtils::ParseArgument(args, "-n", maxevents, false, false); ParserUtils::ParseArgument(args, "-f", fStrategy, false, false); ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false); ParserUtils::ParseArgument(args, "-i", xmlcmds); ParserUtils::ParseArgument(args, "-q", configargs); ParserUtils::ParseArgument(args, "-s", skipevents); ParserUtils::ParseCounter(args, "e", errorcount); ParserUtils::ParseCounter(args, "v", verbocount); ParserUtils::CheckBadArguments(args); // Add extra defaults if none given if (fCardFile.empty() and xmlcmds.empty()) { NUIS_ABORT("No input supplied!"); } if (fOutputFile.empty() and !fCardFile.empty()) { fOutputFile = fCardFile + ".root"; NUIS_ERR(WRN, "No output supplied so saving it to: " << fOutputFile); } else if (fOutputFile.empty()) { NUIS_ABORT("No output file or cardfile supplied!"); } // Configuration Setup ============================= // Check no comp key is available nuiskey fCompKey; if (Config::Get().GetNodes("nuiscomp").empty()) { fCompKey = Config::Get().CreateNode("nuiscomp"); } else { fCompKey = Config::Get().GetNodes("nuiscomp")[0]; } if (!fCardFile.empty()) fCompKey.Set("cardfile", fCardFile); if (!fOutputFile.empty()) fCompKey.Set("outputfile", fOutputFile); if (!fStrategy.empty()) fCompKey.Set("strategy", fStrategy); // Load XML Cardfile configuration.LoadSettings(fCompKey.GetS("cardfile"), ""); // Add Config Args for (size_t i = 0; i < configargs.size(); i++) { configuration.OverrideConfig(configargs[i]); } if (maxevents.compare("-1")) { configuration.OverrideConfig("MAXEVENTS=" + maxevents); } configuration.OverrideConfig("NSKIPEVENTS=" + skipevents); // Finish configuration XML configuration.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml"); // Add Error Verbo Lines verbocount += Config::GetParI("VERBOSITY"); errorcount += Config::GetParI("ERROR"); bool trace = Config::GetParB("TRACE"); std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl; std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl; SETVERBOSITY(verbocount); SETTRACE(trace); // Comparison Setup ======================================== // Proper Setup fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); SetupComparisonsFromXML(); SetupRWEngine(); SetupFCN(); return; }; //************************************* void ComparisonRoutines::SetupComparisonsFromXML() { //************************************* NUIS_LOG(FIT, "Setting up nuiscomp"); // Setup Parameters ------------------------------------------ std::vector parkeys = Config::QueryKeys("parameter"); if (!parkeys.empty()) { NUIS_LOG(FIT, "Number of parameters : " << parkeys.size()); } for (size_t i = 0; i < parkeys.size(); i++) { nuiskey key = parkeys.at(i); // Check for type,name,nom if (!key.Has("type")) { NUIS_ERR(FTL, "No type given for parameter " << i); NUIS_ABORT("type='PARAMETER_TYPE'"); } else if (!key.Has("name")) { NUIS_ERR(FTL, "No name given for parameter " << i); NUIS_ABORT("name='SAMPLE_NAME'"); } else if (!key.Has("nominal")) { NUIS_ERR(FTL, "No nominal given for parameter " << i); NUIS_ABORT("nominal='NOMINAL_VALUE'"); } // Get Inputs std::string partype = key.GetS("type"); std::string parname = key.GetS("name"); double parnom = key.GetD("nominal"); double parlow = parnom - 1; double parhigh = parnom + 1; double parstep = 1; // override if state not given if (!key.Has("state")) { key.SetS("state", "FIX"); } std::string parstate = key.GetS("state"); // Check for incomplete limtis int limdef = ((int)key.Has("low") + (int)key.Has("high") + (int)key.Has("step")); if (limdef > 0 and limdef < 3) { NUIS_ERR(FTL, "Incomplete limit set given for parameter : " << parname); NUIS_ABORT( "Requires: low='LOWER_LIMIT' high='UPPER_LIMIT' step='STEP_SIZE' "); } // Extra limits if (key.Has("low")) { parlow = key.GetD("low"); parhigh = key.GetD("high"); parstep = key.GetD("step"); NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom << " : " << parlow << " < p < " << parhigh << " : " << parstate); } else { NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom << " : " << parstate); } + bool ismirr = false; + if (key.Has("mirror_point")) { + ismirr = true; + mirror_param mir; + mir.mirror_value = key.GetD("mirror_point"); + mir.mirror_above = key.GetB("mirror_above"); + fMirroredParams[parname] = mir; + NUIS_LOG(FIT, + "\t\t" << parname << " is mirrored at " << mir.mirror_value + << " " + << (mir.mirror_above ? "from above" : "from below")); + } + // Convert if required if (parstate.find("ABS") != std::string::npos) { + if (ismirr) { + NUIS_ABORT("Cannot mirror parameters with ABS state!"); + } parnom = FitBase::RWAbsToSigma(partype, parname, parnom); parlow = FitBase::RWAbsToSigma(partype, parname, parlow); parhigh = FitBase::RWAbsToSigma(partype, parname, parhigh); parstep = FitBase::RWAbsToSigma(partype, parname, parstep); } else if (parstate.find("FRAC") != std::string::npos) { + if (ismirr) { + NUIS_ABORT("Cannot mirror parameters with FRAC state!"); + } parnom = FitBase::RWFracToSigma(partype, parname, parnom); parlow = FitBase::RWFracToSigma(partype, parname, parlow); parhigh = FitBase::RWFracToSigma(partype, parname, parhigh); parstep = FitBase::RWFracToSigma(partype, parname, parstep); } // Push into vectors fParams.push_back(parname); fTypeVals[parname] = FitBase::ConvDialType(partype); fCurVals[parname] = parnom; fStateVals[parname] = parstate; } // Setup Samples ---------------------------------------------- std::vector samplekeys = Config::QueryKeys("sample"); if (!samplekeys.empty()) { NUIS_LOG(FIT, "Number of samples : " << samplekeys.size()); } for (size_t i = 0; i < samplekeys.size(); i++) { nuiskey key = samplekeys.at(i); // Get Sample Options std::string samplename = key.GetS("name"); std::string samplefile = key.GetS("input"); std::string sampletype = key.Has("type") ? key.GetS("type") : "DEFAULT"; double samplenorm = key.Has("norm") ? key.GetD("norm") : 1.0; // Print out NUIS_LOG(FIT, "Read Sample " << i << ". : " << samplename << " (" << sampletype << ") [Norm=" << samplenorm << "]" << std::endl << " -> input='" << samplefile << "'"); // If FREE add to parameters otherwise continue if (sampletype.find("FREE") == std::string::npos) { if (samplenorm != 1.0) { NUIS_ERR(FTL, "You provided a sample normalisation but did not specify " "that the sample is free"); NUIS_ABORT("Change so sample contains type=\"FREE\" and re-run"); } continue; } // Form norm dial from samplename + sampletype + "_norm"; std::string normname = samplename + "_norm"; // Check normname not already present if (fTypeVals.find("normname") != fTypeVals.end()) { continue; } // Add new norm dial to list if its passed above checks fParams.push_back(normname); fTypeVals[normname] = kNORM; fStateVals[normname] = sampletype; fCurVals[normname] = samplenorm; } // Setup Fake Parameters ----------------------------- std::vector fakekeys = Config::QueryKeys("fakeparameter"); if (!fakekeys.empty()) { NUIS_LOG(FIT, "Number of fake parameters : " << fakekeys.size()); } for (size_t i = 0; i < fakekeys.size(); i++) { nuiskey key = fakekeys.at(i); // Check for type,name,nom if (!key.Has("name")) { NUIS_ABORT("No name given for fakeparameter " << i); } else if (!key.Has("nominal")) { NUIS_ABORT("No nominal given for fakeparameter " << i); } // Get Inputs std::string parname = key.GetS("name"); double parnom = key.GetD("nominal"); // Push into vectors fFakeVals[parname] = parnom; } } //************************************* void ComparisonRoutines::SetupRWEngine() { //************************************* NUIS_LOG(FIT, "Setting up FitWeight Engine"); for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams[i]; FitBase::GetRW()->IncludeDial(name, fTypeVals.at(name)); } return; } //************************************* void ComparisonRoutines::SetupFCN() { //************************************* NUIS_LOG(FIT, "Building the SampleFCN"); if (fSampleFCN) delete fSampleFCN; Config::Get().out = fOutputRootFile; fOutputRootFile->cd(); fSampleFCN = new JointFCN(fOutputRootFile); SetFakeData(); return; } //************************************* void ComparisonRoutines::SetFakeData() { //************************************* if (fFakeDataInput.empty()) return; if (fFakeDataInput.compare("MC") == 0) { NUIS_LOG(FIT, "Setting fake data from MC starting prediction."); UpdateRWEngine(fFakeVals); FitBase::GetRW()->Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->SetFakeData("MC"); - UpdateRWEngine(fCurVals); + std::map CurVals_wmirr; + for (size_t i = 0; i < fParams.size(); ++i) { + std::string const &pname = fParams[i]; + if (fMirroredParams.count(pname)) { + if (!fMirroredParams[pname].mirror_above && + (fCurVals[pname] < fMirroredParams[pname].mirror_value)) { + double xabove = fMirroredParams[pname].mirror_value - fCurVals[pname]; + CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value + xabove; + std::cout << "\t--Parameter " << pname << " mirrored from " + << fCurVals[pname] << " -> " << CurVals_wmirr[pname] + << std::endl; + } else if (fMirroredParams[pname].mirror_above && + (fCurVals[pname] >= fMirroredParams[pname].mirror_value)) { + double xabove = fCurVals[pname] - fMirroredParams[pname].mirror_value; + CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value - xabove; + std::cout << "\t--Parameter " << pname << " mirrored from " + << fCurVals[pname] << " -> " << CurVals_wmirr[pname] + << std::endl; + } else { + CurVals_wmirr[pname] = fCurVals[pname]; + } + } else { + CurVals_wmirr[pname] = fCurVals[pname]; + } + } + + UpdateRWEngine(CurVals_wmirr); NUIS_LOG(FIT, "Set all data to fake MC predictions."); } else { NUIS_LOG(FIT, "Setting fake data from: " << fFakeDataInput); fSampleFCN->SetFakeData(fFakeDataInput); } return; } /* Fitting Functions */ //************************************* void ComparisonRoutines::UpdateRWEngine( std::map &updateVals) { //************************************* for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams[i]; if (updateVals.find(name) == updateVals.end()) continue; FitBase::GetRW()->SetDialValue(name, updateVals.at(name)); } FitBase::GetRW()->Reconfigure(); return; } //************************************* void ComparisonRoutines::Run() { //************************************* NUIS_LOG(FIT, "Running ComparisonRoutines : " << fStrategy); if (FitPar::Config().GetParB("save_nominal")) { SaveNominal(); } // Parse given routines fRoutines = GeneralUtils::ParseToStr(fStrategy, ","); if (fRoutines.empty()) { NUIS_ABORT("Trying to run ComparisonRoutines with no routines given!"); } for (UInt_t i = 0; i < fRoutines.size(); i++) { std::string routine = fRoutines.at(i); NUIS_LOG(FIT, "Routine: " << routine); if (!routine.compare("Compare")) { - UpdateRWEngine(fCurVals); + + std::map CurVals_wmirr; + for (size_t i = 0; i < fParams.size(); ++i) { + std::string const &pname = fParams[i]; + if (fMirroredParams.count(pname)) { + std::cout << "MA? " << fMirroredParams[pname].mirror_above + << ", OVal: " << fCurVals[pname] + << ", MPoint: " << fMirroredParams[pname].mirror_value + << std::endl; + if (!fMirroredParams[pname].mirror_above && + (fCurVals[pname] < fMirroredParams[pname].mirror_value)) { + double xabove = + fMirroredParams[pname].mirror_value - fCurVals[pname]; + CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value + xabove; + std::cout << "\t--Parameter " << pname << " mirrored from " + << fCurVals[pname] << " -> " << CurVals_wmirr[pname] + << std::endl; + } else if (fMirroredParams[pname].mirror_above && + (fCurVals[pname] >= fMirroredParams[pname].mirror_value)) { + double xabove = + fCurVals[pname] - fMirroredParams[pname].mirror_value; + CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value - xabove; + std::cout << "\t--Parameter " << pname << " mirrored from " + << fCurVals[pname] << " -> " << CurVals_wmirr[pname] + << std::endl; + } else { + CurVals_wmirr[pname] = fCurVals[pname]; + } + } else { + CurVals_wmirr[pname] = fCurVals[pname]; + } + std::cout << "~~~~~~~" << pname << " : " << fCurVals[pname] << " -> " + << CurVals_wmirr[pname] << std::endl; + } + + UpdateRWEngine(CurVals_wmirr); GenerateComparison(); PrintState(); SaveCurrentState(); } } return; } //************************************* void ComparisonRoutines::GenerateComparison() { //************************************* NUIS_LOG(FIT, "Generating Comparison."); // Main Event Loop from event Manager fSampleFCN->ReconfigureAllEvents(); return; } //************************************* void ComparisonRoutines::PrintState() { //************************************* NUIS_LOG(FIT, "------------"); // Count max size int maxcount = 0; for (UInt_t i = 0; i < fParams.size(); i++) { maxcount = max(int(fParams[i].size()), maxcount); } // Header NUIS_LOG(FIT, " # " << left << setw(maxcount) << "Parameter " << " = " << setw(10) << "Value" << " +- " << setw(10) << "Error" << " " << setw(8) << "(Units)" << " " << setw(10) << "Conv. Val" << " +- " << setw(10) << "Conv. Err" << " " << setw(8) << "(Units)"); + std::map CurVals_wmirr; + for (size_t i = 0; i < fParams.size(); ++i) { + std::string const &pname = fParams[i]; + if (fMirroredParams.count(pname)) { + if (!fMirroredParams[pname].mirror_above && + (fCurVals[pname] < fMirroredParams[pname].mirror_value)) { + double xabove = fMirroredParams[pname].mirror_value - fCurVals[pname]; + CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value + xabove; + std::cout << "\t--Parameter " << pname << " mirrored from " + << fCurVals[pname] << " -> " << CurVals_wmirr[pname] + << std::endl; + } else if (fMirroredParams[pname].mirror_above && + (fCurVals[pname] >= fMirroredParams[pname].mirror_value)) { + double xabove = fCurVals[pname] - fMirroredParams[pname].mirror_value; + CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value - xabove; + std::cout << "\t--Parameter " << pname << " mirrored from " + << fCurVals[pname] << " -> " << CurVals_wmirr[pname] + << std::endl; + } else { + CurVals_wmirr[pname] = fCurVals[pname]; + } + } else { + CurVals_wmirr[pname] = fCurVals[pname]; + } + } + // Parameters for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); std::string typestr = FitBase::ConvDialType(fTypeVals[syst]); std::string curunits = "(sig.)"; double curval = fCurVals[syst]; double curerr = 0.0; if (fStateVals[syst].find("ABS") != std::string::npos) { curval = FitBase::RWSigmaToAbs(typestr, syst, curval); curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); curunits = "(Abs.)"; } else if (fStateVals[syst].find("FRAC") != std::string::npos) { curval = FitBase::RWSigmaToFrac(typestr, syst, curval); curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) - FitBase::RWSigmaToFrac(typestr, syst, 0.0)); curunits = "(Frac)"; } std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")"; double convval = FitBase::RWSigmaToAbs(typestr, syst, curval); double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); std::ostringstream curparstring; curparstring << " " << setw(3) << left << i << ". " << setw(maxcount) << syst << " = " << setw(10) << curval << " +- " << setw(10) << curerr << " " << setw(8) << curunits << " " << setw(10) << convval << " +- " << setw(10) << converr << " " << setw(8) << convunits; - NUIS_LOG(FIT, curparstring.str()); + + if (fMirroredParams.count(syst)) { + NUIS_LOG(FIT, "\t\t--> Mirrored at " << fMirroredParams[syst].mirror_value + << " to effective value " + << CurVals_wmirr[syst]); + } } NUIS_LOG(FIT, "------------"); double like = fSampleFCN->GetLikelihood(); NUIS_LOG(FIT, std::left << std::setw(46) << "Likelihood for JointFCN: " << like); NUIS_LOG(FIT, "------------"); } /* Write Functions */ //************************************* void ComparisonRoutines::SaveCurrentState(std::string subdir) { //************************************* NUIS_LOG(FIT, "Saving current full FCN predictions"); // Setup DIRS TDirectory *curdir = gDirectory; if (!subdir.empty()) { TDirectory *newdir = (TDirectory *)gDirectory->mkdir(subdir.c_str()); newdir->cd(); } fSampleFCN->Write(); // Change back to current DIR curdir->cd(); return; } //************************************* void ComparisonRoutines::SaveNominal() { //************************************* fOutputRootFile->cd(); NUIS_LOG(FIT, "Saving Nominal Predictions (be cautious with this)"); FitBase::GetRW()->Reconfigure(); GenerateComparison(); SaveCurrentState("nominal"); }; diff --git a/src/Routines/ComparisonRoutines.h b/src/Routines/ComparisonRoutines.h index 1d967ac..9111ef8 100755 --- a/src/Routines/ComparisonRoutines.h +++ b/src/Routines/ComparisonRoutines.h @@ -1,168 +1,175 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef COMPARISON_ROUTINES_H #define COMPARISON_ROUTINES_H /*! \addtogroup Routines @{ */ #include "TH1.h" #include "TF1.h" #include "TMatrixD.h" #include "TVectorD.h" #include "TSystem.h" #include "TFile.h" #include "TProfile.h" #include #include #include #include #include #include "FitEvent.h" #include "JointFCN.h" #include "GeneralUtils.h" #include "NuisConfig.h" #include "NuisKey.h" #include "FitLogger.h" #include "ParserUtils.h" enum minstate { kErrorStatus = -1, kGoodStatus, kFitError, kNoChange, kFitFinished, kFitUnfinished, kStateChange, }; //************************************* /// Collects all possible fit routines into a single class to avoid repeated code class ComparisonRoutines { //************************************* public: /* Constructor/Destructor */ /// Constructor reads in arguments given at the command line for the fit here. ComparisonRoutines(int argc, char* argv[]); /// Default destructor ~ComparisonRoutines(); /// Reset everything to default/NULL void Init(); /* Input Functions */ /// Queries configuration keys to setup Parameters/Samples/FakeParameters void SetupComparisonsFromXML(); /* Setup Functions */ /// Setups up our custom RW engine with all the parameters passed in the card file void SetupRWEngine(); /// Setups up the jointFCN. void SetupFCN(); /// Set the current data histograms in each sample to the fake data. void SetFakeData(); /* Fitting Functions */ /// Main function to actually start iterating over the different required fit routines void Run(); /// Creates a comparison from FCN void GenerateComparison(); /// Given a new map change the values that the RW engine is currently set to void UpdateRWEngine(std::map& updateVals); /// Print current value void PrintState(); /* Write Functions */ /// Save the sample plots for current MC /// dir if not empty forces plots to be saved in a subdirectory of outputfile void SaveCurrentState(std::string subdir=""); /// Save starting predictions into a seperate folder void SaveNominal(); /* MISC Functions */ /// Get previous fit status from a file Int_t GetStatus(); protected: //! Our Custom ReWeight Object FitWeight* rw; std::string fOutputFile; ///< Output file name // std::string fInputFile; ///< Input file name // TFile* fInputRootFile; ///< TFile* fOutputRootFile; ///< Output ROOT TFile JointFCN* fSampleFCN; ///< Joint Samples Container that handles reconfigures. std::string fCardFile; ///< Input card/XML file. std::string fStrategy; ///< Comparison routine selection. std::vector fRoutines; ///< Split vector of comparison routine selection. std::string fAllowedRoutines; ///< Hard coded list of allowed routines. /// Fake data flag. Can be 'MC' to use 'fake_parameter' /// or 'path_to_file.root' to use previous NUISANCE MC predictions. std::string fFakeDataInput; // Input Dial Vals std::vector fParams; ///< Vector of dial names. std::map fStateVals; ///< Map of dial states std::map fCurVals; ///< Map of dial values std::map fTypeVals; ///< Map of dial type enums. + + struct mirror_param { + double mirror_value; + bool mirror_above; + }; + std::map fMirroredParams; + // Fake Dial Vals std::map fFakeVals; ///< Map of fake data settings. // Configuration nuiskey fCompKey; ///< Configuration Key for this Comparison Instance }; /*! @} */ #endif diff --git a/src/Routines/MinimizerRoutines.cxx b/src/Routines/MinimizerRoutines.cxx index 4e20a4a..9293e51 100755 --- a/src/Routines/MinimizerRoutines.cxx +++ b/src/Routines/MinimizerRoutines.cxx @@ -1,1535 +1,1561 @@ // 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 "MinimizerRoutines.h" #include "Simple_MH_Sampler.h" /* Constructor/Destructor */ //************************ void MinimizerRoutines::Init() { //************************ fInputFile = ""; fInputRootFile = NULL; fOutputFile = ""; fOutputRootFile = NULL; fCovar = NULL; fCovFree = NULL; fCorrel = NULL; fCorFree = NULL; fDecomp = NULL; fDecFree = NULL; fStrategy = "Migrad,FixAtLimBreak,Migrad"; fRoutines.clear(); fCardFile = ""; fFakeDataInput = ""; fSampleFCN = NULL; fMinimizer = NULL; fMinimizerFCN = NULL; fCallFunctor = NULL; fAllowedRoutines = ("Migrad,Simplex,Combined," "Brute,Fumili,ConjugateFR," "ConjugatePR,BFGS,BFGS2," "SteepDesc,GSLSimAn,FixAtLim,FixAtLimBreak," "Chi2Scan1D,Chi2Scan2D,Contours,ErrorBands," "DataToys,MCMC"); }; //************************************* MinimizerRoutines::~MinimizerRoutines(){ //************************************* }; /* Input Functions */ //************************************* MinimizerRoutines::MinimizerRoutines() { //************************************* Init(); } //************************************* MinimizerRoutines::MinimizerRoutines(int argc, char *argv[]) { //************************************* // Initialise Defaults Init(); nuisconfig configuration = Config::Get(); // Default containers std::string cardfile = ""; std::string maxevents = "-1"; int errorcount = 0; int verbocount = 0; std::vector xmlcmds; std::vector configargs; // Make easier to handle arguments. std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv); ParserUtils::ParseArgument(args, "-c", fCardFile, true); ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false); ParserUtils::ParseArgument(args, "-n", maxevents, false, false); ParserUtils::ParseArgument(args, "-f", fStrategy, false, false); ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false); ParserUtils::ParseArgument(args, "-i", xmlcmds); ParserUtils::ParseArgument(args, "-q", configargs); ParserUtils::ParseCounter(args, "e", errorcount); ParserUtils::ParseCounter(args, "v", verbocount); ParserUtils::CheckBadArguments(args); // Add extra defaults if none given if (fCardFile.empty() and xmlcmds.empty()) { NUIS_ABORT("No input supplied!"); } if (fOutputFile.empty() and !fCardFile.empty()) { fOutputFile = fCardFile + ".root"; NUIS_ERR(WRN, "No output supplied so saving it to: " << fOutputFile); } else if (fOutputFile.empty()) { NUIS_ABORT("No output file or cardfile supplied!"); } // Configuration Setup ============================= // Check no comp key is available nuiskey fCompKey; if (Config::Get().GetNodes("nuiscomp").empty()) { fCompKey = Config::Get().CreateNode("nuiscomp"); } else { fCompKey = Config::Get().GetNodes("nuiscomp")[0]; } if (!fCardFile.empty()) fCompKey.Set("cardfile", fCardFile); if (!fOutputFile.empty()) fCompKey.Set("outputfile", fOutputFile); if (!fStrategy.empty()) fCompKey.Set("strategy", fStrategy); // Load XML Cardfile configuration.LoadSettings(fCompKey.GetS("cardfile"), ""); // Add Config Args for (size_t i = 0; i < configargs.size(); i++) { configuration.OverrideConfig(configargs[i]); } if (maxevents.compare("-1")) { configuration.OverrideConfig("MAXEVENTS=" + maxevents); } // Finish configuration XML configuration.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml"); // Add Error Verbo Lines verbocount += Config::GetParI("VERBOSITY"); errorcount += Config::GetParI("ERROR"); NUIS_LOG(FIT, "[ NUISANCE ]: Setting VERBOSITY=" << verbocount); NUIS_LOG(FIT, "[ NUISANCE ]: Setting ERROR=" << errorcount); // FitPar::log_verb = verbocount; SETVERBOSITY(verbocount); // ERR_VERB(errorcount); // Minimizer Setup ======================================== fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); SetupMinimizerFromXML(); SetupCovariance(); SetupRWEngine(); SetupFCN(); return; }; //************************************* void MinimizerRoutines::SetupMinimizerFromXML() { //************************************* NUIS_LOG(FIT, "Setting up nuismin"); // Setup Parameters ------------------------------------------ std::vector parkeys = Config::QueryKeys("parameter"); if (!parkeys.empty()) { NUIS_LOG(FIT, "Number of parameters : " << parkeys.size()); } for (size_t i = 0; i < parkeys.size(); i++) { nuiskey key = parkeys.at(i); // Check for type,name,nom if (!key.Has("type")) { NUIS_ABORT("No type given for parameter " << i); } else if (!key.Has("name")) { NUIS_ABORT("No name given for parameter " << i); } else if (!key.Has("nominal")) { NUIS_ABORT("No nominal given for parameter " << i); } // Get Inputs std::string partype = key.GetS("type"); std::string parname = key.GetS("name"); double parnom = key.GetD("nominal"); double parlow = parnom - 1; double parhigh = parnom + 1; double parstep = 1; // Override state if none given if (!key.Has("state")) { key.SetS("state", "FIX"); } std::string parstate = key.GetS("state"); // Extra limits if (key.Has("low")) { parlow = key.GetD("low"); parhigh = key.GetD("high"); parstep = key.GetD("step"); NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom - << " : " << parlow << " < p < " << parhigh << " : " - << parstate); + << " : " << parlow << " < p < " << parhigh << " : " + << parstate); } else { NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom - << " : " << parstate); + << " : " << parstate); + } + + bool ismirr = false; + if (key.Has("mirror_point")) { + ismirr=true; + mirror_param mir; + mir.mirror_value = key.GetD("mirror_point"); + mir.mirror_above = key.GetB("mirror_above"); + fMirroredParams[parname] = mir; + NUIS_LOG(FIT, + "\t\t" << parname << " is mirrored at " << mir.mirror_value + << " " + << (mir.mirror_above ? "from above" : "from below")); } // Run Parameter Conversion if needed if (parstate.find("ABS") != std::string::npos) { + if(ismirr){ + NUIS_ABORT("Cannot mirror parameters with ABS state!"); + } double opnom = parnom; double oparstep = parstep; parnom = FitBase::RWAbsToSigma(partype, parname, parnom); parlow = FitBase::RWAbsToSigma(partype, parname, parlow); parhigh = FitBase::RWAbsToSigma(partype, parname, parhigh); parstep = FitBase::RWAbsToSigma(partype, parname, opnom + parstep) - parnom; NUIS_LOG(FIT, "ParStep: " << parstep << " (" << oparstep << ")."); } else if (parstate.find("FRAC") != std::string::npos) { + if(ismirr){ + NUIS_ABORT("Cannot mirror parameters with FRAC state!"); + } parnom = FitBase::RWFracToSigma(partype, parname, parnom); parlow = FitBase::RWFracToSigma(partype, parname, parlow); parhigh = FitBase::RWFracToSigma(partype, parname, parhigh); parstep = FitBase::RWFracToSigma(partype, parname, parstep); } // Push into vectors fParams.push_back(parname); fTypeVals[parname] = FitBase::ConvDialType(partype); fStartVals[parname] = parnom; fCurVals[parname] = parnom; fErrorVals[parname] = 0.0; fStateVals[parname] = parstate; bool fixstate = parstate.find("FIX") != std::string::npos; fFixVals[parname] = fixstate; fStartFixVals[parname] = fFixVals[parname]; fMinVals[parname] = parlow; fMaxVals[parname] = parhigh; fStepVals[parname] = parstep; } // Setup Samples ---------------------------------------------- std::vector samplekeys = Config::QueryKeys("sample"); if (!samplekeys.empty()) { NUIS_LOG(FIT, "Number of samples : " << samplekeys.size()); } for (size_t i = 0; i < samplekeys.size(); i++) { nuiskey key = samplekeys.at(i); // Get Sample Options std::string samplename = key.GetS("name"); std::string samplefile = key.GetS("input"); std::string sampletype = key.Has("type") ? key.GetS("type") : "DEFAULT"; double samplenorm = key.Has("norm") ? key.GetD("norm") : 1.0; // Print out - NUIS_LOG(FIT, "Read sample info " << i << " : " << samplename << std::endl - << "\t\t input -> " << samplefile << std::endl - << "\t\t state -> " << sampletype << std::endl - << "\t\t norm -> " << samplenorm); + NUIS_LOG(FIT, "Read sample info " + << i << " : " << samplename << std::endl + << "\t\t input -> " << samplefile << std::endl + << "\t\t state -> " << sampletype << std::endl + << "\t\t norm -> " << samplenorm); // If FREE add to parameters otherwise continue if (sampletype.find("FREE") == std::string::npos) { if (samplenorm != 1.0) { NUIS_ERR(FTL, "You provided a sample normalisation but did not specify " - "that the sample is free"); + "that the sample is free"); NUIS_ABORT("Change so sample contains type=\"FREE\" and re-run"); } continue; } // Form norm dial from samplename + sampletype + "_norm"; std::string normname = samplename + "_norm"; // Check normname not already present if (fTypeVals.find(normname) != fTypeVals.end()) { continue; } // Add new norm dial to list if its passed above checks fParams.push_back(normname); fTypeVals[normname] = kNORM; fStateVals[normname] = sampletype; fStartVals[normname] = samplenorm; fCurVals[normname] = samplenorm; fErrorVals[normname] = 0.0; fMinVals[normname] = 0.1; fMaxVals[normname] = 10.0; fStepVals[normname] = 0.5; bool state = sampletype.find("FREE") == std::string::npos; fFixVals[normname] = state; fStartFixVals[normname] = state; } // Setup Fake Parameters ----------------------------- std::vector fakekeys = Config::QueryKeys("fakeparameter"); if (!fakekeys.empty()) { NUIS_LOG(FIT, "Number of fake parameters : " << fakekeys.size()); } for (size_t i = 0; i < fakekeys.size(); i++) { nuiskey key = fakekeys.at(i); // Check for type,name,nom if (!key.Has("name")) { NUIS_ABORT("No name given for fakeparameter " << i); } else if (!key.Has("nom")) { NUIS_ABORT("No nominal given for fakeparameter " << i); } // Get Inputs std::string parname = key.GetS("name"); double parnom = key.GetD("nom"); // Push into vectors fFakeVals[parname] = parnom; } } /* Setup Functions */ //************************************* void MinimizerRoutines::SetupRWEngine() { //************************************* for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams[i]; FitBase::GetRW()->IncludeDial(name, fTypeVals.at(name)); } UpdateRWEngine(fStartVals); return; } //************************************* void MinimizerRoutines::SetupFCN() { //************************************* NUIS_LOG(FIT, "Making the jointFCN"); if (fSampleFCN) delete fSampleFCN; // fSampleFCN = new JointFCN(fCardFile, fOutputRootFile); fSampleFCN = new JointFCN(fOutputRootFile); SetFakeData(); fMinimizerFCN = new MinimizerFCN(fSampleFCN); fCallFunctor = new ROOT::Math::Functor(*fMinimizerFCN, fParams.size()); fSampleFCN->CreateIterationTree("fit_iterations", FitBase::GetRW()); return; } //****************************************** void MinimizerRoutines::SetupFitter(std::string routine) { //****************************************** // Make the fitter std::string fitclass = ""; std::string fittype = ""; bool UseMCMC = false; // Get correct types if (!routine.compare("Migrad")) { fitclass = "Minuit2"; fittype = "Migrad"; } else if (!routine.compare("Simplex")) { fitclass = "Minuit2"; fittype = "Simplex"; } else if (!routine.compare("Combined")) { fitclass = "Minuit2"; fittype = "Combined"; } else if (!routine.compare("Brute")) { fitclass = "Minuit2"; fittype = "Scan"; } else if (!routine.compare("Fumili")) { fitclass = "Minuit2"; fittype = "Fumili"; } else if (!routine.compare("ConjugateFR")) { fitclass = "GSLMultiMin"; fittype = "ConjugateFR"; } else if (!routine.compare("ConjugatePR")) { fitclass = "GSLMultiMin"; fittype = "ConjugatePR"; } else if (!routine.compare("BFGS")) { fitclass = "GSLMultiMin"; fittype = "BFGS"; } else if (!routine.compare("BFGS2")) { fitclass = "GSLMultiMin"; fittype = "BFGS2"; } else if (!routine.compare("SteepDesc")) { fitclass = "GSLMultiMin"; fittype = "SteepestDescent"; // } else if (!routine.compare("GSLMulti")) { fitclass = "GSLMultiFit"; // fittype = ""; // Doesn't work out of the box } else if (!routine.compare("GSLSimAn")) { fitclass = "GSLSimAn"; fittype = ""; } else if (!routine.compare("MCMC")) { UseMCMC = true; } // make minimizer if (fMinimizer) delete fMinimizer; if (UseMCMC) { fMinimizer = new Simple_MH_Sampler(); } else { fMinimizer = ROOT::Math::Factory::CreateMinimizer(fitclass, fittype); } fMinimizer->SetMaxFunctionCalls(FitPar::Config().GetParI("MAXCALLS")); if (!routine.compare("Brute")) { fMinimizer->SetMaxFunctionCalls(fParams.size() * fParams.size() * 4); fMinimizer->SetMaxIterations(fParams.size() * fParams.size() * 4); } fMinimizer->SetMaxIterations(FitPar::Config().GetParI("MAXITERATIONS")); fMinimizer->SetTolerance(FitPar::Config().GetParD("TOLERANCE")); fMinimizer->SetStrategy(FitPar::Config().GetParI("STRATEGY")); fMinimizer->SetFunction(*fCallFunctor); int ipar = 0; // Add Fit Parameters for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); bool fixed = true; double vstart, vstep, vlow, vhigh; vstart = vstep = vlow = vhigh = 0.0; if (fCurVals.find(syst) != fCurVals.end()) vstart = fCurVals.at(syst); if (fMinVals.find(syst) != fMinVals.end()) vlow = fMinVals.at(syst); if (fMaxVals.find(syst) != fMaxVals.end()) vhigh = fMaxVals.at(syst); if (fStepVals.find(syst) != fStepVals.end()) vstep = fStepVals.at(syst); if (fFixVals.find(syst) != fFixVals.end()) fixed = fFixVals.at(syst); // fix for errors if (vhigh == vlow) vhigh += 1.0; fMinimizer->SetVariable(ipar, syst, vstart, vstep); fMinimizer->SetVariableLimits(ipar, vlow, vhigh); + if (fMirroredParams.count(syst)) { + fSampleFCN->SetVariableMirrored(ipar, fMirroredParams[syst].mirror_value, + fMirroredParams[syst].mirror_above); + } if (fixed) { fMinimizer->FixVariable(ipar); NUIS_LOG(FIT, "Fixed Param: " << syst); } else { NUIS_LOG(FIT, "Free Param: " << syst << " Start:" << vstart - << " Range:" << vlow << " to " << vhigh - << " Step:" << vstep); + << " Range:" << vlow << " to " << vhigh + << " Step:" << vstep); } ipar++; } + fSampleFCN->SetNParams(ipar); NUIS_LOG(FIT, "Setup Minimizer: " << fMinimizer->NDim() << "(NDim) " - << fMinimizer->NFree() << "(NFree)"); + << fMinimizer->NFree() << "(NFree)"); return; } //************************************* // Set fake data from user input void MinimizerRoutines::SetFakeData() { //************************************* // If the fake data input field (-d) isn't provided, return to caller if (fFakeDataInput.empty()) return; // If user specifies -d MC we set the data to the MC // User can also specify fake data parameters to reweight by doing // "fake_parameter" in input card file // "fake_parameter" gets read in ReadCard function (reads to fFakeVals) if (fFakeDataInput.compare("MC") == 0) { NUIS_LOG(FIT, "Setting fake data from MC starting prediction."); // fFakeVals get read in in ReadCard UpdateRWEngine(fFakeVals); // Reconfigure the reweight engine FitBase::GetRW()->Reconfigure(); // Reconfigure all the samples to the new reweight fSampleFCN->ReconfigureAllEvents(); // Feed on and set the fake-data in each measurement class fSampleFCN->SetFakeData("MC"); // Changed the reweight engine values back to the current values // So we start the fit at a different value than what we set the fake-data // to UpdateRWEngine(fCurVals); NUIS_LOG(FIT, "Set all data to fake MC predictions."); } else { fSampleFCN->SetFakeData(fFakeDataInput); } return; } /* Fitting Functions */ //************************************* void MinimizerRoutines::UpdateRWEngine( std::map &updateVals) { //************************************* for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams[i]; if (updateVals.find(name) == updateVals.end()) continue; FitBase::GetRW()->SetDialValue(name, updateVals.at(name)); } FitBase::GetRW()->Reconfigure(); return; } //************************************* void MinimizerRoutines::Run() { //************************************* NUIS_LOG(FIT, "Running MinimizerRoutines : " << fStrategy); if (FitPar::Config().GetParB("save_nominal")) { SaveNominal(); } // Parse given routines fRoutines = GeneralUtils::ParseToStr(fStrategy, ","); if (fRoutines.empty()) { NUIS_ABORT("Trying to run MinimizerRoutines with no routines given!"); } for (UInt_t i = 0; i < fRoutines.size(); i++) { std::string routine = fRoutines.at(i); int fitstate = kFitUnfinished; NUIS_LOG(FIT, "Running Routine: " << routine); // Try Routines if (routine.find("LowStat") != std::string::npos) LowStatRoutine(routine); else if (routine == "FixAtLim") FixAtLimit(); else if (routine == "FixAtLimBreak") fitstate = FixAtLimit(); else if (routine.find("ErrorBands") != std::string::npos) GenerateErrorBands(); else if (routine.find("DataToys") != std::string::npos) ThrowDataToys(); else if (!routine.compare("Chi2Scan1D")) Create1DScans(); else if (!routine.compare("Chi2Scan2D")) Chi2Scan2D(); else fitstate = RunFitRoutine(routine); // If ending early break here if (fitstate == kFitFinished || fitstate == kNoChange) { NUIS_LOG(FIT, "Ending fit routines loop."); break; } } return; } //************************************* int MinimizerRoutines::RunFitRoutine(std::string routine) { //************************************* int endfits = kFitUnfinished; // set fitter at the current start values fOutputRootFile->cd(); SetupFitter(routine); // choose what to do with the minimizer depending on routine. if (!routine.compare("Migrad") or !routine.compare("Simplex") or !routine.compare("Combined") or !routine.compare("Brute") or !routine.compare("Fumili") or !routine.compare("ConjugateFR") or !routine.compare("ConjugatePR") or !routine.compare("BFGS") or !routine.compare("BFGS2") or !routine.compare("SteepDesc") or // !routine.compare("GSLMulti") or !routine.compare("GSLSimAn") or !routine.compare("MCMC")) { if (fMinimizer->NFree() > 0) { NUIS_LOG(FIT, fMinimizer->Minimize()); GetMinimizerState(); } } // other otptions else if (!routine.compare("Contour")) { CreateContours(); } return endfits; } //************************************* void MinimizerRoutines::PrintState() { //************************************* NUIS_LOG(FIT, "------------"); // Count max size int maxcount = 0; for (UInt_t i = 0; i < fParams.size(); i++) { maxcount = max(int(fParams[i].size()), maxcount); } // Header NUIS_LOG(FIT, " # " << left << setw(maxcount) << "Parameter " - << " = " << setw(10) << "Value" - << " +- " << setw(10) << "Error" - << " " << setw(8) << "(Units)" - << " " << setw(10) << "Conv. Val" - << " +- " << setw(10) << "Conv. Err" - << " " << setw(8) << "(Units)"); + << " = " << setw(10) << "Value" + << " +- " << setw(10) << "Error" + << " " << setw(8) << "(Units)" + << " " << setw(10) << "Conv. Val" + << " +- " << setw(10) << "Conv. Err" + << " " << setw(8) << "(Units)"); // Parameters for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); std::string typestr = FitBase::ConvDialType(fTypeVals[syst]); std::string curunits = "(sig.)"; double curval = fCurVals[syst]; double curerr = fErrorVals[syst]; if (fStateVals[syst].find("ABS") != std::string::npos) { curval = FitBase::RWSigmaToAbs(typestr, syst, curval); curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); curunits = "(Abs.)"; } else if (fStateVals[syst].find("FRAC") != std::string::npos) { curval = FitBase::RWSigmaToFrac(typestr, syst, curval); curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) - FitBase::RWSigmaToFrac(typestr, syst, 0.0)); curunits = "(Frac)"; } std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")"; double convval = FitBase::RWSigmaToAbs(typestr, syst, curval); double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); std::ostringstream curparstring; curparstring << " " << setw(3) << left << i << ". " << setw(maxcount) << syst << " = " << setw(10) << curval << " +- " << setw(10) << curerr << " " << setw(8) << curunits << " " << setw(10) << convval << " +- " << setw(10) << converr << " " << setw(8) << convunits; NUIS_LOG(FIT, curparstring.str()); } NUIS_LOG(FIT, "------------"); double like = fSampleFCN->GetLikelihood(); - NUIS_LOG(FIT, std::left << std::setw(46) << "Likelihood for JointFCN: " << like); + NUIS_LOG(FIT, + std::left << std::setw(46) << "Likelihood for JointFCN: " << like); NUIS_LOG(FIT, "------------"); } //************************************* void MinimizerRoutines::GetMinimizerState() { //************************************* NUIS_LOG(FIT, "Minimizer State: "); // Get X and Err const double *values = fMinimizer->X(); const double *errors = fMinimizer->Errors(); // int ipar = 0; for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); fCurVals[syst] = values[i]; fErrorVals[syst] = errors[i]; } PrintState(); // Covar SetupCovariance(); if (fMinimizer->CovMatrixStatus() > 0) { // Fill Full Covar NUIS_LOG(FIT, "Filling covariance"); for (int i = 0; i < fCovar->GetNbinsX(); i++) { for (int j = 0; j < fCovar->GetNbinsY(); j++) { fCovar->SetBinContent(i + 1, j + 1, fMinimizer->CovMatrix(i, j)); } } int freex = 0; int freey = 0; for (int i = 0; i < fCovar->GetNbinsX(); i++) { if (fMinimizer->IsFixedVariable(i)) continue; freey = 0; for (int j = 0; j < fCovar->GetNbinsY(); j++) { if (fMinimizer->IsFixedVariable(j)) continue; fCovFree->SetBinContent(freex + 1, freey + 1, fMinimizer->CovMatrix(i, j)); freey++; } freex++; } fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation"); fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition"); if (fMinimizer->NFree() > 0) { fCorFree = PlotUtils::GetCorrelationPlot(fCovFree, "correlation_free"); fDecFree = PlotUtils::GetDecompPlot(fCovFree, "decomposition_free"); } } return; }; //************************************* void MinimizerRoutines::LowStatRoutine(std::string routine) { //************************************* NUIS_LOG(FIT, "Running Low Statistics Routine: " << routine); int lowstatsevents = FitPar::Config().GetParI("LOWSTATEVENTS"); int maxevents = FitPar::Config().GetParI("MAXEVENTS"); int verbosity = FitPar::Config().GetParI("VERBOSITY"); std::string trueroutine = routine; std::string substring = "LowStat"; trueroutine.erase(trueroutine.find(substring), substring.length()); // Set MAX EVENTS=1000 Config::SetPar("MAXEVENTS", lowstatsevents); Config::SetPar("VERBOSITY", 3); SetupFCN(); RunFitRoutine(trueroutine); Config::SetPar("MAXEVENTS", maxevents); SetupFCN(); Config::SetPar("VERBOSITY", verbosity); return; } //************************************* void MinimizerRoutines::Create1DScans() { //************************************* // 1D Scan Routine // Steps through all free parameters about nominal using the step size // Creates a graph for each free parameter // At the current point create a 1D Scan for all parametes (Uncorrelated) for (UInt_t i = 0; i < fParams.size(); i++) { if (fFixVals[fParams[i]]) continue; NUIS_LOG(FIT, "Running 1D Scan for " << fParams[i]); fSampleFCN->CreateIterationTree(fParams[i] + "_scan1D_iterations", FitBase::GetRW()); double scanmiddlepoint = fCurVals[fParams[i]]; // Determine N points needed double limlow = fMinVals[fParams[i]]; double limhigh = fMaxVals[fParams[i]]; double step = fStepVals[fParams[i]]; int npoints = int(fabs(limhigh - limlow) / (step + 0.)); TH1D *contour = new TH1D(("Chi2Scan1D_" + fParams[i]).c_str(), ("Chi2Scan1D_" + fParams[i] + ";" + fParams[i]).c_str(), npoints, limlow, limhigh); // Fill bins for (int x = 0; x < contour->GetNbinsX(); x++) { // Set X Val fCurVals[fParams[i]] = contour->GetXaxis()->GetBinCenter(x + 1); // Run Eval double *vals = FitUtils::GetArrayFromMap(fParams, fCurVals); double chi2 = fSampleFCN->DoEval(vals); delete vals; // Fill Contour contour->SetBinContent(x + 1, chi2); } // Save contour contour->Write(); // Reset Parameter fCurVals[fParams[i]] = scanmiddlepoint; // Save TTree fSampleFCN->WriteIterationTree(); } return; } //************************************* void MinimizerRoutines::Chi2Scan2D() { //************************************* // Chi2 Scan 2D // Creates a 2D chi2 scan by stepping through all free parameters // Works for all pairwise combos of free parameters // Scan I for (UInt_t i = 0; i < fParams.size(); i++) { if (fFixVals[fParams[i]]) continue; // Scan J for (UInt_t j = 0; j < i; j++) { if (fFixVals[fParams[j]]) continue; fSampleFCN->CreateIterationTree(fParams[i] + "_" + fParams[j] + "_" + "scan2D_iterations", FitBase::GetRW()); double scanmid_i = fCurVals[fParams[i]]; double scanmid_j = fCurVals[fParams[j]]; double limlow_i = fMinVals[fParams[i]]; double limhigh_i = fMaxVals[fParams[i]]; double step_i = fStepVals[fParams[i]]; double limlow_j = fMinVals[fParams[j]]; double limhigh_j = fMaxVals[fParams[j]]; double step_j = fStepVals[fParams[j]]; int npoints_i = int(fabs(limhigh_i - limlow_i) / (step_i + 0.)) + 1; int npoints_j = int(fabs(limhigh_j - limlow_j) / (step_j + 0.)) + 1; TH2D *contour = new TH2D( ("Chi2Scan2D_" + fParams[i] + "_" + fParams[j]).c_str(), ("Chi2Scan2D_" + fParams[i] + "_" + fParams[j] + ";" + fParams[i] + ";" + fParams[j]) .c_str(), npoints_i, limlow_i, limhigh_i, npoints_j, limlow_j, limhigh_j); // Begin Scan NUIS_LOG(FIT, "Running scan for " << fParams[i] << " " << fParams[j]); // Fill bins for (int x = 0; x < contour->GetNbinsX(); x++) { // Set X Val fCurVals[fParams[i]] = contour->GetXaxis()->GetBinCenter(x + 1); // Loop Y for (int y = 0; y < contour->GetNbinsY(); y++) { // Set Y Val fCurVals[fParams[j]] = contour->GetYaxis()->GetBinCenter(y + 1); // Run Eval double *vals = FitUtils::GetArrayFromMap(fParams, fCurVals); double chi2 = fSampleFCN->DoEval(vals); delete vals; // Fill Contour contour->SetBinContent(x + 1, y + 1, chi2); fCurVals[fParams[j]] = scanmid_j; } fCurVals[fParams[i]] = scanmid_i; fCurVals[fParams[j]] = scanmid_j; } // Save contour contour->Write(); // Save Iterations fSampleFCN->WriteIterationTree(); } } return; } //************************************* void MinimizerRoutines::CreateContours() { //************************************* // Use MINUIT for this if possible NUIS_ABORT("Contours not yet implemented as it is really slow!"); return; } //************************************* int MinimizerRoutines::FixAtLimit() { //************************************* bool fixedparam = false; for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); if (fFixVals[syst]) continue; double curVal = fCurVals.at(syst); double minVal = fMinVals.at(syst); double maxVal = fMinVals.at(syst); if (fabs(curVal - minVal) < 0.0001) { fCurVals[syst] = minVal; fFixVals[syst] = true; fixedparam = true; } if (fabs(maxVal - curVal) < 0.0001) { fCurVals[syst] = maxVal; fFixVals[syst] = true; fixedparam = true; } } if (!fixedparam) { NUIS_LOG(FIT, "No dials needed fixing!"); return kNoChange; } else return kStateChange; } /* Write Functions */ //************************************* void MinimizerRoutines::SaveResults() { //************************************* fOutputRootFile->cd(); if (fMinimizer) { SetupCovariance(); SaveMinimizerState(); } SaveCurrentState(); } //************************************* void MinimizerRoutines::SaveMinimizerState() { //************************************* NUIS_LOG(FIT, "Saving Minimizer State"); if (!fMinimizer) { NUIS_ABORT("Can't save minimizer state without min object"); } // Save main fit tree fSampleFCN->WriteIterationTree(); // Get Vals and Errors GetMinimizerState(); // Save tree with fit status std::vector nameVect; std::vector valVect; std::vector errVect; std::vector minVect; std::vector maxVect; std::vector startVect; std::vector endfixVect; std::vector startfixVect; // int NFREEPARS = fMinimizer->NFree(); int NPARS = fMinimizer->NDim(); int ipar = 0; // Dial Vals for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams.at(i); nameVect.push_back(name); valVect.push_back(fCurVals.at(name)); errVect.push_back(fErrorVals.at(name)); minVect.push_back(fMinVals.at(name)); maxVect.push_back(fMaxVals.at(name)); startVect.push_back(fStartVals.at(name)); endfixVect.push_back(fFixVals.at(name)); startfixVect.push_back(fStartFixVals.at(name)); ipar++; } int NFREE = fMinimizer->NFree(); int NDIM = fMinimizer->NDim(); double CHI2 = fSampleFCN->GetLikelihood(); int NBINS = fSampleFCN->GetNDOF(); int NDOF = NBINS - NFREE; // Write fit results TTree *fit_tree = new TTree("fit_result", "fit_result"); fit_tree->Branch("parameter_names", &nameVect); fit_tree->Branch("parameter_values", &valVect); fit_tree->Branch("parameter_errors", &errVect); fit_tree->Branch("parameter_min", &minVect); fit_tree->Branch("parameter_max", &maxVect); fit_tree->Branch("parameter_start", &startVect); fit_tree->Branch("parameter_fix", &endfixVect); fit_tree->Branch("parameter_startfix", &startfixVect); fit_tree->Branch("CHI2", &CHI2, "CHI2/D"); fit_tree->Branch("NDOF", &NDOF, "NDOF/I"); fit_tree->Branch("NBINS", &NBINS, "NBINS/I"); fit_tree->Branch("NDIM", &NDIM, "NDIM/I"); fit_tree->Branch("NFREE", &NFREE, "NFREE/I"); fit_tree->Fill(); fit_tree->Write(); // Make dial variables TH1D dialvar = TH1D("fit_dials", "fit_dials", NPARS, 0, NPARS); TH1D startvar = TH1D("start_dials", "start_dials", NPARS, 0, NPARS); TH1D minvar = TH1D("min_dials", "min_dials", NPARS, 0, NPARS); TH1D maxvar = TH1D("max_dials", "max_dials", NPARS, 0, NPARS); TH1D dialvarfree = TH1D("fit_dials_free", "fit_dials_free", NFREE, 0, NFREE); TH1D startvarfree = TH1D("start_dials_free", "start_dials_free", NFREE, 0, NFREE); TH1D minvarfree = TH1D("min_dials_free", "min_dials_free", NFREE, 0, NFREE); TH1D maxvarfree = TH1D("max_dials_free", "max_dials_free", NFREE, 0, NFREE); int freecount = 0; for (UInt_t i = 0; i < nameVect.size(); i++) { std::string name = nameVect.at(i); dialvar.SetBinContent(i + 1, valVect.at(i)); dialvar.SetBinError(i + 1, errVect.at(i)); dialvar.GetXaxis()->SetBinLabel(i + 1, name.c_str()); startvar.SetBinContent(i + 1, startVect.at(i)); startvar.GetXaxis()->SetBinLabel(i + 1, name.c_str()); minvar.SetBinContent(i + 1, minVect.at(i)); minvar.GetXaxis()->SetBinLabel(i + 1, name.c_str()); maxvar.SetBinContent(i + 1, maxVect.at(i)); maxvar.GetXaxis()->SetBinLabel(i + 1, name.c_str()); if (NFREE > 0) { if (!startfixVect.at(i)) { freecount++; dialvarfree.SetBinContent(freecount, valVect.at(i)); dialvarfree.SetBinError(freecount, errVect.at(i)); dialvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str()); startvarfree.SetBinContent(freecount, startVect.at(i)); startvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str()); minvarfree.SetBinContent(freecount, minVect.at(i)); minvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str()); maxvarfree.SetBinContent(freecount, maxVect.at(i)); maxvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str()); } } } // Save Dial Plots dialvar.Write(); startvar.Write(); minvar.Write(); maxvar.Write(); if (NFREE > 0) { dialvarfree.Write(); startvarfree.Write(); minvarfree.Write(); maxvarfree.Write(); } // Save fit_status plot TH1D statusplot = TH1D("fit_status", "fit_status", 8, 0, 8); std::string fit_labels[8] = {"status", "cov_status", "maxiter", "maxfunc", "iter", "func", "precision", "tolerance"}; double fit_vals[8]; fit_vals[0] = fMinimizer->Status() + 0.; fit_vals[1] = fMinimizer->CovMatrixStatus() + 0.; fit_vals[2] = fMinimizer->MaxIterations() + 0.; fit_vals[3] = fMinimizer->MaxFunctionCalls() + 0.; fit_vals[4] = fMinimizer->NIterations() + 0.; fit_vals[5] = fMinimizer->NCalls() + 0.; fit_vals[6] = fMinimizer->Precision() + 0.; fit_vals[7] = fMinimizer->Tolerance() + 0.; for (int i = 0; i < 8; i++) { statusplot.SetBinContent(i + 1, fit_vals[i]); statusplot.GetXaxis()->SetBinLabel(i + 1, fit_labels[i].c_str()); } statusplot.Write(); // Save Covars if (fCovar) fCovar->Write(); if (fCovFree) fCovFree->Write(); if (fCorrel) fCorrel->Write(); if (fCorFree) fCorFree->Write(); if (fDecomp) fDecomp->Write(); if (fDecFree) fDecFree->Write(); return; } //************************************* void MinimizerRoutines::SaveCurrentState(std::string subdir) { //************************************* NUIS_LOG(FIT, "Saving current full FCN predictions"); // Setup DIRS TDirectory *curdir = gDirectory; if (!subdir.empty()) { TDirectory *newdir = (TDirectory *)gDirectory->mkdir(subdir.c_str()); newdir->cd(); } FitBase::GetRW()->Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->Write(); // Change back to current DIR curdir->cd(); return; } //************************************* void MinimizerRoutines::SaveNominal() { //************************************* fOutputRootFile->cd(); NUIS_LOG(FIT, "Saving Nominal Predictions (be cautious with this)"); FitBase::GetRW()->Reconfigure(); SaveCurrentState("nominal"); }; //************************************* void MinimizerRoutines::SavePrefit() { //************************************* fOutputRootFile->cd(); NUIS_LOG(FIT, "Saving Prefit Predictions"); UpdateRWEngine(fStartVals); SaveCurrentState("prefit"); UpdateRWEngine(fCurVals); }; /* MISC Functions */ //************************************* int MinimizerRoutines::GetStatus() { //************************************* return 0; } //************************************* void MinimizerRoutines::SetupCovariance() { //************************************* // Remove covares if they exist if (fCovar) delete fCovar; if (fCovFree) delete fCovFree; if (fCorrel) delete fCorrel; if (fCorFree) delete fCorFree; if (fDecomp) delete fDecomp; if (fDecFree) delete fDecFree; NUIS_LOG(FIT, "Building covariance matrix.."); int NFREE = 0; int NDIM = 0; // Get NFREE from min or from vals (for cases when doing throws) if (fMinimizer) { NFREE = fMinimizer->NFree(); NDIM = fMinimizer->NDim(); } else { NDIM = fParams.size(); for (UInt_t i = 0; i < fParams.size(); i++) { NUIS_LOG(FIT, "Getting Param " << fParams[i]); if (!fFixVals[fParams[i]]) NFREE++; } } if (NDIM == 0) return; NUIS_LOG(FIT, "NFREE == " << NFREE); fCovar = new TH2D("covariance", "covariance", NDIM, 0, NDIM, NDIM, 0, NDIM); if (NFREE > 0) { fCovFree = new TH2D("covariance_free", "covariance_free", NFREE, 0, NFREE, NFREE, 0, NFREE); } else { fCovFree = NULL; } // Set Bin Labels int countall = 0; int countfree = 0; for (UInt_t i = 0; i < fParams.size(); i++) { fCovar->GetXaxis()->SetBinLabel(countall + 1, fParams[i].c_str()); fCovar->GetYaxis()->SetBinLabel(countall + 1, fParams[i].c_str()); countall++; if (!fFixVals[fParams[i]] and NFREE > 0) { fCovFree->GetXaxis()->SetBinLabel(countfree + 1, fParams[i].c_str()); fCovFree->GetYaxis()->SetBinLabel(countfree + 1, fParams[i].c_str()); countfree++; } } fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation"); fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition"); if (NFREE > 0) { fCorFree = PlotUtils::GetCorrelationPlot(fCovFree, "correlation_free"); fDecFree = PlotUtils::GetDecompPlot(fCovFree, "decomposition_free"); } else { fCorFree = NULL; fDecFree = NULL; } return; }; //************************************* void MinimizerRoutines::ThrowCovariance(bool uniformly) { //************************************* std::vector rands; if (!fDecFree) { NUIS_ERR(WRN, "Trying to throw 0 free parameters"); return; } // Generate Random Gaussians for (Int_t i = 0; i < fDecFree->GetNbinsX(); i++) { rands.push_back(gRandom->Gaus(0.0, 1.0)); } // Reset Thrown Values for (UInt_t i = 0; i < fParams.size(); i++) { fThrownVals[fParams[i]] = fCurVals[fParams[i]]; } // Loop and get decomp for (Int_t i = 0; i < fDecFree->GetNbinsX(); i++) { std::string parname = std::string(fDecFree->GetXaxis()->GetBinLabel(i + 1)); double mod = 0.0; if (!uniformly) { for (Int_t j = 0; j < fDecFree->GetNbinsY(); j++) { mod += rands[j] * fDecFree->GetBinContent(j + 1, i + 1); } } if (fCurVals.find(parname) != fCurVals.end()) { if (uniformly) fThrownVals[parname] = gRandom->Uniform(fMinVals[parname], fMaxVals[parname]); else { fThrownVals[parname] = fCurVals[parname] + mod; } } } // Check Limits for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams[i]; if (fFixVals[syst]) continue; if (fThrownVals[syst] < fMinVals[syst]) fThrownVals[syst] = fMinVals[syst]; if (fThrownVals[syst] > fMaxVals[syst]) fThrownVals[syst] = fMaxVals[syst]; } return; }; //************************************* void MinimizerRoutines::GenerateErrorBands() { //************************************* TDirectory *errorDIR = (TDirectory *)fOutputRootFile->mkdir("error_bands"); errorDIR->cd(); // Make a second file to store throws std::string tempFileName = fOutputFile; if (tempFileName.find(".root") != std::string::npos) tempFileName.erase(tempFileName.find(".root"), 5); tempFileName += ".throws.root"; TFile *tempfile = new TFile(tempFileName.c_str(), "RECREATE"); tempfile->cd(); int nthrows = FitPar::Config().GetParI("error_throws"); UpdateRWEngine(fCurVals); fSampleFCN->ReconfigureAllEvents(); TDirectory *nominal = (TDirectory *)tempfile->mkdir("nominal"); nominal->cd(); fSampleFCN->Write(); TDirectory *outnominal = (TDirectory *)fOutputRootFile->mkdir("nominal_throw"); outnominal->cd(); fSampleFCN->Write(); errorDIR->cd(); TTree *parameterTree = new TTree("throws", "throws"); double chi2; for (UInt_t i = 0; i < fParams.size(); i++) parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]], (fParams[i] + "/D").c_str()); parameterTree->Branch("chi2", &chi2, "chi2/D"); bool uniformly = FitPar::Config().GetParB("error_uniform"); // Run Throws and save for (Int_t i = 0; i < nthrows; i++) { TDirectory *throwfolder = (TDirectory *)tempfile->mkdir(Form("throw_%i", i)); throwfolder->cd(); // Generate Random Parameter Throw ThrowCovariance(uniformly); // Run Eval double *vals = FitUtils::GetArrayFromMap(fParams, fThrownVals); chi2 = fSampleFCN->DoEval(vals); delete vals; // Save the FCN fSampleFCN->Write(); parameterTree->Fill(); } errorDIR->cd(); fDecFree->Write(); fCovFree->Write(); parameterTree->Write(); delete parameterTree; // Now go through the keys in the temporary file and look for TH1D, and TH2D // plots TIter next(nominal->GetListOfKeys()); TKey *key; while ((key = (TKey *)next())) { TClass *cl = gROOT->GetClass(key->GetClassName()); if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue; TH1D *baseplot = (TH1D *)key->ReadObj(); std::string plotname = std::string(baseplot->GetName()); int nbins = baseplot->GetNbinsX() * baseplot->GetNbinsY(); // Setup TProfile with RMS option TProfile *tprof = new TProfile((plotname + "_prof").c_str(), (plotname + "_prof").c_str(), nbins, 0, nbins, "S"); // Setup The TTREE double *bincontents; bincontents = new double[nbins]; double *binlowest; binlowest = new double[nbins]; double *binhighest; binhighest = new double[nbins]; errorDIR->cd(); TTree *bintree = new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str()); for (Int_t i = 0; i < nbins; i++) { bincontents[i] = 0.0; binhighest[i] = 0.0; binlowest[i] = 0.0; bintree->Branch(Form("content_%i", i), &bincontents[i], Form("content_%i/D", i)); } for (Int_t i = 0; i < nthrows; i++) { TH1 *newplot = (TH1 *)tempfile->Get(Form(("throw_%i/" + plotname).c_str(), i)); for (Int_t j = 0; j < nbins; j++) { tprof->Fill(j + 0.5, newplot->GetBinContent(j + 1)); bincontents[j] = newplot->GetBinContent(j + 1); if (bincontents[j] < binlowest[j] or i == 0) binlowest[j] = bincontents[j]; if (bincontents[j] > binhighest[j] or i == 0) binhighest[j] = bincontents[j]; } errorDIR->cd(); bintree->Fill(); delete newplot; } errorDIR->cd(); for (Int_t j = 0; j < nbins; j++) { if (!uniformly) { baseplot->SetBinError(j + 1, tprof->GetBinError(j + 1)); } else { baseplot->SetBinContent(j + 1, (binlowest[j] + binhighest[j]) / 2.0); baseplot->SetBinError(j + 1, (binhighest[j] - binlowest[j]) / 2.0); } } errorDIR->cd(); baseplot->Write(); tprof->Write(); bintree->Write(); delete baseplot; delete tprof; delete bintree; delete[] bincontents; } return; }; void MinimizerRoutines::ThrowDataToys() { NUIS_LOG(FIT, "Generating Toy Data Throws"); int verb = Config::GetParI("VERBOSITY"); SETVERBOSITY(FIT); int nthrows = FitPar::Config().GetParI("NToyThrows"); double maxlike = -1.0; double minlike = -1.0; std::vector values; for (int i = 0; i < 1.E4; i++) { fSampleFCN->ThrowDataToy(); double like = fSampleFCN->GetLikelihood(); values.push_back(like); if (maxlike == -1.0 or like > maxlike) maxlike = like; if (minlike == -1.0 or like < minlike) minlike = like; } SETVERBOSITY(verb); // Fill Histogram TH1D *likes = new TH1D("toydatalikelihood", "toydatalikelihood", int(sqrt(nthrows)), minlike, maxlike); for (size_t i = 0; i < values.size(); i++) { likes->Fill(values[i]); } // Save to file NUIS_LOG(FIT, "Writing toy data throws"); fOutputRootFile->cd(); likes->Write(); } diff --git a/src/Routines/MinimizerRoutines.h b/src/Routines/MinimizerRoutines.h index 3918e7a..dcff46b 100755 --- a/src/Routines/MinimizerRoutines.h +++ b/src/Routines/MinimizerRoutines.h @@ -1,278 +1,284 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef MINIMIZER_ROUTINES_H #define MINIMIZER_ROUTINES_H /*! * \addtogroup Minimizer * @{ */ #include "TH1.h" #include "TF1.h" #include "TMatrixD.h" #include "TVectorD.h" #ifdef ROOT6_USE_FIT_FITTER_INTERFACE #include "Fit/Fitter.h" #else #include "Minuit2/FCNBase.h" #include "TFitterMinuit.h" #endif #include "TSystem.h" #include "TFile.h" #include "TProfile.h" #include #include #include #include #include #include "FitEvent.h" #include "JointFCN.h" #include "MinimizerFCN.h" #include "Math/Minimizer.h" #include "Math/Factory.h" #include "Math/Functor.h" #include "FitLogger.h" #include "ParserUtils.h" enum minstate { kErrorStatus = -1, kGoodStatus, kFitError, kNoChange, kFitFinished, kFitUnfinished, kStateChange, }; //************************************* //! Collects all possible fit routines into a single class to avoid repeated code class MinimizerRoutines{ //************************************* public: /* Constructor/Destructor */ MinimizerRoutines(); //! Constructor reads in arguments given at the command line for the fit here. MinimizerRoutines(int argc, char* argv[]); //! Default destructor ~MinimizerRoutines(); //! Reset everything to default/NULL void Init(); /* Input Functions */ //! Splits the arguments ready for initial setup void ParseArgs(int argc, char* argv[]); //! Sorts out configuration and verbosity right at the very start. //! Calls readCard to set everything else up. void InitialSetup(); //! Loops through each line of the card file and passes it to other read functions void ReadCard(std::string cardfile); //! Check for parameter string in the line and assign the correct type. //! Fills maps for each of the parameters int ReadParameters(std::string parstring); //! Reads in fake parameters and assigns them (Requires the parameter to be included as a normal parameter as well) int ReadFakeDataPars(std::string parstring); //! Read in the samples so we can set up the free normalisation dials if required int ReadSamples(std::string sampleString); void SetupMinimizerFromXML(); /* Setup Functions */ //! Setup the configuration given the arguments passed at the commandline and card file void SetupConfig(); //! Setups up our custom RW engine with all the parameters passed in the card file void SetupRWEngine(); //! Setups up the jointFCN. void SetupFCN(); //! Sets up the minimizerObj for ROOT. there are cases where this is called repeatedly, e.g. If you are using a brute force scan before using Migrad. void SetupFitter(std::string routine); //! Set the current data histograms in each sample to the fake data. void SetFakeData(); //! Setup the covariances with the correct dimensions. At the start this is either uncorrelated or merged given all the input covariances. //! At the end of the fit this produces the blank covariances which can then be filled by the minimizerObj with best fit covariances. void SetupCovariance(); /* Fitting Functions */ //! Main function to actually start iterating over the different required fit routines void Run(); //! Given a new map change the values that the RW engine is currently set to void UpdateRWEngine(std::map& updateVals); //! Given a single routine (see tutorial for options) run that fit routine now. int RunFitRoutine(std::string routine); //! Get the current state of minimizerObj and fill it into currentVals and currentNorms void GetMinimizerState(); //! Print current value void PrintState(); //! Performs a fit routine where the input.maxevents is set to a much lower value to try and move closer to the best fit minimum. void LowStatRoutine(std::string routine); //! Perform a chi2 scan in 1D around the current point void Create1DScans(); //! Perform a chi2 scan in 2D around the current point void Chi2Scan2D(); //! Currently a placeholder NEEDS UPDATING void CreateContours(); //! If any currentVals are close to the limits set them to the limit and fix them int FixAtLimit(); //! Throw the current covariance of dial values we have, and fill the thrownVals and thrownNorms maps. //! If uniformly is true parameters will be thrown uniformly between their upper and lower limits. void ThrowCovariance(bool uniformly); //! Given the covariance we currently have generate error bands by throwing the covariance. //! The FitPar config "error_uniform" defines whether to throw using the covariance or uniformly. //! The FitPar config "error_throws" defines how many throws are needed. //! Currently only supports TH1D plots. void GenerateErrorBands(); /* Write Functions */ //! Write plots and TTrees listing the minimizerObj result of the fit to file void SaveMinimizerState(); //! Save the sample plots for current MC //! dir if not empty forces plots to be saved in a subdirectory of outputfile void SaveCurrentState(std::string subdir=""); //! Save starting predictions into a seperate folder void SaveNominal(); //! Save predictions before the fit is ran into a seperate folder void SavePrefit(); void SaveResults(); /* MISC Functions */ //! Get previous fit status from a file Int_t GetStatus(); /// Makes a histogram of likelihoods when throwing the data according to its statistics void ThrowDataToys(); protected: //! Our Custom ReWeight Object FitWeight* rw; std::string fOutputFile; std::string fInputFile; TFile* fInputRootFile; TFile* fOutputRootFile; //! Flag for whether the fit should be continued if an output file is already found. bool fitContinue; //! Minimizer Object for handling roots different minimizer methods ROOT::Math::Minimizer* fMinimizer; JointFCN* fSampleFCN; MinimizerFCN* fMinimizerFCN; ROOT::Math::Functor* fCallFunctor; int nfreepars; std::string fCardFile; std::string fStrategy; std::vector fRoutines; std::string fAllowedRoutines; std::string fFakeDataInput; // Input Dial Vals //! Vector of dial names std::vector fParams; std::map fStateVals; std::map fStartVals; std::map fCurVals; std::map fErrorVals; std::map fMinVals; std::map fMaxVals; std::map fStepVals; std::map fTypeVals; std::map fFixVals; std::map fStartFixVals; + struct mirror_param { + double mirror_value; + bool mirror_above; + }; + std::map fMirroredParams; + //! Vector of fake parameter names std::map fFakeVals; //! Map of thrown parameter names and values (After ThrowCovariance) std::map fThrownVals; TH2D* fCorrel; TH2D* fDecomp; TH2D* fCovar; TH2D* fCorFree; TH2D* fDecFree; TH2D* fCovFree; nuiskey fCompKey; }; /*! @} */ #endif