diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx index f3eb4e6..0c4609d 100755 --- a/src/FCN/JointFCN.cxx +++ b/src/FCN/JointFCN.cxx @@ -1,1136 +1,1119 @@ #include "JointFCN.h" #include #include "FitUtils.h" - //*************************************************** 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) { -//*************************************************** + //*************************************************** LOG(FIT) << " Creating new iteration container! " << std::endl; DestroyIterationTree(); fIterationTreeName = name; // Add sample likelihoods and ndof - for (MeasListConstIter iter = fSamples.begin(); - iter != fSamples.end(); + 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++) { - + 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(); + fSampleN = fSamples.size() + fPulls.size(); fSampleLikes = new double[fSampleN]; - fSampleNDOF = new int[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 ); + fNameValues.push_back(dials[i]); + fCurrentValues.push_back(0.0); } - fNDials = dials.size(); + fNDials = dials.size(); fDialVals = new double[fNDials]; // Set IterationTree Flag fIterationTree = true; - } //*************************************************** void JointFCN::DestroyIterationTree() { -//*************************************************** + //*************************************************** fIterationCount.clear(); fCurrentValues.clear(); fNameValues.clear(); fIterationValues.clear(); - } //*************************************************** void JointFCN::WriteIterationTree() { -//*************************************************** + //*************************************************** LOG(FIT) << "Writing iteration tree" << std::endl; // Make a new TTree - TTree* itree = new TTree(fIterationTreeName.c_str(), - fIterationTreeName.c_str()); + TTree* itree = + new TTree(fIterationTreeName.c_str(), fIterationTreeName.c_str()); double* vals = new double[fNameValues.size()]; int count = 0; itree->Branch("iteration", &count, "Iteration/I"); for (int i = 0; i < fNameValues.size(); i++) { - itree->Branch( fNameValues[i].c_str(), - &vals[i], - (fNameValues[i] + "/D").c_str() ); + 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 ); + fIterationCount.push_back(fCurIter); fIterationValues.push_back(fCurrentValues); } //*************************************************** double JointFCN::DoEval(const double* x) { //*************************************************** // WEIGHT ENGINE fDialChanged = FitBase::GetRW()->HasRWDialChanged(x); FitBase::GetRW()->UpdateWeightEngine(x); if (fDialChanged) { FitBase::GetRW()->Reconfigure(); FitBase::EvtManager().ResetWeightFlags(); } if (LOG_LEVEL(REC)) { FitBase::GetRW()->Print(); } // SORT SAMPLES ReconfigureSamples(); // GET TEST STAT fLikelihood = GetLikelihood(); - fNDOF = GetNDOF(); + fNDOF = GetNDOF(); // PRINT PROGRESS LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood << std::endl; // UPDATE TREE if (fIterationTree) FillIterationTree(FitBase::GetRW()); return fLikelihood; } //*************************************************** int JointFCN::GetNDOF() { //*************************************************** int totaldof = 0; int count = 0; // Total number of Free bins in each MC prediction for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; int dof = exp->GetNDOF(); // Save Seperate DOF if (fIterationTree) { fSampleNDOF[count] = dof; } // Add to total totaldof += dof; count++; } // Loop over pulls for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; double dof = pull->GetLikelihood(); // Save seperate DOF if (fIterationTree) { fSampleNDOF[count] = dof; } // Add to total totaldof += dof; count++; } // Set Data Variable if (fIterationTree) { fSampleNDOF[count] = totaldof; } return totaldof; } //*************************************************** double JointFCN::GetLikelihood() { //*************************************************** LOG(MIN) << std::left << std::setw(43) << "Getting likelihoods..." << " : " << "-2logL" << std::endl; // Loop and add up likelihoods in an uncorrelated way double like = 0.0; int count = 0; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; double newlike = exp->GetLikelihood(); int ndof = exp->GetNDOF(); // Save seperate likelihoods if (fIterationTree) { fSampleLikes[count] = newlike; } LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : " << newlike << "/" << ndof << std::endl; // Add Weight Scaling // like *= FitBase::GetRW()->GetSampleLikelihoodWeight(exp->GetName()); // Add to total like += newlike; count++; } // Loop over pulls for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; double newlike = pull->GetLikelihood(); // Save seperate likelihoods if (fIterationTree) { fSampleLikes[count] = newlike; } // Add to total like += newlike; count++; } // Set Data Variable fLikelihood = like; if (fIterationTree) { fSampleLikes[count] = fLikelihood; } return like; }; void JointFCN::LoadSamples(std::vector samplekeys) { LOG(MIN) << "Loading Samples : " << samplekeys.size() << std::endl; for (size_t i = 0; i < samplekeys.size(); i++) { nuiskey key = samplekeys[i]; // Get Sample Options std::string samplename = key.GetS("name"); std::string samplefile = key.GetS("input"); std::string sampletype = key.GetS("type"); std::string fakeData = ""; LOG(MIN) << "Loading Sample : " << samplename << std::endl; fOutputDir->cd(); MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(key); if (!NewLoadedSample) { ERR(FTL) << "Could not load sample provided: " << samplename << std::endl; ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx" << std::endl; throw; } else { fSamples.push_back(NewLoadedSample); } } } //*************************************************** void JointFCN::LoadPulls(std::vector pullkeys) { -//*************************************************** + //*************************************************** for (size_t i = 0; i < pullkeys.size(); i++) { nuiskey key = pullkeys[i]; std::string pullname = key.GetS("name"); std::string pullfile = key.GetS("input"); std::string pulltype = key.GetS("type"); fOutputDir->cd(); fPulls.push_back(new ParamPull(pullname, pullfile, pulltype)); } } //*************************************************** void JointFCN::ReconfigureSamples(bool fullconfig) { -//*************************************************** + //*************************************************** int starttime = time(NULL); LOG(REC) << "Starting Reconfigure iter. " << this->fCurIter << std::endl; // std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled << // std::endl; // Event Manager Reconf if (fUsingEventManager) { if (!fullconfig and fMCFilled) ReconfigureFastUsingManager(); else ReconfigureUsingManager(); } else { // Loop over all Measurement Classes for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; // If RW Either do signal or full reconfigure. if (fDialChanged or !fMCFilled or fullconfig) { if (!fullconfig and fMCFilled) exp->ReconfigureFast(); else exp->Reconfigure(); // If RW Not needed just do normalisation } else { exp->Renormalise(); } } } // Loop over pulls and update for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; pull->Reconfigure(); } fMCFilled = true; LOG(MIN) << "Finished Reconfigure iter. " << fCurIter << " in " << time(NULL) - starttime << "s" << std::endl; fCurIter++; } //*************************************************** void JointFCN::ReconfigureSignal() { -//*************************************************** + //*************************************************** ReconfigureSamples(false); } //*************************************************** void JointFCN::ReconfigureAllEvents() { //*************************************************** FitBase::GetRW()->Reconfigure(); FitBase::EvtManager().ResetWeightFlags(); ReconfigureSamples(true); } std::vector JointFCN::GetInputList() { std::vector InputList; fIsAllSplines = true; MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); std::vector subsamples = exp->GetSubSamples(); for (size_t i = 0; i < subsamples.size(); i++) { InputHandlerBase* inp = subsamples[i]->GetInput(); if (std::find(InputList.begin(), InputList.end(), inp) == InputList.end()) { if (subsamples[i]->GetInput()->GetType() != kSPLINEPARAMETER) fIsAllSplines = false; InputList.push_back(subsamples[i]->GetInput()); } } } return InputList; } std::vector JointFCN::GetSubSampleList() { std::vector SampleList; MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); std::vector subsamples = exp->GetSubSamples(); for (size_t i = 0; i < subsamples.size(); i++) { SampleList.push_back(subsamples[i]); } } return SampleList; } //*************************************************** void JointFCN::ReconfigureUsingManager() { -//*************************************************** + //*************************************************** // 'Slow' Event Manager Reconfigure LOG(REC) << "Event Manager Reconfigure" << std::endl; int timestart = time(NULL); // Reset all samples MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ResetAll(); } // If we are siving signal, reset all containers. bool savesignal = (FitPar::Config().GetParB("SignalReconfigures")); if (savesignal) { // Reset all of our event signal vectors fSignalEventBoxes.clear(); fSignalEventFlags.clear(); fSampleSignalFlags.clear(); fSignalEventSplines.clear(); } // Make sure we have a list of inputs if (fInputList.empty()) { fInputList = GetInputList(); fSubSampleList = GetSubSampleList(); } // If all inputs are splines make sure the readers are told // they need to be reconfigured. std::vector::iterator inp_iter = fInputList.begin(); if (fIsAllSplines) { for (; inp_iter != fInputList.end(); inp_iter++) { InputHandlerBase* curinput = (*inp_iter); // Tell reader in each BaseEvent it needs a Reconfigure next weight calc. BaseFitEvt* curevent = curinput->FirstBaseEvent(); if (curevent->fSplineRead) { curevent->fSplineRead->SetNeedsReconfigure(true); } } } // MAIN INPUT LOOP ==================== int fillcount = 0; int inputcount = 0; inp_iter = fInputList.begin(); // Loop over each input in manager for (; inp_iter != fInputList.end(); inp_iter++) { InputHandlerBase* curinput = (*inp_iter); // Get event information FitEvent* curevent = curinput->FirstNuisanceEvent(); curinput->CreateCache(); int i = 0; int nevents = curinput->GetNEvents(); int countwidth = nevents / 5; // Start event loop iterating until we get a NULL pointer. while (curevent) { // Get Event Weight curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent); curevent->Weight = curevent->RWWeight * curevent->InputWeight; double rwweight = curevent->Weight; // std::cout << "RWWeight = " << curevent->RWWeight << " " << // curevent->InputWeight << std::endl; // Logging // std::cout << CHECKLOG(1) << std::endl; if (LOGGING(REC)) { - if (i % countwidth == 0) { + if (countwidth && (i % countwidth == 0)) { QLOG(REC, curinput->GetName() - << " : Processed " << i << " events. [M, W] = [" - << curevent->Mode << ", " << rwweight << "]"); + << " : Processed " << i << " events. [M, W] = [" + << curevent->Mode << ", " << rwweight << "]"); } } // Setup flag for if signal found in at least one sample bool foundsignal = false; // Create a new signal bitset for this event std::vector signalbitset(fSubSampleList.size()); // Create a new signal box vector for this event std::vector signalboxes; // Start measurement iterator size_t measitercount = 0; std::vector::iterator meas_iter = - fSubSampleList.begin(); + fSubSampleList.begin(); // Loop over all subsamples (sub in JointMeas) for (; meas_iter != fSubSampleList.end(); meas_iter++) { MeasurementBase* curmeas = (*meas_iter); // Compare input pointers, to current input, skip if not. // Pointer tells us if it matches without doing ID checks. if (curinput != curmeas->GetInput()) { if (savesignal) { // Set bit to 0 as definitely not signal signalbitset[measitercount] = 0; } // Count up what measurement we are on. measitercount++; // Skip sample as input not signal. continue; } // Fill events for matching inputs. MeasurementVariableBox* box = curmeas->FillVariableBox(curevent); bool signal = curmeas->isSignal(curevent); curmeas->SetSignal(signal); curmeas->FillHistograms(curevent->Weight); // If its Signal tally up fills if (signal) { fillcount++; } // If we are saving signal/splines fill the bitset if (savesignal) { signalbitset[measitercount] = signal; } // If signal save a clone of the event box for use later. if (savesignal and signal) { foundsignal = true; signalboxes.push_back(box->CloneSignalBox()); } // Keep track of Measurement we are on. measitercount++; } // Once we've filled the measurements, if saving signal // push back if any sample flagged this event as signal if (savesignal) { fSignalEventFlags.push_back(foundsignal); } // Save the vector of signal boxes for this event if (savesignal and foundsignal) { fSignalEventBoxes.push_back(signalboxes); fSampleSignalFlags.push_back(signalbitset); } // If all inputs are splines we can save the spline coefficients // for fast in memory reconfigures later. if (fIsAllSplines and savesignal and foundsignal) { // Make temp vector to push back with std::vector coeff; for (size_t l = 0; l < (UInt_t)curevent->fSplineRead->GetNPar(); l++) { coeff.push_back(curevent->fSplineCoeff[l]); } // Push back to signal event splines. Kept in sync with // fSignalEventBoxes size. // int splinecount = fSignalEventSplines.size(); fSignalEventSplines.push_back(coeff); // if (splinecount % 1000 == 0) { // std::cout << "Pushed Back Coeff " << splinecount << " : "; // for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++) // { // std::cout << " " << fSignalEventSplines[splinecount][l]; // } // std::cout << std::endl; // } } // Clean up vectors once done with this event signalboxes.clear(); signalbitset.clear(); // Iterate to the next event. curevent = curinput->NextNuisanceEvent(); i++; } // curinput->RemoveCache(); // Keep track of what input we are on. inputcount++; } // End of Event Loop =============================== // Now event loop is finished loop over all Measurements // Converting Binned events to XSec Distributions iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ConvertEventRates(); } // Print out statements on approximate memory usage for profiling. LOG(REC) << "Filled " << fillcount << " signal events." << std::endl; if (savesignal) { int mem = - ( // sizeof(fSignalEventBoxes) + - // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) + - sizeof(MeasurementVariableBox1D) * fillcount) * - 1E-6; + ( // sizeof(fSignalEventBoxes) + + // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) + + sizeof(MeasurementVariableBox1D) * fillcount) * + 1E-6; LOG(REC) << " -> Saved " << fillcount << " signal boxes for faster access. (~" << mem << " MB)" << std::endl; if (fIsAllSplines and !fSignalEventSplines.empty()) { int splmem = sizeof(float) * fSignalEventSplines.size() * fSignalEventSplines[0].size() * 1E-6; LOG(REC) << " -> Saved " << fillcount << " " << fSignalEventSplines.size() << " spline sets into memory. (~" << splmem << " MB)" << std::endl; } } LOG(REC) << "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart << std::endl; // Check SignalReconfigures works for all samples if (savesignal) { double likefull = GetLikelihood(); ReconfigureFastUsingManager(); double likefast = GetLikelihood(); - if (fabs(likefull - likefast) > 0.0001) - { - ERROR(FTL, "Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast); - ERROR(FTL, "This means some samples you are using are not setup to use SignalReconfigures=1"); + if (fabs(likefull - likefast) > 0.0001) { + ERROR(FTL, "Fast and Full Likelihoods DIFFER! : " << likefull << " : " + << likefast); + ERROR(FTL, + "This means some samples you are using are not setup to use " + "SignalReconfigures=1"); ERROR(FTL, "Please turn OFF signal reconfigures."); throw; } else { - LOG(FIT) << "Likelihoods for FULL and FAST match. Will use FAST next time." << std::endl; + LOG(FIT) + << "Likelihoods for FULL and FAST match. Will use FAST next time." + << std::endl; } } // End of reconfigure return; }; //*************************************************** void JointFCN::ReconfigureFastUsingManager() { -//*************************************************** + //*************************************************** LOG(FIT) << " -> Doing FAST using manager" << std::endl; // Get Start time for profilling int timestart = time(NULL); // Reset all samples MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ResetAll(); } // Check for saved variables if not do a full reconfigure. if (fSignalEventFlags.empty()) { ERR(WRN) << "Signal Flags Empty! Using normal manager." << std::endl; ReconfigureUsingManager(); return; } bool fFillNuisanceEvent = - FitPar::Config().GetParB("FullEventOnSignalReconfigure"); + FitPar::Config().GetParB("FullEventOnSignalReconfigure"); // Setup fast vector iterators. std::vector::iterator inpsig_iter = fSignalEventFlags.begin(); std::vector >::iterator box_iter = - fSignalEventBoxes.begin(); + fSignalEventBoxes.begin(); std::vector >::iterator spline_iter = - fSignalEventSplines.begin(); + fSignalEventSplines.begin(); std::vector >::iterator samsig_iter = - fSampleSignalFlags.begin(); + fSampleSignalFlags.begin(); int splinecount = 0; // Setup stuff for logging int fillcount = 0; int nevents = fSignalEventFlags.size(); int countwidth = nevents / 20; // If All Splines tell splines they need a reconfigure. std::vector::iterator inp_iter = fInputList.begin(); if (fIsAllSplines) { LOG(REC) << "All Spline Inputs so using fast spline loop." << std::endl; for (; inp_iter != fInputList.end(); inp_iter++) { InputHandlerBase* curinput = (*inp_iter); // Tell each fSplineRead in BaseFitEvent to reconf next weight calc BaseFitEvt* curevent = curinput->FirstBaseEvent(); if (curevent->fSplineRead) curevent->fSplineRead->SetNeedsReconfigure(true); } } // Loop over all possible spline inputs double* coreeventweights = new double[fSignalEventBoxes.size()]; splinecount = 0; inp_iter = fInputList.begin(); inpsig_iter = fSignalEventFlags.begin(); spline_iter = fSignalEventSplines.begin(); // Loop over all signal flags // For each valid signal flag add one to splinecount // Get Splines from that count and add to weight // Add splinecount int sigcount = 0; splinecount = 0; // #pragma omp parallel for shared(splinecount,sigcount) for (uint iinput = 0; iinput < fInputList.size(); iinput++) { InputHandlerBase* curinput = fInputList[iinput]; BaseFitEvt* curevent = curinput->FirstBaseEvent(); for (int i = 0; i < curinput->GetNEvents(); i++) { double rwweight = 0.0; if (fSignalEventFlags[sigcount]) { // Get Event Info if (!fIsAllSplines) { if (fFillNuisanceEvent) curinput->GetNuisanceEvent(i); else curevent = curinput->GetBaseEvent(i); } else { curevent->fSplineCoeff = &fSignalEventSplines[splinecount][0]; } curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent); curevent->Weight = curevent->RWWeight * curevent->InputWeight; rwweight = curevent->Weight; coreeventweights[splinecount] = rwweight; if (countwidth && ((splinecount % countwidth) == 0)) { LOG(REC) << "Processed " << splinecount << " event weights. W = " << rwweight << std::endl; } // #pragma omp atomic splinecount++; } // #pragma omp atomic sigcount++; } } LOG(SAM) << "Processed event weights." << std::endl; // #pragma omp barrier // Reset Iterators inpsig_iter = fSignalEventFlags.begin(); spline_iter = fSignalEventSplines.begin(); box_iter = fSignalEventBoxes.begin(); samsig_iter = fSampleSignalFlags.begin(); int nsplineweights = splinecount; splinecount = 0; // Start of Fast Event Loop ============================ // Start input iterators // Loop over number of inputs for (int ispline = 0; ispline < nsplineweights; ispline++) { double rwweight = coreeventweights[ispline]; // Get iterators for this event std::vector::iterator subsamsig_iter = (*samsig_iter).begin(); std::vector::iterator subbox_iter = - (*box_iter).begin(); + (*box_iter).begin(); // Loop over all sub measurements. std::vector::iterator meas_iter = fSubSampleList.begin(); for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) { MeasurementBase* curmeas = (*meas_iter); // If event flagged as signal for this sample fill from the box. if (*subsamsig_iter) { curmeas->SetSignal(true); curmeas->FillHistogramsFromBox((*subbox_iter), rwweight); // Move onto next box if there is one. subbox_iter++; fillcount++; } } if (ispline % countwidth == 0) { LOG(REC) << "Filled " << ispline << " sample weights." << std::endl; } // Iterate over the main signal event containers. samsig_iter++; box_iter++; spline_iter++; splinecount++; } // End of Fast Event Loop =================== LOG(SAM) << "Filled sample distributions." << std::endl; // Now loop over all Measurements // Convert Binned events iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ConvertEventRates(); } // Cleanup coreeventweights if (fIsAllSplines) { delete coreeventweights; } // Print some reconfigure profiling. LOG(REC) << "Filled " << fillcount << " signal events." << std::endl; LOG(REC) << "Time taken ReconfigureFastUsingManager() : " << time(NULL) - timestart << std::endl; } //*************************************************** void JointFCN::Write() { -//*************************************************** + //*************************************************** // Save a likelihood/ndof plot LOG(MIN) << "Writing likelihood plot.." << std::endl; std::vector likes; std::vector ndofs; std::vector names; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; double like = exp->GetLikelihood(); double ndof = exp->GetNDOF(); std::string name = exp->GetName(); likes.push_back(like); ndofs.push_back(ndof); names.push_back(name); } - TH1D likehist = TH1D("likelihood_hist", "likelihood_hist", - likes.size(), 0.0, double(likes.size())); - TH1D ndofhist = TH1D("ndof_hist", "ndof_hist", - ndofs.size(), 0.0, double(ndofs.size())); - TH1D divhist = TH1D("likedivndof_hist", "likedivndof_hist", - likes.size(), 0.0, double(likes.size())); + TH1D likehist = TH1D("likelihood_hist", "likelihood_hist", likes.size(), 0.0, + double(likes.size())); + TH1D ndofhist = + TH1D("ndof_hist", "ndof_hist", ndofs.size(), 0.0, double(ndofs.size())); + TH1D divhist = TH1D("likedivndof_hist", "likedivndof_hist", likes.size(), 0.0, + double(likes.size())); for (size_t i = 0; i < likehist.GetNbinsX(); i++) { likehist.SetBinContent(i + 1, likes[i]); ndofhist.SetBinContent(i + 1, ndofs[i]); if (ndofs[i] != 0.0) { divhist.SetBinContent(i + 1, likes[i] / ndofs[i]); } likehist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); ndofhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); divhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); } likehist.Write(); ndofhist.Write(); divhist.Write(); // Loop over individual experiments and call Write LOG(MIN) << "Writing each of the data classes..." << std::endl; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; exp->Write(); } // Save Pull Terms for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; pull->Write(); } if (FitPar::Config().GetParB("EventManager")) { // Get list of inputs std::map fInputs = - FitBase::EvtManager().GetInputs(); + FitBase::EvtManager().GetInputs(); std::map::const_iterator iterInp; for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) { InputHandlerBase* input = (iterInp->second); input->GetFluxHistogram()->Write(); input->GetXSecHistogram()->Write(); input->GetEventHistogram()->Write(); } } }; //*************************************************** void JointFCN::SetFakeData(std::string fakeinput) { -//*************************************************** + //*************************************************** LOG(MIN) << "Setting fake data from " << fakeinput << std::endl; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; exp->SetFakeDataValues(fakeinput); } return; } //*************************************************** void JointFCN::ThrowDataToy() { -//*************************************************** + //*************************************************** for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; exp->ThrowDataToy(); } return; } //*************************************************** 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(); + 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++) { + 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; LOG(MIN) << "Likelihoods : " << std::endl; // Loop over samples first - for (MeasListConstIter iter = fSamples.begin(); - iter != fSamples.end(); + 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 LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : " << singlelike << std::endl; } - // Loop over pulls second - for (PullListConstIter iter = fPulls.begin(); - iter != fPulls.end(); - iter++) { + 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 LOG(MIN) << "-> " << std::left << std::setw(40) << pull->GetName() << " : " << singlelike << std::endl; - } // 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(); + 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++) { + 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/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx index 5f88e1d..3306167 100644 --- a/src/FitBase/Measurement1D.cxx +++ b/src/FitBase/Measurement1D.cxx @@ -1,1903 +1,1905 @@ // Copyright 2016 L. Pickering, P. Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This ile 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 "Measurement1D.h" //******************************************************************** Measurement1D::Measurement1D(void) { //******************************************************************** // XSec Scalings fScaleFactor = -1.0; fCurrentNorm = 1.0; // Histograms fDataHist = NULL; fDataTrue = NULL; fMCHist = NULL; fMCFine = NULL; fMCWeighted = NULL; fMaskHist = NULL; // Covar covar = NULL; fFullCovar = NULL; fShapeCovar = NULL; fCovar = NULL; fInvert = NULL; fDecomp = NULL; // Fake Data fFakeDataInput = ""; fFakeDataFile = NULL; // Options fDefaultTypes = "FIX/FULL/CHI2"; fAllowedTypes = "FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK/NOWIDTH"; fIsFix = false; fIsShape = false; fIsFree = false; fIsDiag = false; fIsFull = false; fAddNormPen = false; fIsMask = false; fIsChi2SVD = false; fIsRawEvents = false; fIsNoWidth = false; fIsDifXSec = false; fIsEnu1D = false; // Inputs fInput = NULL; fRW = NULL; // Extra Histograms fMCHist_Modes = NULL; } //******************************************************************** Measurement1D::~Measurement1D(void) { //******************************************************************** if (fDataHist) delete fDataHist; if (fDataTrue) delete fDataTrue; if (fMCHist) delete fMCHist; if (fMCFine) delete fMCFine; if (fMCWeighted) delete fMCWeighted; if (fMaskHist) delete fMaskHist; if (covar) delete covar; if (fFullCovar) delete fFullCovar; if (fShapeCovar) delete fShapeCovar; if (fCovar) delete fCovar; if (fInvert) delete fInvert; if (fDecomp) delete fDecomp; } //******************************************************************** void Measurement1D::FinaliseSampleSettings() { //******************************************************************** MeasurementBase::FinaliseSampleSettings(); // Setup naming + renaming fName = fSettings.GetName(); fSettings.SetS("originalname", fName); if (fSettings.Has("rename")) { fName = fSettings.GetS("rename"); fSettings.SetS("name", fName); } // Setup all other options LOG(SAM) << "Finalising Sample Settings: " << fName << std::endl; if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) { fIsRawEvents = true; LOG(SAM) << "Found event rate measurement but using poisson likelihoods." << std::endl; } if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) { fIsEnu1D = true; LOG(SAM) << "::" << fName << "::" << std::endl; LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, " << "not flux averaged!" << std::endl; } if (fIsEnu1D && fIsRawEvents) { LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!" << std::endl; LOG(SAM) << "Check experiment constructor for " << fName << " and correct this!" << std::endl; LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl; exit(-1); } if (!fRW) fRW = FitBase::GetRW(); if (!fInput and !fIsJoint) SetupInputs(fSettings.GetS("input")); // Setup options SetFitOptions(fDefaultTypes); // defaults SetFitOptions(fSettings.GetS("type")); // user specified EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min")); EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max")); if (fAddNormPen) { if (fNormError <= 0.0) { ERR(WRN) << "Norm error for class " << fName << " is 0.0!" << std::endl; ERR(WRN) << "If you want to use it please add fNormError=VAL" << std::endl; throw; } } } //******************************************************************** void Measurement1D::CreateDataHistogram(int dimx, double* binx) { //******************************************************************** if (fDataHist) delete fDataHist; fDataHist = new TH1D( (fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str(), dimx, binx) ; } //******************************************************************** void Measurement1D::SetDataFromTextFile(std::string datafile) { //******************************************************************** LOG(SAM) << "Reading data from text file: " << datafile << std::endl; fDataHist = PlotUtils::GetTH1DFromFile(datafile, fSettings.GetName() + "_data", fSettings.GetFullTitles()); } //******************************************************************** void Measurement1D::SetDataFromRootFile(std::string datafile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading data from root file: " << datafile << ";" << histname << std::endl; fDataHist = PlotUtils::GetTH1DFromRootFile(datafile, histname); fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(), (fSettings.GetFullTitles()).c_str()); return; }; //******************************************************************** void Measurement1D::SetEmptyData(){ //******************************************************************** fDataHist = new TH1D("EMPTY_DATA","EMPTY_DATA",1,0.0,1.0); } //******************************************************************** void Measurement1D::SetPoissonErrors() { //******************************************************************** if (!fDataHist) { ERR(FTL) << "Need a data hist to setup possion errors! " << std::endl; ERR(FTL) << "Setup Data First!" << std::endl; throw; } for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) { fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1))); } } //******************************************************************** void Measurement1D::SetCovarFromDiagonal(TH1D* data) { //******************************************************************** if (!data and fDataHist) { data = fDataHist; } if (data) { LOG(SAM) << "Setting diagonal covariance for: " << data->GetName() << std::endl; fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } else { ERR(FTL) << "No data input provided to set diagonal covar from!" << std::endl; } // if (!fIsDiag) { // ERR(FTL) << "SetCovarMatrixFromDiag called for measurement " // << "that is not set as diagonal." << std::endl; // throw; // } } //******************************************************************** void Measurement1D::SetCovarFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } LOG(SAM) << "Reading covariance from text file: " << covfile << std::endl; fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarFromMultipleTextFiles(std::string covfiles, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } std::vector covList = GeneralUtils::ParseToStr(covfiles, ";"); fFullCovar = new TMatrixDSym(dim); for (uint i = 0; i < covList.size(); ++i){ LOG(SAM) << "Reading covariance from text file: " << covList[i] << std::endl; TMatrixDSym* temp_cov = StatUtils::GetCovarFromTextFile(covList[i], dim); (*fFullCovar) += (*temp_cov); delete temp_cov; } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarFromRootFile(std::string covfile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading covariance from text file: " << covfile << ";" << histname << std::endl; fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarInvertFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } LOG(SAM) << "Reading inverted covariance from text file: " << covfile << std::endl; covar = StatUtils::GetCovarFromTextFile(covfile, dim); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCovarInvertFromRootFile(std::string covfile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading inverted covariance from text file: " << covfile << ";" << histname << std::endl; covar = StatUtils::GetCovarFromRootFile(covfile, histname); fFullCovar = StatUtils::GetInvert(covar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCorrelationFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) dim = fDataHist->GetNbinsX(); LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << dim << std::endl; TMatrixDSym* correlation = StatUtils::GetCovarFromTextFile(covfile, dim); if (!fDataHist) { ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n" << "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl; throw; } // Fill covar from data errors and correlations fFullCovar = new TMatrixDSym(dim); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } // Fill other covars. covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete correlation; } //******************************************************************** void Measurement1D::SetCorrelationFromMultipleTextFiles(std::string corrfiles, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } std::vector corrList = GeneralUtils::ParseToStr(corrfiles, ";"); fFullCovar = new TMatrixDSym(dim); for (uint i = 0; i < corrList.size(); ++i){ LOG(SAM) << "Reading covariance from text file: " << corrList[i] << std::endl; TMatrixDSym* temp_cov = StatUtils::GetCovarFromTextFile(corrList[i], dim); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*temp_cov)(i, j) = (*temp_cov)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } (*fFullCovar) += (*temp_cov); delete temp_cov; } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void Measurement1D::SetCorrelationFromRootFile(std::string covfile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << histname << std::endl; TMatrixDSym* correlation = StatUtils::GetCovarFromRootFile(covfile, histname); if (!fDataHist) { ERR(FTL) << "Trying to set correlations from text file but there is no data to build it from. \n" << "In constructor make sure data is set before SetCorrelationFromTextFile is called. \n" << std::endl; throw; } // Fill covar from data errors and correlations fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX()); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { for (int j = 0; j < fDataHist->GetNbinsX(); j++) { (*fFullCovar)(i, j) = (*correlation)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1.E76; } } // Fill other covars. covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete correlation; } //******************************************************************** void Measurement1D::SetCholDecompFromTextFile(std::string covfile, int dim) { //******************************************************************** if (dim == -1) { dim = fDataHist->GetNbinsX(); } LOG(SAM) << "Reading cholesky from text file: " << covfile << std::endl; TMatrixD* temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim); TMatrixD* trans = (TMatrixD*)temp->Clone(); trans->T(); (*trans) *= (*temp); fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), ""); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } //******************************************************************** void Measurement1D::SetCholDecompFromRootFile(std::string covfile, std::string histname) { //******************************************************************** LOG(SAM) << "Reading cholesky decomp from root file: " << covfile << ";" << histname << std::endl; TMatrixD* temp = StatUtils::GetMatrixFromRootFile(covfile, histname); TMatrixD* trans = (TMatrixD*)temp->Clone(); trans->T(); (*trans) *= (*temp); fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), ""); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); delete temp; delete trans; } void Measurement1D::SetShapeCovar(){ // Return if this is missing any pre-requisites if (!fFullCovar) return; if (!fDataHist) return; // Also return if it's bloody stupid under the circumstances if (fIsDiag) return; fShapeCovar = StatUtils::ExtractShapeOnlyCovar(fFullCovar, fDataHist); return; } //******************************************************************** void Measurement1D::ScaleData(double scale) { //******************************************************************** fDataHist->Scale(scale); } //******************************************************************** void Measurement1D::ScaleDataErrors(double scale) { //******************************************************************** for (int i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinError(i + 1, fDataHist->GetBinError(i + 1) * scale); } } //******************************************************************** void Measurement1D::ScaleCovar(double scale) { //******************************************************************** (*fFullCovar) *= scale; (*covar) *= 1.0 / scale; (*fDecomp) *= sqrt(scale); } //******************************************************************** void Measurement1D::SetBinMask(std::string maskfile) { //******************************************************************** if (!fIsMask) return; LOG(SAM) << "Reading bin mask from file: " << maskfile << std::endl; // Create a mask histogram with dim of data int nbins = fDataHist->GetNbinsX(); fMaskHist = new TH1I((fSettings.GetName() + "_BINMASK").c_str(), (fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(), nbins, 0, nbins); std::string line; std::ifstream mask(maskfile.c_str(), std::ifstream::in); if (!mask.is_open()) { LOG(FTL) << " Cannot find mask file." << std::endl; throw; } while (std::getline(mask >> std::ws, line, '\n')) { std::vector entries = GeneralUtils::ParseToInt(line, " "); // Skip lines with poorly formatted lines if (entries.size() < 2) { LOG(WRN) << "Measurement1D::SetBinMask(), couldn't parse line: " << line << std::endl; continue; } // The first index should be the bin number, the second should be the mask // value. int val = 0; if (entries[1] > 0) val = 1; fMaskHist->SetBinContent(entries[0], val); } // Apply masking by setting masked data bins to zero PlotUtils::MaskBins(fDataHist, fMaskHist); return; } //******************************************************************** void Measurement1D::FinaliseMeasurement() { //******************************************************************** LOG(SAM) << "Finalising Measurement: " << fName << std::endl; if (fSettings.GetB("onlymc")){ if (fDataHist) delete fDataHist; fDataHist = new TH1D("empty_data","empty_data",1,0.0,1.0); } // Make sure data is setup if (!fDataHist) { ERR(FTL) << "No data has been setup inside " << fName << " constructor!" << std::endl; throw; } // Make sure covariances are setup if (!fFullCovar) { fIsDiag = true; SetCovarFromDiagonal(fDataHist); } if (!covar) { covar = StatUtils::GetInvert(fFullCovar); } if (!fDecomp) { fDecomp = StatUtils::GetDecomp(fFullCovar); } // Push the diagonals of fFullCovar onto the data histogram // Comment this out until the covariance/data scaling is consistent! StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1E-38); // If shape only, set covar and fDecomp using the shape-only matrix (if set) if (fIsShape && fShapeCovar and FitPar::Config().GetParB("UseShapeCovar")){ if (covar) delete covar; covar = StatUtils::GetInvert(fShapeCovar); if (fDecomp) delete fDecomp; fDecomp = StatUtils::GetDecomp(fFullCovar); } // Setup fMCHist from data fMCHist = (TH1D*)fDataHist->Clone(); fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(), (fSettings.GetFullTitles()).c_str()); fMCHist->Reset(); // Setup fMCFine fMCFine = new TH1D("mcfine", "mcfine", fDataHist->GetNbinsX() * 8, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(fDataHist->GetNbinsX() + 1)); fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(), (fSettings.GetFullTitles()).c_str()); fMCFine->Reset(); // Setup MC Stat fMCStat = (TH1D*)fMCHist->Clone(); fMCStat->Reset(); // Search drawopts for possible types to include by default std::string drawopts = FitPar::Config().GetParS("drawopts"); if (drawopts.find("MODES") != std::string::npos) { fMCHist_Modes = new TrueModeStack( (fSettings.GetName() + "_MODES").c_str(), ("True Channels"), fMCHist); SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write); } // Setup bin masks using sample name if (fIsMask) { std::string curname = fName; std::string origname = fSettings.GetS("originalname"); // Check rename.mask std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask"); // Check origname.mask if (maskloc.empty()) maskloc = FitPar::Config().GetParDIR(origname + ".mask"); // Check database if (maskloc.empty()) { maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask"; } // Setup Bin Mask SetBinMask(maskloc); } if (fScaleFactor < 0) { ERR(FTL) << "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__ << std::endl; ERR(FTL) << "fScaleFactor = " << fScaleFactor << std::endl; ERR(FTL) << "EXITING" << std::endl; throw; } // Create and fill Weighted Histogram if (!fMCWeighted) { fMCWeighted = (TH1D*)fMCHist->Clone(); fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(), (fName + "_MCWGHTS" + fPlotTitles).c_str()); fMCWeighted->GetYaxis()->SetTitle("Weighted Events"); } } //******************************************************************** void Measurement1D::SetFitOptions(std::string opt) { //******************************************************************** // Do nothing if default given if (opt == "DEFAULT") return; // CHECK Conflicting Fit Options std::vector fit_option_allow = GeneralUtils::ParseToStr(fAllowedTypes, "/"); for (UInt_t i = 0; i < fit_option_allow.size(); i++) { std::vector fit_option_section = GeneralUtils::ParseToStr(fit_option_allow.at(i), ","); bool found_option = false; for (UInt_t j = 0; j < fit_option_section.size(); j++) { std::string av_opt = fit_option_section.at(j); if (!found_option and opt.find(av_opt) != std::string::npos) { found_option = true; } else if (found_option and opt.find(av_opt) != std::string::npos) { ERR(FTL) << "ERROR: Conflicting fit options provided: " << opt << std::endl << "Conflicting group = " << fit_option_section.at(i) << std::endl << "You should only supply one of these options in card file." << std::endl; throw; } } } // Check all options are allowed std::vector fit_options_input = GeneralUtils::ParseToStr(opt, "/"); for (UInt_t i = 0; i < fit_options_input.size(); i++) { if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) { ERR(FTL) << "ERROR: Fit Option '" << fit_options_input.at(i) << "' Provided is not allowed for this measurement." << std::endl; ERR(FTL) << "Fit Options should be provided as a '/' seperated list " "(e.g. FREE/DIAG/NORM)" << std::endl; ERR(FTL) << "Available options for " << fName << " are '" << fAllowedTypes << "'" << std::endl; throw; } } // Set TYPE fFitType = opt; // FIX,SHAPE,FREE if (opt.find("FIX") != std::string::npos) { fIsFree = fIsShape = false; fIsFix = true; } else if (opt.find("SHAPE") != std::string::npos) { fIsFree = fIsFix = false; fIsShape = true; } else if (opt.find("FREE") != std::string::npos) { fIsFix = fIsShape = false; fIsFree = true; } // DIAG,FULL (or default to full) if (opt.find("DIAG") != std::string::npos) { fIsDiag = true; fIsFull = false; } else if (opt.find("FULL") != std::string::npos) { fIsDiag = false; fIsFull = true; } // CHI2/LL (OTHERS?) if (opt.find("LOG") != std::string::npos) { fIsChi2 = false; ERR(FTL) << "No other LIKELIHOODS properly supported!" << std::endl; ERR(FTL) << "Try to use a chi2!" << std::endl; throw; } else { fIsChi2 = true; } // EXTRAS if (opt.find("RAW") != std::string::npos) fIsRawEvents = true; if (opt.find("NOWIDTH") != std::string::npos) fIsNoWidth = true; if (opt.find("DIF") != std::string::npos) fIsDifXSec = true; if (opt.find("ENU1D") != std::string::npos) fIsEnu1D = true; if (opt.find("NORM") != std::string::npos) fAddNormPen = true; if (opt.find("MASK") != std::string::npos) fIsMask = true; return; }; //******************************************************************** void Measurement1D::SetSmearingMatrix(std::string smearfile, int truedim, int recodim) { //******************************************************************** // The smearing matrix describes the migration from true bins (rows) to reco // bins (columns) // Counter over the true bins! int row = 0; std::string line; std::ifstream smear(smearfile.c_str(), std::ifstream::in); // Note that the smearing matrix may be rectangular. fSmearMatrix = new TMatrixD(truedim, recodim); if (smear.is_open()) LOG(SAM) << "Reading smearing matrix from file: " << smearfile << std::endl; else ERR(FTL) << "Smearing matrix provided is incorrect: " << smearfile << std::endl; while (std::getline(smear >> std::ws, line, '\n')) { int column = 0; std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { (*fSmearMatrix)(row, column) = (*iter) / 100.; // Convert to fraction from // percentage (this may not be // general enough) column++; } row++; } return; } //******************************************************************** void Measurement1D::ApplySmearingMatrix() { //******************************************************************** if (!fSmearMatrix) { ERR(WRN) << fName << ": attempted to apply smearing matrix, but none was set" << std::endl; return; } TH1D* unsmeared = (TH1D*)fMCHist->Clone(); TH1D* smeared = (TH1D*)fMCHist->Clone(); smeared->Reset(); // Loop over reconstructed bins // true = row; reco = column for (int rbin = 0; rbin < fSmearMatrix->GetNcols(); ++rbin) { // Sum up the constributions from all true bins double rBinVal = 0; // Loop over true bins for (int tbin = 0; tbin < fSmearMatrix->GetNrows(); ++tbin) { rBinVal += (*fSmearMatrix)(tbin, rbin) * unsmeared->GetBinContent(tbin + 1); } smeared->SetBinContent(rbin + 1, rBinVal); } fMCHist = (TH1D*)smeared->Clone(); return; } /* Reconfigure LOOP */ //******************************************************************** void Measurement1D::ResetAll() { //******************************************************************** fMCHist->Reset(); fMCFine->Reset(); fMCStat->Reset(); return; }; //******************************************************************** void Measurement1D::FillHistograms() { //******************************************************************** if (Signal) { + QLOG(DEB, "Fill MCHist: " << fXVar << ", " << Weight); + fMCHist->Fill(fXVar, Weight); fMCFine->Fill(fXVar, Weight); fMCStat->Fill(fXVar, 1.0); if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, Weight); } return; }; //******************************************************************** void Measurement1D::ScaleEvents() { //******************************************************************** // Fill MCWeighted; // for (int i = 0; i < fMCHist->GetNbinsX(); i++) { // fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1)); // fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1)); // } // Setup Stat ratios for MC and MC Fine double* statratio = new double[fMCHist->GetNbinsX()]; for (int i = 0; i < fMCHist->GetNbinsX(); i++) { if (fMCHist->GetBinContent(i + 1) != 0) { statratio[i] = fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1); } else { statratio[i] = 0.0; } } double* statratiofine = new double[fMCFine->GetNbinsX()]; for (int i = 0; i < fMCFine->GetNbinsX(); i++) { if (fMCFine->GetBinContent(i + 1) != 0) { statratiofine[i] = fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1); } else { statratiofine[i] = 0.0; } } // Scaling for raw event rates if (fIsRawEvents) { double datamcratio = fDataHist->Integral() / fMCHist->Integral(); fMCHist->Scale(datamcratio); fMCFine->Scale(datamcratio); if (fMCHist_Modes) fMCHist_Modes->Scale(datamcratio); // Scaling for XSec as function of Enu } else if (fIsEnu1D) { PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(), GetEventHistogram(), fScaleFactor, fNEvents); PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(), GetEventHistogram(), fScaleFactor, fNEvents); // if (fMCHist_Modes) { // PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(), // GetEventHistogram(), fScaleFactor, // fNEvents); // } } else if (fIsNoWidth) { fMCHist->Scale(fScaleFactor); fMCFine->Scale(fScaleFactor); if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor); // Any other differential scaling } else { fMCHist->Scale(fScaleFactor, "width"); fMCFine->Scale(fScaleFactor, "width"); if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width"); } // Proper error scaling - ROOT Freaks out with xsec weights sometimes for (int i = 0; i < fMCStat->GetNbinsX(); i++) { fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]); } for (int i = 0; i < fMCFine->GetNbinsX(); i++) { fMCFine->SetBinError(i + 1, fMCFine->GetBinContent(i + 1) * statratiofine[i]); } // Clean up delete statratio; delete statratiofine; return; }; //******************************************************************** void Measurement1D::ApplyNormScale(double norm) { //******************************************************************** fCurrentNorm = norm; fMCHist->Scale(1.0 / norm); fMCFine->Scale(1.0 / norm); return; }; /* Statistic Functions - Outsources to StatUtils */ //******************************************************************** int Measurement1D::GetNDOF() { //******************************************************************** int ndof = fDataHist->GetNbinsX(); if (fMaskHist and fIsMask) ndof -= fMaskHist->Integral(); return ndof; } //******************************************************************** double Measurement1D::GetLikelihood() { //******************************************************************** // If this is for a ratio, there is no data histogram to compare to! if (fNoData || !fDataHist) return 0.; // Apply Masking to MC if Required. if (fIsMask and fMaskHist) { PlotUtils::MaskBins(fMCHist, fMaskHist); } // Sort Shape Scaling double scaleF = 0.0; // TODO Include !fIsRawEvents if (fIsShape) { if (fMCHist->Integral(1, fMCHist->GetNbinsX(), "width")) { scaleF = fDataHist->Integral(1, fDataHist->GetNbinsX(), "width") / fMCHist->Integral(1, fMCHist->GetNbinsX(), "width"); fMCHist->Scale(scaleF); fMCFine->Scale(scaleF); } } // Likelihood Calculation double stat = 0.; if (fIsChi2) { if (fIsRawEvents) { stat = StatUtils::GetChi2FromEventRate(fDataHist, fMCHist, fMaskHist); } else if (fIsDiag) { stat = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMaskHist); } else if (!fIsDiag and !fIsRawEvents) { stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist); } } // Sort Penalty Terms if (fAddNormPen) { double penalty = (1. - fCurrentNorm) * (1. - fCurrentNorm) / (fNormError * fNormError); stat += penalty; } // Return to normal scaling if (fIsShape) { // and !FitPar::Config().GetParB("saveshapescaling")) { fMCHist->Scale(1. / scaleF); fMCFine->Scale(1. / scaleF); } fLikelihood = stat; return stat; } /* Fake Data Functions */ //******************************************************************** void Measurement1D::SetFakeDataValues(std::string fakeOption) { //******************************************************************** // Setup original/datatrue TH1D* tempdata = (TH1D*) fDataHist->Clone(); if (!fIsFakeData) { fIsFakeData = true; // Make a copy of the original data histogram. if (!fDataOrig) fDataOrig = (TH1D*)fDataHist->Clone((fName + "_data_original").c_str()); } else { ResetFakeData(); } // Setup Inputs fFakeDataInput = fakeOption; LOG(SAM) << "Setting fake data from : " << fFakeDataInput << std::endl; // From MC if (fFakeDataInput.compare("MC") == 0) { fDataHist = (TH1D*)fMCHist->Clone((fName + "_MC").c_str()); // Fake File } else { if (!fFakeDataFile) fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ"); fDataHist = (TH1D*)fFakeDataFile->Get((fName + "_MC").c_str()); } // Setup Data Hist fDataHist->SetNameTitle((fName + "_FAKE").c_str(), (fName + fPlotTitles).c_str()); // Replace Data True if (fDataTrue) delete fDataTrue; fDataTrue = (TH1D*)fDataHist->Clone(); fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(), (fName + fPlotTitles).c_str()); // Make a new covariance for fake data hist. int nbins = fDataHist->GetNbinsX(); double alpha_i = 0.0; double alpha_j = 0.0; for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { alpha_i = fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1); alpha_j = fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1); (*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j); } } // Setup Covariances if (covar) delete covar; covar = StatUtils::GetInvert(fFullCovar); if (fDecomp) delete fDecomp; fDecomp = StatUtils::GetInvert(fFullCovar); delete tempdata; return; }; //******************************************************************** void Measurement1D::ResetFakeData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH1D*)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str()); } } //******************************************************************** void Measurement1D::ResetData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH1D*)fDataOrig->Clone((fSettings.GetName() + "_data").c_str()); } fIsFakeData = false; } //******************************************************************** void Measurement1D::ThrowCovariance() { //******************************************************************** // Take a fDecomposition and use it to throw the current dataset. // Requires fDataTrue also be set incase used repeatedly. if (!fDataTrue) fDataTrue = (TH1D*) fDataHist->Clone(); if (fDataHist) delete fDataHist; fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); return; }; //******************************************************************** void Measurement1D::ThrowDataToy(){ //******************************************************************** if (!fDataTrue) fDataTrue = (TH1D*) fDataHist->Clone(); if (fMCHist) delete fMCHist; fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); } /* Access Functions */ //******************************************************************** TH1D* Measurement1D::GetMCHistogram() { //******************************************************************** if (!fMCHist) return fMCHist; std::ostringstream chi2; chi2 << std::setprecision(5) << this->GetLikelihood(); int linecolor = kRed; int linestyle = 1; int linewidth = 1; int fillcolor = 0; int fillstyle = 1001; // if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor"); // if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle"); // if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth"); // if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor"); // if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle"); fMCHist->SetTitle(chi2.str().c_str()); fMCHist->SetLineColor(linecolor); fMCHist->SetLineStyle(linestyle); fMCHist->SetLineWidth(linewidth); fMCHist->SetFillColor(fillcolor); fMCHist->SetFillStyle(fillstyle); return fMCHist; }; //******************************************************************** TH1D* Measurement1D::GetDataHistogram() { //******************************************************************** if (!fDataHist) return fDataHist; int datacolor = kBlack; int datastyle = 1; int datawidth = 1; // if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor"); // if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle"); // if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth"); fDataHist->SetLineColor(datacolor); fDataHist->SetLineWidth(datawidth); fDataHist->SetMarkerStyle(datastyle); return fDataHist; }; /* Write Functions */ // Save all the histograms at once //******************************************************************** void Measurement1D::Write(std::string drawOpt) { //******************************************************************** // Get Draw Options drawOpt = FitPar::Config().GetParS("drawopts"); // Write Settigns if (drawOpt.find("SETTINGS") != std::string::npos){ fSettings.Set("#chi^{2}",fLikelihood); fSettings.Set("NDOF", this->GetNDOF() ); fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF() ); fSettings.Write(); } // Write Data/MC GetDataList().at(0)->Write(); GetMCList().at(0)->Write(); if((fEvtRateScaleFactor != 0xdeadbeef) && GetMCList().at(0)){ TH1D * PredictedEvtRate = static_cast(GetMCList().at(0)->Clone()); PredictedEvtRate->Scale(fEvtRateScaleFactor); PredictedEvtRate->GetYaxis()->SetTitle("Predicted event rate"); PredictedEvtRate->Write(); } // Write Fine Histogram if (drawOpt.find("FINE") != std::string::npos) GetFineList().at(0)->Write(); // Write Weighted Histogram if (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted) fMCWeighted->Write(); // Save Flux/Evt if no event manager if (!FitPar::Config().GetParB("EventManager")) { if (drawOpt.find("FLUX") != std::string::npos && GetFluxHistogram()) GetFluxHistogram()->Write(); if (drawOpt.find("EVT") != std::string::npos && GetEventHistogram()) GetEventHistogram()->Write(); if (drawOpt.find("XSEC") != std::string::npos && GetEventHistogram()) GetXSecHistogram()->Write(); } // Write Mask if (fIsMask && (drawOpt.find("MASK") != std::string::npos)) { fMaskHist->Write(); } // Write Covariances if (drawOpt.find("COV") != std::string::npos && fFullCovar) { PlotUtils::GetFullCovarPlot(fFullCovar, fSettings.GetName()); } if (drawOpt.find("INVCOV") != std::string::npos && covar) { PlotUtils::GetInvCovarPlot(covar, fSettings.GetName()); } if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) { PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName()); } // // Likelihood residual plots // if (drawOpt.find("RESIDUAL") != std::string::npos) { // WriteResidualPlots(); // } // Ratio and Shape Plots if (drawOpt.find("RATIO") != std::string::npos) { WriteRatioPlot(); } if (drawOpt.find("SHAPE") != std::string::npos) { WriteShapePlot(); if (drawOpt.find("RATIO") != std::string::npos) WriteShapeRatioPlot(); } // // RATIO // if (drawOpt.find("CANVMC") != std::string::npos) { // TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist); // c1->Write(); // delete c1; // } // // PDG // if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) { // TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes); // c2->Write(); // delete c2; // } // Write Extra Histograms AutoWriteExtraTH1(); WriteExtraHistograms(); // Returning LOG(SAM) << "Written Histograms: " << fName << std::endl; return; } //******************************************************************** void Measurement1D::WriteRatioPlot() { //******************************************************************** // Setup mc data ratios TH1D* dataRatio = (TH1D*)fDataHist->Clone((fName + "_data_RATIO").c_str()); TH1D* mcRatio = (TH1D*)fMCHist->Clone((fName + "_MC_RATIO").c_str()); // Extra MC Data Ratios for (int i = 0; i < mcRatio->GetNbinsX(); i++) { dataRatio->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1)); dataRatio->SetBinError(i + 1, fDataHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1)); mcRatio->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1) / fMCHist->GetBinContent(i + 1)); mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1)); } // Write ratios mcRatio->Write(); dataRatio->Write(); delete mcRatio; delete dataRatio; } //******************************************************************** void Measurement1D::WriteShapePlot() { //******************************************************************** TH1D* mcShape = (TH1D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); TH1D* dataShape = (TH1D*)fDataHist->Clone((fName + "_data_SHAPE").c_str()); if (fShapeCovar) StatUtils::SetDataErrorFromCov(dataShape, fShapeCovar, 1E-38); double shapeScale = 1.0; if (fIsRawEvents) { shapeScale = fDataHist->Integral() / fMCHist->Integral(); } else { shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); } mcShape->Scale(shapeScale); std::stringstream ss; ss << shapeScale; mcShape->SetTitle(ss.str().c_str()); mcShape->SetLineWidth(3); mcShape->SetLineStyle(7); mcShape->Write(); dataShape->Write(); delete mcShape; } //******************************************************************** void Measurement1D::WriteShapeRatioPlot() { //******************************************************************** // Get a mcshape histogram TH1D* mcShape = (TH1D*)fMCHist->Clone((fName + "_MC_SHAPE").c_str()); double shapeScale = 1.0; if (fIsRawEvents) { shapeScale = fDataHist->Integral() / fMCHist->Integral(); } else { shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width"); } mcShape->Scale(shapeScale); // Create shape ratio histograms TH1D* mcShapeRatio = (TH1D*)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str()); TH1D* dataShapeRatio = (TH1D*)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str()); // Divide the histograms mcShapeRatio->Divide(mcShape); dataShapeRatio->Divide(mcShape); // Colour the shape ratio plots mcShapeRatio->SetLineWidth(3); mcShapeRatio->SetLineStyle(7); mcShapeRatio->Write(); dataShapeRatio->Write(); delete mcShapeRatio; delete dataShapeRatio; } //// CRAP TO BE REMOVED //******************************************************************** void Measurement1D::SetupMeasurement(std::string inputfile, std::string type, FitWeight * rw, std::string fkdt) { //******************************************************************** nuiskey samplekey = Config::CreateKey("sample"); samplekey.Set("name", fName); samplekey.Set("type",type); samplekey.Set("input",inputfile); fSettings = LoadSampleSettings(samplekey); // Reset everything to NULL // Init(); // Check if name contains Evt, indicating that it is a raw number of events // measurements and should thus be treated as once fIsRawEvents = false; if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) { fIsRawEvents = true; LOG(SAM) << "Found event rate measurement but fIsRawEvents == false!" << std::endl; LOG(SAM) << "Overriding this and setting fIsRawEvents == true!" << std::endl; } fIsEnu1D = false; if (fName.find("XSec_1DEnu") != std::string::npos) { fIsEnu1D = true; LOG(SAM) << "::" << fName << "::" << std::endl; LOG(SAM) << "Found XSec Enu measurement, applying flux integrated scaling, " "not flux averaged!" << std::endl; } if (fIsEnu1D && fIsRawEvents) { LOG(SAM) << "Found 1D Enu XSec distribution AND fIsRawEvents, is this " "really correct?!" << std::endl; LOG(SAM) << "Check experiment constructor for " << fName << " and correct this!" << std::endl; LOG(SAM) << "I live in " << __FILE__ << ":" << __LINE__ << std::endl; exit(-1); } fRW = rw; if (!fInput and !fIsJoint) SetupInputs(inputfile); // Set Default Options SetFitOptions(fDefaultTypes); // Set Passed Options SetFitOptions(type); // Still adding support for flat flux inputs // // Set Enu Flux Scaling // if (isFlatFluxFolding) this->Input()->ApplyFluxFolding( // this->defaultFluxHist ); // FinaliseMeasurement(); } //******************************************************************** void Measurement1D::SetupDefaultHist() { //******************************************************************** // Setup fMCHist fMCHist = (TH1D*)fDataHist->Clone(); fMCHist->SetNameTitle((fName + "_MC").c_str(), (fName + "_MC" + fPlotTitles).c_str()); // Setup fMCFine Int_t nBins = fMCHist->GetNbinsX(); fMCFine = new TH1D( (fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(), nBins * 6, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(nBins + 1)); fMCStat = (TH1D*)fMCHist->Clone(); fMCStat->Reset(); fMCHist->Reset(); fMCFine->Reset(); // Setup the NEUT Mode Array PlotUtils::CreateNeutModeArray((TH1D*)fMCHist, (TH1**)fMCHist_PDG); PlotUtils::ResetNeutModeArray((TH1**)fMCHist_PDG); // Setup bin masks using sample name if (fIsMask) { std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask"); if (maskloc.empty()) { maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask"; } SetBinMask(maskloc); } fMCHist_Modes = new TrueModeStack( (fName + "_MODES").c_str(), ("True Channels"), fMCHist); SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write); return; } //******************************************************************** void Measurement1D::SetDataValues(std::string dataFile) { //******************************************************************** // Override this function if the input file isn't in a suitable format LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl; fDataHist = PlotUtils::GetTH1DFromFile(dataFile, (fName + "_data"), fPlotTitles); fDataTrue = (TH1D*)fDataHist->Clone(); // Number of data points is number of bins fNDataPointsX = fDataHist->GetXaxis()->GetNbins(); return; }; //******************************************************************** void Measurement1D::SetDataFromDatabase(std::string inhistfile, std::string histname) { //******************************************************************** LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname << std::endl; fDataHist = PlotUtils::GetTH1DFromRootFile( (GeneralUtils::GetTopLevelDir() + "/data/" + inhistfile), histname); fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str()); return; }; //******************************************************************** void Measurement1D::SetDataFromFile(std::string inhistfile, std::string histname) { //******************************************************************** LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname << std::endl; fDataHist = PlotUtils::GetTH1DFromRootFile((inhistfile), histname); fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str()); return; }; //******************************************************************** void Measurement1D::SetCovarMatrix(std::string covarFile) { //******************************************************************** // Covariance function, only really used when reading in the MB Covariances. TFile* tempFile = new TFile(covarFile.c_str(), "READ"); TH2D* covarPlot = new TH2D(); // TH2D* decmpPlot = new TH2D(); TH2D* covarInvPlot = new TH2D(); TH2D* fFullCovarPlot = new TH2D(); std::string covName = ""; std::string covOption = FitPar::Config().GetParS("thrown_covariance"); if (fIsShape || fIsFree) covName = "shp_"; if (fIsDiag) covName += "diag"; else covName += "full"; covarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str()); covarInvPlot = (TH2D*)tempFile->Get((covName + "covinv").c_str()); if (!covOption.compare("SUB")) fFullCovarPlot = (TH2D*)tempFile->Get((covName + "cov").c_str()); else if (!covOption.compare("FULL")) fFullCovarPlot = (TH2D*)tempFile->Get("fullcov"); else ERR(WRN) << "Incorrect thrown_covariance option in parameters." << std::endl; int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral()); int covdim = int(fDataHist->GetNbinsX()); this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); fDecomp = new TMatrixDSym(dim); int row, column = 0; row = 0; column = 0; for (Int_t i = 0; i < covdim; i++) { // if (this->masked->GetBinContent(i+1) > 0) continue; for (Int_t j = 0; j < covdim; j++) { // if (this->masked->GetBinContent(j+1) > 0) continue; (*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1); (*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1); column++; } column = 0; row++; } // Set bin errors on data if (!fIsDiag) { StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar); } // Get Deteriminant and inverse matrix // fCovDet = this->covar->Determinant(); TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); return; }; //******************************************************************** // Sets the covariance matrix from a provided file in a text format // scale is a multiplicative pre-factor to apply in the case where the // covariance is given in some unit (e.g. 1E-38) void Measurement1D::SetCovarMatrixFromText(std::string covarFile, int dim, double scale) { //******************************************************************** // Make a counter to track the line number int row = 0; std::string line; std::ifstream covarread(covarFile.c_str(), std::ifstream::in); this->covar = new TMatrixDSym(dim); fFullCovar = new TMatrixDSym(dim); if (covarread.is_open()) LOG(SAM) << "Reading covariance matrix from file: " << covarFile << std::endl; else ERR(FTL) << "Covariance matrix provided is incorrect: " << covarFile << std::endl; // Loop over the lines in the file while (std::getline(covarread >> std::ws, line, '\n')) { int column = 0; // Loop over entries and insert them into matrix std::vector entries = GeneralUtils::ParseToDbl(line, " "); if (entries.size() <= 1) { ERR(WRN) << "SetCovarMatrixFromText -> Covariance matrix only has <= 1 " "entries on this line: " << row << std::endl; } for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { (*covar)(row, column) = *iter; (*fFullCovar)(row, column) = *iter; column++; } row++; } covarread.close(); // Scale the actualy covariance matrix by some multiplicative factor (*fFullCovar) *= scale; // Robust matrix inversion method TDecompSVD LU = TDecompSVD(*this->covar); // THIS IS ACTUALLY THE INVERSE COVARIANCE MATRIXA AAAAARGH delete this->covar; this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); // Now need to multiply by the scaling factor // If the covariance (*this->covar) *= 1. / (scale); return; }; //******************************************************************** void Measurement1D::SetCovarMatrixFromCorrText(std::string corrFile, int dim) { //******************************************************************** // Make a counter to track the line number int row = 0; std::string line; std::ifstream corr(corrFile.c_str(), std::ifstream::in); this->covar = new TMatrixDSym(dim); this->fFullCovar = new TMatrixDSym(dim); if (corr.is_open()) LOG(SAM) << "Reading and converting correlation matrix from file: " << corrFile << std::endl; else { ERR(FTL) << "Correlation matrix provided is incorrect: " << corrFile << std::endl; exit(-1); } while (std::getline(corr >> std::ws, line, '\n')) { int column = 0; // Loop over entries and insert them into matrix // Multiply by the errors to get the covariance, rather than the correlation // matrix std::vector entries = GeneralUtils::ParseToDbl(line, " "); for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { double val = (*iter) * this->fDataHist->GetBinError(row + 1) * 1E38 * this->fDataHist->GetBinError(column + 1) * 1E38; if (val == 0) { ERR(FTL) << "Found a zero value in the covariance matrix, assuming " "this is an error!" << std::endl; exit(-1); } (*this->covar)(row, column) = val; (*this->fFullCovar)(row, column) = val; column++; } row++; } // Robust matrix inversion method TDecompSVD LU = TDecompSVD(*this->covar); delete this->covar; this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); return; }; //******************************************************************** // FullUnits refers to if we have "real" unscaled units in the covariance matrix, e.g. 1E-76. // If this is the case we need to scale it so that the chi2 contribution is correct // NUISANCE internally assumes the covariance matrix has units of 1E76 void Measurement1D::SetCovarFromDataFile(std::string covarFile, std::string covName, bool FullUnits) { //******************************************************************** LOG(SAM) << "Getting covariance from " << covarFile << "->" << covName << std::endl; TFile* tempFile = new TFile(covarFile.c_str(), "READ"); TH2D* covPlot = (TH2D*)tempFile->Get(covName.c_str()); covPlot->SetDirectory(0); // Scale the covariance matrix if it comes in normal units if (FullUnits) { covPlot->Scale(1.E76); } int dim = covPlot->GetNbinsX(); fFullCovar = new TMatrixDSym(dim); for (int i = 0; i < dim; i++) { for (int j = 0; j < dim; j++) { (*fFullCovar)(i, j) = covPlot->GetBinContent(i + 1, j + 1); } } this->covar = (TMatrixDSym*)fFullCovar->Clone(); fDecomp = (TMatrixDSym*)fFullCovar->Clone(); TDecompSVD LU = TDecompSVD(*this->covar); this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), ""); TDecompChol LUChol = TDecompChol(*fDecomp); LUChol.Decompose(); fDecomp = new TMatrixDSym(dim, LU.GetU().GetMatrixArray(), ""); return; }; // //******************************************************************** // void Measurement1D::SetBinMask(std::string maskFile) { // //******************************************************************** // // Create a mask histogram. // int nbins = fDataHist->GetNbinsX(); // fMaskHist = // new TH1I((fName + "_fMaskHist").c_str(), // (fName + "_fMaskHist; Bin; Mask?").c_str(), nbins, 0, nbins); // std::string line; // std::ifstream mask(maskFile.c_str(), std::ifstream::in); // if (mask.is_open()) // LOG(SAM) << "Reading bin mask from file: " << maskFile << std::endl; // else // LOG(FTL) << " Cannot find mask file." << std::endl; // while (std::getline(mask >> std::ws, line, '\n')) { // std::vector entries = GeneralUtils::ParseToInt(line, " "); // // Skip lines with poorly formatted lines // if (entries.size() < 2) { // LOG(WRN) << "Measurement1D::SetBinMask(), couldn't parse line: " << line // << std::endl; // continue; // } // // The first index should be the bin number, the second should be the mask // // value. // fMaskHist->SetBinContent(entries[0], entries[1]); // } // // Set masked data bins to zero // PlotUtils::MaskBins(fDataHist, fMaskHist); // return; // } // //******************************************************************** // void Measurement1D::GetBinContents(std::vector& cont, // std::vector& err) { // //******************************************************************** // // Return a vector of the main bin contents // for (int i = 0; i < fMCHist->GetNbinsX(); i++) { // cont.push_back(fMCHist->GetBinContent(i + 1)); // err.push_back(fMCHist->GetBinError(i + 1)); // } // return; // }; /* XSec Functions */ // //******************************************************************** // void Measurement1D::SetFluxHistogram(std::string fluxFile, int minE, int // maxE, // double fluxNorm) { // //******************************************************************** // // Note this expects the flux bins to be given in terms of MeV // LOG(SAM) << "Reading flux from file: " << fluxFile << std::endl; // TGraph f(fluxFile.c_str(), "%lg %lg"); // fFluxHist = // new TH1D((fName + "_flux").c_str(), (fName + "; E_{#nu} (GeV)").c_str(), // f.GetN() - 1, minE, maxE); // Double_t* yVal = f.GetY(); // for (int i = 0; i < fFluxHist->GetNbinsX(); ++i) // fFluxHist->SetBinContent(i + 1, yVal[i] * fluxNorm); // }; // //******************************************************************** // double Measurement1D::TotalIntegratedFlux(std::string intOpt, double low, // double high) { // //******************************************************************** // if (fInput->GetType() == kGiBUU) { // return 1.0; // } // // The default case of low = -9999.9 and high = -9999.9 // if (low == -9999.9) low = this->EnuMin; // if (high == -9999.9) high = this->EnuMax; // int minBin = fFluxHist->GetXaxis()->FindBin(low); // int maxBin = fFluxHist->GetXaxis()->FindBin(high); // // Get integral over custom range // double integral = fFluxHist->Integral(minBin, maxBin + 1, intOpt.c_str()); // return integral; // }; diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx index 61751f8..50ef1bc 100644 --- a/src/InputHandler/FitEvent.cxx +++ b/src/InputHandler/FitEvent.cxx @@ -1,429 +1,429 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is pddrt 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 "FitEvent.h" #include #include "TObjArray.h" FitEvent::FitEvent() { fGenInfo = NULL; kRemoveFSIParticles = true; kRemoveUndefParticles = true; AllocateParticleStack(400); }; void FitEvent::AddGeneratorInfo(GeneratorInfoBase* gen) { fGenInfo = gen; gen->AllocateParticleStack(kMaxParticles); } void FitEvent::AllocateParticleStack(int stacksize) { LOG(DEB) << "Allocating particle stack of size: " << stacksize << std::endl; kMaxParticles = stacksize; fParticleList = new FitParticle*[kMaxParticles]; fParticleMom = new double*[kMaxParticles]; fParticleState = new UInt_t[kMaxParticles]; fParticlePDG = new int[kMaxParticles]; fOrigParticleMom = new double*[kMaxParticles]; fOrigParticleState = new UInt_t[kMaxParticles]; fOrigParticlePDG = new int[kMaxParticles]; for (size_t i = 0; i < kMaxParticles; i++) { fParticleList[i] = NULL; fParticleMom[i] = new double[4]; fOrigParticleMom[i] = new double[4]; } if (fGenInfo) fGenInfo->AllocateParticleStack(kMaxParticles); } void FitEvent::ExpandParticleStack(int stacksize) { DeallocateParticleStack(); AllocateParticleStack(stacksize); } void FitEvent::DeallocateParticleStack() { for (size_t i = 0; i < kMaxParticles; i++) { if (fParticleList[i]) delete fParticleList[i]; delete fParticleMom[i]; delete fOrigParticleMom[i]; } delete fParticleMom; delete fOrigParticleMom; delete fParticleList; delete fParticleState; delete fParticlePDG; delete fOrigParticleState; delete fOrigParticlePDG; if (fGenInfo) fGenInfo->DeallocateParticleStack(); kMaxParticles = 0; } void FitEvent::ClearFitParticles() { for (size_t i = 0; i < kMaxParticles; i++) { fParticleList[i] = NULL; } } void FitEvent::FreeFitParticles() { for (size_t i = 0; i < kMaxParticles; i++) { FitParticle* fp = fParticleList[i]; if (fp) delete fp; fParticleList[i] = NULL; } } void FitEvent::ResetParticleList() { for (unsigned int i = 0; i < kMaxParticles; i++) { FitParticle* fp = fParticleList[i]; if (fp) delete fp; fParticleList[i] = NULL; } } void FitEvent::HardReset() { for (unsigned int i = 0; i < kMaxParticles; i++) { fParticleList[i] = NULL; } } void FitEvent::ResetEvent() { Mode = 9999; fEventNo = -1; fTotCrs = -1.0; fTargetA = -1; fTargetZ = -1; fTargetH = -1; fBound = false; fNParticles = 0; if (fGenInfo) fGenInfo->Reset(); for (unsigned int i = 0; i < kMaxParticles; i++) { if (fParticleList[i]) delete fParticleList[i]; fParticleList[i] = NULL; continue; fParticlePDG[i] = 0; fParticleState[i] = kUndefinedState; fParticleMom[i][0] = 0.0; fParticleMom[i][1] = 0.0; fParticleMom[i][2] = 0.0; fParticleMom[i][3] = 0.0; fOrigParticlePDG[i] = 0; fOrigParticleState[i] = kUndefinedState; fOrigParticleMom[i][0] = 0.0; fOrigParticleMom[i][1] = 0.0; fOrigParticleMom[i][2] = 0.0; fOrigParticleMom[i][3] = 0.0; } } void FitEvent::OrderStack() { // Copy current stack int npart = fNParticles; for (int i = 0; i < npart; i++) { fOrigParticlePDG[i] = fParticlePDG[i]; fOrigParticleState[i] = fParticleState[i]; fOrigParticleMom[i][0] = fParticleMom[i][0]; fOrigParticleMom[i][1] = fParticleMom[i][1]; fOrigParticleMom[i][2] = fParticleMom[i][2]; fOrigParticleMom[i][3] = fParticleMom[i][3]; } // Now run loops for each particle fNParticles = 0; int stateorder[6] = {kInitialState, kFinalState, kFSIState, kNuclearInitial, kNuclearRemnant, kUndefinedState}; for (int s = 0; s < 6; s++) { for (int i = 0; i < npart; i++) { if ((UInt_t)fOrigParticleState[i] != (UInt_t)stateorder[s]) continue; fParticlePDG[fNParticles] = fOrigParticlePDG[i]; fParticleState[fNParticles] = fOrigParticleState[i]; fParticleMom[fNParticles][0] = fOrigParticleMom[i][0]; fParticleMom[fNParticles][1] = fOrigParticleMom[i][1]; fParticleMom[fNParticles][2] = fOrigParticleMom[i][2]; fParticleMom[fNParticles][3] = fOrigParticleMom[i][3]; fNParticles++; } } if (LOG_LEVEL(DEB)) { LOG(DEB) << "Ordered stack" << std::endl; for (int i = 0; i < fNParticles; i++) { LOG(DEB) << "Particle " << i << ". " << fParticlePDG[i] << " " << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " << fParticleMom[i][2] << " " << fParticleMom[i][3] << " " << fParticleState[i] << std::endl; } } if (fNParticles != npart) { ERR(FTL) << "Dropped some particles when ordering the stack!" << std::endl; } return; } void FitEvent::Print() { if (LOG_LEVEL(FIT)) { LOG(FIT) << "FITEvent print" << std::endl; - LOG(FIT) << "Mode: " << Mode << std::endl; + LOG(FIT) << "Mode: " << Mode << ", Weight: " << InputWeight << std::endl; LOG(FIT) << "Particles: " << fNParticles << std::endl; LOG(FIT) << " -> Particle Stack " << std::endl; for (int i = 0; i < fNParticles; i++) { LOG(FIT) << " -> -> " << i << ". " << fParticlePDG[i] << " " << fParticleState[i] << " " << " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1] << ", " << fParticleMom[i][2] << ", " << fParticleMom[i][3] << ")." << std::endl; } } return; } /* Read/Write own event class */ void FitEvent::SetBranchAddress(TChain* tn) { tn->SetBranchAddress("Mode", &Mode); tn->SetBranchAddress("EventNo", &fEventNo); tn->SetBranchAddress("TotCrs", &fTotCrs); tn->SetBranchAddress("TargetA", &fTargetA); tn->SetBranchAddress("TargetH", &fTargetH); tn->SetBranchAddress("Bound", &fBound); tn->SetBranchAddress("RWWeight", &SavedRWWeight); tn->SetBranchAddress("InputWeight", &InputWeight); } void FitEvent::AddBranchesToTree(TTree* tn) { tn->Branch("Mode", &Mode, "Mode/I"); tn->Branch("EventNo", &fEventNo, "EventNo/i"); tn->Branch("TotCrs", &fTotCrs, "TotCrs/D"); tn->Branch("TargetA", &fTargetA, "TargetA/I"); tn->Branch("TargetH", &fTargetH, "TargetH/I"); tn->Branch("Bound", &fBound, "Bound/O"); tn->Branch("RWWeight", &RWWeight, "RWWeight/D"); tn->Branch("InputWeight", &InputWeight, "InputWeight/D"); tn->Branch("NParticles", &fNParticles, "NParticles/I"); tn->Branch("ParticleState", fOrigParticleState, "ParticleState[NParticles]/i"); tn->Branch("ParticlePDG", fOrigParticlePDG, "ParticlePDG[NParticles]/I"); tn->Branch("ParticleMom", fOrigParticleMom, "ParticleMom[NParticles][4]/D"); } // ------- EVENT ACCESS FUNCTION --------- // TLorentzVector FitEvent::GetParticleP4(int index) const { if (index == -1 or index >= fNParticles) return TLorentzVector(); return TLorentzVector(fParticleMom[index][0], fParticleMom[index][1], fParticleMom[index][2], fParticleMom[index][3]); } TVector3 FitEvent::GetParticleP3(int index) const { if (index == -1 or index >= fNParticles) return TVector3(); return TVector3(fParticleMom[index][0], fParticleMom[index][1], fParticleMom[index][2]); } double FitEvent::GetParticleMom(int index) const { if (index == -1 or index >= fNParticles) return 0.0; return sqrt(fParticleMom[index][0] * fParticleMom[index][0] + fParticleMom[index][1] * fParticleMom[index][1] + fParticleMom[index][2] * fParticleMom[index][2]); } double FitEvent::GetParticleMom2(int index) const { if (index == -1 or index >= fNParticles) return 0.0; return fabs((fParticleMom[index][0] * fParticleMom[index][0] + fParticleMom[index][1] * fParticleMom[index][1] + fParticleMom[index][2] * fParticleMom[index][2])); } double FitEvent::GetParticleE(int index) const { if (index == -1 or index >= fNParticles) return 0.0; return fParticleMom[index][3]; } int FitEvent::GetParticleState(int index) const { if (index == -1 or index >= fNParticles) return kUndefinedState; return (fParticleState[index]); } int FitEvent::GetParticlePDG(int index) const { if (index == -1 or index >= fNParticles) return 0; return (fParticlePDG[index]); } FitParticle* FitEvent::GetParticle(int const i) { // Check Valid Index if (i == -1) { return NULL; } // Check Valid if (i > fNParticles) { ERR(FTL) << "Requesting particle beyond stack!" << std::endl << "i = " << i << " N = " << fNParticles << std::endl << "Mode = " << Mode << std::endl; throw; } if (!fParticleList[i]) { /* std::cout << "Creating particle with values i " << i << " "; std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " << fParticleMom[i][2] << " " << fParticleMom[i][3] << " "; std::cout << fParticlePDG[i] << " " << fParticleState[i] << std::endl; */ fParticleList[i] = new FitParticle(fParticleMom[i][0], fParticleMom[i][1], fParticleMom[i][2], fParticleMom[i][3], fParticlePDG[i], fParticleState[i]); } else { /* std::cout << "Filling particle with values i " << i << " "; std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " << fParticleMom[i][2] << " " << fParticleMom[i][3] << " "; std::cout << fParticlePDG[i] << " "<< fParticleState[i] <SetValues(fParticleMom[i][0], fParticleMom[i][1], fParticleMom[i][2], fParticleMom[i][3], fParticlePDG[i], fParticleState[i]); } return fParticleList[i]; } bool FitEvent::HasParticle(int const pdg, int const state) const { bool found = false; for (int i = 0; i < fNParticles; i++) { if (state != -1 && fParticleState[i] != (uint)state) continue; if (fParticlePDG[i] == pdg) found = true; } return found; } int FitEvent::NumParticle(int const pdg, int const state) const { int nfound = 0; for (int i = 0; i < fNParticles; i++) { if (state != -1 and fParticleState[i] != (uint)state) continue; if (pdg == 0 or fParticlePDG[i] == pdg) nfound += 1; } return nfound; } std::vector FitEvent::GetAllParticleIndices(int const pdg, int const state) const { std::vector indexlist; for (int i = 0; i < fNParticles; i++) { if (state != -1 and fParticleState[i] != (uint)state) continue; if (pdg == 0 or fParticlePDG[i] == pdg) { indexlist.push_back(i); } } return indexlist; } std::vector FitEvent::GetAllParticle(int const pdg, int const state) { std::vector indexlist = GetAllParticleIndices(pdg, state); std::vector plist; for (std::vector::iterator iter = indexlist.begin(); iter != indexlist.end(); iter++) { plist.push_back(GetParticle((*iter))); } return plist; } int FitEvent::GetHMParticleIndex(int const pdg, int const state) const { double maxmom2 = -9999999.9; int maxind = -1; for (int i = 0; i < fNParticles; i++) { if (state != -1 and fParticleState[i] != (uint)state) continue; if (pdg == 0 or fParticlePDG[i] == pdg) { double newmom2 = GetParticleMom2(i); if (newmom2 > maxmom2) { maxind = i; maxmom2 = newmom2; } } } return maxind; } int FitEvent::GetBeamNeutrinoIndex(void) const { for (int i = 0; i < fNParticles; i++) { if (fParticleState[i] != kInitialState) continue; int pdg = abs(fParticlePDG[i]); if (pdg == 12 or pdg == 14 or pdg == 16) { return i; } } return 0; } int FitEvent::GetBeamElectronIndex(void) const { return GetHMISParticleIndex(11); } int FitEvent::GetBeamPionIndex(void) const { return GetHMISParticleIndex(PhysConst::pdg_pions); } int FitEvent::NumFSMesons() { int nMesons = 0; for (int i = 0; i < fNParticles; i++) { if (fParticleState[i] != kFinalState) continue; if (abs(fParticlePDG[i]) >= 111 && abs(fParticlePDG[i]) <= 557) nMesons += 1; } return nMesons; } int FitEvent::NumFSLeptons(void) const { int nLeptons = 0; for (int i = 0; i < fNParticles; i++) { if (fParticleState[i] != kFinalState) continue; if (abs(fParticlePDG[i]) == 11 || abs(fParticlePDG[i]) == 13 || abs(fParticlePDG[i]) == 15) nLeptons += 1; } return nLeptons; } diff --git a/src/Logger/FitLogger.h b/src/Logger/FitLogger.h index ce46f6c..2fede17 100644 --- a/src/Logger/FitLogger.h +++ b/src/Logger/FitLogger.h @@ -1,221 +1,221 @@ // 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 FITLOGGER_HPP #define FITLOGGER_HPP /*! * \addtogroup FitBase * @{ */ #include #include #include #include #include "TRandom3.h" #include "Initialiser.h" #include "NuisConfig.h" #define RESET "\033[0m" #define BLACK "\033[30m" /* Black */ #define RED "\033[31m" /* Red */ #define GREEN "\033[32m" /* Green */ #define YELLOW "\033[33m" /* Yellow */ #define BLUE "\033[34m" /* Blue */ #define MAGENTA "\033[35m" /* Magenta */ #define CYAN "\033[36m" /* Cyan */ #define WHITE "\033[37m" /* White */ #define BOLDBLACK "\033[1m\033[30m" /* Bold Black */ #define BOLDRED "\033[1m\033[31m" /* Bold Red */ #define BOLDGREEN "\033[1m\033[32m" /* Bold Green */ #define BOLDYELLOW "\033[1m\033[33m" /* Bold Yellow */ #define BOLDBLUE "\033[1m\033[34m" /* Bold Blue */ #define BOLDMAGENTA "\033[1m\033[35m" /* Bold Magenta */ #define BOLDCYAN "\033[1m\033[36m" /* Bold Cyan */ #define BOLDWHITE "\033[1m\033[37m" /* Bold White */ namespace Logger { extern int log_verb; //!< Current VERBOSITY extern int err_verb; //!< Current ERROR VERBOSITY extern bool external_verb; extern bool use_colors; //!< Use BASH Terminal Colors Flag extern bool super_rainbow_mode; //!< For when fitting gets boring. extern unsigned int super_rainbow_mode_colour; extern bool showtrace; // Quick Tracing for debugging extern int nloggercalls; extern int timelastlog; extern std::streambuf* default_cout; //!< Where the STDOUT stream is currently directed extern std::streambuf* default_cerr; //!< Where the STDERR stream is currently directed extern std::ofstream redirect_stream; //!< Where should unwanted messages be thrown } /// Returns full path to file currently in #define __FILENAME__ \ (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__) // ------ LOGGER FUNCTIONS ------------ // namespace Logger { /// NULL Output Stream extern std::ofstream __LOG_nullstream; /// Logging Stream extern std::ostream* __LOG_outstream; } /// Fitter VERBOSITY Enumerations /// These go through the different depths of the fitter. /// /// 0 QUIET - Little output. /// 1 FIT - Top Level Minimizer Status /// 2 MIN - Output from the FCN Minimizer Functions /// 3 SAM - Output from each of the samples during setup etc /// 4 REC - Output during each reconfigure. Percentage progress etc. /// 5 SIG - Output during every signal event that is found. /// 6 EVT - Output during every event. /// -1 DEB - Will print only debugging info wherever a LOG(DEB) statement was /// made -enum __LOG_levels { DEB = -1, QUIET, FIT, MIN, SAM, REC, SIG, EVT }; +enum __LOG_levels { QUIET = 0, FIT, MIN, SAM, REC, SIG, EVT, DEB }; /// Returns log level for a given file/function int __GETLOG_LEVEL(int level, const char* filename, const char* funct); /// Actually runs the logger std::ostream& __OUTLOG(int level, const char* filename, const char* funct, int line); /// Global Logging Definitions #define QLOG(level, stream) \ { \ if (Logger::log_verb >= \ __GETLOG_LEVEL(level, __FILENAME__, __FUNCTION__)) { \ __OUTLOG(level, __FILENAME__, __FUNCTION__, __LINE__) << stream \ << std::endl; \ } \ }; #define BREAK(level) \ { \ \ if (Logger::log_verb >= \ __GETLOG_LEVEL(level, __FILENAME__, __FUNCTION__)) { \ __OUTLOG(level, __FILENAME__, __FUNCTION__, __LINE__) << std::endl; \ } \ }; /// Return whether logging level is valid bool LOGGING(int level); /// Set Global Verbosity void SETVERBOSITY(int level); /// Set Global Verbosity from String void SETVERBOSITY(std::string verb); /// Set Trace Option void SETTRACE(bool val); // ----------- ERROR FUNCTIONS ---------- // /// Error Stream extern std::ostream* __ERR_outstream; /// Fitter ERROR VERBOSITY Enumerations /// /// 0 QUIET - No Error Output /// 1 FTL - Show errors only if fatal /// 2 WRN - Show Warning messages enum __ERR_levels { ERRQUIET = 0, FTL, WRN }; /// Actually runs the error messager std::ostream& __OUTERR(int level, const char* filename, const char* funct, int line); /// Error Logging Function #define ERROR(level, stream) \ { \ __OUTERR(level, __FILENAME__, __FUNCTION__, __LINE__) << stream \ << std::endl; \ }; // ----------- ERROR HANDLING ------------- // /// Exit the program with given error message stream #define THROW(stream) \ { \ __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) << stream \ << std::endl; \ __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) \ << "Attempting to save output file." << std::endl; \ - if (Config::Get().out && Config::Get().out->IsOpen()) { \ - Config::Get().out->Write(); \ - Config::Get().out->Close(); \ + if (Config::Get().out && Config::Get().out->IsOpen()) { \ + Config::Get().out->Write(); \ + Config::Get().out->Close(); \ __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) << "Done." \ << std::endl; \ } else { \ __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) \ << "No output file set." << std::endl; \ } \ __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) << "Exiting!" \ << std::endl; \ std::abort(); \ } // ----------- External Logging ----------- // void SETEXTERNALVERBOSITY(int level); void StopTalking(); void StartTalking(); extern "C" { void shhnuisancepythiaitokay_(void); void canihaznuisancepythia_(void); } // ---------- LEGACY FUNCTIONS -------------- // bool LOG_LEVEL(int level); //! Set LOG VERBOSITY from a string void LOG_VERB(std::string verb); inline void LOG_VERB(int verb) { Logger::log_verb = verb; }; void SET_TRACE(bool val); //! Set ERROR VERBOSITY from a string void ERR_VERB(std::string verb); inline void ERR_VERB(int verb) { Logger::err_verb = verb; }; /// Logging Function. Use as a string stream. e.g. LOG(SAM) << "This sample is /// dope." << std::endl; std::ostream& _LOG(int level, const char* filename, const char* funct, int line); #define LOG(level) _LOG(level, __FILENAME__, __FUNCTION__, __LINE__) //! Error Function. Use as a string stream. e.g. ERR(FTL) << "The fit is //! completely buggered." << std::endl; std::ostream& _ERR(int level, const char* filename, const char* funct, int line); #define ERR(level) _ERR(level, __FILENAME__, __FUNCTION__, __LINE__) /*! @} */ #endif diff --git a/src/Reweight/ModeNormEngine.h b/src/Reweight/ModeNormEngine.h index 4915298..6b5da6f 100644 --- a/src/Reweight/ModeNormEngine.h +++ b/src/Reweight/ModeNormEngine.h @@ -1,82 +1,87 @@ #ifndef ModeNormEngine_H #define ModeNormEngine_H #include "FitLogger.h" #include "GeneratorUtils.h" #include "WeightEngineBase.h" class ModeNormEngine : public WeightEngineBase { public: ModeNormEngine(std::string name = "ModeNormEngine") : fName(name){}; ~ModeNormEngine(){}; void IncludeDial(std::string name, double startval) { int rwenum = Reweight::ConvDial(name, kMODENORM); int mode = Reweight::RemoveDialType(rwenum); if (fDialEnumIndex.count(mode)) { THROW("Mode dial: " << mode << " already included. Cannot include twice."); } fDialEnumIndex[mode] = fDialValues.size(); fDialValues.push_back(startval); QLOG(FIT, "Added mode dial for mode: " << mode); } void SetDialValue(int rwenum, double val) { int mode = Reweight::RemoveDialType(rwenum); if (!fDialEnumIndex.count(mode)) { THROW("Mode dial: " << mode << " has not been included. Cannot set value."); } + QLOG(DEB, "[INFO]: ModeNormEngine ObsMode: " << mode << " weight " << val + << ", rwenum = " << rwenum); fDialValues[fDialEnumIndex[mode]] = val; } void SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kMODENORM), val); } void Reconfigure(bool silent = false) { (void)silent; } static int ModeToDial(int mode) { return 60 + mode; } double CalcWeight(BaseFitEvt* evt) { int mode = ModeToDial(abs(evt->Mode)); if (!fDialEnumIndex.count(mode)) { return 1; } + QLOG(DEB, "[INFO]: Ev mode " + << evt->Mode << ", ObsMode: " << mode + << ", weight = " << fDialValues[fDialEnumIndex[mode]]); return fDialValues[fDialEnumIndex[mode]]; }; bool NeedsEventReWeight() { return false; }; double GetDialValue(std::string name) { int rwenum = Reweight::ConvDial(name, kMODENORM); int mode = Reweight::RemoveDialType(rwenum); if (fDialEnumIndex.count(mode)) { return fDialValues[fDialEnumIndex[mode]]; } else { return 0xdeadbeef; } } static int SystEnumFromString(std::string const& name) { std::vector splits = GeneralUtils::ParseToStr(name, "_"); if (splits.size() != 2) { ERR(FTL) << "Attempting to parse dial name: \"" << name << "\" as a mode norm dial but failed. Expect e.g. \"mode_2\"." << std::endl; } int mode_num = GeneralUtils::StrToInt(splits[1]); if (!mode_num) { ERR(FTL) << "Attempting to parse dial name: \"" << name << "\" as a mode norm dial but failed." << std::endl; throw; } return ModeToDial(mode_num); } std::map fDialEnumIndex; std::vector fDialValues; std::string fName; }; #endif diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx index 7204c64..1a1bd77 100644 --- a/src/Reweight/WeightUtils.cxx +++ b/src/Reweight/WeightUtils.cxx @@ -1,615 +1,614 @@ #include "WeightUtils.h" #include "FitLogger.h" #ifdef __T2KREW_ENABLED__ #include "T2KGenieReWeight.h" #include "T2KNIWGReWeight.h" #include "T2KNIWGUtils.h" #include "T2KNeutReWeight.h" #include "T2KNeutUtils.h" #include "T2KReWeight.h" using namespace t2krew; #endif #ifdef __NIWG_ENABLED__ #include "NIWGReWeight.h" #include "NIWGReWeight1piAngle.h" #include "NIWGReWeight2010a.h" #include "NIWGReWeight2012a.h" #include "NIWGReWeight2014a.h" #include "NIWGReWeightDeltaMass.h" #include "NIWGReWeightEffectiveRPA.h" #include "NIWGReWeightHadronMultSwitch.h" #include "NIWGReWeightMEC.h" #include "NIWGReWeightPiMult.h" #include "NIWGReWeightProtonFSIbug.h" #include "NIWGReWeightRPA.h" #include "NIWGReWeightSpectralFunc.h" #include "NIWGReWeightSplineEnu.h" #include "NIWGSyst.h" #include "NIWGSystUncertainty.h" #endif #ifdef __NEUT_ENABLED__ #include "NReWeight.h" #include "NReWeightCasc.h" #include "NReWeightNuXSecCCQE.h" #include "NReWeightNuXSecCCRES.h" #include "NReWeightNuXSecCOH.h" #include "NReWeightNuXSecDIS.h" #include "NReWeightNuXSecNC.h" #include "NReWeightNuXSecNCEL.h" #include "NReWeightNuXSecNCRES.h" #include "NReWeightNuXSecRES.h" #include "NReWeightNuclPiless.h" #include "NSyst.h" #include "NSystUncertainty.h" #include "neutpart.h" #include "neutvect.h" #endif #ifdef __NUWRO_ENABLED__ #include "event1.h" #endif #ifdef __NUWRO_REWEIGHT_ENABLED__ #include "NuwroReWeight.h" #include "NuwroReWeight_FlagNorm.h" #include "NuwroReWeight_QEL.h" #include "NuwroReWeight_SPP.h" #include "NuwroSyst.h" #include "NuwroSystUncertainty.h" #endif #ifdef __GENIE_ENABLED__ #include "EVGCore/EventRecord.h" #include "EVGCore/EventRecord.h" #include "GHEP/GHepRecord.h" #include "GSyst.h" #include "GSystUncertainty.h" #include "Ntuple/NtpMCEventRecord.h" #include "ReWeight/GReWeight.h" #include "ReWeight/GReWeightAGKY.h" #include "ReWeight/GReWeightDISNuclMod.h" #include "ReWeight/GReWeightFGM.h" #include "ReWeight/GReWeightFZone.h" #include "ReWeight/GReWeightINuke.h" #include "ReWeight/GReWeightNonResonanceBkg.h" #include "ReWeight/GReWeightNuXSecCCQE.h" #include "ReWeight/GReWeightNuXSecCCQEvec.h" #include "ReWeight/GReWeightNuXSecCCRES.h" #include "ReWeight/GReWeightNuXSecCOH.h" #include "ReWeight/GReWeightNuXSecDIS.h" #include "ReWeight/GReWeightNuXSecNC.h" #include "ReWeight/GReWeightNuXSecNCEL.h" #include "ReWeight/GReWeightNuXSecNCRES.h" #include "ReWeight/GReWeightResonanceDecay.h" using namespace genie; using namespace genie::rew; #endif #include "GlobalDialList.h" #include "ModeNormEngine.h" #include "NUISANCESyst.h" #include "OscWeightEngine.h" //******************************************************************** TF1 FitBase::GetRWConvFunction(std::string const &type, std::string const &name) { //******************************************************************** std::string dialfunc = "x"; std::string parType = type; double low = -10000.0; double high = 10000.0; if (parType.find("parameter") == std::string::npos) parType += "_parameter"; std::string line; std::ifstream card( (GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card") .c_str(), std::ifstream::in); while (std::getline(card >> std::ws, line, '\n')) { std::vector inputlist = GeneralUtils::ParseToStr(line, " "); // Check the line length if (inputlist.size() < 4) continue; // Check whether this is a comment if (inputlist[0].c_str()[0] == '#') continue; // Check whether this is the correct parameter type if (inputlist[0].compare(parType) != 0) continue; // Check the parameter name if (inputlist[1].compare(name) != 0) continue; // inputlist[2] should be the units... ignore for now dialfunc = inputlist[3]; // High and low are optional, check whether they exist if (inputlist.size() > 4) low = GeneralUtils::StrToDbl(inputlist[4]); if (inputlist.size() > 5) high = GeneralUtils::StrToDbl(inputlist[5]); } TF1 convfunc = TF1((name + "_convfunc").c_str(), dialfunc.c_str(), low, high); return convfunc; } //******************************************************************** std::string FitBase::GetRWUnits(std::string const &type, std::string const &name) { //******************************************************************** std::string unit = "sig."; std::string parType = type; if (parType.find("parameter") == std::string::npos) { parType += "_parameter"; } std::string line; std::ifstream card( (GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card") .c_str(), std::ifstream::in); while (std::getline(card >> std::ws, line, '\n')) { std::vector inputlist = GeneralUtils::ParseToStr(line, " "); // Check the line length if (inputlist.size() < 3) continue; // Check whether this is a comment if (inputlist[0].c_str()[0] == '#') continue; // Check whether this is the correct parameter type if (inputlist[0].compare(parType) != 0) continue; // Check the parameter name if (inputlist[1].compare(name) != 0) continue; unit = inputlist[2]; break; } return unit; } //******************************************************************** double FitBase::RWAbsToSigma(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.GetX(val); if (fabs(conv_val) < 1E-10) conv_val = 0.0; - std::cout << "AbsToSigma(" << name << ") = " << val << " -> " << conv_val - << std::endl; + QLOG(FIT, "AbsToSigma(" << name << ") = " << val << " -> " << conv_val); return conv_val; } //******************************************************************** double FitBase::RWSigmaToAbs(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.Eval(val); return conv_val; } //******************************************************************** double FitBase::RWFracToSigma(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.GetX((val * f1.Eval(0.0))); if (fabs(conv_val) < 1E-10) conv_val = 0.0; return conv_val; } //******************************************************************** double FitBase::RWSigmaToFrac(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.Eval(val) / f1.Eval(0.0); return conv_val; } int FitBase::ConvDialType(std::string const &type) { if (!type.compare("neut_parameter")) return kNEUT; else if (!type.compare("niwg_parameter")) return kNIWG; else if (!type.compare("nuwro_parameter")) return kNUWRO; else if (!type.compare("t2k_parameter")) return kT2K; else if (!type.compare("genie_parameter")) return kGENIE; else if (!type.compare("custom_parameter")) return kCUSTOM; else if (!type.compare("norm_parameter")) return kNORM; else if (!type.compare("likeweight_parameter")) return kLIKEWEIGHT; else if (!type.compare("spline_parameter")) return kSPLINEPARAMETER; else if (!type.compare("osc_parameter")) return kOSCILLATION; else if (!type.compare("modenorm_parameter")) return kMODENORM; else return kUNKNOWN; } std::string FitBase::ConvDialType(int type) { switch (type) { case kNEUT: { return "neut_parameter"; } case kNIWG: { return "niwg_parameter"; } case kNUWRO: { return "nuwro_parameter"; } case kT2K: { return "t2k_parameter"; } case kGENIE: { return "genie_parameter"; } case kNORM: { return "norm_parameter"; } case kCUSTOM: { return "custom_parameter"; } case kLIKEWEIGHT: { return "likeweight_parameter"; } case kSPLINEPARAMETER: { return "spline_parameter"; } case kOSCILLATION: { return "osc_parameter"; } case kMODENORM: { return "modenorm_parameter"; } default: return "unknown_parameter"; } } int FitBase::GetDialEnum(std::string const &type, std::string const &name) { return FitBase::GetDialEnum(FitBase::ConvDialType(type), name); } int FitBase::GetDialEnum(int type, std::string const &name) { int offset = type * 1000; int this_enum = Reweight::kNoDialFound; // Not Found - std::cout << "Getting dial enum " << type << " " << name << std::endl; + QLOG(DEB, "Getting dial enum " << type << " " << name); // Select Types switch (type) { // NEUT DIAL TYPE case kNEUT: { #ifdef __NEUT_ENABLED__ int neut_enum = (int)neut::rew::NSyst::FromString(name); if (neut_enum != 0) { this_enum = neut_enum + offset; } #else this_enum = Reweight::kNoTypeFound; // Not enabled #endif break; } // NIWG DIAL TYPE case kNIWG: { #ifdef __NIWG_ENABLED__ int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name); if (niwg_enum != 0) { this_enum = niwg_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } // NUWRO DIAL TYPE case kNUWRO: { #ifdef __NUWRO_REWEIGHT_ENABLED__ int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name); if (nuwro_enum > 0) { this_enum = nuwro_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif } // GENIE DIAL TYPE case kGENIE: { #ifdef __GENIE_ENABLED__ int genie_enum = (int)genie::rew::GSyst::FromString(name); if (genie_enum > 0) { this_enum = genie_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } case kCUSTOM: { int custom_enum = 0; // PLACEHOLDER this_enum = custom_enum + offset; break; } // T2K DIAL TYPE case kT2K: { #ifdef __T2KREW_ENABLED__ int t2k_enum = (int)t2krew::T2KSyst::FromString(name); if (t2k_enum > 0) { this_enum = t2k_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } case kNORM: { if (gNormEnums.find(name) == gNormEnums.end()) { gNormEnums[name] = gNormEnums.size() + 1 + offset; } this_enum = gNormEnums[name]; break; } case kLIKEWEIGHT: { if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) { gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset; } this_enum = gLikeWeightEnums[name]; break; } case kSPLINEPARAMETER: { if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) { gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset; } this_enum = gSplineParameterEnums[name]; break; } case kOSCILLATION: { #ifdef __PROB3PP_ENABLED__ int oscEnum = OscWeightEngine::SystEnumFromString(name); if (oscEnum != 0) { this_enum = oscEnum + offset; } #else this_enum = Reweight::kNoTypeFound; // Not enabled #endif } case kMODENORM: { size_t us_pos = name.find_first_of('_'); std::string numstr = name.substr(us_pos + 1); int mode_num = std::atoi(numstr.c_str()); LOG(FTL) << "Getting mode num " << mode_num << std::endl; if (!mode_num) { ERR(FTL) << "Attempting to parse dial name: \"" << name << "\" as a mode norm dial but failed." << std::endl; throw; } this_enum = 60 + mode_num + offset; break; } } // If Not Enabled if (this_enum == Reweight::kNoTypeFound) { ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type) << std::endl; ERR(FTL) << "Check dial " << name << std::endl; } // If Not Found if (this_enum == Reweight::kNoDialFound) { ERR(FTL) << "Dial " << name << " not found." << std::endl; } return this_enum; } int Reweight::ConvDialType(std::string const &type) { return FitBase::ConvDialType(type); } std::string Reweight::ConvDialType(int type) { return FitBase::ConvDialType(type); } int Reweight::GetDialType(int type) { int t = (type / 1000); return t > kMODENORM ? Reweight::kNoDialFound : t; } int Reweight::RemoveDialType(int type) { return (type % 1000); } int Reweight::NEUTEnumFromName(std::string const &name) { #ifdef __NEUT_ENABLED__ int neutenum = (int)neut::rew::NSyst::FromString(name); return (neutenum > 0) ? neutenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::NIWGEnumFromName(std::string const &name) { #ifdef __NIWG_ENABLED__ int niwgenum = (int)niwg::rew::NIWGSyst::FromString(name); return (niwgenum != 0) ? niwgenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::NUWROEnumFromName(std::string const &name) { #ifdef __NUWRO_REWEIGHT_ENABLED__ int nuwroenum = (int)nuwro::rew::NuwroSyst::FromString(name); return (nuwroenum > 0) ? nuwroenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::GENIEEnumFromName(std::string const &name) { #ifdef __GENIE_ENABLED__ int genieenum = (int)genie::rew::GSyst::FromString(name); return (genieenum > 0) ? genieenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::T2KEnumFromName(std::string const &name) { #ifdef __T2KREW_ENABLED__ int t2kenum = (int)t2krew::T2KSyst::FromString(name); return (t2kenum > 0) ? t2kenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::OscillationEnumFromName(std::string const &name) { #ifdef __PROB3PP_ENABLED__ int oscEnum = OscWeightEngine::SystEnumFromString(name); return (oscEnum > 0) ? oscEnum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::NUISANCEEnumFromName(std::string const &name, int type) { int nuisenum = Reweight::DialList().EnumFromNameAndType(name, type); return nuisenum; } int Reweight::CustomEnumFromName(std::string const &name) { int custenum = Reweight::ConvertNUISANCEDial(name); return custenum; } int Reweight::ConvDial(std::string const &name, std::string const &type, bool exceptions) { return Reweight::ConvDial(name, Reweight::ConvDialType(type), exceptions); } int Reweight::ConvDial(std::string const &fullname, int type, bool exceptions) { std::string name = GeneralUtils::ParseToStr(fullname, ",")[0]; // Only use first dial given // Produce offset seperating each type. int offset = type * 1000; int genenum = Reweight::kNoDialFound; switch (type) { case kNEUT: genenum = NEUTEnumFromName(name); break; case kNIWG: genenum = NIWGEnumFromName(name); break; case kNUWRO: genenum = NUWROEnumFromName(name); break; case kGENIE: genenum = GENIEEnumFromName(name); break; case kT2K: genenum = T2KEnumFromName(name); break; case kCUSTOM: genenum = CustomEnumFromName(name); break; case kNORM: case kLIKEWEIGHT: case kSPLINEPARAMETER: case kNEWSPLINE: genenum = NUISANCEEnumFromName(name, type); break; case kOSCILLATION: genenum = OscillationEnumFromName(name); break; case kMODENORM: genenum = ModeNormEngine::SystEnumFromString(name); break; default: genenum = Reweight::kNoTypeFound; break; } // Throw if required. if (exceptions) { // If Not Enabled if (genenum == Reweight::kGeneratorNotBuilt) { ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type) << std::endl; ERR(FTL) << "Check dial " << name << std::endl; throw; } // If no type enabled if (genenum == Reweight::kNoTypeFound) { ERR(FTL) << "Type mismatch inside ConvDialEnum" << std::endl; throw; } // If Not Found if (genenum == Reweight::kNoDialFound) { ERR(FTL) << "Dial " << name << " not found." << std::endl; throw; } } // Add offset if no issue int nuisenum = genenum; if ((genenum != Reweight::kGeneratorNotBuilt) && (genenum != Reweight::kNoTypeFound) && (genenum != Reweight::kNoDialFound)) { nuisenum += offset; } // Now register dial Reweight::DialList().RegisterDialEnum(name, type, nuisenum); return nuisenum; } std::string Reweight::ConvDial(int nuisenum) { for (size_t i = 0; i < Reweight::DialList().fAllDialEnums.size(); i++) { if (Reweight::DialList().fAllDialEnums[i] == nuisenum) { return Reweight::DialList().fAllDialNames[i]; } } LOG(FIT) << "Cannot find dial with enum = " << nuisenum << std::endl; return ""; } diff --git a/src/Routines/MinimizerRoutines.cxx b/src/Routines/MinimizerRoutines.cxx index 369cf45..9e30dd4 100755 --- a/src/Routines/MinimizerRoutines.cxx +++ b/src/Routines/MinimizerRoutines.cxx @@ -1,1511 +1,1517 @@ // 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()) { ERR(FTL) << "No input supplied!" << std::endl; throw; } if (fOutputFile.empty() and !fCardFile.empty()) { fOutputFile = fCardFile + ".root"; ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl; } else if (fOutputFile.empty()) { ERR(FTL) << "No output file or cardfile supplied!" << std::endl; throw; } // 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"); std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl; std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl; // 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() { //************************************* LOG(FIT) << "Setting up nuismin" << std::endl; // Setup Parameters ------------------------------------------ std::vector parkeys = Config::QueryKeys("parameter"); if (!parkeys.empty()) { LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl; } for (size_t i = 0; i < parkeys.size(); i++) { nuiskey key = parkeys.at(i); // Check for type,name,nom if (!key.Has("type")) { ERR(FTL) << "No type given for parameter " << i << std::endl; throw; } else if (!key.Has("name")) { ERR(FTL) << "No name given for parameter " << i << std::endl; throw; } else if (!key.Has("nominal")) { ERR(FTL) << "No nominal given for parameter " << i << std::endl; throw; } // 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"); LOG(FIT) << "Read " << partype << " : " << parname << " = " << parnom << " : " << parlow << " < p < " << parhigh << " : " << parstate << std::endl; } else { LOG(FIT) << "Read " << partype << " : " << parname << " = " << parnom << " : " << parstate << std::endl; } // Run Parameter Conversion if needed if (parstate.find("ABS") != std::string::npos) { 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; std::cout << "ParStep: " << parstep << " (" << oparstep << ")." << std::endl; } else if (parstate.find("FRAC") != std::string::npos) { 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()) { LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl; } 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 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 << std::endl; // If FREE add to parameters otherwise continue if (sampletype.find("FREE") == std::string::npos) { 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()) { LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl; } for (size_t i = 0; i < fakekeys.size(); i++) { nuiskey key = fakekeys.at(i); // Check for type,name,nom if (!key.Has("name")) { ERR(FTL) << "No name given for fakeparameter " << i << std::endl; throw; } else if (!key.Has("nom")) { ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl; throw; } // 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() { //************************************* LOG(FIT) << "Making the jointFCN" << std::endl; 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("minimizer.maxcalls")); if (!routine.compare("Brute")) { fMinimizer->SetMaxFunctionCalls(fParams.size() * fParams.size() * 4); fMinimizer->SetMaxIterations(fParams.size() * fParams.size() * 4); } fMinimizer->SetMaxIterations( FitPar::Config().GetParI("minimizer.maxiterations")); fMinimizer->SetTolerance(FitPar::Config().GetParD("minimizer.tolerance")); fMinimizer->SetStrategy(FitPar::Config().GetParI("minimizer.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 (fixed) { fMinimizer->FixVariable(ipar); LOG(FIT) << "Fixed Param: " << syst << std::endl; } else { LOG(FIT) << "Free Param: " << syst << " Start:" << vstart << " Range:" << vlow << " to " << vhigh << " Step:" << vstep << std::endl; } ipar++; } LOG(FIT) << "Setup Minimizer: " << fMinimizer->NDim() << "(NDim) " << fMinimizer->NFree() << "(NFree)" << std::endl; 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) { LOG(FIT) << "Setting fake data from MC starting prediction." << std::endl; // 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); LOG(FIT) << "Set all data to fake MC predictions." << std::endl; } 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() { //************************************* LOG(FIT) << "Running MinimizerRoutines : " << fStrategy << std::endl; if (FitPar::Config().GetParB("save_nominal")) { SaveNominal(); } // Parse given routines fRoutines = GeneralUtils::ParseToStr(fStrategy, ","); if (fRoutines.empty()) { ERR(FTL) << "Trying to run MinimizerRoutines with no routines given!" << std::endl; throw; } for (UInt_t i = 0; i < fRoutines.size(); i++) { std::string routine = fRoutines.at(i); int fitstate = kFitUnfinished; LOG(FIT) << "Running Routine: " << routine << std::endl; // 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) { LOG(FIT) << "Ending fit routines loop." << std::endl; 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) { LOG(FIT) << fMinimizer->Minimize() << std::endl; GetMinimizerState(); } } // other otptions else if (!routine.compare("Contour")) { CreateContours(); } return endfits; } //************************************* void MinimizerRoutines::PrintState() { //************************************* LOG(FIT) << "------------" << std::endl; // Count max size int maxcount = 0; for (UInt_t i = 0; i < fParams.size(); i++) { maxcount = max(int(fParams[i].size()), maxcount); } // Header 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::endl; // 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; LOG(FIT) << curparstring.str() << std::endl; } LOG(FIT) << "------------" << std::endl; double like = fSampleFCN->GetLikelihood(); LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl; LOG(FIT) << "------------" << std::endl; } //************************************* void MinimizerRoutines::GetMinimizerState() { //************************************* LOG(FIT) << "Minimizer State: " << std::endl; // 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 std::cout << "Filling covariance" << std::endl; 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"); } } std::cout << "Got STATE" << std::endl; return; }; //************************************* void MinimizerRoutines::LowStatRoutine(std::string routine) { //************************************* LOG(FIT) << "Running Low Statistics Routine: " << routine << std::endl; int lowstatsevents = FitPar::Config().GetParI("minimizer.lowstatevents"); int maxevents = FitPar::Config().GetParI("input.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("input.maxevents", lowstatsevents); Config::SetPar("VERBOSITY", 3); SetupFCN(); RunFitRoutine(trueroutine); Config::SetPar("input.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; LOG(FIT) << "Running 1D Scan for " << fParams[i] << std::endl; 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 LOG(FIT) << "Running scan for " << fParams[i] << " " << fParams[j] << std::endl; // 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 ERR(FTL) << " Contours not yet implemented as it is really slow!" << std::endl; throw; 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) { LOG(FIT) << "No dials needed fixing!" << std::endl; return kNoChange; } else return kStateChange; } /* Write Functions */ //************************************* void MinimizerRoutines::SaveResults() { //************************************* fOutputRootFile->cd(); if (fMinimizer) { SetupCovariance(); SaveMinimizerState(); } SaveCurrentState(); } //************************************* void MinimizerRoutines::SaveMinimizerState() { //************************************* std::cout << "Saving Minimizer State" << std::endl; if (!fMinimizer) { ERR(FTL) << "Can't save minimizer state without min object" << std::endl; throw; } // 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) { //************************************* LOG(FIT) << "Saving current full FCN predictions" << std::endl; // 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(); LOG(FIT) << "Saving Nominal Predictions (be cautious with this)" << std::endl; FitBase::GetRW()->Reconfigure(); SaveCurrentState("nominal"); }; //************************************* void MinimizerRoutines::SavePrefit() { //************************************* fOutputRootFile->cd(); LOG(FIT) << "Saving Prefit Predictions" << std::endl; 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; LOG(FIT) << "Building covariance matrix.." << std::endl; int NFREE = 0; int NDIM = 0; // Get NFREE from min or from vals (for cases when doing throws) if (fMinimizer) { std::cout << "NFREE FROM MINIMIZER" << std::endl; NFREE = fMinimizer->NFree(); NDIM = fMinimizer->NDim(); } else { NDIM = fParams.size(); for (UInt_t i = 0; i < fParams.size(); i++) { std::cout << "Getting Param " << fParams[i] << std::endl; if (!fFixVals[fParams[i]]) NFREE++; } } if (NDIM == 0) return; LOG(FIT) << "NFREE == " << NFREE << std::endl; 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++) { std::cout << "Getting Param " << i << std::endl; std::cout << "ParamI = " << fParams[i] << std::endl; 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++; } } std::cout << "Filling Matrices" << std::endl; 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; } std::cout << " Set the covariance" << std::endl; return; }; //************************************* void MinimizerRoutines::ThrowCovariance(bool uniformly) { //************************************* std::vector rands; if (!fDecFree) { ERR(WRN) << "Trying to throw 0 free parameters" << std::endl; 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() { LOG(FIT) << "Generating Toy Data Throws" << std::endl; 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 LOG(FIT) << "Writing toy data throws" << std::endl; fOutputRootFile->cd(); likes->Write(); } diff --git a/src/Routines/MinimizerRoutines.h b/src/Routines/MinimizerRoutines.h index 86b2d1d..3918e7a 100755 --- a/src/Routines/MinimizerRoutines.h +++ b/src/Routines/MinimizerRoutines.h @@ -1,276 +1,278 @@ // 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; //! 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 diff --git a/src/Statistical/StatUtils.cxx b/src/Statistical/StatUtils.cxx index 2648470..e8d23b1 100644 --- a/src/Statistical/StatUtils.cxx +++ b/src/Statistical/StatUtils.cxx @@ -1,1303 +1,1304 @@ // 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 "TH1D.h" #include "StatUtils.h" -#include "NuisConfig.h" #include "GeneralUtils.h" +#include "NuisConfig.h" +#include "TH1D.h" //******************************************************************* Double_t StatUtils::GetChi2FromDiag(TH1D* data, TH1D* mc, TH1I* mask) { -//******************************************************************* + //******************************************************************* Double_t Chi2 = 0.0; TH1D* calc_data = (TH1D*)data->Clone(); - TH1D* calc_mc = (TH1D*)mc->Clone(); + TH1D* calc_mc = (TH1D*)mc->Clone(); // Add MC Error to data if required if (FitPar::Config().GetParB("statutils.addmcerror")) { for (int i = 0; i < calc_data->GetNbinsX(); i++) { - double dterr = calc_data->GetBinError(i + 1); double mcerr = calc_mc->GetBinError(i + 1); if (dterr > 0.0) { calc_data->SetBinError(i + 1, sqrt(dterr * dterr + mcerr * mcerr)); } } } // Apply masking if required if (mask) { calc_data = ApplyHistogramMasking(data, mask); - calc_mc = ApplyHistogramMasking(mc, mask); + calc_mc = ApplyHistogramMasking(mc, mask); } // Iterate over bins in X for (int i = 0; i < calc_data->GetNbinsX(); i++) { - // Ignore bins with zero data or zero bin error if (calc_data->GetBinError(i + 1) <= 0.0 || - calc_data->GetBinContent(i + 1) == 0.0) continue; + calc_data->GetBinContent(i + 1) == 0.0) + continue; // Take mc data difference - double diff = calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1); + double diff = + calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1); double err = calc_data->GetBinError(i + 1); Chi2 += (diff * diff) / (err * err); - } // cleanup delete calc_data; delete calc_mc; return Chi2; }; //******************************************************************* -Double_t StatUtils::GetChi2FromDiag(TH2D* data, TH2D* mc, - TH2I* map, TH2I* mask) { -//******************************************************************* +Double_t StatUtils::GetChi2FromDiag(TH2D* data, TH2D* mc, TH2I* map, + TH2I* mask) { + //******************************************************************* // Generate a simple map - if (!map) map = GenerateMap(data); + if (!map) map = GenerateMap(data); // Convert to 1D Histograms TH1D* data_1D = MapToTH1D(data, map); - TH1D* mc_1D = MapToTH1D(mc, map); + TH1D* mc_1D = MapToTH1D(mc, map); TH1I* mask_1D = MapToMask(mask, map); // Calculate 1D chi2 from 1D Plots - Double_t Chi2 = StatUtils:: GetChi2FromDiag(data_1D, mc_1D, mask_1D); + Double_t Chi2 = StatUtils::GetChi2FromDiag(data_1D, mc_1D, mask_1D); // CleanUp delete data_1D; delete mc_1D; delete mask_1D; return Chi2; }; //******************************************************************* -Double_t StatUtils::GetChi2FromCov(TH1D* data, TH1D* mc, - TMatrixDSym* invcov, TH1I* mask, - double data_scale, double covar_scale) { -//******************************************************************* +Double_t StatUtils::GetChi2FromCov(TH1D* data, TH1D* mc, TMatrixDSym* invcov, + TH1I* mask, double data_scale, + double covar_scale) { + //******************************************************************* Double_t Chi2 = 0.0; - TMatrixDSym* calc_cov = (TMatrixDSym*) invcov->Clone(); - TH1D* calc_data = (TH1D*) data->Clone(); - TH1D* calc_mc = (TH1D*) mc->Clone(); + TMatrixDSym* calc_cov = (TMatrixDSym*)invcov->Clone(); + TH1D* calc_data = (TH1D*)data->Clone(); + TH1D* calc_mc = (TH1D*)mc->Clone(); // If a mask if applied we need to apply it before the matrix is inverted if (mask) { - calc_cov = ApplyInvertedMatrixMasking(invcov, mask); + calc_cov = ApplyInvertedMatrixMasking(invcov, mask); calc_data = ApplyHistogramMasking(data, mask); - calc_mc = ApplyHistogramMasking(mc, mask); + calc_mc = ApplyHistogramMasking(mc, mask); } // Add MC Error to data if required if (FitPar::Config().GetParB("statutils.addmcerror")) { - // Make temp cov TMatrixDSym* newcov = StatUtils::GetInvert(calc_cov); // Add MC err to diag for (int i = 0; i < calc_data->GetNbinsX(); i++) { - double mcerr = calc_mc->GetBinError(i + 1) * sqrt(covar_scale); double oldval = (*newcov)(i, i); - std::cout << "Adding cov stat " << mcerr*mcerr << " to " << (*newcov)(i,i) << std::endl; + std::cout << "Adding cov stat " << mcerr * mcerr << " to " + << (*newcov)(i, i) << std::endl; (*newcov)(i, i) = oldval + mcerr * mcerr; } // Reset the calc_cov to new invert delete calc_cov; calc_cov = GetInvert(newcov); // Delete the tempcov delete newcov; } calc_data->Scale(data_scale); - calc_mc ->Scale(data_scale); + calc_mc->Scale(data_scale); (*calc_cov) *= covar_scale; // iterate over bins in X (i,j) + QLOG(DEB, "START Chi2 Calculation================="); for (int i = 0; i < calc_data->GetNbinsX(); i++) { - + QLOG(DEB, + "[CHI2] i = " << i << " [" + << calc_data->GetXaxis()->GetBinLowEdge(i + 1) << " -- " + << calc_data->GetXaxis()->GetBinUpEdge(i + 1) << "]."); for (int j = 0; j < calc_data->GetNbinsX(); j++) { - - if (calc_data->GetBinContent(i + 1) != 0 || calc_mc->GetBinContent(i + 1) != 0) { - - LOG(DEB) << "i j = " << i << " " << j << std::endl; - LOG(DEB) << "Calc_data mc i = " << calc_data->GetBinContent(i + 1) - << " " << calc_mc->GetBinContent(i + 1) - << " Dif = " - << ( calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1) ) - << std::endl; - - LOG(DEB) << "Calc_data mc i = " << calc_data->GetBinContent(j + 1) - << " " << calc_mc->GetBinContent(j + 1) - << " Dif = " - << ( calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1) ) - << std::endl; - - LOG(DEB) << "Covar = " << (*calc_cov)(i, j) << std::endl; - - LOG(DEB) << "Cont chi2 = " \ - << ( ( calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1) ) \ - * (*calc_cov)(i, j) \ - * ( calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1))) - << " " << Chi2 << std::endl; - - Chi2 += ( ( calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1) ) \ - * (*calc_cov)(i, j) \ - * ( calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1) ) ); + QLOG(DEB, "[CHI2]\t j = " + << i << " [" << calc_data->GetXaxis()->GetBinLowEdge(j + 1) + << " -- " << calc_data->GetXaxis()->GetBinUpEdge(j + 1) + << "]."); + if ((calc_data->GetBinContent(i + 1) != 0 || + calc_mc->GetBinContent(i + 1) != 0) && + ((*calc_cov)(i, j) != 0)) { + QLOG(DEB, "[CHI2]\t\t Chi2 contribution (i,j) = (" << i << "," << j + << ")"); + QLOG(DEB, "[CHI2]\t\t Data - MC(i) = " + << calc_data->GetBinContent(i + 1) << " - " + << calc_mc->GetBinContent(i + 1) << " = " + << (calc_data->GetBinContent(i + 1) - + calc_mc->GetBinContent(i + 1))); + + QLOG(DEB, "[CHI2]\t\t Data - MC(j) = " + << calc_data->GetBinContent(j + 1) << " - " + << calc_mc->GetBinContent(j + 1) << " = " + << (calc_data->GetBinContent(j + 1) - + calc_mc->GetBinContent(j + 1))); + + QLOG(DEB, "[CHI2]\t\t Covar = " << (*calc_cov)(i, j)); + + QLOG(DEB, "[CHI2]\t\t Cont chi2 = " + << ((calc_data->GetBinContent(i + 1) - + calc_mc->GetBinContent(i + 1)) * + (*calc_cov)(i, j) * (calc_data->GetBinContent(j + 1) - + calc_mc->GetBinContent(j + 1))) + << " " << Chi2); + + Chi2 += + ((calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1)) * + (*calc_cov)(i, j) * + (calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1))); } else { - - LOG(DEB) << "Error on bin (i,j) = (" << i << "," << j << ")" << std::endl; - LOG(DEB) << "data->GetBinContent(i+1) = " << calc_data->GetBinContent(i + 1) << std::endl; - LOG(DEB) << "mc->GetBinContent(i+1) = " << calc_mc->GetBinContent(i + 1) << std::endl; - LOG(DEB) << "Adding zero to chi2 instead of dying horrifically " << std::endl; - + QLOG(DEB, "Skipping chi2 contribution (i,j) = (" + << i << "," << j + << "), Data = " << calc_data->GetBinContent(i + 1) + << ", MC = " << calc_mc->GetBinContent(i + 1) + << ", Cov = " << (*calc_cov)(i, j)); Chi2 += 0.; } } } // Cleanup delete calc_cov; delete calc_data; delete calc_mc; return Chi2; } - -//******************************************************************* -Double_t StatUtils::GetChi2FromCov( TH2D* data, TH2D* mc, - TMatrixDSym* invcov, TH2I* map, TH2I* mask) { //******************************************************************* +Double_t StatUtils::GetChi2FromCov(TH2D* data, TH2D* mc, TMatrixDSym* invcov, + TH2I* map, TH2I* mask) { + //******************************************************************* // Generate a simple map if (!map) { map = StatUtils::GenerateMap(data); } // Convert to 1D Histograms TH1D* data_1D = MapToTH1D(data, map); - TH1D* mc_1D = MapToTH1D(mc, map); + TH1D* mc_1D = MapToTH1D(mc, map); TH1I* mask_1D = MapToMask(mask, map); // Calculate 1D chi2 from 1D Plots - Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D); + Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D); // CleanUp delete data_1D; delete mc_1D; delete mask_1D; return Chi2; } //******************************************************************* -Double_t StatUtils::GetChi2FromSVD( TH1D* data, TH1D* mc, - TMatrixDSym* cov, TH1I* mask) { -//******************************************************************* +Double_t StatUtils::GetChi2FromSVD(TH1D* data, TH1D* mc, TMatrixDSym* cov, + TH1I* mask) { + //******************************************************************* Double_t Chi2 = 0.0; - TMatrixDSym* calc_cov = (TMatrixDSym*) cov->Clone(); - TH1D* calc_data = (TH1D*) data->Clone(); - TH1D* calc_mc = (TH1D*) mc->Clone(); + TMatrixDSym* calc_cov = (TMatrixDSym*)cov->Clone(); + TH1D* calc_data = (TH1D*)data->Clone(); + TH1D* calc_mc = (TH1D*)mc->Clone(); // If a mask if applied we need to apply it before the matrix is inverted if (mask) { - calc_cov = StatUtils::ApplyMatrixMasking(cov, mask); + calc_cov = StatUtils::ApplyMatrixMasking(cov, mask); calc_data = StatUtils::ApplyHistogramMasking(data, mask); - calc_mc = StatUtils::ApplyHistogramMasking(mc, mask); + calc_mc = StatUtils::ApplyHistogramMasking(mc, mask); } // Decompose matrix TDecompSVD LU = TDecompSVD((*calc_cov)); LU.Decompose(); - TMatrixDSym* cov_U = new TMatrixDSym(calc_data->GetNbinsX(), LU .GetU().GetMatrixArray(), ""); - TVectorD* cov_S = new TVectorD( LU.GetSig() ); + TMatrixDSym* cov_U = + new TMatrixDSym(calc_data->GetNbinsX(), LU.GetU().GetMatrixArray(), ""); + TVectorD* cov_S = new TVectorD(LU.GetSig()); // Apply basis rotation before adding up chi2 Double_t rotated_difference = 0.0; for (int i = 0; i < calc_data->GetNbinsX(); i++) { rotated_difference = 0.0; // Rotate basis of Data - MC for (int j = 0; j < calc_data->GetNbinsY(); j++) - rotated_difference += ( calc_data->GetBinContent(j + 1) - calc_mc ->GetBinContent(j + 1) ) * (*cov_U)(j, i) ; + rotated_difference += + (calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1)) * + (*cov_U)(j, i); // Divide by rotated error cov_S Chi2 += rotated_difference * rotated_difference * 1E76 / (*cov_S)(i); - } // Cleanup delete calc_cov; delete calc_data; delete calc_mc; delete cov_U; delete cov_S; return Chi2; } - -//******************************************************************* -Double_t StatUtils::GetChi2FromSVD( TH2D* data, TH2D* mc, - TMatrixDSym* cov, TH2I* map, TH2I* mask) { //******************************************************************* +Double_t StatUtils::GetChi2FromSVD(TH2D* data, TH2D* mc, TMatrixDSym* cov, + TH2I* map, TH2I* mask) { + //******************************************************************* // Generate a simple map - if (!map) - map = StatUtils::GenerateMap(data); + if (!map) map = StatUtils::GenerateMap(data); // Convert to 1D Histograms TH1D* data_1D = MapToTH1D(data, map); - TH1D* mc_1D = MapToTH1D(mc, map); + TH1D* mc_1D = MapToTH1D(mc, map); TH1I* mask_1D = MapToMask(mask, map); // Calculate from 1D Double_t Chi2 = StatUtils::GetChi2FromSVD(data_1D, mc_1D, cov, mask_1D); // CleanUp delete data_1D; delete mc_1D; delete mask_1D; return Chi2; } //******************************************************************* double StatUtils::GetChi2FromEventRate(TH1D* data, TH1D* mc, TH1I* mask) { -//******************************************************************* + //******************************************************************* - // If just an event rate, for chi2 just use Poission Likelihood to calculate the chi2 component + // If just an event rate, for chi2 just use Poission Likelihood to calculate + // the chi2 component double chi2 = 0.0; TH1D* calc_data = (TH1D*)data->Clone(); - TH1D* calc_mc = (TH1D*)mc->Clone(); + TH1D* calc_mc = (TH1D*)mc->Clone(); // Apply masking if required if (mask) { calc_data = ApplyHistogramMasking(data, mask); - calc_mc = ApplyHistogramMasking(mc, mask); + calc_mc = ApplyHistogramMasking(mc, mask); } // Iterate over bins in X for (int i = 0; i < calc_data->GetNbinsX(); i++) { - double dt = calc_data->GetBinContent(i + 1); double mc = calc_mc->GetBinContent(i + 1); if (mc <= 0) continue; if (dt <= 0) { // Only add difference chi2 += 2 * (mc - dt); } else { // Do the chi2 for Poisson distributions - chi2 += 2 * (mc - dt + (dt * log(dt / mc))); + chi2 += 2 * (mc - dt + (dt * log(dt / mc))); } /* LOG(REC)<<"Evt Chi2 cont = "<Clone(); // If a mask is provided we need to apply it before getting NDOF if (mask) { calc_hist = StatUtils::ApplyHistogramMasking(hist, mask); } // NDOF is defined as total number of bins with non-zero errors Int_t NDOF = 0; for (int i = 0; i < calc_hist->GetNbinsX(); i++) { if (calc_hist->GetBinError(i + 1) > 0.0) NDOF++; } delete calc_hist; return NDOF; }; - //******************************************************************* Int_t StatUtils::GetNDOF(TH2D* hist, TH2I* map, TH2I* mask) { -//******************************************************************* + //******************************************************************* Int_t NDOF = 0; if (!map) map = StatUtils::GenerateMap(hist); for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { - if (mask->GetBinContent(i + 1, j + 1)) continue; if (map->GetBinContent(i + 1, j + 1) <= 0) continue; NDOF++; - } } return NDOF; }; - - -//******************************************************************* -TH1D* StatUtils::ThrowHistogram(TH1D* hist, TMatrixDSym* cov, bool throwdiag, TH1I* mask) { //******************************************************************* +TH1D* StatUtils::ThrowHistogram(TH1D* hist, TMatrixDSym* cov, bool throwdiag, + TH1I* mask) { + //******************************************************************* - TH1D* calc_hist = (TH1D*) hist->Clone( (std::string(hist->GetName()) + "_THROW" ).c_str() ); - TMatrixDSym* calc_cov = (TMatrixDSym*) cov->Clone(); + TH1D* calc_hist = + (TH1D*)hist->Clone((std::string(hist->GetName()) + "_THROW").c_str()); + TMatrixDSym* calc_cov = (TMatrixDSym*)cov->Clone(); Double_t correl_val = 0.0; // If a mask if applied we need to apply it before the matrix is decomposed if (mask) { - calc_cov = ApplyMatrixMasking(cov, mask); + calc_cov = ApplyMatrixMasking(cov, mask); calc_hist = ApplyHistogramMasking(calc_hist, mask); } // If a covariance is provided we need a preset random vector and a decomp std::vector rand_val; TMatrixDSym* decomp_cov; if (cov) { for (int i = 0; i < hist->GetNbinsX(); i++) { rand_val.push_back(gRandom->Gaus(0.0, 1.0)); } // Decomp the matrix decomp_cov = StatUtils::GetDecomp(calc_cov); } // iterate over bins for (int i = 0; i < hist->GetNbinsX(); i++) { - - // By Default the errors on the histogram are thrown uncorrelated to the other errors + // By Default the errors on the histogram are thrown uncorrelated to the + // other errors // if (throwdiag) { // calc_hist->SetBinContent(i + 1, (calc_hist->GetBinContent(i + 1) + \ // gRandom->Gaus(0.0, 1.0) * calc_hist->GetBinError(i + 1)) ); // } // If a covariance is provided that is also thrown if (cov) { correl_val = 0.0; for (int j = 0; j < hist->GetNbinsX(); j++) { - correl_val += rand_val[j] * (*decomp_cov)(j, i) ; + correl_val += rand_val[j] * (*decomp_cov)(j, i); } - calc_hist->SetBinContent(i + 1, (calc_hist->GetBinContent(i + 1) + correl_val * 1E-38)); + calc_hist->SetBinContent( + i + 1, (calc_hist->GetBinContent(i + 1) + correl_val * 1E-38)); } } delete calc_cov; delete decomp_cov; // return this new thrown data return calc_hist; }; //******************************************************************* -TH2D* StatUtils::ThrowHistogram(TH2D* hist, TMatrixDSym* cov, TH2I* map, bool throwdiag, TH2I* mask) { -//******************************************************************* +TH2D* StatUtils::ThrowHistogram(TH2D* hist, TMatrixDSym* cov, TH2I* map, + bool throwdiag, TH2I* mask) { + //******************************************************************* // PLACEHOLDER!!!!!!!!! // Currently no support for throwing 2D Histograms from a covariance - (void) hist; - (void) cov; - (void) map; - (void) throwdiag; - (void) mask; + (void)hist; + (void)cov; + (void)map; + (void)throwdiag; + (void)mask; // /todo // Sort maps if required // Throw the covariance for a 1D plot // Unmap back to 2D Histogram return hist; } - - - - - //******************************************************************* TH1D* StatUtils::ApplyHistogramMasking(TH1D* hist, TH1I* mask) { -//******************************************************************* + //******************************************************************* - if (!mask) return ( (TH1D*)hist->Clone() ); + if (!mask) return ((TH1D*)hist->Clone()); - // This masking is only sufficient for chi2 calculations, and will have dodgy bin edges. + // This masking is only sufficient for chi2 calculations, and will have dodgy + // bin edges. // Get New Bin Count Int_t NBins = 0; for (int i = 0; i < hist->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1)) continue; NBins++; } // Make new hist std::string newmaskname = std::string(hist->GetName()) + "_MSKD"; - TH1D* calc_hist = new TH1D( newmaskname.c_str(), newmaskname.c_str(), NBins, 0, NBins); + TH1D* calc_hist = + new TH1D(newmaskname.c_str(), newmaskname.c_str(), NBins, 0, NBins); // fill new hist int binindex = 0; for (int i = 0; i < hist->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1)) { - LOG(REC) << "Applying mask to bin " << i + 1 << " " << hist->GetName() << std::endl; + LOG(REC) << "Applying mask to bin " << i + 1 << " " << hist->GetName() + << std::endl; continue; } calc_hist->SetBinContent(binindex + 1, hist->GetBinContent(i + 1)); calc_hist->SetBinError(binindex + 1, hist->GetBinError(i + 1)); binindex++; } return calc_hist; }; //******************************************************************* TH2D* StatUtils::ApplyHistogramMasking(TH2D* hist, TH2I* mask) { -//******************************************************************* + //******************************************************************* - TH2D* newhist = (TH2D*) hist->Clone(); + TH2D* newhist = (TH2D*)hist->Clone(); if (!mask) return newhist; for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { - if (mask->GetBinContent(i + 1, j + 1) > 0) { newhist->SetBinContent(i + 1, j + 1, 0.0); newhist->SetBinContent(i + 1, j + 1, 0.0); } } } return newhist; } //******************************************************************* TMatrixDSym* StatUtils::ApplyMatrixMasking(TMatrixDSym* mat, TH1I* mask) { -//******************************************************************* + //******************************************************************* if (!mask) return (TMatrixDSym*)(mat->Clone()); // Get New Bin Count Int_t NBins = 0; for (int i = 0; i < mask->GetNbinsX(); i++) { if (mask->GetBinContent(i + 1)) continue; NBins++; } // make new matrix TMatrixDSym* calc_mat = new TMatrixDSym(NBins); int col, row; // Need to mask out bins in the current matrix row = 0; for (int i = 0; i < mask->GetNbinsX(); i++) { col = 0; // skip if masked if (mask->GetBinContent(i + 1) > 0.5) continue; for (int j = 0; j < mask->GetNbinsX(); j++) { - // skip if masked if (mask->GetBinContent(j + 1) > 0.5) continue; (*calc_mat)(row, col) = (*mat)(i, j); col++; } row++; } return calc_mat; }; //******************************************************************* -TMatrixDSym* StatUtils::ApplyMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map) { -//******************************************************************* +TMatrixDSym* StatUtils::ApplyMatrixMasking(TMatrixDSym* mat, TH2D* data, + TH2I* mask, TH2I* map) { + //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1I* mask_1D = StatUtils::MapToMask(mask, map); - TMatrixDSym* newmat = StatUtils::ApplyMatrixMasking(mat, mask_1D); + TMatrixDSym* newmat = StatUtils::ApplyMatrixMasking(mat, mask_1D); delete mask_1D; return newmat; } - - -//******************************************************************* -TMatrixDSym* StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH1I* mask) { //******************************************************************* +TMatrixDSym* StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym* mat, + TH1I* mask) { + //******************************************************************* TMatrixDSym* new_mat = GetInvert(mat); TMatrixDSym* masked_mat = ApplyMatrixMasking(new_mat, mask); TMatrixDSym* inverted_mat = GetInvert(masked_mat); delete masked_mat; delete new_mat; return inverted_mat; }; - - -//******************************************************************* -TMatrixDSym* StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH2D* data, TH2I* mask, TH2I* map) { //******************************************************************* +TMatrixDSym* StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym* mat, TH2D* data, + TH2I* mask, TH2I* map) { + //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1I* mask_1D = StatUtils::MapToMask(mask, map); TMatrixDSym* newmat = ApplyInvertedMatrixMasking(mat, mask_1D); delete mask_1D; return newmat; } - - - //******************************************************************* TMatrixDSym* StatUtils::GetInvert(TMatrixDSym* mat) { -//******************************************************************* + //******************************************************************* TMatrixDSym* new_mat = (TMatrixDSym*)mat->Clone(); // Check for diagonal bool non_diagonal = false; for (int i = 0; i < new_mat->GetNrows(); i++) { for (int j = 0; j < new_mat->GetNrows(); j++) { if (i == j) continue; if ((*new_mat)(i, j) != 0.0) { non_diagonal = true; break; } } } // If diag, just flip the diag if (!non_diagonal or new_mat->GetNrows() == 1) { for (int i = 0; i < new_mat->GetNrows(); i++) { if ((*new_mat)(i, i) != 0.0) (*new_mat)(i, i) = 1.0 / (*new_mat)(i, i); else (*new_mat)(i, i) = 0.0; } return new_mat; } - // Invert full matrix TDecompSVD LU = TDecompSVD((*new_mat)); - new_mat = new TMatrixDSym(new_mat->GetNrows(), LU.Invert().GetMatrixArray(), ""); + new_mat = + new TMatrixDSym(new_mat->GetNrows(), LU.Invert().GetMatrixArray(), ""); return new_mat; } //******************************************************************* TMatrixDSym* StatUtils::GetDecomp(TMatrixDSym* mat) { -//******************************************************************* + //******************************************************************* TMatrixDSym* new_mat = (TMatrixDSym*)mat->Clone(); int nrows = new_mat->GetNrows(); // Check for diagonal bool diagonal = true; for (int i = 0; i < nrows; i++) { for (int j = 0; j < nrows; j++) { if (i == j) continue; if ((*new_mat)(i, j) != 0.0) { diagonal = false; break; } } } // If diag, just flip the diag if (diagonal or nrows == 1) { for (int i = 0; i < nrows; i++) { if ((*new_mat)(i, i) > 0.0) (*new_mat)(i, i) = sqrt((*new_mat)(i, i)); else (*new_mat)(i, i) = 0.0; } return new_mat; } TDecompChol LU = TDecompChol(*new_mat); LU.Decompose(); delete new_mat; TMatrixDSym* dec_mat = new TMatrixDSym(nrows, LU.GetU().GetMatrixArray(), ""); return dec_mat; } - - - - //******************************************************************* void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH1D* hist, double norm) { -//******************************************************************* + //******************************************************************* if (!mat) mat = MakeDiagonalCovarMatrix(hist); int nbins = mat->GetNrows(); TMatrixDSym* new_mat = new TMatrixDSym(nbins); for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { - double valx = hist->GetBinContent(i + 1) * 1E38; double valy = hist->GetBinContent(j + 1) * 1E38; (*new_mat)(i, j) = (*mat)(i, j) + norm * norm * valx * valy; - } } // Swap the two delete mat; mat = new_mat; return; }; //******************************************************************* -void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH2D* data, double norm, TH2I* map ) { -//******************************************************************* +void StatUtils::ForceNormIntoCovar(TMatrixDSym* mat, TH2D* data, double norm, + TH2I* map) { + //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1D* data_1D = MapToTH1D(data, map); StatUtils::ForceNormIntoCovar(mat, data_1D, norm); delete data_1D; return; } //******************************************************************* TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH1D* data, double scaleF) { -//******************************************************************* + //******************************************************************* TMatrixDSym* newmat = new TMatrixDSym(data->GetNbinsX()); for (int i = 0; i < data->GetNbinsX(); i++) { - (*newmat)(i, i) = data->GetBinError(i + 1) * data->GetBinError(i + 1) * scaleF * scaleF; + (*newmat)(i, i) = + data->GetBinError(i + 1) * data->GetBinError(i + 1) * scaleF * scaleF; } return newmat; } - - -//******************************************************************* -TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH2D* data, TH2I* map, double scaleF) { //******************************************************************* +TMatrixDSym* StatUtils::MakeDiagonalCovarMatrix(TH2D* data, TH2I* map, + double scaleF) { + //******************************************************************* if (!map) map = StatUtils::GenerateMap(data); TH1D* data_1D = MapToTH1D(data, map); return StatUtils::MakeDiagonalCovarMatrix(data_1D, scaleF); }; - -//******************************************************************* -void StatUtils::SetDataErrorFromCov(TH1D* data, TMatrixDSym* cov, double scale) { //******************************************************************* +void StatUtils::SetDataErrorFromCov(TH1D* data, TMatrixDSym* cov, + double scale) { + //******************************************************************* // Check if (cov->GetNrows() != data->GetNbinsX()) { - ERR(WRN) << "Nrows in cov don't match nbins in data for SetDataErrorFromCov" << std::endl; + ERR(WRN) << "Nrows in cov don't match nbins in data for SetDataErrorFromCov" + << std::endl; } // Set bin errors form cov diag for (int i = 0; i < data->GetNbinsX(); i++) { - data->SetBinError(i + 1, sqrt((*cov)(i, i)) * scale ); + data->SetBinError(i + 1, sqrt((*cov)(i, i)) * scale); } return; } //******************************************************************* -void StatUtils::SetDataErrorFromCov(TH2D* data, TMatrixDSym* cov, TH2I* map, double scale) { -//******************************************************************* +void StatUtils::SetDataErrorFromCov(TH2D* data, TMatrixDSym* cov, TH2I* map, + double scale) { + //******************************************************************* // Create map if required if (!map) map = StatUtils::GenerateMap(data); // Set Bin Errors from cov diag int count = 0; for (int i = 0; i < data->GetNbinsX(); i++) { for (int j = 0; j < data->GetNbinsY(); j++) { - if (data->GetBinContent(i + 1, j + 1) == 0.0) continue; count = map->GetBinContent(i + 1, j + 1) - 1; - data->SetBinError(i + 1, j + 1, sqrt((*cov)(count, count)) * scale ); - + data->SetBinError(i + 1, j + 1, sqrt((*cov)(count, count)) * scale); } } return; } - -TMatrixDSym* StatUtils::ExtractShapeOnlyCovar(TMatrixDSym* full_covar, TH1* data_hist, double data_scale){ - +TMatrixDSym* StatUtils::ExtractShapeOnlyCovar(TMatrixDSym* full_covar, + TH1* data_hist, + double data_scale) { int nbins = full_covar->GetNrows(); TMatrixDSym* shape_covar = new TMatrixDSym(nbins); // Check nobody is being silly - if (data_hist->GetNbinsX() != nbins){ - ERR(WRN) << "Inconsistent matrix and data histogram passed to StatUtils::ExtractShapeOnlyCovar!" << std::endl; - ERR(WRN) << "data_hist has " << data_hist->GetNbinsX() << " matrix has " << nbins << std::endl; + if (data_hist->GetNbinsX() != nbins) { + ERR(WRN) << "Inconsistent matrix and data histogram passed to " + "StatUtils::ExtractShapeOnlyCovar!" + << std::endl; + ERR(WRN) << "data_hist has " << data_hist->GetNbinsX() << " matrix has " + << nbins << std::endl; int err_bins = data_hist->GetNbinsX(); if (nbins > err_bins) err_bins = nbins; - for (int i = 0; i < err_bins; ++i){ - ERR(WRN) << "Matrix diag. = " << (*full_covar)(i, i) << " data = " << data_hist->GetBinContent(i+1) << std::endl; + for (int i = 0; i < err_bins; ++i) { + ERR(WRN) << "Matrix diag. = " << (*full_covar)(i, i) + << " data = " << data_hist->GetBinContent(i + 1) << std::endl; } return NULL; } - double total_data = 0; + double total_data = 0; double total_covar = 0; - + // Initial loop to calculate some constants for (int i = 0; i < nbins; ++i) { - total_data += data_hist->GetBinContent(i+1)*data_scale; + total_data += data_hist->GetBinContent(i + 1) * data_scale; for (int j = 0; j < nbins; ++j) { - total_covar += (*full_covar)(i,j); + total_covar += (*full_covar)(i, j); } } - - if (total_data == 0 || total_covar == 0){ - ERR(WRN) << "Stupid matrix or data histogram passed to StatUtils::ExtractShapeOnlyCovar! Ignoring..." << std::endl; + + if (total_data == 0 || total_covar == 0) { + ERR(WRN) << "Stupid matrix or data histogram passed to " + "StatUtils::ExtractShapeOnlyCovar! Ignoring..." + << std::endl; return NULL; } - LOG(SAM) << "Norm error = " << sqrt(total_covar)/total_data << std::endl; - + LOG(SAM) << "Norm error = " << sqrt(total_covar) / total_data << std::endl; + // Now loop over and calculate the shape-only matrix for (int i = 0; i < nbins; ++i) { - double data_i = data_hist->GetBinContent(i+1)*data_scale; + double data_i = data_hist->GetBinContent(i + 1) * data_scale; for (int j = 0; j < nbins; ++j) { - double data_j = data_hist->GetBinContent(j+1)*data_scale; - - double norm_term = data_i*data_j*total_covar/total_data/total_data; + double data_j = data_hist->GetBinContent(j + 1) * data_scale; + + double norm_term = + data_i * data_j * total_covar / total_data / total_data; double mix_sum1 = 0; double mix_sum2 = 0; - - for (int k = 0; k < nbins; ++k){ - mix_sum1 += (*full_covar)(k,j); - mix_sum2 += (*full_covar)(i,k); + + for (int k = 0; k < nbins; ++k) { + mix_sum1 += (*full_covar)(k, j); + mix_sum2 += (*full_covar)(i, k); } - double mix_term1 = data_i*(mix_sum1/total_data - total_covar*data_j/total_data/total_data); - double mix_term2 = data_j*(mix_sum2/total_data - total_covar*data_i/total_data/total_data); - - (*shape_covar)(i, j) = (*full_covar)(i, j) - mix_term1 - mix_term2 - norm_term; + double mix_term1 = + data_i * (mix_sum1 / total_data - + total_covar * data_j / total_data / total_data); + double mix_term2 = + data_j * (mix_sum2 / total_data - + total_covar * data_i / total_data / total_data); + + (*shape_covar)(i, j) = + (*full_covar)(i, j) - mix_term1 - mix_term2 - norm_term; } } return shape_covar; } - //******************************************************************* TH2I* StatUtils::GenerateMap(TH2D* hist) { -//******************************************************************* + //******************************************************************* std::string maptitle = std::string(hist->GetName()) + "_MAP"; - TH2I* map = new TH2I( maptitle.c_str(), maptitle.c_str(), - hist->GetNbinsX(), 0, hist->GetNbinsX(), - hist->GetNbinsY(), 0, hist->GetNbinsY()); + TH2I* map = + new TH2I(maptitle.c_str(), maptitle.c_str(), hist->GetNbinsX(), 0, + hist->GetNbinsX(), hist->GetNbinsY(), 0, hist->GetNbinsY()); Int_t index = 1; for (int i = 0; i < hist->GetNbinsX(); i++) { for (int j = 0; j < hist->GetNbinsY(); j++) { - - if (hist->GetBinContent(i + 1, j + 1) > 0 && hist->GetBinError(i + 1, j + 1) > 0) { - + if (hist->GetBinContent(i + 1, j + 1) > 0 && + hist->GetBinError(i + 1, j + 1) > 0) { map->SetBinContent(i + 1, j + 1, index); index++; } else { map->SetBinContent(i + 1, j + 1, 0); } } } return map; } //******************************************************************* TH1D* StatUtils::MapToTH1D(TH2D* hist, TH2I* map) { -//******************************************************************* + //******************************************************************* if (!hist) return NULL; // Get N bins for 1D plot Int_t Nbins = map->GetMaximum(); std::string name1D = std::string(hist->GetName()) + "_1D"; // Make new 1D Hist TH1D* newhist = new TH1D(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins); // map bin contents for (int i = 0; i < map->GetNbinsX(); i++) { for (int j = 0; j < map->GetNbinsY(); j++) { - if (map->GetBinContent(i + 1, j + 1) == 0) continue; - newhist->SetBinContent(map->GetBinContent(i + 1, j + 1), hist->GetBinContent(i + 1, j + 1)); - newhist->SetBinError(map->GetBinContent(i + 1, j + 1), hist->GetBinError(i + 1, j + 1)); + newhist->SetBinContent(map->GetBinContent(i + 1, j + 1), + hist->GetBinContent(i + 1, j + 1)); + newhist->SetBinError(map->GetBinContent(i + 1, j + 1), + hist->GetBinError(i + 1, j + 1)); } } // return return newhist; } - //******************************************************************* TH1I* StatUtils::MapToMask(TH2I* hist, TH2I* map) { -//******************************************************************* + //******************************************************************* TH1I* newhist = NULL; if (!hist) return newhist; // Get N bins for 1D plot Int_t Nbins = map->GetMaximum(); std::string name1D = std::string(hist->GetName()) + "_1D"; // Make new 1D Hist newhist = new TH1I(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins); // map bin contents for (int i = 0; i < map->GetNbinsX(); i++) { for (int j = 0; j < map->GetNbinsY(); j++) { - if (map->GetBinContent(i + 1, j + 1) == 0) continue; - newhist->SetBinContent(map->GetBinContent(i + 1, j + 1), hist->GetBinContent(i + 1, j + 1)); + newhist->SetBinContent(map->GetBinContent(i + 1, j + 1), + hist->GetBinContent(i + 1, j + 1)); } } // return return newhist; } - TMatrixDSym* StatUtils::GetCovarFromCorrel(TMatrixDSym* correl, TH1D* data) { - int nbins = correl->GetNrows(); TMatrixDSym* covar = new TMatrixDSym(nbins); for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { - (*covar)(i, j) = (*correl)(i, j) * data->GetBinError(i + 1) * data->GetBinError(j + 1); + (*covar)(i, j) = + (*correl)(i, j) * data->GetBinError(i + 1) * data->GetBinError(j + 1); } } return covar; } - - -//******************************************************************* -TMatrixD* StatUtils::GetMatrixFromTextFile(std::string covfile, int dimx, int dimy) { //******************************************************************* +TMatrixD* StatUtils::GetMatrixFromTextFile(std::string covfile, int dimx, + int dimy) { + //******************************************************************* // Determine dim if (dimx == -1 and dimy == -1) { std::string line; std::ifstream covar(covfile.c_str(), std::ifstream::in); int row = 0; while (std::getline(covar >> std::ws, line, '\n')) { int column = 0; std::vector entries = GeneralUtils::ParseToDbl(line, " "); if (entries.size() <= 1) { - ERR(WRN) << "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 " - "entries on this line: " << row << std::endl; + ERR(WRN) << "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 " + "entries on this line: " + << row << std::endl; } for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { column++; if (column > dimx) dimx = column; } row++; if (row > dimy) dimy = row; } } // Or assume symmetric if (dimx != -1 and dimy == -1) { dimy = dimx; } assert(dimy != -1 && " matrix dimy not set."); // Make new matrix TMatrixD* mat = new TMatrixD(dimx, dimy); std::string line; std::ifstream covar(covfile.c_str(), std::ifstream::in); int row = 0; while (std::getline(covar >> std::ws, line, '\n')) { int column = 0; std::vector entries = GeneralUtils::ParseToDbl(line, " "); if (entries.size() <= 1) { ERR(WRN) << "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 " - "entries on this line: " << row << std::endl; + "entries on this line: " + << row << std::endl; } for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { - // Check Rows - //assert(row > mat->GetNrows() && " covar rows doesn't match matrix rows."); - //assert(column > mat->GetNcols() && " covar cols doesn't match matrix cols."); + // assert(row > mat->GetNrows() && " covar rows doesn't match matrix + // rows."); + // assert(column > mat->GetNcols() && " covar cols doesn't match matrix + // cols."); // Fill Matrix (*mat)(row, column) = (*iter); column++; } row++; } return mat; } //******************************************************************* -TMatrixD* StatUtils::GetMatrixFromRootFile(std::string covfile, std::string histname) { -//******************************************************************* +TMatrixD* StatUtils::GetMatrixFromRootFile(std::string covfile, + std::string histname) { + //******************************************************************* std::string inputfile = covfile + ";" + histname; std::vector splitfile = GeneralUtils::ParseToStr(inputfile, ";"); if (splitfile.size() < 2) { ERR(FTL) << "No object name given!" << std::endl; throw; } // Get file TFile* tempfile = new TFile(splitfile[0].c_str(), "READ"); // Get Object TObject* obj = tempfile->Get(splitfile[1].c_str()); if (!obj) { ERR(FTL) << "Object " << splitfile[1] << " doesn't exist!" << std::endl; throw; } // Try casting TMatrixD* mat = dynamic_cast(obj); if (mat) { - TMatrixD* newmat = (TMatrixD*)mat->Clone(); delete mat; tempfile->Close(); return newmat; } TMatrixDSym* matsym = dynamic_cast(obj); if (matsym) { - TMatrixD* newmat = new TMatrixD(matsym->GetNrows(), matsym->GetNrows()); for (int i = 0; i < matsym->GetNrows(); i++) { for (int j = 0; j < matsym->GetNrows(); j++) { (*newmat)(i, j) = (*matsym)(i, j); } } delete matsym; tempfile->Close(); return newmat; } TH2D* mathist = dynamic_cast(obj); if (mathist) { TMatrixD* newmat = new TMatrixD(mathist->GetNbinsX(), mathist->GetNbinsX()); for (int i = 0; i < mathist->GetNbinsX(); i++) { for (int j = 0; j < mathist->GetNbinsX(); j++) { (*newmat)(i, j) = mathist->GetBinContent(i + 1, j + 1); } } delete mathist; tempfile->Close(); return newmat; } return NULL; } //******************************************************************* -TMatrixDSym* StatUtils::GetCovarFromTextFile(std::string covfile, int dim){ -//******************************************************************* +TMatrixDSym* StatUtils::GetCovarFromTextFile(std::string covfile, int dim) { + //******************************************************************* // Delete TempMat TMatrixD* tempmat = GetMatrixFromTextFile(covfile, dim, dim); // Make a symmetric covariance TMatrixDSym* newmat = new TMatrixDSym(tempmat->GetNrows()); - for (int i = 0; i < tempmat->GetNrows(); i++){ - for (int j = 0; j < tempmat->GetNrows(); j++){ - (*newmat)(i,j) = (*tempmat)(i,j); + for (int i = 0; i < tempmat->GetNrows(); i++) { + for (int j = 0; j < tempmat->GetNrows(); j++) { + (*newmat)(i, j) = (*tempmat)(i, j); } } - delete tempmat; return newmat; } //******************************************************************* -TMatrixDSym* StatUtils::GetCovarFromRootFile(std::string covfile, std::string histname){ -//******************************************************************* +TMatrixDSym* StatUtils::GetCovarFromRootFile(std::string covfile, + std::string histname) { + //******************************************************************* TMatrixD* tempmat = GetMatrixFromRootFile(covfile, histname); TMatrixDSym* newmat = new TMatrixDSym(tempmat->GetNrows()); - for (int i = 0; i < tempmat->GetNrows(); i++){ - for (int j = 0; j < tempmat->GetNrows(); j++){ - (*newmat)(i,j) = (*tempmat)(i,j); + for (int i = 0; i < tempmat->GetNrows(); i++) { + for (int j = 0; j < tempmat->GetNrows(); j++) { + (*newmat)(i, j) = (*tempmat)(i, j); } } delete tempmat; return newmat; } diff --git a/src/Tests/CMakeLists.txt b/src/Tests/CMakeLists.txt index 6fe0313..913f041 100644 --- a/src/Tests/CMakeLists.txt +++ b/src/Tests/CMakeLists.txt @@ -1,40 +1,51 @@ # 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_directories(${MINIMUM_INCLUDE_DIRECTORIES}) include_directories(${CMAKE_SOURCE_DIR}/src/Routines) include_directories(${CMAKE_SOURCE_DIR}/src/FCN) include_directories(${CMAKE_SOURCE_DIR}/src/MCStudies) include_directories(${CMAKE_SOURCE_DIR}/src/Smearceptance) include_directories(${EXP_INCLUDE_DIRECTORIES}) -SET(TESTAPPS SignalDefTests ParserTests SmearceptanceTests) +SET(TESTAPPS SignalDefTests ParserTests SmearceptanceTests FitMechanicsTest) foreach(appimpl ${TESTAPPS}) add_executable(${appimpl} ${appimpl}.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};${appimpl}) target_link_libraries(${appimpl} ${MODULETargets}) target_link_libraries(${appimpl} ${CMAKE_DEPENDLIB_FLAGS}) target_link_libraries(${appimpl} ${ROOT_LIBS}) if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") set_target_properties(${appimpl} PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) endif() install(TARGETS ${appimpl} DESTINATION tests) add_test(${appimpl} ${appimpl} 1) endforeach() + +add_library(DummySample SHARED DummySample.cxx) +target_link_libraries(DummySample ${MODULETargets}) +target_link_libraries(DummySample ${CMAKE_DEPENDLIB_FLAGS}) +target_link_libraries(DummySample ${ROOT_LIBS}) + +if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") + set_target_properties(DummySample PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) +endif() + +install(TARGETS DummySample DESTINATION tests) diff --git a/src/Tests/ConstructibleFitEvent.h b/src/Tests/ConstructibleFitEvent.h new file mode 100644 index 0000000..a8b8759 --- /dev/null +++ b/src/Tests/ConstructibleFitEvent.h @@ -0,0 +1,59 @@ +#include "TLorentzVector.h" +#include "TRandom3.h" +#include "FitEvent.h" + +struct ConstructibleFitEvent : public FitEvent { + ConstructibleFitEvent() : FitEvent() { fNParticles = 0; } + void AddPart(double Mom[4], size_t State, int PDG) { + fParticleMom[fNParticles][0] = Mom[0]; + fParticleMom[fNParticles][1] = Mom[1]; + fParticleMom[fNParticles][2] = Mom[2]; + fParticleMom[fNParticles][3] = Mom[3]; + fParticleState[fNParticles] = State; + fParticlePDG[fNParticles] = PDG; + fNParticles++; + } + void SetMode(int mode) { Mode = mode; } + std::string ToString() { + std::stringstream ss(""); + ss << "Mode: " << Mode << std::endl; + ss << "Particles: " << fNParticles << std::endl; + ss << " -> Particle Stack " << std::endl; + for (int i = 0; i < fNParticles; i++) { + ss << " -> -> " << i << ". " << fParticlePDG[i] << " " + << fParticleState[i] << " " + << " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1] << ", " + << fParticleMom[i][2] << ", " << fParticleMom[i][3] << ")." + << std::endl; + } + return ss.str(); + } +}; + +template +ConstructibleFitEvent MakePDGStackEvent(int (&ISpdgs)[N], int (&FSpdgs)[M], + int Mode = 1) { + ConstructibleFitEvent fe; + double MomHolder[4] = {0, 0, 1E3, 1E3}; + for (size_t p_it = 0; p_it < N; ++p_it) { + fe.AddPart(MomHolder, kInitialState, ISpdgs[p_it]); + } + + TRandom3 rnd; + TLorentzVector rnd4M; + TVector3 rn3M; + for (size_t p_it = 0; p_it < M; ++p_it) { + /// Could do better and actually get the correct masses... + rn3M.SetMagThetaPhi(fabs(rnd.Gaus(200, 50)), rnd.Uniform(M_PI), + 2 * rnd.Uniform(M_PI)); + rnd4M.SetVectM(rn3M, 105); + MomHolder[0] = rnd4M[0]; + MomHolder[1] = rnd4M[1]; + MomHolder[2] = rnd4M[2]; + MomHolder[3] = rnd4M[3]; + fe.AddPart(MomHolder, kFinalState, FSpdgs[p_it]); + } + fe.SetMode(Mode); + fe.OrderStack(); + return fe; +} diff --git a/src/Tests/ConstructibleInputHandler.h b/src/Tests/ConstructibleInputHandler.h new file mode 100644 index 0000000..3c8ee1f --- /dev/null +++ b/src/Tests/ConstructibleInputHandler.h @@ -0,0 +1,25 @@ +#include "InputHandler.h" + +struct ConstructibleInputHandler : public InputHandlerBase { + std::vector FitEvents; + + ConstructibleInputHandler( + std::string const& name) { + fName = name; + } + + void AddFitEvent(FitEvent *fe) { + FitEvents.push_back(fe); + + fNEvents = FitEvents.size(); + + std::cout << "[INFO]: Added event " << std::endl; + fe->Print(); + } + + FitEvent* GetNuisanceEvent( + const UInt_t entry, const bool lightweight) { + if (entry >= (UInt_t)fNEvents) return NULL; + return FitEvents[entry]; + } +}; diff --git a/src/Tests/DummyMinimizer.h b/src/Tests/DummyMinimizer.h new file mode 100644 index 0000000..5ef73b3 --- /dev/null +++ b/src/Tests/DummyMinimizer.h @@ -0,0 +1,61 @@ +#include "MinimizerRoutines.h" + +struct DummyMinimizer : public MinimizerRoutines { + DummyMinimizer() : MinimizerRoutines() { + nuiskey p_m1 = Config::CreateKey("parameter"); + + p_m1.SetS("type", "modenorm_parameter"); + p_m1.SetS("name", "mode_1"); + p_m1.SetD("nominal", 1); + p_m1.SetS("state", "FREE"); + p_m1.SetD("low", 0); + p_m1.SetD("high", 3); + p_m1.SetD("step", 0.1); + + nuiskey p_m2 = Config::CreateKey("parameter"); + + p_m2.SetS("type", "modenorm_parameter"); + p_m2.SetS("name", "mode_2"); + p_m2.SetD("nominal", 1); + p_m2.SetS("state", "FREE"); + p_m2.SetD("low", 0); + p_m2.SetD("high", 3); + p_m2.SetD("step", 0.1); + + nuiskey p_m3 = Config::CreateKey("parameter"); + + p_m3.SetS("type", "modenorm_parameter"); + p_m3.SetS("name", "mode_3"); + p_m3.SetD("nominal", 1); + p_m3.SetS("state", "FREE"); + p_m3.SetD("low", 0); + p_m3.SetD("high", 3); + p_m3.SetD("step", 0.1); + + nuiskey dummysample = Config::CreateKey("sample"); + dummysample.SetS("name", "DummySample"); + + Config::SetPar("dynamic_sample.path", + std::string(getenv("NUISANCE")) + "/build/Linux/tests"); + + Config::SetPar("EventManager", false); + // SETVERBOSITY(DEB); + + + + std::cout << "[INFO]: Parameter config:" << std::endl; + nuisconfig::GetConfig().PrintXML(NULL); + std::cout << "================" << std::endl; + + fOutputRootFile = new TFile("TestOutput.root", "RECREATE"); + + SetupMinimizerFromXML(); + SetupRWEngine(); + SetupFCN(); + } + + std::string GetParamName (int i) { return fParams[i]; } + double GetParamVal (int i) { return fCurVals[GetParamName(i)]; } + double GetParamError (int i) { return fErrorVals[GetParamName(i)]; } + +}; diff --git a/src/Tests/DummySample.cxx b/src/Tests/DummySample.cxx new file mode 100644 index 0000000..0734dc9 --- /dev/null +++ b/src/Tests/DummySample.cxx @@ -0,0 +1,89 @@ +#include "ConstructibleFitEvent.h" +#include "ConstructibleInputHandler.h" +#include "Measurement1D.h" + +struct DummySample : public Measurement1D { + + std::vector FitEvents; + DummySample(nuiskey samplekey) { + ConstructibleInputHandler *cih = new ConstructibleInputHandler("DummyIHandler"); + fInput = cih; + + fSettings = SampleSettings(samplekey); + fSettings.SetTitle("DummySample"); + FinaliseSampleSettings(); + + fDataHist = new TH1D("data", "", 3, 1, 4); + fMCHist = new TH1D("data", "", 3, 1, 4); + + double FakeDataError = Config::GetParD("FakeDataError"); + + fScaleFactor = 1; + + int IS[] = {14}; + int FS_CC0pi_1[] = {13, 2212}; + int FS_CC1pi_1[] = {13, 2212, 211}; + int FS_CC2pi_1[] = {13, 2212, 211, -211}; + FitEvents.push_back(new ConstructibleFitEvent(MakePDGStackEvent(IS, FS_CC0pi_1))); + FitEvents.back()->SetMode(1); + FitEvents.push_back(new ConstructibleFitEvent(MakePDGStackEvent(IS, FS_CC1pi_1))); + FitEvents.back()->SetMode(2); + FitEvents.push_back(new ConstructibleFitEvent(MakePDGStackEvent(IS, FS_CC2pi_1))); + FitEvents.back()->SetMode(3); + + for(size_t fe_it = 0; fe_it < FitEvents.size(); ++fe_it){ + cih->AddFitEvent(FitEvents[fe_it]); + fDataHist->SetBinContent(FitEvents[fe_it]->Mode, 1.0 + float(fe_it)*0.5); + fDataHist->SetBinError(FitEvents[fe_it]->Mode, FakeDataError); + } + + SetupDefaultHist(); + SetCovarFromDiagonal(); + FinaliseMeasurement(); + } + void FillEventVariables(FitEvent* nvect) { fXVar = nvect->Mode; } + + bool isSignal(FitEvent* nvect) { return true; } + + virtual ~DummySample(){ + for(size_t fe_it = 0; fe_it < FitEvents.size(); ++fe_it){ + delete FitEvents[fe_it]; + } + } +}; + +static char const* SampleNames[] = {"DummySample"}; +static int const NSamples = 1; + +extern "C" { +int DSF_NSamples() { return NSamples; } +char const* DSF_GetSampleName(int i) { + if (i < NSamples) { + return SampleNames[i]; + } + return 0; +} +MeasurementBase* DSF_GetSample(int i, void* samplekey) { + nuiskey* sk = reinterpret_cast(samplekey); + if (!sk) { + return 0; + } + + if (sk->GetS("name") != DSF_GetSampleName(i)) { + std::cout + << "[ERROR]: When instantiating dynamic sample. Samplekey named: " + << sk->GetS("name") + << ", but requested sample named: " << DSF_GetSampleName(i) + << ". It is possible that the nuiskey object is lost in translation. " + "Was NUISANCE and this dynamic sample manifest built with the same " + "environment and compiler?" + << std::endl; + } + + if (i == 0) { + return new DummySample(*sk); + } + return 0; +} +void DSF_DestroySample(MeasurementBase* mb) { delete mb; } +} diff --git a/src/Tests/FitMechanicsTest.cxx b/src/Tests/FitMechanicsTest.cxx new file mode 100644 index 0000000..842c995 --- /dev/null +++ b/src/Tests/FitMechanicsTest.cxx @@ -0,0 +1,36 @@ +#include "DummyMinimizer.h" + +int main(){ + + double FakeDataError = 1E-4; + Config::SetPar("FakeDataError",FakeDataError); + + DummyMinimizer min; + + std::cout << "======Running simple minimizer test======" << std::endl; + + min.Run(); + + std::cout << "======Done minimizing, checking output======" << std::endl; + + std::vector pnames; + pnames.push_back("mode_1"); + pnames.push_back("mode_2"); + pnames.push_back("mode_3"); + + std::vector postfitvals; + postfitvals.push_back(1); + postfitvals.push_back(1.5); + postfitvals.push_back(2); + + for(size_t i = 0; i < 3; ++i){ + std::cout << "\tTest param name: " << pnames[i] << " == " << min.GetParamName(i) << std::endl; + assert(pnames[i]==min.GetParamName(i)); + std::cout << "\tTest param post-fit value : fabs(" << postfitvals[i] << " - " << min.GetParamVal(i) << ") < " << min.GetParamError(i) << std::endl; + assert(fabs(postfitvals[i] - min.GetParamVal(i)) < min.GetParamError(i)); + std::cout << "\tTest post-fit error " << min.GetParamError(i) << " < " << FakeDataError << std::endl; + assert(min.GetParamError(i) < FakeDataError); + } + + return 0; +} diff --git a/src/Tests/SignalDefTests.cxx b/src/Tests/SignalDefTests.cxx index 551f66b..031526f 100644 --- a/src/Tests/SignalDefTests.cxx +++ b/src/Tests/SignalDefTests.cxx @@ -1,447 +1,389 @@ #include #include -#include "TLorentzVector.h" -#include "TRandom3.h" - -#include "FitEvent.h" #include "SignalDef.h" +#include "ConstructibleFitEvent.h" -struct ConstructibleFitEvent : public FitEvent { - ConstructibleFitEvent() : FitEvent() { fNParticles = 0; } - void AddPart(double Mom[4], size_t State, int PDG) { - fParticleMom[fNParticles][0] = Mom[0]; - fParticleMom[fNParticles][1] = Mom[1]; - fParticleMom[fNParticles][2] = Mom[2]; - fParticleMom[fNParticles][3] = Mom[3]; - fParticleState[fNParticles] = State; - fParticlePDG[fNParticles] = PDG; - fNParticles++; - } - void SetMode(int mode) { Mode = mode; } - std::string ToString() { - std::stringstream ss(""); - ss << "Mode: " << Mode << std::endl; - ss << "Particles: " << fNParticles << std::endl; - ss << " -> Particle Stack " << std::endl; - for (int i = 0; i < fNParticles; i++) { - ss << " -> -> " << i << ". " << fParticlePDG[i] << " " - << fParticleState[i] << " " - << " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1] << ", " - << fParticleMom[i][2] << ", " << fParticleMom[i][3] << ")." - << std::endl; - } - return ss.str(); - } -}; - -template -ConstructibleFitEvent MakePDGStackEvent(int (&ISpdgs)[N], int (&FSpdgs)[M], - int Mode = 1) { - ConstructibleFitEvent fe; - double MomHolder[4] = {0, 0, 1E3, 1E3}; - for (size_t p_it = 0; p_it < N; ++p_it) { - fe.AddPart(MomHolder, kInitialState, ISpdgs[p_it]); - } - - TRandom3 rnd; - TLorentzVector rnd4M; - TVector3 rn3M; - for (size_t p_it = 0; p_it < M; ++p_it) { - /// Could do better and actually get the correct masses... - rn3M.SetMagThetaPhi(fabs(rnd.Gaus(200, 50)), rnd.Uniform(M_PI), - 2 * rnd.Uniform(M_PI)); - rnd4M.SetVectM(rn3M, 105); - MomHolder[0] = rnd4M[0]; - MomHolder[1] = rnd4M[1]; - MomHolder[2] = rnd4M[2]; - MomHolder[3] = rnd4M[3]; - fe.AddPart(MomHolder, kFinalState, FSpdgs[p_it]); - } - fe.SetMode(Mode); - fe.OrderStack(); - return fe; -} int main(int argc, char const *argv[]) { bool FailOnFail = (argc > 1); LOG_VERB(SAM); LOG(FIT) << "* Running SignalDef Tests" << std::endl; LOG(FIT) << "***************************************************" << std::endl; int IS[] = {14}; int FS_CC0pi_1[] = {13, 2112, 2212, 2112, 2212}; ConstructibleFitEvent fe_CC0pi_1 = MakePDGStackEvent(IS, FS_CC0pi_1); int FS_CC0pi_2[] = {13, 2112, 2212, 2112, 2212}; ConstructibleFitEvent fe_CC0pi_2 = MakePDGStackEvent(IS, FS_CC0pi_2, 2); int FS_CC0pi_3[] = {-13, 2112, 2212, 2112, 2212}; ConstructibleFitEvent fe_CC0pi_3 = MakePDGStackEvent(IS, FS_CC0pi_3); int FS_CC0pi_4[] = {13, 2112, 2212}; ConstructibleFitEvent fe_CC0pi_4 = MakePDGStackEvent(IS, FS_CC0pi_4, 12); int FS_CC1pip_1[] = {13, 2212, 2112, 211}; ConstructibleFitEvent fe_CC1pip_1 = MakePDGStackEvent(IS, FS_CC1pip_1, 2); int FS_CC1pim_1[] = {13, -211, 2212, 2112}; ConstructibleFitEvent fe_CC1pim_1 = MakePDGStackEvent(IS, FS_CC1pim_1, 11); int FS_CC1pi0_1[] = {13, 2212, 111, 2112}; ConstructibleFitEvent fe_CC1pi0_1 = MakePDGStackEvent(IS, FS_CC1pi0_1, 12); int FS_CC1pi0_2[] = {11, 2212, 111, 2112}; ConstructibleFitEvent fe_CC1pi0_2 = MakePDGStackEvent(IS, FS_CC1pi0_2, 12); int FS_CCNpi_1[] = {13, 2212, 111, 2112, 211}; ConstructibleFitEvent fe_CCNpi_1 = MakePDGStackEvent(IS, FS_CCNpi_1, 21); int FS_CCNpi_2[] = {13, -211, 211, 2112, 211}; ConstructibleFitEvent fe_CCNpi_2 = MakePDGStackEvent(IS, FS_CCNpi_2, 21); int FS_CCNpi_3[] = {13, 2212, 111, 211, -211}; ConstructibleFitEvent fe_CCNpi_3 = MakePDGStackEvent(IS, FS_CCNpi_3, 26); int FS_CCNpi_4[] = {13, 2212, 111, 111, 111}; ConstructibleFitEvent fe_CCNpi_4 = MakePDGStackEvent(IS, FS_CCNpi_4, 26); int FS_CCCOH_1[] = {13, 211}; ConstructibleFitEvent fe_CCCOH_1 = MakePDGStackEvent(IS, FS_CCCOH_1, 16); int FS_NCel_1[] = {14, 2112}; ConstructibleFitEvent fe_NCel_1 = MakePDGStackEvent(IS, FS_NCel_1, 52); int FS_NCel_2[] = {12, 2212}; ConstructibleFitEvent fe_NCel_2 = MakePDGStackEvent(IS, FS_NCel_2, 51); int FS_NC1pi_1[] = {14, 2112, -211}; ConstructibleFitEvent fe_NC1pi_1 = MakePDGStackEvent(IS, FS_NC1pi_1, 31); int FS_NCNpi_1[] = {14, 2212, 211}; ConstructibleFitEvent fe_NCNpi_1 = MakePDGStackEvent(IS, FS_NCNpi_1, 32); LOG(FIT) << "* Testing: SignalDef::isCCINC" << std::endl; std::map isCCINC_PassExpectations; isCCINC_PassExpectations[&fe_CC0pi_1] = true; // numu CC0pi isCCINC_PassExpectations[&fe_CC0pi_2] = true; // numu CC0pi (2p2h) isCCINC_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi isCCINC_PassExpectations[&fe_CC0pi_4] = true; // numu CC0pi (RES) isCCINC_PassExpectations[&fe_CC1pip_1] = true; // numu CC1pip (2p2h) isCCINC_PassExpectations[&fe_CC1pim_1] = true; // numu CC1pim isCCINC_PassExpectations[&fe_CC1pi0_1] = true; // numu CC1pi0 isCCINC_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0 isCCINC_PassExpectations[&fe_CCNpi_1] = true; // numu CC multi pi isCCINC_PassExpectations[&fe_CCNpi_2] = true; // numu CC multi pi isCCINC_PassExpectations[&fe_CCNpi_3] = true; // numu CC multi pi isCCINC_PassExpectations[&fe_CCNpi_4] = true; // numu CC multi pi isCCINC_PassExpectations[&fe_CCCOH_1] = true; // numu CC COH pi isCCINC_PassExpectations[&fe_NCel_1] = false; // numu NCEl isCCINC_PassExpectations[&fe_NCel_2] = false; // nue NCEl isCCINC_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi isCCINC_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi size_t ctr = 0; for (std::map::iterator fe_it = isCCINC_PassExpectations.begin(); fe_it != isCCINC_PassExpectations.end(); ++fe_it, ++ctr) { bool res = SignalDef::isCCINC(fe_it->first, 14); if (res != fe_it->second) { ERR(FTL) << "Event: (" << ctr << ")\n" << fe_it->first->ToString() << std::endl; ERR(FTL) << (res ? "passed" : "failed") << " SignalDef::isCCINC unexpectedly." << std::endl; } else { LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed") << " as expected." << std::endl; } if (FailOnFail) { assert(res == fe_it->second); } } LOG(FIT) << "* Testing: SignalDef::isNCINC" << std::endl; std::map isNCINC_PassExpectations; isNCINC_PassExpectations[&fe_CC0pi_1] = false; // numu CC0pi isNCINC_PassExpectations[&fe_CC0pi_2] = false; // numu CC0pi (2p2h) isNCINC_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi isNCINC_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES) isNCINC_PassExpectations[&fe_CC1pip_1] = false; // numu CC1pip (2p2h) isNCINC_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim isNCINC_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0 isNCINC_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0 isNCINC_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi isNCINC_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi isNCINC_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi isNCINC_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi isNCINC_PassExpectations[&fe_CCCOH_1] = false; // numu CC COH pi isNCINC_PassExpectations[&fe_NCel_1] = true; // numu NCEl isNCINC_PassExpectations[&fe_NCel_2] = false; // nue NCEl isNCINC_PassExpectations[&fe_NC1pi_1] = true; // numu NC1pi isNCINC_PassExpectations[&fe_NCNpi_1] = true; // numu NC multi pi ctr = 0; for (std::map::iterator fe_it = isNCINC_PassExpectations.begin(); fe_it != isNCINC_PassExpectations.end(); ++fe_it, ++ctr) { bool res = SignalDef::isNCINC(fe_it->first, 14); if (res != fe_it->second) { ERR(FTL) << "Event: (" << ctr << ")\n" << fe_it->first->ToString() << std::endl; ERR(FTL) << (res ? "passed" : "failed") << " SignalDef::isNCINC unexpectedly." << std::endl; } else { LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed") << " as expected." << std::endl; } if (FailOnFail) { assert(res == fe_it->second); } } LOG(FIT) << "* Testing: SignalDef::isCC0pi" << std::endl; std::map isCC0pi_PassExpectations; isCC0pi_PassExpectations[&fe_CC0pi_1] = true; // numu CC0pi isCC0pi_PassExpectations[&fe_CC0pi_2] = true; // numu CC0pi (2p2h) isCC0pi_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi isCC0pi_PassExpectations[&fe_CC0pi_4] = true; // numu CC0pi (RES) isCC0pi_PassExpectations[&fe_CC1pip_1] = false; // numu CC1pip (2p2h) isCC0pi_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim isCC0pi_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0 isCC0pi_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0 isCC0pi_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi isCC0pi_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi isCC0pi_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi isCC0pi_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi isCC0pi_PassExpectations[&fe_CCCOH_1] = false; // numu CC COH pi isCC0pi_PassExpectations[&fe_NCel_1] = false; // numu NCEl isCC0pi_PassExpectations[&fe_NCel_2] = false; // nue NCEl isCC0pi_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi isCC0pi_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi ctr = 0; for (std::map::iterator fe_it = isCC0pi_PassExpectations.begin(); fe_it != isCC0pi_PassExpectations.end(); ++fe_it, ++ctr) { bool res = SignalDef::isCC0pi(fe_it->first, 14); if (res != fe_it->second) { ERR(FTL) << "Event: (" << ctr << ")\n" << fe_it->first->ToString() << " " << std::endl; ERR(FTL) << (res ? "passed" : "failed") << " SignalDef::isCC0pi unexpectedly." << std::endl; } else { LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed") << " as expected." << std::endl; } if (FailOnFail) { assert(res == fe_it->second); } } LOG(FIT) << "* Testing: SignalDef::isCCQELike" << std::endl; std::map isCCQELike_PassExpectations; isCCQELike_PassExpectations[&fe_CC0pi_1] = true; // numu CC0pi isCCQELike_PassExpectations[&fe_CC0pi_2] = true; // numu CC0pi (2p2h) isCCQELike_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi isCCQELike_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES) isCCQELike_PassExpectations[&fe_CC1pip_1] = true; // numu CC1pip (2p2h) isCCQELike_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim isCCQELike_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0 isCCQELike_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0 isCCQELike_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi isCCQELike_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi isCCQELike_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi isCCQELike_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi isCCQELike_PassExpectations[&fe_CCCOH_1] = false; // numu CC COH pi isCCQELike_PassExpectations[&fe_NCel_1] = false; // numu NCEl isCCQELike_PassExpectations[&fe_NCel_2] = false; // nue NCEl isCCQELike_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi isCCQELike_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi ctr = 0; for (std::map::iterator fe_it = isCCQELike_PassExpectations.begin(); fe_it != isCCQELike_PassExpectations.end(); ++fe_it, ++ctr) { bool res = SignalDef::isCCQELike(fe_it->first, 14); if (res != fe_it->second) { ERR(FTL) << "Event: (" << ctr << ")\n" << fe_it->first->ToString() << std::endl; ERR(FTL) << (res ? "passed" : "failed") << " SignalDef::isCCQELike unexpectedly." << std::endl; } else { LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed") << " as expected." << std::endl; } if (FailOnFail) { assert(res == fe_it->second); } } LOG(FIT) << "* Testing: SignalDef::isCCQE" << std::endl; std::map isCCQE_PassExpectations; isCCQE_PassExpectations[&fe_CC0pi_1] = true; // numu CC0pi isCCQE_PassExpectations[&fe_CC0pi_2] = false; // numu CC0pi (2p2h) isCCQE_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi isCCQE_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES) isCCQE_PassExpectations[&fe_CC1pip_1] = false; // numu CC1pip (2p2h) isCCQE_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim isCCQE_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0 isCCQE_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0 isCCQE_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi isCCQE_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi isCCQE_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi isCCQE_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi isCCQE_PassExpectations[&fe_CCCOH_1] = false; // numu CC COH pi isCCQE_PassExpectations[&fe_NCel_1] = false; // numu NCEl isCCQE_PassExpectations[&fe_NCel_2] = false; // nue NCEl isCCQE_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi isCCQE_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi ctr = 0; for (std::map::iterator fe_it = isCCQE_PassExpectations.begin(); fe_it != isCCQE_PassExpectations.end(); ++fe_it, ++ctr) { bool res = SignalDef::isCCQE(fe_it->first, 14); if (res != fe_it->second) { ERR(FTL) << "Event: (" << ctr << ")\n" << fe_it->first->ToString() << std::endl; ERR(FTL) << (res ? "passed" : "failed") << " SignalDef::isCCQE unexpectedly." << std::endl; } else { LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed") << " as expected." << std::endl; } if (FailOnFail) { assert(res == fe_it->second); } } LOG(FIT) << "* Testing: SignalDef::isCCCOH" << std::endl; std::map isCCCOH_PassExpectations; isCCCOH_PassExpectations[&fe_CC0pi_1] = false; // numu CC0pi isCCCOH_PassExpectations[&fe_CC0pi_2] = false; // numu CC0pi (2p2h) isCCCOH_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi isCCCOH_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES) isCCCOH_PassExpectations[&fe_CC1pip_1] = false; // numu CC1pip (2p2h) isCCCOH_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim isCCCOH_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0 isCCCOH_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0 isCCCOH_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi isCCCOH_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi isCCCOH_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi isCCCOH_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi isCCCOH_PassExpectations[&fe_CCCOH_1] = true; // numu CC COH pi isCCCOH_PassExpectations[&fe_NCel_1] = false; // numu NCEl isCCCOH_PassExpectations[&fe_NCel_2] = false; // nue NCEl isCCCOH_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi isCCCOH_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi ctr = 0; for (std::map::iterator fe_it = isCCCOH_PassExpectations.begin(); fe_it != isCCCOH_PassExpectations.end(); ++fe_it, ++ctr) { bool res = SignalDef::isCCCOH(fe_it->first, 14, 211); if (res != fe_it->second) { ERR(FTL) << "Event: (" << ctr << ")\n" << fe_it->first->ToString() << std::endl; ERR(FTL) << (res ? "passed" : "failed") << " SignalDef::isCCCOH unexpectedly." << std::endl; } else { LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed") << " as expected." << std::endl; } if (FailOnFail) { assert(res == fe_it->second); } } LOG(FIT) << "* Testing: SignalDef::isCC1pi" << std::endl; std::map isCC1pi_PassExpectations; isCC1pi_PassExpectations[&fe_CC0pi_1] = false; // numu CC0pi isCC1pi_PassExpectations[&fe_CC0pi_2] = false; // numu CC0pi (2p2h) isCC1pi_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi isCC1pi_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES) isCC1pi_PassExpectations[&fe_CC1pip_1] = true; // numu CC1pip (2p2h) isCC1pi_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim isCC1pi_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0 isCC1pi_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0 isCC1pi_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi isCC1pi_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi isCC1pi_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi isCC1pi_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi isCC1pi_PassExpectations[&fe_CCCOH_1] = true; // numu CC COH pi isCC1pi_PassExpectations[&fe_NCel_1] = false; // numu NCEl isCC1pi_PassExpectations[&fe_NCel_2] = false; // nue NCEl isCC1pi_PassExpectations[&fe_NC1pi_1] = false; // numu NC1pi isCC1pi_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi ctr = 0; for (std::map::iterator fe_it = isCC1pi_PassExpectations.begin(); fe_it != isCC1pi_PassExpectations.end(); ++fe_it, ++ctr) { bool res = SignalDef::isCC1pi(fe_it->first, 14, 211); if (res != fe_it->second) { ERR(FTL) << "Event: (" << ctr << ")\n" << fe_it->first->ToString() << std::endl; ERR(FTL) << (res ? "passed" : "failed") << " SignalDef::isCC1pi unexpectedly." << std::endl; } else { LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed") << " as expected." << std::endl; } if (FailOnFail) { assert(res == fe_it->second); } } LOG(FIT) << "* Testing: SignalDef::isNC1pi" << std::endl; std::map isNC1pi_PassExpectations; isNC1pi_PassExpectations[&fe_CC0pi_1] = false; // numu CC0pi isNC1pi_PassExpectations[&fe_CC0pi_2] = false; // numu CC0pi (2p2h) isNC1pi_PassExpectations[&fe_CC0pi_3] = false; // numub CC0pi isCCINC_PassExpectations[&fe_CC0pi_4] = false; // numu CC0pi (RES) isNC1pi_PassExpectations[&fe_CC1pip_1] = false; // numu CC1pip (2p2h) isNC1pi_PassExpectations[&fe_CC1pim_1] = false; // numu CC1pim isNC1pi_PassExpectations[&fe_CC1pi0_1] = false; // numu CC1pi0 isNC1pi_PassExpectations[&fe_CC1pi0_2] = false; // nue CC1pi0 isNC1pi_PassExpectations[&fe_CCNpi_1] = false; // numu CC multi pi isNC1pi_PassExpectations[&fe_CCNpi_2] = false; // numu CC multi pi isNC1pi_PassExpectations[&fe_CCNpi_3] = false; // numu CC multi pi isNC1pi_PassExpectations[&fe_CCNpi_4] = false; // numu CC multi pi isCCINC_PassExpectations[&fe_CCCOH_1] = false; // numu CC COH pi isNC1pi_PassExpectations[&fe_NCel_1] = false; // numu NCEl isNC1pi_PassExpectations[&fe_NCel_2] = false; // nue NCEl isCCINC_PassExpectations[&fe_NC1pi_1] = true; // numu NC1pi isCCINC_PassExpectations[&fe_NCNpi_1] = false; // numu NC multi pi ctr = 0; for (std::map::iterator fe_it = isNC1pi_PassExpectations.begin(); fe_it != isNC1pi_PassExpectations.end(); ++fe_it, ++ctr) { bool res = SignalDef::isNC1pi(fe_it->first, 14, -211); if (res != fe_it->second) { ERR(FTL) << "Event: (" << ctr << ")\n" << fe_it->first->ToString() << std::endl; ERR(FTL) << (res ? "passed" : "failed") << " SignalDef::isNC1pi unexpectedly." << std::endl; } else { LOG(SAM) << "Event: (" << ctr << ") " << (res ? "passed" : "failed") << " as expected." << std::endl; } if (FailOnFail) { assert(res == fe_it->second); } } // SignalDef::isCCWithFS(&fe,14); }