diff --git a/app/neut_NUISANCE.cxx b/app/neut_NUISANCE.cxx index b319240..a2e1dbc 100644 --- a/app/neut_NUISANCE.cxx +++ b/app/neut_NUISANCE.cxx @@ -1,839 +1,839 @@ #ifdef __NEUT_ENABLED__ #include "ComparisonRoutines.h" #include "ParserUtils.h" #include "TargetUtils.h" #ifdef WINDOWS #include #define GetCurrentDir _getcwd #else #include #define GetCurrentDir getcwd #endif // All possible inputs std::string gOptEnergyDef; std::vector gOptEnergyRange; int gOptNumberEvents = -1; int gOptNumberTestEvents = 5E6; std::string gOptGeneratorList = "Default"; std::string gOptCrossSections = "Default"; // If default this will look in $NUISANCE/data/neut/Default_params.txt int gOptSeed = time(NULL); std::string gOptTargetDef = ""; std::string gOptFluxDef = ""; std::string gOptOutputFile = ""; int gOptRunNumber = -1; using namespace TargetUtils; void GetCommandLineArgs (int argc, char ** argv); void PrintSyntax (void); std::string GETCWD(){ char cCurrentPath[FILENAME_MAX]; if (!GetCurrentDir(cCurrentPath, sizeof(cCurrentPath))){ THROW("CANT FIND CURRENT DIRECTORY!"); } std::string curdir = std::string(cCurrentPath); return curdir; } std::string ExpandPath(std::string name){ // Get Current char cCurrentPath[FILENAME_MAX]; if (!GetCurrentDir(cCurrentPath, sizeof(cCurrentPath))){ THROW("CANT FIND CURRENT DIRECTORY!"); } std::string curdir = std::string(cCurrentPath); // If first entry is not / then add the current working directory if (!name.empty() and name.at(0) != '/'){ name = curdir + "/" + name; } return name; } std::string GetBaseName(std::string name){ std::vector splitfile = GeneralUtils::ParseToStr(name,"/"); std::string filename = ""; if (splitfile.size() == 1){ filename = splitfile[0]; } else if (splitfile.size() > 1){ filename = splitfile[splitfile.size()-1]; } else { THROW("Cannot split filename: " << name); } return filename; } std::string GetDynamicModes(std::string list, bool neutrino){ LOG(FIT) << "Using " << list << " to define interaction modes for Neutrino=" << neutrino << std::endl; std::map modes; std::vector ids; // Create vector of ids for the print out and future reference /* C C nu nub C 1: CC Q.E. CC Q.E.( Free ) C 2-4: CC 1pi CC 1pi C 5: CC DIS 1320 CC DIS 1.3 < W < 2.0 C 6-9: NC 1pi NC 1pi C 10: NC DIS 1320 NC DIS 1.3 < W < 2.0 C 11: NC els CC Q.E.( Bound ) C 12: NC els NC els C 13: NC els NC els C 14: coherent NC els C 15: coherent coherent C 16: CC eta coherent C 17 NC eta CC eta C 18: NC eta NC eta C 19: CC K NC eta C 20 NC K CC K C 21: NC K NC K C 22: N/A NC K C 23: CC DIS CC DIS (W > 2.0) C 24: NC DIS NC DIS (W > 2.0) C 25: CC 1 gamma CC 1 gamma C 26: NC 1 gamma NC 1 gamma C 27: NC 1 gamma NC 1 gamma C 28: 2p2h 2p2h */ ids.push_back("crsmode_CCQE" ); ids.push_back("crsmode_CC2P2H" ); ids.push_back("crsmode_CC1pi"); ids.push_back("crsmode_CCDIS_lowW" ); ids.push_back("crsmode_NC1pi"); ids.push_back("crsmode_NCDIS_lowW" ); ids.push_back("crsmode_NCEL"); ids.push_back("crsmode_CCCOH"); ids.push_back("crsmode_NCCOH"); ids.push_back("crsmode_CCETA"); ids.push_back("crsmode_NCETA"); ids.push_back("crsmode_CCKAON"); ids.push_back("crsmode_NCKAON"); ids.push_back("crsmode_CCDIS_highW"); ids.push_back("crsmode_NCDIS_highW"); ids.push_back("crsmode_CCGAMMA"); ids.push_back("crsmode_NCGAMMA"); // Now define possible models if (!list.compare("Default")){ // Everything but MEC modes["crsmode_CCQE"] = 1; modes["crsmode_CC2P2H"] = 0; modes["crsmode_CC1pi"] = 1; modes["crsmode_CCDIS_lowW"] = 1; modes["crsmode_CCCOH"] = 1; modes["crsmode_CCETA"] = 1; modes["crsmode_CCKAON"] = 1; modes["crsmode_CCDIS_highW"] = 1; modes["crsmode_CCGAMMA"] = 1; modes["crsmode_NC1pi"] = 1; modes["crsmode_NCDIS_lowW"] = 1; modes["crsmode_NCEL"] = 1; modes["crsmode_NCCOH"] = 1; modes["crsmode_NCETA"] = 1; modes["crsmode_NCKAON"] = 1; modes["crsmode_NCDIS_highW"] = 1; modes["crsmode_NCGAMMA"] = 1; } else if (!list.compare("DefaultFree")){ modes["crsmode_CCQE"] = 1; modes["crsmode_CC2P2H"] = 0; modes["crsmode_CC1pi"] = 1; modes["crsmode_CCDIS_lowW"] = 1; modes["crsmode_CCCOH"] = 0; modes["crsmode_CCETA"] = 1; modes["crsmode_CCKAON"] = 1; modes["crsmode_CCDIS_highW"] = 1; modes["crsmode_CCGAMMA"] = 1; modes["crsmode_NC1pi"] = 1; modes["crsmode_NCDIS_lowW"] = 1; modes["crsmode_NCEL"] = 1; modes["crsmode_NCCOH"] = 0; modes["crsmode_NCETA"] = 1; modes["crsmode_NCKAON"] = 1; modes["crsmode_NCDIS_highW"] = 1; modes["crsmode_NCGAMMA"] = 1; } else if (!list.compare("Default+MEC")){ modes["crsmode_CCQE"] = 1; modes["crsmode_CC2P2H"] = 1; modes["crsmode_CC1pi"] = 1; modes["crsmode_CCDIS_lowW"] = 1; modes["crsmode_CCCOH"] = 1; modes["crsmode_CCETA"] = 1; modes["crsmode_CCKAON"] = 1; modes["crsmode_CCDIS_highW"] = 1; modes["crsmode_CCGAMMA"] = 1; modes["crsmode_NC1pi"] = 1; modes["crsmode_NCDIS_lowW"] = 1; modes["crsmode_NCEL"] = 1; modes["crsmode_NCCOH"] = 1; modes["crsmode_NCETA"] = 1; modes["crsmode_NCKAON"] = 1; modes["crsmode_NCDIS_highW"] = 1; modes["crsmode_NCGAMMA"] = 1; } else { THROW("Event generator list " << list << " not found!"); } // Now we actually have to make the conversion because NEUTS modes organisation are a mess. /* C C nu nub C 1: CC Q.E. CC Q.E.( Free ) C 2-4: CC 1pi CC 1pi C 5: CC DIS 1320 CC DIS 1.3 < W < 2.0 C 6-9: NC 1pi NC 1pi C 10: NC DIS 1320 NC DIS 1.3 < W < 2.0 C 11: NC els CC Q.E.( Bound ) C 12: NC els NC els C 13: NC els NC els C 14: coherent NC els C 15: coherent coherent C 16: CC eta coherent C 17 NC eta CC eta C 18: NC eta NC eta C 19: CC K NC eta C 20 NC K CC K C 21: NC K NC K C 22: N/A NC K C 23: CC DIS CC DIS (W > 2.0) C 24: NC DIS NC DIS (W > 2.0) C 25: CC 1 gamma CC 1 gamma C 26,27: NC 1 gamma NC 1 gamma */ std::string modestring_neutrino = "NEUT-CRS "; std::string modestring_antineutrino = "NEUT-CRSB "; // Neutrino First if (neutrino){ // Fill empty NEUT-CRSB for (size_t i = 0; i < 27; i++){ modestring_antineutrino += " 0"; } modestring_neutrino += (modes["crsmode_CCQE"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_CC1pi"]? " 1 1 1" : " 0 0 0"); modestring_neutrino += (modes["crsmode_CCDIS_lowW"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_NC1pi"]? " 1 1 1 1" : " 0 0 0 0"); modestring_neutrino += (modes["crsmode_NCDIS_lowW"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_NCEL"]? " 1 1 1" : " 0 0 0"); modestring_neutrino += (modes["crsmode_CCCOH"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_NCCOH"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_CCETA"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_NCETA"]? " 1 1" : " 0 0"); modestring_neutrino += (modes["crsmode_CCKAON"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_NCKAON"]? " 1 1" : " 0 0"); modestring_neutrino += " 1"; // /NA modestring_neutrino += (modes["crsmode_CCDIS_highW"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_NCDIS_highW"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_CCGAMMA"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_NCGAMMA"]? " 1" : " 0"); modestring_neutrino += (modes["crsmode_CC2P2H"]? " 1" : " 0"); } else { // Fill Empty NEUT CRS for (size_t i = 0; i < 27; i++){ modestring_neutrino += " 0"; } modestring_antineutrino += (modes["crsmode_CCQE"]? " 1" : " 0"); modestring_antineutrino += (modes["crsmode_CC1pi"]? " 1 1 1" : " 0 0 0"); modestring_antineutrino += (modes["crsmode_CCDIS_lowW"]? " 1" : " 0"); modestring_antineutrino += (modes["crsmode_NC1pi"]? " 1 1 1 1" : " 0 0 0 0"); modestring_antineutrino += (modes["crsmode_NCDIS_lowW"]? " 1" : " 0"); modestring_antineutrino += (modes["crsmode_CCQE"]? " 1" : " 0"); modestring_antineutrino += (modes["crsmode_NCEL"]? " 1 1 1" : " 0 0 0"); modestring_antineutrino += (modes["crsmode_CCCOH"]? " 1" : " 0"); modestring_antineutrino += (modes["crsmode_NCCOH"]? " 1" : " 0"); modestring_antineutrino += (modes["crsmode_CCETA"]? " 1" : " 0"); modestring_antineutrino += (modes["crsmode_NCETA"]? " 1 1" : " 0 0"); modestring_antineutrino += (modes["crsmode_CCKAON"]? " 1" : " 0"); modestring_antineutrino += (modes["crsmode_NCKAON"]? " 1 1" : " 0 0"); modestring_antineutrino += (modes["crsmode_CCDIS_highW"]? " 1" : " 0"); modestring_antineutrino+= (modes["crsmode_NCDIS_highW"]? " 1" : " 0"); modestring_antineutrino+= (modes["crsmode_CCGAMMA"]? " 1" : " 0"); modestring_antineutrino+= (modes["crsmode_NCGAMMA"]? " 1" : " 0"); modestring_antineutrino+= (modes["crsmode_CC2P2H"]? " 1" : " 0"); } return "NEUT-MODE -1 \n" + modestring_neutrino + "\n" + modestring_antineutrino; } std::map MakeNewFluxFile(std::string flux, std::string out){ LOG(FIT) << "Using " << flux << " to define NEUT beam." << std::endl; std::vector fluxargs = GeneralUtils::ParseToStr(flux,","); if (fluxargs.size() < 2){ THROW("Expected flux in the format: file.root,hist_name1[pdg1],... : reveived : " << flux); } // Build Map std::map fluxmap; fluxmap["beam_type"] = "6"; fluxmap["beam_inputroot"] = fluxargs[0]; fluxmap["beam_inputroot_nue"] = ""; fluxmap["beam_inputroot_nueb"] = ""; fluxmap["beam_inputroot_numu"] = ""; fluxmap["beam_inputroot_numub"] = ""; fluxmap["beam_inputroot_nutau"] = ""; fluxmap["beam_inputroot_nutaub"] = ""; // Split by beam bdgs for (uint i = 1; i < fluxargs.size(); i++){ std::string histdef = fluxargs[i]; string::size_type open_bracket = histdef.find("["); string::size_type close_bracket = histdef.find("]"); string::size_type ibeg = 0; string::size_type iend = open_bracket; string::size_type jbeg = open_bracket+1; string::size_type jend = close_bracket-1; std::string name = std::string(histdef.substr(ibeg,iend).c_str()); int pdg = atoi(histdef.substr(jbeg,jend).c_str()); if (pdg == 12) fluxmap["beam_inputroot_nue"] = name; else if (pdg ==-12) fluxmap["beam_inputroot_nueb"] = name; else if (pdg == 14) fluxmap["beam_inputroot_numu"] = name; else if (pdg ==-14) fluxmap["beam_inputroot_numub"] = name; else if (pdg == 16) fluxmap["beam_inputroot_tau"] = name; else if (pdg ==-16) fluxmap["beam_inputroot_taub"] = name; } // Now create a new flux file matching the output file std::cout << " -> Moving flux from '" + fluxmap["beam_inputroot"] + "' to current directory to keep everything organised." << std::endl; TFile* fluxread = new TFile(fluxmap["beam_inputroot"].c_str(),"READ"); TFile* fluxwrite = new TFile((out + ".flux.root").c_str(),"RECREATE"); for(std::map::iterator iter = fluxmap.begin(); iter != fluxmap.end(); iter++){ TH1* temp = (TH1*)fluxread->Get(iter->second.c_str()); if (!temp) continue; TH1D* cuthist = (TH1D*)temp->Clone(); // Restrict energy range if required if (gOptEnergyRange.size() == 2){ for (int i = 0; i < cuthist->GetNbinsX(); i++){ if (cuthist->GetXaxis()->GetBinCenter(i+1) < gOptEnergyRange[0] or cuthist->GetXaxis()->GetBinCenter(i+1) > gOptEnergyRange[1]){ cuthist->SetBinContent(i+1, 0.0); } } } // Check Flux if (cuthist->Integral() <= 0.0){ THROW("Flux histogram " << iter->second << " has integral <= 0.0"); } // Save fluxwrite->cd(); cuthist->Write(); } std::cout << " ->-> Saved to : " << (out + ".flux.root") << std::endl; fluxmap["beam_inputroot"] = (out + ".flux.root"); fluxwrite->Close(); return fluxmap; } std::string GetFluxDefinition( std::string fluxfile, std::string fluxhist, std::string fluxid ){ // Get base name of flux file as its being copied to NEUT Run Directory std::vector splitfluxfile = GeneralUtils::ParseToStr(fluxfile,"/"); std::string filename = ""; if (splitfluxfile.size() == 1){ filename = splitfluxfile[0]; } else if (splitfluxfile.size() > 1){ filename = splitfluxfile[splitfluxfile.size()-1]; } else { THROW("NO FILENAME FOR FLUX DEFINITION FOUND!"); } // Build string std::string fluxparams = ""; fluxparams += "EVCT-FILENM \'" + filename + "\' \n"; fluxparams += "EVCT-HISTNM \'" + fluxhist + "\' \n"; fluxparams += "EVCT-INMEV 0 \n"; fluxparams += "EVCT-MPV 3 \n"; // Set PDG Code if (!fluxid.compare("nue")) fluxparams += "EVCT-IDPT 12"; else if (!fluxid.compare("nueb")) fluxparams += "EVCT-IDPT -12"; else if (!fluxid.compare("numu")) fluxparams += "EVCT-IDPT 14"; else if (!fluxid.compare("numub")) fluxparams += "EVCT-IDPT -14"; else if (!fluxid.compare("nutau")) fluxparams += "EVCT-IDPT 16"; else if (!fluxid.compare("nutaub")) fluxparams += "EVCT-IDPT -16"; else { THROW("UNKNOWN FLUX ID GIVEN!"); } return fluxparams; } std::string GetTargetDefinition(std::string target){ LOG(FIT) << "Defining NEUT Target from : " << target << std::endl; // Target is given as either a single PDG, or a combo with the total number of nucleons std::vector trgts = GeneralUtils::ParseToStr(target,","); std::string targetstring = ""; // NEUT only lets us pass C and CH type inputs. // Single Target if (trgts.size() == 1){ int PDG = GeneralUtils::StrToInt(trgts[0]); int Z = TargetUtils::GetTargetZFromPDG(PDG); int N = TargetUtils::GetTargetAFromPDG(PDG) - Z; targetstring += "NEUT-NUMBNDP " + GeneralUtils::IntToStr(Z) + "\n"; targetstring += "NEUT-NUMBNDN " + GeneralUtils::IntToStr(N) + "\n"; targetstring += "NEUT-NUMFREP 0\n"; targetstring += "NEUT-NUMATOM " + GeneralUtils::IntToStr(Z+N) + "\n"; // Combined target } else if (trgts.size() == 3){ int NUCLEONS = GeneralUtils::StrToInt(trgts[0]); std::string target1 = trgts[1]; std::string target2 = trgts[2]; // Parse target strings string::size_type open_bracket = target1.find("["); string::size_type close_bracket = target1.find("]"); string::size_type ibeg = 0; string::size_type iend = open_bracket; string::size_type jbeg = open_bracket+1; string::size_type jend = close_bracket-1; int PDG1 = atoi(target1.substr(ibeg,iend).c_str()); double W1 = atof(target1.substr(jbeg,jend).c_str()); open_bracket = target2.find("["); close_bracket = target2.find("]"); ibeg = 0; iend = open_bracket; jbeg = open_bracket+1; jend = close_bracket-1; int PDG2 = atoi(target2.substr(ibeg,iend).c_str()); double W2 = atof(target2.substr(jbeg,jend).c_str()); // Can only have H as a secondary target! if (PDG1 != 1000010010 && PDG2 != 1000010010){ THROW("NEUT Doesn't support composite targets apart fromn Target+H. E.g. CH"); } // Switch so H is PDG2 if given if (PDG1 == 1000010010 && PDG2 != 1000010010){ PDG1 = PDG2; PDG2 = 1000010010; double temp1 = W1; W1 = W2; W2 = temp1; } // Now build string int Z = TargetUtils::GetTargetZFromPDG(PDG1); int N = TargetUtils::GetTargetAFromPDG(PDG1) - Z; int NHydrogen = int(W2 * double(NUCLEONS)); if (double(W2*double(NUCLEONS)) - (double(NHydrogen)) > 0.5){ NHydrogen++; // awkward rounding bug fix } targetstring += "NEUT-NUMBNDP " + GeneralUtils::IntToStr(Z) + "\n"; targetstring += "NEUT-NUMBNDN " + GeneralUtils::IntToStr(N) + "\n"; targetstring += "NEUT-NUMFREP " + GeneralUtils::IntToStr(NHydrogen) + "\n"; targetstring += "NEUT-NUMATOM " + GeneralUtils::IntToStr(Z+N) + "\n"; } else { THROW("NEUT only supports single targets or ones with a secondary H!"); } return targetstring; } std::string GetEventAndSeedDefinition(int nevents, int seed){ std::string eventdef = ""; eventdef += "EVCT-NEVT " + GeneralUtils::IntToStr(nevents) + "\n"; LOG(FIT) << "Event Definition: " << std::endl; std::cout << " -> EVCT-NEVT : " << nevents << std::endl; return eventdef; } //____________________________________________________________________________ int main(int argc, char ** argv) { LOG(FIT) << "==== RUNNING neut_nuisance Event Generator =====" << std::endl; GetCommandLineArgs(argc,argv); std::string neutroot = std::string(getenv("NEUT_ROOT")) + "/src/neutsmpl/"; // Calculate the dynamic modes definition // Read Target string std::string targetparamsdef = GetTargetDefinition(gOptTargetDef); //____________________________ // NEUT doesn't let us do combined flux inputs so have to loop over each flux. std::map newfluxdef = MakeNewFluxFile(gOptFluxDef,gOptOutputFile); // Quick fix to make flux also save to pid_time.root.flux.root std::stringstream ss; ss << getpid() << "_" << time(NULL) << ".root"; newfluxdef = MakeNewFluxFile(gOptFluxDef,ss.str()); // Copy this file to the NEUT working directory LOG(FIT) << "Copying flux to NEUT working directory" << std::endl; system(("cp -v " + newfluxdef["beam_inputroot"] + " " + neutroot + "/").c_str()); TFile* fluxrootfile = new TFile( newfluxdef["beam_inputroot"].c_str(), "READ"); // Setup possible beams and get relative fractions std::vector possiblefluxids; std::vector fluxfractions; possiblefluxids.push_back("nue"); possiblefluxids.push_back("nueb"); possiblefluxids.push_back("numu"); possiblefluxids.push_back("numub"); possiblefluxids.push_back("tau"); possiblefluxids.push_back("taub"); // Total up integrals double totintflux = 0.0; for (size_t i = 0; i < possiblefluxids.size(); i++){ if (newfluxdef["beam_inputroot_" + possiblefluxids[i]].empty()){ fluxfractions.push_back(0.0); } else { TH1D* fluxhist = (TH1D*) fluxrootfile->Get( newfluxdef["beam_inputroot_" + possiblefluxids[i]].c_str() ); if (!fluxhist){ THROW("FLUX HIST : " << newfluxdef["beam_inputroot_" + possiblefluxids[i]] << " not found!"); } fluxfractions.push_back( fluxhist->Integral() ); totintflux += fluxhist->Integral(); } } fluxrootfile->Close(); // Now loop over and actually generate jobs! for (size_t i = 0; i < possiblefluxids.size(); i++){ if (fluxfractions[i] == 0.0) continue; // Get number of events for this subbeam int nevents = int( double(gOptNumberEvents) * fluxfractions[i] / totintflux ); std::cout << "NEVENTS = " << gOptNumberEvents << " " << fluxfractions[i] <<" " << totintflux << " " << nevents << std::endl; std::string eventparamsdef = GetEventAndSeedDefinition(nevents, gOptSeed); bool neutrino = true; if (possiblefluxids[i].find("b") != std::string::npos) neutrino = false; std::string dynparamsdef = GetDynamicModes(gOptGeneratorList, neutrino); std::string fluxparamsdef = GetFluxDefinition( newfluxdef["beam_inputroot"], newfluxdef["beam_inputroot_" + possiblefluxids[i]], possiblefluxids[i] ); LOG(FIT) << "==== Generating CardFiles NEUT! ===" << std::endl; std::cout << dynparamsdef << std::endl; std::cout << targetparamsdef << std::endl; std::cout << eventparamsdef << std::endl; std::cout << fluxparamsdef << std::endl; // Create card file std::ifstream incardfile; std::ofstream outcardfile; std::string line; incardfile.open( gOptCrossSections.c_str(), ios::in ); outcardfile.open( (gOptOutputFile + "." + possiblefluxids[i] + ".par").c_str(), ios::out ); // Copy base card file if (incardfile.is_open()){ while( getline (incardfile, line) ){ outcardfile << line << '\n'; } } else { THROW( "Cannot find card file : " << gOptCrossSections ); } // Now copy our strings outcardfile << eventparamsdef<< '\n'; outcardfile << dynparamsdef << '\n'; outcardfile << targetparamsdef<< '\n'; outcardfile << fluxparamsdef<< '\n'; // Close card and keep name for future use. outcardfile.close(); } LOG(FIT) << "GENERATING" << std::endl; for (size_t i = 0; i < possiblefluxids.size(); i++){ if (fluxfractions[i] == 0.0) continue; int nevents = int( double(gOptNumberEvents) * fluxfractions[i] / totintflux ); if (nevents <= 0) continue; std::string cardfile = ExpandPath(gOptOutputFile + "." + possiblefluxids[i] + ".par"); std::string outputfile = ExpandPath(gOptOutputFile + "." + possiblefluxids[i] + ".root"); std::string basecardfile = GetBaseName(cardfile); std::string baseoutputfile = ss.str(); std::cout << "CARDFILE = " << cardfile << " : " << basecardfile << std::endl; std::cout << "OUTPUTFILE = " << outputfile << " : " << baseoutputfile << std::endl; system(("cp " + cardfile + " " + neutroot).c_str()); std::string cwd = GETCWD(); chdir(neutroot.c_str()); - int attempts = 0; + // int attempts = 0; // while(true){ // Break if too many attempts // attempts++; // if (attempts > 20) continue; // Actually run neutroot2 system(("./neutroot2 " + basecardfile + " " + baseoutputfile).c_str()); // Check the output is valid, sometimes NEUT aborts mid run. // TFile* f = new TFile(baseoutputfile.c_str(),"READ"); // if (!f or f->IsZombie()) continue; // Check neutttree is there and filled correctly. // TTree* tn = (TTree*) f->Get("neuttree"); // if (!tn) continue; // if (tn->GetEntries() < nevents * 0.9) continue; // break; // } // Move the finished file back and clean this directory of card files system(("mv -v " + baseoutputfile + " " + outputfile).c_str()); system(("rm -v " + basecardfile).c_str()); chdir(cwd.c_str()); } return 0; } //____________________________________________________________________________ void GetCommandLineArgs(int argc, char ** argv) { // Check for -h flag. for (int i = 0; i < argc; i++){ if (!std::string(argv[i]).compare("-h")) PrintSyntax(); } // Format is neut -r run_number -n n events std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv); ParserUtils::ParseArgument(args, "-n", gOptNumberEvents, false); if (gOptNumberEvents == -1){ THROW( "No event count passed to neut_NUISANCE!"); } // Flux/Energy Specs ParserUtils::ParseArgument(args, "-e", gOptEnergyDef, false); gOptEnergyRange = GeneralUtils::ParseToDbl(gOptEnergyDef,","); ParserUtils::ParseArgument(args, "-f", gOptFluxDef, false); if (gOptFluxDef.empty() and gOptEnergyRange.size() < 1){ THROW("No flux or energy range given to neut_nuisance!"); } else if (gOptFluxDef.empty() and gOptEnergyRange.size() == 1){ // Fixed energy, make sure -p is given THROW("neut_NUISANCE cannot yet do fixed energy!"); } else if (gOptFluxDef.empty() and gOptEnergyRange.size() == 2){ // Uniform energy range THROW("neut_NUISANCE cannot yet do a uniform energy range!"); } else if (!gOptFluxDef.empty()){ // Try to convert the flux definition if possible. std::string convflux = BeamUtils::ConvertFluxIDs(gOptFluxDef); if (!convflux.empty()) gOptFluxDef = convflux; } else { THROW("Unknown flux energy range combination!"); } ParserUtils::ParseArgument(args, "-t", gOptTargetDef, false); if (gOptTargetDef.empty()){ THROW("No Target passed to neut_nuisance! use the '-t' argument."); } else { std::string convtarget = TargetUtils::ConvertTargetIDs(gOptTargetDef); if (!convtarget.empty()) gOptTargetDef = convtarget; } ParserUtils::ParseArgument(args, "-r", gOptRunNumber, false); ParserUtils::ParseArgument(args, "-o", gOptOutputFile, false); if (gOptOutputFile.empty()){ if (gOptRunNumber == -1) gOptRunNumber = 1; LOG(FIT) << "No output file given! Saving file to : neutgen." << gOptRunNumber << ".neutvect.root" << std::endl; gOptOutputFile = "neutgen." + GeneralUtils::IntToStr(gOptRunNumber) + ".neutvect.root"; } else { // if no run number given leave as is, else add run number. if (gOptRunNumber != -1){ gOptOutputFile += "." + GeneralUtils::IntToStr(gOptRunNumber) + ".root"; } else { gOptRunNumber = 0; } } ParserUtils::ParseArgument(args, "--cross-section", gOptCrossSections, false); if (!gOptCrossSections.compare("Default")){ LOG(FIT) << "No Parameters File passed. Using default neut one." << std::endl; char * const var = getenv("NUISANCE"); if (!var) { std::cout << "Cannot find top level directory! Set the NUISANCE environmental variable" << std::endl; exit(-1); } std::string topnuisancedir = string(var); gOptCrossSections = topnuisancedir + "/neut/Default_params.txt"; } ParserUtils::ParseArgument(args, "--event-generator-list", gOptGeneratorList, false); // Final Check and output if (args.size() > 0){ PrintSyntax(); ParserUtils::CheckBadArguments(args); } LOG(FIT) << "Generating Neut Events with the following properties:" << std::endl << " -> Energy : " << gOptEnergyDef << " (" << gOptEnergyRange.size() << ")" << std::endl << " -> NEvents : " << gOptNumberEvents << std::endl << " -> Generators : " << gOptGeneratorList << std::endl << " -> XSecPars : " << gOptCrossSections << std::endl << " -> Target : " << gOptTargetDef << std::endl << " -> Flux : " << gOptFluxDef << std::endl << " -> Output : " << gOptOutputFile << std::endl << " -> Run : " << gOptRunNumber << std::endl; return; } //____________________________________________________________________________ void PrintSyntax(void) { LOG(FIT) << "\n\n" << "Syntax:" << "\n" << "\n neut_nuisance [-h]" << "\n -n nev" << "\n -f flux_description" << "\n -t target_description" << "\n [ -r run_number ]" << "\n [ -o output_file ]" << "\n [ --cross-section /path/to/params.txt ]" << "\n [ --event-generator-list mode_definition ]" << "\n \n"; LOG(FIT) << "\n\n Arguments:" << "\n" << "\n -n nev" << "\n -> Total number of events to generate (e.g. 2500000)" << "\n" << "\n -f flux_description" << "\n Definition of the flux to be read in from a ROOT file." << "\n" << "\n Multiple histograms can be read in from the same file using" << "\n the format '-f file.root,hist1[pdg1],hist2[pdg2]" << "\n e.g. \'-f ./flux/myfluxfile.root,numu_flux[14],numubar_flux[-14]\'" << "\n" << "\n A flux can also be given according to any of the flux IDs shown" << "\n at the end of this help message." << "\n e.g. \' -f MINERvA_fhc_numu\' " << "\n" << "\n WARNING: NEUT can't actually generate combined fluxes yet" << "\n if you want a composite flux, pass them in as normal, but the app" << "\n will generate you the files seperately, with reduced nevents in each" << "\n so that the statistics are roughly okay." << "\n" << "\n -t target_description" << "\n Definition of the target to be used. Multiple targets can be given." << "\n" << "\n To pass a single target just provide the target PDG" << "\n e.g. \' -t 1000060120 \'" << "\n" << "\n To pass a combined target provide a list containing the following" << "\n \' -t TotalNucleons,Target1[Weight1],Target2[Weight2],.. where the " << "\n TotalNucleons is the total nucleons combined, Target1 is the PDG " << "\n of the first target, and Weight1 is the fractional weight of the " << "\n first target." << "\n e.g. \' -t 13,1000060120[0.9231],1000010010[0.0769] \'" << "\n" << "\n Target can also be specified by the target IDs given at the end of" << "\n this help message." << "\n e.g. \' -t CH2 \'" << "\n" << "\n WARNING: NEUT can only generate A+H targets. E.g. CH or CH2 will work, but " << "\n Fe+Pb will not. You will have to generate each seperately if you want" << "\n something other than A+NH." << "\n" << "\n -r run_number" << "\n run number ID that can be used when generating large samples in small " << "\n jobs. Must be an integer. When given neut_nuisance will update the " << "\n output file from 'output.root' to 'output.root.run_number.root'" << "\n" << "\n -o output_file" << "\n Path to the output_file you want to save events to." << "\n" << "\n If this is not given but '-r' is then events will be saved to " << "\n the file 'neutgen.run_number.neutvect.root'" << "\n" << "\n If a run number is given alongside '-o' then events will be saved " << "\n to 'output.root.run_number.root'" << "\n" << "\n --cross-section /path/to/params.txt" << "\n Path to the neut model definition. If this is not given, then this " << "\n will default to $NUISANCE/data/neut/Default_params.txt" << "\n" << "\n Look in $NUISANCE/data/neut/Default_params.txt for examples when " << "\n writing your own card files." << "\n" << "\n --event-generator-list mode_definition" << "\n Name of modes to run. This sets the CRS and CRSB values in NEUT." << "\n e.g. --event-generator-list Default+MEC" << "\n" << "\n Allowed mode_definitions are given at the end of this help message." << "\n" << "\n\n"; std::cout << "-----------------"<. *******************************************************************************/ #include "MINERvA_SignalDef.h" #include "MINERvA_CCQE_XSec_1DQ2_antinu.h" //******************************************************************** MINERvA_CCQE_XSec_1DQ2_antinu::MINERvA_CCQE_XSec_1DQ2_antinu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCQE_XSec_1DQ2_antinu sample. \n" \ "Target: CH \n" \ "Flux: MINERvA Forward Horn Current Numubar \n" \ "Signal: True CCQE/2p2h defined at the vertex level \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("Q^{2}_{QE} (GeV^{2})"); fSettings.SetYTitle("d#sigma/dQ_{QE}^{2} (cm^{2}/GeV^{2})"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL"); fSettings.SetEnuRange(1.5, 10.0); fSettings.DefineAllowedTargets("C,H"); isFluxFix = !fSettings.Found("name", "_oldflux"); fullphasespace = !fSettings.Found("name", "_20deg"); // CCQELike plot information fSettings.SetTitle("MINERvA_CCQE_XSec_1DQ2_antinu"); fSaveExtra = FitPar::Config().GetParB("MINERvASaveExtraCCQE"); std::string basedir = FitPar::GetDataBase() + "/MINERvA/CCQE/"; std::string datafilename = ""; std::string covarfilename = ""; // Full Phase Space if (fullphasespace) { if (isFluxFix) { if (fIsShape) { ERR(WRN) << "SHAPE likelihood comparison not available for MINERvA " << "datasets with fixed flux information. NUISANCE will scale MC to match " << "data normalization but full covariance will be used. " << std::endl; } datafilename = "Q2QE_numubar_data_fluxfix.txt"; covarfilename = "Q2QE_numubar_covar_fluxfix.txt"; // Correlation Matrix } else { if (fIsShape) { datafilename = "Q2QE_numubar_data_SHAPE-extracted.txt"; covarfilename = "Q2QE_numubar_covar_SHAPE-extracted.txt"; // correlation } else { datafilename = "Q2QE_numubar_data.txt"; covarfilename = "Q2QE_numubar_covar.txt"; // Correlation } } // Restricted Phase Space } else { if (isFluxFix) { if (fIsShape) { ERR(WRN) << "SHAPE likelihood comparison not available for MINERvA " << "datasets with fixed flux information. NUISANCE will scale MC to match " << "data normalization but full covariance will be used. " << std::endl; } datafilename = "20deg_Q2QE_numubar_data_fluxfix.txt"; covarfilename = "20deg_Q2QE_numubar_covar_fluxfix.txt"; // Correlation } else { if (fIsShape) { datafilename = "20deg_Q2QE_numubar_data_SHAPE-extracted.txt"; covarfilename = "20deg_Q2QE_numubar_covar_SHAPE-extracted.txt"; // Correlation } else { datafilename = "20deg_Q2QE_numubar_data.txt"; covarfilename = "20deg_Q2QE_numubar_covar.txt"; // Correlation } } } fSettings.SetDataInput( basedir + datafilename ); fSettings.SetCovarInput( basedir + covarfilename ); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) * 13. / 7. / TotalIntegratedFlux(); // Plot Setup ------------------------------------------------------- SetDataFromTextFile( fSettings.GetDataInput() ); SetCorrelationFromTextFile( fSettings.GetCovarInput() ); if (fSaveExtra){ fExtra_Eav = new TH1D((fSettings.GetName() + "Eav").c_str(),"Eav",30, 0.0, 1.0); fExtra_Eav_MODES = new MINERvAUtils::ModeStack((fSettings.GetName() + "EavMODES").c_str(),"EavMODES", fExtra_Eav); // fExtra_Eav_MODES->SetReverseStack(); SetAutoProcess(fExtra_Eav); SetAutoProcess(fExtra_Eav_MODES); fExtra_EavQ2 = new TH2D((fSettings.GetName() + "EavQ2").c_str(),"EavQ2",12,0.0,2.0,50,0.0,1.0); fExtra_EavQ2_MODES = new MINERvAUtils::ModeStack((fSettings.GetName() + "EavQ2MODES").c_str(),"EavQ2MODES", fExtra_EavQ2); // fExtra_EavQ2_MODES->SetReverseStack(); SetAutoProcess(fExtra_EavQ2); SetAutoProcess(fExtra_EavQ2_MODES); fEavQ2Cut = new TF1((fSettings.GetName() + "f1").c_str(),"0.03 + 0.3*x",0.0,2.0); SetAutoProcess(fEavQ2Cut); } // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CCQE_XSec_1DQ2_antinu::FillEventVariables(FitEvent *event) { //******************************************************************** if (event->NumFSParticle(-13) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(-13)->fP; double ThetaMu = Pnu.Vect().Angle(Pmu.Vect()); double q2qe = FitUtils::Q2QErec(Pmu, cos(ThetaMu), 30., true); fXVar = q2qe; if (fSaveExtra){ double pm = 938.2720813; double pe = pm + 120.0; FitParticle* fakeproton = new FitParticle(0.0,0.0,sqrt(pe*pe-pm*pm), pe, 2212, kFinalState); double range = MINERvAUtils::RangeInScintillator(fakeproton, 100); pm = 139.57018; pe = pm + 65; - FitParticle* fakepion = new FitParticle(0.0,0.0,sqrt(pe*pe-pm*pm), pe, 211, kFinalState); - double pionrange = MINERvAUtils::RangeInScintillator(fakepion, 100); + // FitParticle* fakepion = new FitParticle(0.0,0.0,sqrt(pe*pe-pm*pm), pe, 211, kFinalState); + // double pionrange = MINERvAUtils::RangeInScintillator(fakepion, 100); double Eav = 0.0; - for (int i = 0; i < event->NParticles(); i++){ + for (uint i = 0; i < event->NParticles(); i++){ if (event->GetParticleState(i) != kFinalState) continue; int pid = event->GetParticlePDG(i); double ParticleEav = 0.0; if (abs(pid) == 13) continue; if (pid != 2112 and pid != 22 and pid != 111){ ParticleEav = MINERvAUtils::GetEDepositOutsideRangeInScintillator(event->GetParticle(i), range) / 1.E3; } else if (pid == 22 or pid == 111){ ParticleEav = event->GetParticle(i)->fP.E()/1.E3; } Eav += ParticleEav; } fExtra_Eav->Fill(Eav); fExtra_Eav_MODES->Fill( event->Mode, Eav ); fExtra_EavQ2->Fill(q2qe,Eav); fExtra_EavQ2_MODES->Fill( event->Mode, q2qe, Eav ); } return; } //******************************************************************** bool MINERvA_CCQE_XSec_1DQ2_antinu::isSignal(FitEvent *event) { //******************************************************************* return SignalDef::isCCQEnumubar_MINERvA(event, EnuMin, EnuMax, fullphasespace); } diff --git a/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_nu.cxx b/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_nu.cxx index 14adf03..fef868f 100644 --- a/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_nu.cxx +++ b/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_nu.cxx @@ -1,199 +1,199 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "MINERvA_SignalDef.h" #include "MINERvA_CCQE_XSec_1DQ2_nu.h" //******************************************************************** MINERvA_CCQE_XSec_1DQ2_nu::MINERvA_CCQE_XSec_1DQ2_nu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "MINERvA_CCQE_XSec_1DQ2_nu sample. \n" \ "Target: CH \n" \ "Flux: MINERvA Forward Horn Current Numu \n" \ "Signal: True CCQE/2p2h defined at the vertex level \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("Q^{2}_{QE} (GeV^{2})"); fSettings.SetYTitle("d#sigma/dQ_{QE}^{2} (cm^{2}/GeV^{2})"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL"); fSettings.SetEnuRange(1.5, 10.0); fSettings.DefineAllowedTargets("C,H"); isFluxFix = !fSettings.Found("name", "_oldflux"); fullphasespace = !fSettings.Found("name", "_20deg"); fSaveExtra = FitPar::Config().GetParB("MINERvASaveExtraCCQE"); // CCQELike plot information fSettings.SetTitle("MINERvA_CCQE_XSec_1DQ2_nu"); std::string basedir = FitPar::GetDataBase() + "/MINERvA/CCQE/"; std::string datafilename = ""; std::string covarfilename = ""; // Full Phase Space if (fullphasespace) { if (isFluxFix) { if (fIsShape) { ERR(WRN) << "SHAPE likelihood comparison not available for MINERvA " << "datasets with fixed flux information. NUISANCE will scale MC to match " << "data normalization but full covariance will be used. " << std::endl; } datafilename = "Q2QE_numu_data_fluxfix.txt"; covarfilename = "Q2QE_numu_covar_fluxfix.txt"; } else { if (fIsShape) { datafilename = "Q2QE_numu_data_SHAPE-extracted.txt"; covarfilename = "Q2QE_numu_covar_SHAPE-extracted.txt"; } else { datafilename = "Q2QE_numu_data.txt"; covarfilename = "Q2QE_numu_covar.txt"; } } // Restricted Phase Space } else { if (isFluxFix) { if (fIsShape) { ERR(WRN) << "SHAPE likelihood comparison not available for MINERvA " << "datasets with fixed flux information. NUISANCE will scale MC to match " << "data normalization but full covariance will be used. " << std::endl; } datafilename = "20deg_Q2QE_numu_data_fluxfix.txt"; covarfilename = "20deg_Q2QE_numu_covar_fluxfix.txt"; } else { if (fIsShape) { datafilename = "20deg_Q2QE_numu_data_SHAPE-extracted.txt"; covarfilename = "20deg_Q2QE_numu_covar_SHAPE-extracted.txt"; } else { datafilename = "20deg_Q2QE_numu_data.txt"; covarfilename = "20deg_Q2QE_numu_covar.txt"; } } } fSettings.SetDataInput( basedir + datafilename ); fSettings.SetCovarInput( basedir + covarfilename ); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 * 13.0 / 6.0 / (fNEvents + 0.)) / TotalIntegratedFlux(); // Plot Setup ------------------------------------------------------- SetDataFromTextFile( fSettings.GetDataInput() ); if (!isFluxFix or !fullphasespace){ SetCorrelationFromTextFile( fSettings.GetCovarInput() ); } else { SetCovarFromTextFile( fSettings.GetCovarInput() ); } if (fSaveExtra){ fExtra_Eav = new TH1D((fSettings.GetName() + "Eav").c_str(),"Eav",30, 0.0, 1.0); fExtra_Eav_MODES = new MINERvAUtils::ModeStack((fSettings.GetName() + "EavMODES").c_str(),"EavMODES", fExtra_Eav); // fExtra_Eav_MODES->SetReverseStack(); SetAutoProcess(fExtra_Eav); SetAutoProcess(fExtra_Eav_MODES); fExtra_EavQ2 = new TH2D((fSettings.GetName() + "EavQ2").c_str(),"EavQ2",12,0.0,2.0,50,0.0,1.0); fExtra_EavQ2_MODES = new MINERvAUtils::ModeStack((fSettings.GetName() + "EavQ2MODES").c_str(),"EavQ2MODES", fExtra_EavQ2); // fExtra_EavQ2_MODES->SetReverseStack(); SetAutoProcess(fExtra_EavQ2); SetAutoProcess(fExtra_EavQ2_MODES); fEavQ2Cut = new TF1((fSettings.GetName() + "f1").c_str(),"0.03 + 0.3*x",0.0,2.0); SetAutoProcess(fEavQ2Cut); } // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** void MINERvA_CCQE_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) { //******************************************************************** if (event->NumFSParticle(13) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double ThetaMu = Pnu.Vect().Angle(Pmu.Vect()); double q2qe = FitUtils::Q2QErec(Pmu, cos(ThetaMu), 34., true); // Set binning variable fXVar = q2qe; // Calculate from scratch the vertex size if (fSaveExtra){ double pm = 938.2720813; double pe = pm + 120.0; FitParticle* fakeproton = new FitParticle(0.0,0.0,sqrt(pe*pe-pm*pm), pe, 2212, kFinalState); double range = MINERvAUtils::RangeInScintillator(fakeproton, 100); pm = 139.57018; pe = pm + 65; - FitParticle* fakepion = new FitParticle(0.0,0.0,sqrt(pe*pe-pm*pm), pe, 211, kFinalState); - double pionrange = MINERvAUtils::RangeInScintillator(fakepion, 100); + // FitParticle* fakepion = new FitParticle(0.0,0.0,sqrt(pe*pe-pm*pm), pe, 211, kFinalState); + // double pionrange = MINERvAUtils::RangeInScintillator(fakepion, 100); double Eav = 0.0; - for (int i = 0; i < event->NParticles(); i++){ + for (uint i = 0; i < event->NParticles(); i++){ if (event->GetParticleState(i) != kFinalState) continue; int pid = event->GetParticlePDG(i); double ParticleEav = 0.0; if (pid == 13) continue; if (pid != 2112 and pid != 22 and pid != 111){ ParticleEav = MINERvAUtils::GetEDepositOutsideRangeInScintillator(event->GetParticle(i), range) / 1.E3; } else if (pid == 22 or pid == 111){ ParticleEav = event->GetParticle(i)->fP.E()/1.E3; } Eav += ParticleEav; } fExtra_Eav->Fill(Eav); fExtra_Eav_MODES->Fill( event->Mode, Eav ); fExtra_EavQ2->Fill(q2qe,Eav); fExtra_EavQ2_MODES->Fill( event->Mode, q2qe, Eav ); } return; } //******************************************************************** bool MINERvA_CCQE_XSec_1DQ2_nu::isSignal(FitEvent *event) { //******************************************************************* return SignalDef::isCCQEnumu_MINERvA(event, EnuMin, EnuMax, fullphasespace); } diff --git a/src/Routines/MinimizerRoutines.cxx b/src/Routines/MinimizerRoutines.cxx index c258da8..9d1dc26 100755 --- a/src/Routines/MinimizerRoutines.cxx +++ b/src/Routines/MinimizerRoutines.cxx @@ -1,1510 +1,1510 @@ // 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" /* 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", + "SteepDesc,GSLSimAn,FixAtLim,FixAtLimBreak," + "Chi2Scan1D,Chi2Scan2D,Contours,ErrorBands," "DataToys"); }; //************************************* MinimizerRoutines::~MinimizerRoutines() { //************************************* }; /* Input Functions */ //************************************* 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.AddS("cardfile", fCardFile); if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile); if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy); // Load XML Cardfile configuration.LoadConfig( fCompKey.GetS("cardfile"), ""); // Add CMD XML Structs for (size_t i = 0; i < xmlcmds.size(); i++) { configuration.AddXMLLine(xmlcmds[i]); } // 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.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml"); // Add Error Verbo Lines verbocount += Config::Get().GetParI("VERBOSITY"); errorcount += Config::Get().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) { parnom = FitBase::RWAbsToSigma( partype, parname, parnom ); parlow = FitBase::RWAbsToSigma( partype, parname, parlow ); parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh ); parstep = FitBase::RWAbsToSigma( partype, parname, parstep ); } else if (parstate.find("FRAC") != std::string::npos) { 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 = ""; // 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 = ""; } // make minimizer if (fMinimizer) delete fMinimizer; 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")) { 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 FitPar::Config().SetParI("input.maxevents", lowstatsevents); FitPar::Config().SetParI("VERBOSITY", 3); SetupFCN(); RunFitRoutine(trueroutine); FitPar::Config().SetParI("input.maxevents", maxevents); SetupFCN(); FitPar::Config().SetParI("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::Get().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(); }