diff --git a/MINERvA_NOvA_WS/FlatAnalysis.cc b/MINERvA_NOvA_WS/FlatAnalysis.cc deleted file mode 100644 index e428144..0000000 --- a/MINERvA_NOvA_WS/FlatAnalysis.cc +++ /dev/null @@ -1,518 +0,0 @@ -int MyPalette[255]; - -// Make the TPalette -void SetCustomPalette() { - // Uncomment these to use the blue->white->red palette (good for correlation matrices) - TColor::InitializeColors(); - const Int_t NRGBs = 5; - const Int_t NCont = 255; - - Double_t stops[NRGBs] = { 0.00, 0.25, 0.50, 0.75, 1.00 }; - Double_t red[NRGBs] = { 0.00, 0.25, 1.00, 1.00, 0.50 }; - Double_t green[NRGBs] = { 0.00, 0.25, 1.00, 0.25, 0.00 }; - Double_t blue[NRGBs] = { 0.50, 1.00, 1.00, 0.25, 0.00 }; - - int start = TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); - gStyle->SetNumberContours(NCont); - for (int i = 0; i < NCont; i++) { - MyPalette[i] = start+i; - } -} -void FlatAnalysis(std::string filename, std::string novaname) { - - gStyle->SetOptStat(0); - SetCustomPalette(); - - TFile *file = new TFile(filename.c_str(), "OPEN"); - TTree *tree = (TTree*)file->Get("FlatTree_VARS"); - - TFile *filenova = new TFile(novaname.c_str(), "OPEN"); - TTree *treenova = (TTree*)filenova->Get("FlatTree_VARS"); - - // Mode, cc, PDGnu, tgtA, PDGLep - // ELep, CosLep, Q2, q0, q3, Enu_true, Enu_QE, Q2_QE, W_nuc_rest, W, x, y - // nfsp, px[nfsp], py[nfsp], pz[nfsp], E[nfsp], pdg[nfsp] - // Weight (useless), InputWeight (useless), RWWeight (useless) - // fScaleFactor, CustomWeight, CustomWeightArray[6] - // - // flagCCINC, flagNCINC - // flagCCQE, flagCC0pi, flagCCQELike, flagNCEL, flagNC0pi - // flagCCcoh, flagNCcoh, - // flagCC1pip, flagNC1pip - // flagCC1pim, flagNC1pim - // flagCC1pi0, flagNC1pi0 - // - double scalefactor = tree->GetMinimum("fScaleFactor"); - double scalefactornova = treenova->GetMinimum("fScaleFactor"); - if (scalefactor < 1E-42) { - scalefactor *= 1E42; - scalefactornova *= 1E42; - } - std::cout << scalefactor << std::endl; - std::cout << scalefactornova << std::endl; - - //const int WeightNumbers = 7; - const int WeightNumbers = 1; - std::string WeightType[WeightNumbers]; - // The total weights - //WeightType[0] = "1"; - WeightType[0] = "Weight"; - //WeightType[1] = "CustomWeight"; - //WeightType[2] = "CustomWeightArray[0]"; - //WeightType[3] = "CustomWeightArray[1]"; - //WeightType[4] = "CustomWeightArray[2]"; - //WeightType[5] = "CustomWeightArray[3]"; - //WeightType[6] = "CustomWeightArray[4]"; - //WeightType[7] = "CustomWeightArray[5]"; - - std::string Names[WeightNumbers]; - //Names[0] = "No Weight"; - Names[0] = "Total Weight"; - //Names[2] = "M_{A}^{QE} Weight"; - //Names[3] = "Non-res Weight"; - //Names[4] = "RPA QE Weight"; - //Names[5] = "RPA RES Weight"; - //Names[6] = "MEC Weight"; - //Names[7] = "Total Weight"; - - // What we want to draw - std::vector Draws2D; - //Draws2D.push_back("ELep:CosLep"); - Draws2D.push_back("Q2:W"); - Draws2D.push_back("q0:q3"); - Draws2D.push_back("Q2:Enu_true"); - Draws2D.push_back("Eav:Enu_true"); - - const int nModes = 8; - std::string Modes[nModes]; - Modes[0]="Mode<30"; - Modes[1]="Mode==1"; - Modes[2]="Mode==2"; - Modes[3]="Mode==11"; - Modes[4]="Mode==12"; - Modes[5]="Mode==13"; - Modes[6]="Mode==21"; - Modes[7]="Mode==26"; - - std::string ModeNames[nModes]; - ModeNames[0]="CC-Inclusive"; - ModeNames[1]="CCQE"; - ModeNames[2]="2p2h"; - ModeNames[3]="CC1#pi^{+}1p"; - ModeNames[4]="CC1#pi^{0}"; - ModeNames[5]="CC1#pi^{+}1n"; - ModeNames[5]="CC Multi-#pi"; - ModeNames[6]="CC DIS"; - - std::string LimX2D[Draws2D.size()]; - //LimX2D[0]="0.0, 1.0"; - LimX2D[0]="1.0, 2.0"; - LimX2D[1]="0.0, 1.5"; - LimX2D[2]="1.0, 3.0"; - LimX2D[3]="1.0, 3.0"; - - std::string LimY2D[Draws2D.size()]; - //LimY2D[0]="0.0, 3.0"; - LimY2D[0]="0.0, 2.0"; - LimY2D[1]="0.0, 1.5"; - LimY2D[2]="0.0, 2.0"; - LimY2D[3]="0.0, 0.5"; - - std::vector Draws1D; - //Draws1D.push_back("ELep"); - //Draws1D.push_back("CosLep"); - Draws1D.push_back("Q2"); - Draws1D.push_back("W"); - Draws1D.push_back("x"); - Draws1D.push_back("y"); - Draws1D.push_back("Eav"); - - std::string LimX1D[Draws1D.size()]; - LimX1D[0] = "0.0, 2.0"; - LimX1D[1] = "1.0, 3.0"; - LimX1D[2] = "0.0, 1.0"; - LimX1D[3] = "0.0, 1.0"; - LimX1D[4] = "0.0, 0.5"; - - TCanvas *canv = new TCanvas("canv", "canv", 1024, 1024); - canv->SetLeftMargin(canv->GetLeftMargin()*1.1); - canv->SetBottomMargin(canv->GetBottomMargin()*0.8); - canv->SetRightMargin(canv->GetRightMargin()*1.5); - canv->SetTopMargin(canv->GetTopMargin()*0.8); - std::string outname = filename; - outname = outname.substr(0, outname.find(".root")); - outname += "_output"; - canv->Print((outname+".pdf[").c_str()); - - // The number of bins - int nBins = 70; - - TFile *output = new TFile((outname+".root").c_str(), "RECREATE"); - - // Draw the 2D - for (int i = 0; i < Draws2D.size(); ++i) { - - // Split the 2D into two parts - std::string ypart = Draws2D[i].substr(0, Draws2D[i].find(":")); - std::string xpart = Draws2D[i].substr(Draws2D[i].find(":")+1, Draws2D[i].length()); - - for (int k = 0; k < nModes; ++k) { - std::string ModeString = Modes[k]; - while (ModeString.find("(") != std::string::npos) { - ModeString.replace(ModeString.find("("), 1, "_"); - } - while (ModeString.find(")") != std::string::npos) { - ModeString.replace(ModeString.find(")"), 1, "_"); - } - while (ModeString.find("=") != std::string::npos) { - ModeString.replace(ModeString.find("="), 1, "_"); - } - - std::stringstream refname; - refname << ypart << ":" << xpart << ">>" << ypart << "_" << xpart << "_" << "noscale" << "_" << ModeString << "(" << nBins << ", "; - std::stringstream ss2; - ss2 << ", " << nBins << ", "; - std::string DrawCmd = refname.str() + LimX2D[i] + ss2.str() + LimY2D[i] + std::string(")"); - std::string WeightDraw = Modes[k]; - file->cd(); - //std::cout << DrawCmd << std::endl; - //std::cout << WeightDraw << std::endl; - tree->Draw(DrawCmd.c_str(), WeightDraw.c_str(), "colz"); - std::string refname_str = ypart+"_"+xpart+"_noscale_"+ModeString; - - TH2D *refplot = (TH2D*)gDirectory->Get(refname_str.c_str())->Clone(); - //std::cout << refname_str << std::endl; - refplot->Scale(scalefactor, "width"); - refplot->GetXaxis()->SetTitle(xpart.c_str()); - refplot->GetYaxis()->SetTitle(ypart.c_str()); - - gStyle->SetPalette(51); - // Make the string to draw - std::stringstream ssnovaref; - ssnovaref << ypart << ":" << xpart << ">>" << ypart << "_" << xpart << "_noscale_" << "CustomWeight" << "_" << ModeString << "_nova(" << nBins << ", "; - - filenova->cd(); - std::string DrawCmdnova = ssnovaref.str() + LimX2D[i] + ss2.str() + LimY2D[i] + std::string(")"); - std::string WeightDrawnova = Modes[k]; - //std::cout << DrawCmdnova << std::endl; - //std::cout << WeightDrawnova << std::endl; - treenova->Draw(DrawCmdnova.c_str(), WeightDrawnova.c_str(), "colz"); - std::string refname_strnova = ypart+"_"+xpart+"_noscale_CustomWeight_"+ModeString+"_nova"; - //std::cout << refname_strnova << std::endl; - TH2D *refplotnova = (TH2D*)gDirectory->Get(refname_strnova.c_str())->Clone(); - refplotnova->Scale(scalefactornova, "width"); - refplotnova->GetXaxis()->SetTitle(xpart.c_str()); - refplotnova->GetYaxis()->SetTitle(ypart.c_str()); - - // Loop over the different weights - for (int j = 0; j < WeightNumbers; ++j) { - - std::string WeightString = WeightType[j]; - while (WeightString.find("[")!=std::string::npos) { - WeightString.replace(WeightString.find("["), 1, "_"); - } - while (WeightString.find("]")!=std::string::npos) { - WeightString.replace(WeightString.find("]"), 1, "_"); - } - - // Make the string to draw - std::stringstream ss; - ss << ypart << ":" << xpart << ">>" << ypart << "_" << xpart << "_" << WeightString << "_" << ModeString << "(" << nBins << ", "; - gStyle->SetPalette(51); - - DrawCmd = ss.str() + LimX2D[i] + ss2.str() + LimY2D[i] + std::string(")"); - WeightDraw = WeightType[j] + "*(" + Modes[k] + ")"; - canv->Clear(); - file->cd(); - tree->Draw(DrawCmd.c_str(), WeightDraw.c_str(), "colz"); - - std::string getname = ypart+"_"+xpart+"_"+WeightString+"_"+ModeString; - TH2D *plot = (TH2D*)gDirectory->Get(getname.c_str())->Clone(); - plot->Scale(scalefactor, "width"); - - // Make the string to draw - std::stringstream ssnova; - ssnova << ypart << ":" << xpart << ">>" << ypart << "_" << xpart << "_" << "CustomWeight" << "_" << ModeString << "_nova(" << nBins << ", "; - - DrawCmdnova = ssnova.str() + LimX2D[i] + ss2.str() + LimY2D[i] + std::string(")"); - WeightDrawnova = std::string("CustomWeight") + "*(" + Modes[k] + ")"; - treenova->Draw(DrawCmdnova.c_str(), WeightDrawnova.c_str(), "colz"); - std::string getnamenova = ypart+"_"+xpart+"_"+"CustomWeight"+"_"+ModeString+"_nova"; - TH2D *plotnova = (TH2D*)gDirectory->Get(getnamenova.c_str())->Clone(); - - plotnova->Scale(scalefactornova, "width"); - - if (plot->Integral() == refplot->Integral()) { - delete plot; - continue; - } - - double maximum = plot->GetMaximum(); - if (maximum < refplot->GetMaximum()) maximum = refplot->GetMaximum(); - plot->GetZaxis()->SetRangeUser(-1*maximum/1000.0, maximum); - refplot->GetZaxis()->SetRangeUser(-1*maximum/1000.0, maximum); - - plotnova->GetZaxis()->SetRangeUser(-1*maximum/1000.0, maximum); - refplotnova->GetZaxis()->SetRangeUser(-1*maximum/1000.0, maximum); - - plot->GetXaxis()->SetTitle(xpart.c_str()); - plot->GetYaxis()->SetTitle(ypart.c_str()); - canv->Clear(); - plot->Draw("colz"); - canv->Print((outname+".pdf").c_str()); - - canv->Clear(); - refplot->Draw("colz"); - canv->Print((outname+".pdf").c_str()); - - plotnova->GetXaxis()->SetTitle(xpart.c_str()); - plotnova->GetYaxis()->SetTitle(ypart.c_str()); - canv->Clear(); - plotnova->Draw("colz"); - canv->Print((outname+".pdf").c_str()); - - canv->Clear(); - refplotnova->Draw("colz"); - canv->Print((outname+".pdf").c_str()); - - gStyle->SetPalette(255, MyPalette); - canv->Update(); - // Make a difference plot - TH2D *diffplot = (TH2D*)plot->Clone(); - plot->Divide(refplot); - //if (j == 0) maximum = plot->GetMaximum(); - plot->GetZaxis()->SetRangeUser(0, 2); - - std::string bla = (Names[j]+std::string(" MINERvA ")+ModeNames[k]); - plot->SetTitle(bla.c_str()); - canv->Clear(); - plot->Draw("colz"); - canv->Print((outname+".pdf").c_str()); - - output->cd(); - plot->Write(); - - TH2D *diffplotnova = (TH2D*)plotnova->Clone(); - plotnova->Divide(refplotnova); - //if (j == 0) maximum = plot->GetMaximum(); - plotnova->GetZaxis()->SetRangeUser(0, 2); - - std::string bla = (Names[j]+std::string(" NOvA ")+ModeNames[k]); - plotnova->SetTitle(bla.c_str()); - canv->Clear(); - plotnova->Draw("colz"); - canv->Print((outname+".pdf").c_str()); - - output->cd(); - plotnova->Write(); - - // Make the difference plot - for (int a = 0; a < plot->GetXaxis()->GetNbins(); ++a) { - for (int b = 0; b < plot->GetYaxis()->GetNbins(); ++b) { - if (diffplot->GetBinContent(a+1,b+1) != 0) { - diffplot->SetBinContent(a+1, b+1, (diffplot->GetBinContent(a+1,b+1)-diffplotnova->GetBinContent(a+1, b+1))/diffplot->GetBinContent(a+1,b+1)); - } else { - diffplot->SetBinContent(a+1, b+1, 0); - } - } - } - diffplot->GetZaxis()->SetRangeUser(-1.0, 1.0); - diffplot->SetTitle((std::string("(MINERvA tuned-NOvA tuned)/MINERvA tuned ")+ModeNames[k]).c_str()); - diffplot->SetName((std::string(diffplot->GetName())+"_diff").c_str()); - - canv->Clear(); - diffplot->Draw("colz"); - canv->Print((outname+".pdf").c_str()); - - delete plot; - delete plotnova; - delete diffplot; - delete diffplotnova; - } - delete refplot; - delete refplotnova; - } - } - - nBins = 40; - /// ///////////// - // Now do the 1D - for (int i = 0; i < Draws1D.size(); ++i) { - - for (int k = 0; k < nModes; ++k) { - std::string ModeString = Modes[k]; - while (ModeString.find("(") != std::string::npos) { - ModeString.replace(ModeString.find("("), 1, "_"); - } - while (ModeString.find(")") != std::string::npos) { - ModeString.replace(ModeString.find(")"), 1, "_"); - } - while (ModeString.find("=") != std::string::npos) { - ModeString.replace(ModeString.find("="), 1, "_"); - } - - // Make the string to draw - std::stringstream refname; - refname << Draws1D[i] << ">>" << Draws1D[i] << "_" << "noscale" << "_" << ModeString << "(" << nBins << ", "; - std::string DrawCmd = refname.str() + LimX1D[i] + std::string(")"); - - std::string WeightDraw = std::string("(") + Modes[k] + ")"; - file->cd(); - tree->Draw(DrawCmd.c_str(), WeightDraw.c_str(), "colz"); - std::string refname_str = Draws1D[i]+"_noscale_"+ModeString; - TH1D *refplot = (TH1D*)gDirectory->Get(refname_str.c_str())->Clone(); - refplot->Scale(scalefactor, "width"); - refplot->GetXaxis()->SetTitle(Draws1D[i].c_str()); - - // Make the string to draw - std::stringstream ssnovaref; - ssnovaref << Draws1D[i] << ">>" << Draws1D[i] << "_noscale_" << "CustomWeight" << "_" << ModeString << "_nova(" << nBins << ", "; - - filenova->cd(); - std::string DrawCmdnova = ssnovaref.str() + LimX1D[i] + std::string(")"); - std::string WeightDrawnova = std::string("(") + Modes[k] + ")"; - treenova->Draw(DrawCmdnova.c_str(), WeightDrawnova.c_str(), "colz"); - std::string refname_strnova = Draws1D[i]+"_noscale_CustomWeight_"+ModeString+"_nova"; - //std::cout << refname_strnova << std::endl; - TH1D *refplotnova = (TH1D*)gDirectory->Get(refname_strnova.c_str())->Clone(); - refplotnova->Scale(scalefactornova, "width"); - refplotnova->GetXaxis()->SetTitle(Draws1D[i].c_str()); - - // Loop over the different weights - for (int j = 0; j < WeightNumbers; ++j) { - - std::string WeightString = WeightType[j]; - while (WeightString.find("[")!=std::string::npos) { - WeightString.replace(WeightString.find("["), 1, "_"); - } - while (WeightString.find("]")!=std::string::npos) { - WeightString.replace(WeightString.find("]"), 1, "_"); - } - - // Make the string to draw - std::stringstream ss; - ss << Draws1D[i] << ">>" << Draws1D[i] << "_" << WeightString << "_" << ModeString << "(" << nBins << ", "; - - DrawCmd = ss.str() + LimX1D[i] + std::string(")"); - - WeightDraw = WeightType[j] + "*(" + Modes[k] + ")"; - canv->Clear(); - file->cd(); - tree->Draw(DrawCmd.c_str(), WeightDraw.c_str(), "colz"); - std::string getname = Draws1D[i]+"_"+WeightString+"_"+ModeString; - TH1D *plot = (TH1D*)gDirectory->Get(getname.c_str())->Clone(); - plot->Scale(scalefactor, "width"); - - // Make the string to draw - std::stringstream ssnova; - ssnova << Draws1D[i] << ">>" << Draws1D[i] << "_" << "CustomWeight" << "_" << ModeString << "_nova(" << nBins << ", "; - DrawCmdnova = ssnova.str() + LimX1D[i] + std::string(")"); - WeightDrawnova = std::string("CustomWeight") + "*(" + Modes[k] + ")"; - canv->Clear(); - treenova->Draw(DrawCmdnova.c_str(), WeightDrawnova.c_str(), "colz"); - std::string getnamenova = Draws1D[i]+"_CustomWeight"+"_"+ModeString+"_nova"; - TH2D *plotnova = (TH2D*)gDirectory->Get(getnamenova.c_str())->Clone(); - plotnova->Scale(scalefactornova, "width"); - - if (plot->Integral() == refplot->Integral()) { - delete plot; - continue; - } - - double maximum = plot->GetMaximum(); - if (maximum < refplot->GetMaximum()) maximum = refplot->GetMaximum(); - maximum *= 1.5; - plot->GetYaxis()->SetRangeUser(maximum, maximum); - refplot->GetYaxis()->SetRangeUser(0, maximum); - - plotnova->GetYaxis()->SetRangeUser(0, maximum); - refplotnova->GetYaxis()->SetRangeUser(0, maximum); - - plot->GetXaxis()->SetTitle(Draws1D[i].c_str()); - canv->Clear(); - refplot->SetLineColor(kRed); - refplot->SetLineWidth(2); - refplot->SetLineStyle(kDashed); - refplot->Draw(); - - plot->SetLineColor(kRed); - plot->SetLineWidth(2); - plot->Draw("same"); - - plotnova->GetXaxis()->SetTitle(Draws1D[i].c_str()); - refplotnova->SetLineColor(kBlue); - refplotnova->SetLineStyle(kDashed); - refplotnova->SetLineWidth(2); - refplotnova->Draw("same"); - - plotnova->SetLineColor(kBlue); - plotnova->SetLineWidth(2); - plotnova->Draw("same"); - - TLegend *leg = new TLegend(0.5, 0.5, 0.95, 0.95); - leg->AddEntry(plot, "MINERvA tuned", "l"); - leg->AddEntry(refplot, "MINERvA untuned", "l"); - leg->AddEntry(plotnova, "NOvA tuned","l"); - leg->AddEntry(refplotnova, "NOvA untuned", "l"); - leg->SetBorderSize(0); - leg->SetFillStyle(0); - leg->SetFillColor(0); - leg->Draw("same"); - - canv->Print((outname+".pdf").c_str()); - delete leg; - - plot->Divide(refplot); - plot->GetYaxis()->SetRangeUser(0, 2); - TLine *line = new TLine(plot->GetXaxis()->GetBinLowEdge(1), 1, plot->GetXaxis()->GetBinLowEdge(plot->GetXaxis()->GetNbins()+1), 1); - line->SetLineStyle(kDashed); - line->SetLineWidth(2); - line->SetLineColor(kBlack); - - std::string bla = (Names[j]+std::string(" Comp ")+ModeNames[k]); - plot->SetTitle(bla.c_str()); - plot->SetLineColor(kRed); - plot->SetLineWidth(2); - canv->Clear(); - plot->Draw(); - - plotnova->Divide(refplotnova); - plotnova->GetYaxis()->SetRangeUser(0, 2); - TLine *line = new TLine(plot->GetXaxis()->GetBinLowEdge(1), 1, plot->GetXaxis()->GetBinLowEdge(plot->GetXaxis()->GetNbins()+1), 1); - line->SetLineStyle(kDashed); - line->SetLineWidth(2); - line->SetLineColor(kBlack); - - std::string bla = (Names[j]+std::string(" NOvA ")+ModeNames[k]); - plotnova->SetTitle(bla.c_str()); - plotnova->SetLineColor(kBlue); - plotnova->SetLineWidth(2); - plotnova->Draw("same"); - line->Draw("same"); - - TLegend *leg = new TLegend(0.5, 0.5, 0.95, 0.95); - leg->AddEntry(plot, "MINERvA tuned/untuned", "l"); - leg->AddEntry(plotnova, "NOvA tuned/untuned","l"); - leg->SetBorderSize(0); - leg->SetFillStyle(0); - leg->SetFillColor(0); - leg->Draw("same"); - - canv->Print((outname+".pdf").c_str()); - - output->cd(); - plot->Write(); - plotnova->Write(); - - delete plot; - delete plotnova; - delete line; - } - delete refplot; - delete refplotnova; - } - } - - canv->Print((outname+".pdf]").c_str()); - output->Close(); -} diff --git a/MINERvA_NOvA_WS/exp_list.xml b/MINERvA_NOvA_WS/exp_list.xml deleted file mode 100644 index 9f71121..0000000 --- a/MINERvA_NOvA_WS/exp_list.xml +++ /dev/null @@ -1,149 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/MINERvA_NOvA_WS/minerva_tune.xml b/MINERvA_NOvA_WS/minerva_tune.xml deleted file mode 100644 index 6a32c67..0000000 --- a/MINERvA_NOvA_WS/minerva_tune.xml +++ /dev/null @@ -1,38 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/MINERvA_NOvA_WS/plotscript.cc b/MINERvA_NOvA_WS/plotscript.cc deleted file mode 100644 index acdeb05..0000000 --- a/MINERvA_NOvA_WS/plotscript.cc +++ /dev/null @@ -1,245 +0,0 @@ -void plotscript(std::string file1, std::string file2, std::string file3, std::string file4) { - // file1 untuned MINERvA - // file2 untuned NOvA - // file3 tuned MINERvA - // file4 tuned NOvA - TFile *f1 = new TFile(file1.c_str(), "OPEN"); - TFile *f2 = new TFile(file2.c_str(), "OPEN"); - TFile *f3 = new TFile(file3.c_str(), "OPEN"); - TFile *f4 = new TFile(file4.c_str(), "OPEN"); - - // Loop through all - TIter next(f1->GetListOfKeys()); - TKey *key; - std::string name; - - TCanvas *canv = new TCanvas("canv", "canv", 1024, 1024); - canv->SetRightMargin(0.05); - canv->Print((file1+"_"+file2+"_comp.pdf[").c_str()); - - // Loop through all entries - while ((key = (TKey*)next())) { - - // Get name of object - name = std::string(key->GetName()); - - if (name.find("_data") != std::string::npos && name.find("SHAPE") == std::string::npos && name.find("RATIO") == std::string::npos && name.find("sample") == std::string::npos) { - if (name == "MiniBooNE_CCQE_XSec_1DQ2_nu_data" || name == "MiniBooNE_CCQE_XSec_1DQ2_antinu_data" || name == "MiniBooNE_CCQELike_XSec_1DQ2_nu_CCQELIKE_BKG_data" || name == "MiniBooNE_CCQELike_XSec_1DQ2_antinu_data") continue; - if (name.find("_1D") != std::string::npos && name.find("MINERvA_CCinc_XSec") != std::string::npos) continue; - - std::string bla = name.substr(0, name.find("_data")); - std::string dataName = name; - - // Get class of object - TClass *cl = gROOT->GetClass(key->GetClassName()); - - // Check TH1D inheritance and that both files contain the name - if (!f2->GetListOfKeys()->Contains(name.c_str())) continue; - if (!f3->GetListOfKeys()->Contains(name.c_str())) continue; - if (!f4->GetListOfKeys()->Contains(name.c_str())) continue; - - if (std::string((key->ReadObj())->ClassName()).compare("TH1D") == 0) { - std::cout << name << std::endl; - - bool error = false; - TH1D *data = (TH1D*)(f1->Get(name.c_str())); - if (data == NULL) { - std::cerr << "no data" << std::endl; - error = true; - } - - TH1D *MC = (TH1D*)(f1->Get((bla+"_MC").c_str())); - if (MC == NULL) { - std::cerr << "no MC" << std::endl; - error = true; - } - TH1D *MC2 = (TH1D*)(f2->Get((bla+"_MC").c_str())); - if (MC2 == NULL) { - std::cerr << "no MC2" << std::endl; - error = true; - } - - TH1D *MC3 = (TH1D*)(f3->Get((bla+"_MC").c_str())); - if (MC3 == NULL) { - std::cerr << "no MC3" << std::endl; - error = true; - } - - TH1D *MC4 = (TH1D*)(f4->Get((bla+"_MC").c_str())); - if (MC4 == NULL) { - std::cerr << "no MC4" << std::endl; - error = true; - } - - if (error) continue; - - data->SetMinimum(0); - for (int i = 0; i < MC->GetNbinsX() + 1; i++) { - MC->SetBinError(i+1, 0); - MC2->SetBinError(i+1, 0); - } - - canv->cd(); - data->Draw(); - MC->Draw("same,hist"); - MC2->Draw("same,hist"); - MC3->Draw("same,hist"); - MC4->Draw("same,hist"); - - double maximum = data->GetMaximum(); - if (maximum < MC->GetMaximum()) maximum = MC->GetMaximum(); - data->GetYaxis()->SetRangeUser(0, maximum*1.7); - - data->SetLineWidth(2); - MC->SetLineWidth(2); - MC2->SetLineWidth(2); - MC3->SetLineWidth(2); - MC4->SetLineWidth(2); - - MC->SetLineColor(kRed-9); - MC->SetMarkerStyle(0); - MC->SetLineStyle(kDashed); - - MC2->SetLineColor(kRed+2); - MC2->SetMarkerStyle(0); - - MC3->SetLineColor(kBlue-9); - MC3->SetMarkerStyle(0); - MC3->SetLineStyle(kDashed); - - MC4->SetLineColor(kBlue+2); - MC4->SetMarkerStyle(0); - - TLegend *leg = new TLegend(0.6, 0.5, 0.95, 0.90); - leg->AddEntry(data, "Data", "le"); - double chi2ndof = std::atof(MC->GetTitle())/MC->GetNbinsX(); - double chi2ndof2 = std::atof(MC2->GetTitle())/MC->GetNbinsX(); - double chi2ndof3 = std::atof(MC3->GetTitle())/MC->GetNbinsX(); - double chi2ndof4 = std::atof(MC4->GetTitle())/MC->GetNbinsX(); - leg->AddEntry(MC, Form("MINERvA untuned, #chi^{2}/ndof=%.2f", chi2ndof), "l"); - leg->AddEntry(MC2, Form("MINERvA tuned, #chi^{2}/ndof=%.2f", chi2ndof2), "l"); - leg->AddEntry(MC3, Form("NOvA untuned, #chi^{2}/ndof=%.2f", chi2ndof3), "l"); - leg->AddEntry(MC4, Form("NOvA tuned, #chi^{2}/ndof=%.2f", chi2ndof4), "l"); - leg->SetBorderSize(0); - leg->SetFillStyle(0); - leg->SetFillColor(0); - leg->Draw("same"); - - canv->Print((file1+"_"+file2+"_comp.pdf").c_str()); - delete leg; - - // Now draw the THStacks - - } else if (std::string((key->ReadObj())->ClassName()).compare("TH2D") == 0) { - std::cout << name << std::endl; - - // Do projections of the TH2D - bool error = false; - TH2D *data2d = (TH2D*)(f1->Get(name.c_str())); - if (data2d == NULL) { - std::cerr << "no data" << std::endl; - error = true; - } - - TH2D *MC2d = (TH2D*)(f1->Get((bla+"_MC").c_str())); - if (MC2d == NULL) { - std::cerr << "no MC" << std::endl; - error = true; - } - TH2D *MC22d = (TH2D*)(f2->Get((bla+"_MC").c_str())); - if (MC22d == NULL) { - std::cerr << "no MC2" << std::endl; - error = true; - } - - TH2D *MC32d = (TH2D*)(f3->Get((bla+"_MC").c_str())); - if (MC32d == NULL) { - std::cerr << "no MC3" << std::endl; - error = true; - } - - TH2D *MC42d = (TH2D*)(f4->Get((bla+"_MC").c_str())); - if (MC42d == NULL) { - std::cerr << "no MC4" << std::endl; - error = true; - } - - if (error) continue; - - /* - data2d->SetMinimum(0); - for (int i = 0; i < MC2d->GetXaxis()->GetNbins()+1; i++) { - for (int j = 0; i < MC2d->GetYaxis()->GetNbins()+1; j++) { - MC2d->SetBinError(i+1, j+1, 0); - MC22d->SetBinError(i+1, j+1, 0); - MC32d->SetBinError(i+1, j+1, 0); - MC42d->SetBinError(i+1, j+1, 0); - } - } - */ - // Do the projections - int nbinsx = data2d->GetXaxis()->GetNbins(); - for (int i = 1; i < nbinsx+1; ++i) { - TH1D *datax = data2d->ProjectionY(Form("%s_%i,%i", data2d->GetName(), 0, i), i, i); - TH1D *mcx = MC2d ->ProjectionY(Form("%s_%i,%i", MC2d->GetName(), 1, i), i, i); - TH1D *mc2x = MC22d ->ProjectionY(Form("%s_%i,%i", MC22d->GetName(), 2, i), i, i); - TH1D *mc3x = MC32d ->ProjectionY(Form("%s_%i,%i", MC32d->GetName(), 3, i), i, i); - TH1D *mc4x = MC42d ->ProjectionY(Form("%s_%i,%i", MC42d->GetName(), 4, i), i, i); - - datax->Draw(); - mcx->Draw("same,hist"); - mc2x->Draw("same,hist"); - mc3x->Draw("same,hist"); - mc4x->Draw("same,hist"); - - double maximum = datax->GetMaximum(); - if (maximum < MC->GetMaximum()) maximum = mcx->GetMaximum(); - datax->GetYaxis()->SetRangeUser(0, maximum*1.7); - if (name == "MINERvA_CCinc_XSec_2DEavq3_nu_data") { - datax->GetXaxis()->SetRangeUser(0, 0.6); - datax->GetYaxis()->SetRangeUser(0, 6E-42); - } - - mcx->SetLineColor(kRed-9); - mcx->SetMarkerStyle(0); - mcx->SetLineStyle(kDashed); - - mc2x->SetLineColor(kRed+2); - mc2x->SetMarkerStyle(0); - - mc3x->SetLineColor(kBlue-9); - mc3x->SetMarkerStyle(0); - mc3x->SetLineStyle(kDashed); - - mc4x->SetLineColor(kBlue+2); - mc4x->SetMarkerStyle(0); - - TLegend *leg = new TLegend(0.6, 0.5, 0.95, 0.90); - int nbins = MC->GetXaxis()->GetNbins()*MC->GetYaxis()->GetNbins(); - double chi2ndof = std::atof(MC->GetTitle())/nbins; - double chi2ndof2 = std::atof(MC2->GetTitle())/nbins; - double chi2ndof3 = std::atof(MC3->GetTitle())/nbins; - double chi2ndof4 = std::atof(MC4->GetTitle())/nbins; - leg->AddEntry(datax, Form("Data %.2f-%.2f", data2d->GetXaxis()->GetBinLowEdge(i), data2d->GetXaxis()->GetBinLowEdge(i+1)), "le"); - leg->AddEntry(MC, Form("MINERvA untuned, #chi^{2}/ndof=%.2f", chi2ndof), "l"); - leg->AddEntry(MC2, Form("MINERvA tuned, #chi^{2}/ndof=%.2f", chi2ndof2), "l"); - leg->AddEntry(MC3, Form("NOvA untuned, #chi^{2}/ndof=%.2f", chi2ndof3), "l"); - leg->AddEntry(MC4, Form("NOvA tuned, #chi^{2}/ndof=%.2f", chi2ndof4), "l"); - leg->SetBorderSize(0); - leg->SetFillStyle(0); - leg->SetFillColor(0); - leg->Draw("same"); - - canv->Print((file1+"_"+file2+"_comp.pdf").c_str()); - delete leg; - } - - - } - - - } - - } - canv->Print((file1+"_"+file2+"_comp.pdf]").c_str()); -} diff --git a/README b/README index d41546c..356149b 100644 --- a/README +++ b/README @@ -1,149 +1,147 @@ # NUISANCE special release for MINERvA-NOvA workshop(s) [Indico page, 29 Sep 2018](https://minerva-docdb.fnal.gov/cgi-bin/private/DisplayMeeting?conferenceid=6457) ---- ## Description NUISANCE release used for MINERvA-NOvA workshop on 29 Sep 2018. Used to compare the NOvA tuned GENIE 2.12.10 interaction simulation to MINERvA's tuned GENIE 2.8.6/4. Talk is available on [FNAL DocDB](https://minerva-docdb.fnal.gov/cgi-bin/private/ShowDocument?docid=20487) and [here](http://www.hep.ph.ic.ac.uk/~cvw09/files/MINERvA_NOvA_WS_forNOVA.pdf). Requires GENIE 2.8.4 and above and the generated GENIE files. Please read the known issues! ---- ## Instructions Follow the usual instructions to install NUISANCE. In a nutshell: 1. Install GENIE to some directory, export that directory as `${GENIE}` 2. Source the `environment_setup.sh` file for GENIE, setting up all its dependencies 3. Get NUISANCE, get this branch, set the `${NUISANCE}` environment variable to the nuisance folder 4. Source `${NUISANCE}/myenv.sh` which sets up all your environment variables (e.g. ROOT) and generator dependencies (e.g. GENIE) — an example file is provided 5. Source `${NUISANCE}/freshbuild.sh`, making sure it contains `-DUSE_MINERvA_RW=1` (for parts of the MINERvA tune), and `-DUSE_GENIE=1` (for GENIE event libraries and linking), with any other generators you want included 6. NUISANCE should now be built in `${NUISANCE}/build` with executables in `${NUISANCE}/build/app` 7. `${NUISANCE}/build/app/PrepareGENIE` prepares GENIE output files for NUISANCE use. It converts the GHep format into something more manageable and produces cross-sections per target and per process. Its syntax is found by `./PrepareGENIE --help` ---- ## Generated GENIE files The GENIE files used are located [on the box](https://imperialcollegelondon.box.com/v/MINERvA-NOvA-WS).The file naming is `EXPERIMENT_FluxAndTarget_WithTUNEModel*.root`, where `EXPERIMENT` is the experiment whose flux and target generated the events, and `TUNE` is the NOvA or MINERvA tune. The `WithNOvAModel` files should have an additional `TTree` called `nova_wgts`, which are the weights provided by Jeremy Wolcott (Tufts). See below for the meaning of these weights. ---- ## Running flat-tree converters, external data comparisons, and producing plots The folder `${NUISANCE}/MINERvA_NOvA_WS` contains `xml` cardfiles for NUISANCE to produce external data comparison plots. It also contains a file `minerva_tune.xml` which is used to apply the MINERvA tune. The NOvA tune is applied by the code looking for the `nova_wgts` branch in the GENIE file. If it can't find it it moves on, if it does it applies the total weight from that file as a reweight. All the below assumes your raw GENIE file has been converted using `PrepareGENIE`, outlined in step 7 above. ### To run the flat-tree converter: Do `./nuisflat -i "GENIE:${GENIE_PREPARED_FILE}" -o "${GENIE_PREPARED_FILE%%.root}_FlatTree_Tune.root"` Add `-c minerva_tune.xml` if you want to apply the MINERvA tune. The flat-tree variables are intuitively named, and their meaning can be found in `${NUISANCE}/src/MCStudies/GenericFlux_Vectors.cxx`. To apply the tunings use the tree entries `Weight` for the MINERvA tune and `CustomWeight` for the NOvA tune. The `CustomWeightArray[6]` contains the different parts of the NOvA tune and are ordered as ``` fNUISANCEEvent->CustomWeight = NOVAw; fNUISANCEEvent->CustomWeightArray[0] = MAQEw; fNUISANCEEvent->CustomWeightArray[1] = NonResw; fNUISANCEEvent->CustomWeightArray[2] = RPAQEw; fNUISANCEEvent->CustomWeightArray[3] = RPARESw; fNUISANCEEvent->CustomWeightArray[4] = MECw; fNUISANCEEvent->CustomWeightArray[5] = NOVAw; ``` For the MINERvA tune the `CustomWeight` and `CustomWeightArray` are all 1.0. These are all set in `${NUISANCE}/src/InputHandler/GENIEInputHandler.cxx`. To scale from the event rates to cross sections apply `fScaleFactor` to your distributions. All entries have their own `fScaleFactor` but for GENIE they are all the same, so doing `double scale = tree->GetMinimum("fScaleFactor")` and applying that to a given histogram provides the scaling. ### To compare against external data: Do `./nuiscomp -c cardfile.xml -o GENIE_FILE.root` where the `cardfile.xml` provides NUISANCE with what experiments you want to run and where the generated events live. The example cardfile `${NUISANCE}/MINERvA_NOvA_WS/exp_list.xml` provides the general structure of the `xml` cardfiles. To find all supported experiments run `${NUISANCE}/scripts/nuissamples` ### Producing plots: #### Flat-tree An example flat-tree analysis ROOT scripts is provided in `${NUISANCE}/MINERvA_NOvA_WS/FlatAnalysis.cc`. It takes two arguments: ``` root -b -q -l 'FlatAnalysis.cc("MINERvA_FlatTree.root", "NOvA_FlatTree.root")' ``` where `MINERvA_FlatTree.root` is the flat-tree produced with the MINERvA tuning. `NOvA_FlatTree.root` is the equivalent for the NOvA tuning. #### Data comparisons `${NUISANCE}/MINERvA_NOvA_WS/plotscript.cc` takes four arguments: ``` root -b -q -l 'plotscript.cc("MINERvA_untuned.root", "MINERvA_tuned.root", "NOvA_untuned.root", "NOvA_tuned.root")' ``` which are the output files from `nuiscomp` for the MINERvA untuned, tuned, NOvA untuned and tuned generator files. These are far from perfect and exclude some known problems with some data distributions. So please contact us below if you need a hand or have requests! ---- ## Known issues The cross section scaling variable `fScaleFactor` does not play well with the NOvA target. Will be fixed. Some of the data distributions have broken in moving beyond v2r7. E.g. MiniBooNE CCQE Q2 seems to have issues under some ROOT versions (e.g. 5.34/36). Will be fixed. The `${NUISANCE}/data/reweight` directory does not exist. This is because the RikRPA files are very heavy (beyond git recommendations), so have to be added manually. Email us below to sort this out. ---- ## Contact Join our [Slack channel](https://nuisance-xsec.slack.com) Visit our [website](https://nuisance.hepforge.org) Send an email to [us](mailto:nuisance@projects.hepforge.org), [Clarence](mailto:cwret@fnal.gov), [Luke](mailto:picker24@fnal.gov), and [Kevin](mailto:kevin@rochester.edu) ---- ## Old readme Author Luke Pickering, Patrick Stowell, Callum Wilkinson, Clarence Wret Package Manager: (p.stowell@sheffield.ac.uk) NUISANCE Package v2r7 README 20/07/2016 --- Compilation The following instructions should be used to build the fitter after checking out from scratch. 1. Make sure environmental variables required for the generators you wish to build against are set. 2. In the top nuisance directory make a new build directory: "$mkdir build && cd build" 3. Run CMAKE with compiler flags for which generators are required "$ cmake -DUSE_NuWro=1 -DUSE_NEUT=1 -DUSE_GENIE=0 ../" 4. Build the fitter "$ make -j" 5. Build documentation "$ make docs" 6. Install to build location "$ make install" 7. Source the setup script "$ source Linux/setup.sh" If you prefer, most configure variables can be entered through a cmake UI, such as ccmake. e.g. "$ mkdir build && cd build && ccmake ../" --- Adding Classes The fitter is designed to be easily extended by adding new measurement classes whilst keeping the input convertors and tuning functionality the same. The Devel module folder is setup with some examples of how to add new classes into the framework. Feel free to email me if there are difficulties adding new measurements. --- Running Fits Whilst running fits is relatively quick and simple, there are now a large range of possible options. Doxygen Documentation is being added to the $NUISANCE/doc/html folder. Refer there for guidance on how to properly formulate a card file. - - diff --git a/README.md b/README.md deleted file mode 100644 index 5326413..0000000 --- a/README.md +++ /dev/null @@ -1,102 +0,0 @@ -# NUISANCE special release for MINERvA-NOvA workshop(s) -[Indico page, 29 Sep 2018](https://minerva-docdb.fnal.gov/cgi-bin/private/DisplayMeeting?conferenceid=6457) - ----- -## Description -NUISANCE release used for MINERvA-NOvA workshop on 29 Sep 2018. Used to compare the NOvA tuned GENIE 2.12.10 interaction simulation to MINERvA's tuned GENIE 2.8.6/4. - -Talk is available on [FNAL DocDB](https://minerva-docdb.fnal.gov/cgi-bin/private/ShowDocument?docid=20487) and [here](http://www.hep.ph.ic.ac.uk/~cvw09/files/MINERvA_NOvA_WS_forNOVA.pdf). - -Requires GENIE 2.8.4 and above and the generated GENIE files. - -Please read the known issues! - ----- -## Instructions -Follow the usual instructions to install NUISANCE. - -In a nutshell: - -1. Install GENIE to some directory, export that directory as ``${GENIE}`` -2. Source the ``environment_setup.sh`` file for GENIE, setting up all its dependencies -3. Get NUISANCE, get this branch, set the ``${NUISANCE}`` environment variable to the nuisance folder -4. Source ``${NUISANCE}/myenv.sh`` which sets up all your environment variables (e.g. ROOT) and generator dependencies (e.g. GENIE) — an example file is provided -5. Source ``${NUISANCE}/freshbuild.sh``, making sure it contains ``-DUSE_MINERvA_RW=1`` (for parts of the MINERvA tune), and ``-DUSE_GENIE=1`` (for GENIE event libraries and linking), with any other generators you want included -6. NUISANCE should now be built in ``${NUISANCE}/build`` with executables in ``${NUISANCE}/build/app`` -7. ``${NUISANCE}/build/app/PrepareGENIE`` prepares GENIE output files for NUISANCE use. It converts the GHep format into something more manageable and produces cross-sections per target and per process. Its syntax is found by ``./PrepareGENIE --help`` - ----- -## Running flat-tree converters, external data comparisons, and producing plots -The folder ``${NUISANCE}/MINERvA_NOvA_WS`` contains ``xml`` cardfiles for NUISANCE to produce external data comparison plots. It also contains a file ``minerva_tune.xml`` which is used to apply the MINERvA tune. - -The NOvA tune is applied by the code looking for the ``nova_wgts`` branch in the GENIE file. If it can't find it it moves on, if it does it applies the total weight from that file as a reweight. - -All the below assumes your raw GENIE file has been converted using ``PrepareGENIE``, outlined in step 7 above. - -### To run the flat-tree converter: -Do ``./nuisflat -i "GENIE:${GENIE_PREPARED_FILE}" -o "${GENIE_PREPARED_FILE%%.root}_FlatTree_Tune.root" -`` Add ``-c minerva_tune.xml`` if you want to apply the MINERvA tune. - -The flat-tree variables are intuitively named, and their meaning can be found in ``${NUISANCE}/src/MCStudies/GenericFlux_Vectors.cxx``. - -To apply the tunings use the tree entries ``Weight`` for the MINERvA tune and ``CustomWeight`` for the NOvA tune. The ``CustomWeightArray[6]`` contains the different parts of the NOvA tune and are ordered as - -``` -fNUISANCEEvent->CustomWeight = NOVAw; -fNUISANCEEvent->CustomWeightArray[0] = MAQEw; -fNUISANCEEvent->CustomWeightArray[1] = NonResw; -fNUISANCEEvent->CustomWeightArray[2] = RPAQEw; -fNUISANCEEvent->CustomWeightArray[3] = RPARESw; -fNUISANCEEvent->CustomWeightArray[4] = MECw; -fNUISANCEEvent->CustomWeightArray[5] = NOVAw; -``` -For the MINERvA tune the `CustomWeight` and `CustomWeightArray` are all 1.0. - -To scale from the event rates to cross sections apply `fScaleFactor` to your distributions. All entries have their own ``fScaleFactor`` but for GENIE they are all the same, so doing ``double scale = tree->GetMinimum("fScaleFactor")`` and applying that to a given histogram provides the scaling. - - -### To compare against external data: -Do ``./nuiscomp -c cardfile.xml -o GENIE_FILE.root`` where the `cardfile.xml` provides NUISANCE with what experiments you want to run and where the generated events live. - -The example cardfile ``${NUISANCE}/MINERvA_NOvA_WS/exp_list.xml`` provides the general structure of the `xml` cardfiles. - -To find all supported experiments run ``${NUISANCE}/scripts/nuissamples`` - -### Producing plots: -#### Flat-tree -An example flat-tree analysis ROOT scripts is provided in `${NUISANCE}/MINERvA_NOvA_WS/FlatAnalysis.cc`. It takes two arguments: - -``` -root -b -q -l 'FlatAnalysis.cc("MINERvA_FlatTree.root", "NOvA_FlatTree.root")' -``` - -where `MINERvA_FlatTree.root` is the flat-tree produced with the MINERvA tuning. `NOvA_FlatTree.root` is the equivalent for the NOvA tuning. - -#### Data comparisons -``${NUISANCE}/MINERvA_NOvA_WS/plotscript.cc`` takes four arguments: - -``` -root -b -q -l 'plotscript.cc("MINERvA_untuned.root", "MINERvA_tuned.root", "NOvA_untuned.root", "NOvA_tuned.root")' -``` - -which are the output files from `nuiscomp` for the MINERvA untuned, tuned, NOvA untuned and tuned generator files. - -These are far from perfect and exclude some known problems with some data distributions. So please contact us below if you need a hand or have requests! - - ----- -## Known issues -The cross section scaling variable `fScaleFactor` does not play well with the NOvA target. Will be fixed. - -Some of the data distributions have broken in moving beyond v2r7. E.g. MiniBooNE CCQE Q2 seems to have issues under some ROOT versions (e.g. 5.34/36). Will be fixed. - -The ``${NUISANCE}/data/reweight`` directory does not exist. This is because the RikRPA files are very heavy (beyond git recommendations), so have to be added manually. Email us below to sort this out. - ----- -## Contact -Join our [Slack channel](nuisance-xsec.slack.com) - -Visit our [website](nuisance.hepforge.org) - -Send an email to [us](mailto:nuisance@projects.hepforge.org), [Clarence](mailto:cwret@fnal.gov), [Luke](mailto:picker24@fnal.gov), and [Kevin](mailto:kevin@rochester.edu) - diff --git a/app/PrepareGENIE.cxx b/app/PrepareGENIE.cxx index d7ac3d2..9521eaf 100644 --- a/app/PrepareGENIE.cxx +++ b/app/PrepareGENIE.cxx @@ -1,799 +1,860 @@ #include #include #include "FitLogger.h" #include "PlotUtils.h" #include "TFile.h" #include "TH1D.h" #include "TTree.h" #ifdef __GENIE_ENABLED__ #include "Conventions/Units.h" #include "GHEP/GHepParticle.h" #include "PDG/PDGUtils.h" #endif std::string gInputFiles = ""; std::string gOutputFile = ""; std::string gFluxFile = ""; std::string gTarget = ""; double MonoEnergy; int gNEvents = -999; bool IsMonoE = false; void PrintOptions(); void ParseOptions(int argc, char* argv[]); void RunGENIEPrepareMono(std::string input, std::string target, std::string output); void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output); int main(int argc, char* argv[]) { ParseOptions(argc, argv); if (IsMonoE) { RunGENIEPrepareMono(gInputFiles, gTarget, gOutputFile); } else { RunGENIEPrepare(gInputFiles, gFluxFile, gTarget, gOutputFile); } } void RunGENIEPrepareMono(std::string input, std::string target, std::string output) { LOG(FIT) << "Running GENIE Prepare in mono energetic with E = " << MonoEnergy << " GeV" << std::endl; // Setup TTree TChain* tn = new TChain("gtree"); tn->AddFile(input.c_str()); if (tn->GetFile() == NULL) { tn->Print(); ERROR(FTL, "gtree not located in GENIE file: " << input); THROW("Check your inputs, they may need to be completely regenerated!"); throw; } int nevt = tn->GetEntries(); if (gNEvents != -999) { LOG(FIT) << "Overriding number of events by user from " << nevt << " to " << gNEvents << std::endl; nevt = gNEvents; } NtpMCEventRecord* genientpl = NULL; tn->SetBranchAddress("gmcrec", &genientpl); // Have the TH1D go from MonoEnergy/2 to MonoEnergy/2 TH1D* fluxhist = new TH1D("flux", "flux", 1000, MonoEnergy/2., MonoEnergy*2.); fluxhist->Fill(MonoEnergy); fluxhist->Scale(1, "width"); // Make Event Hist TH1D* eventhist = (TH1D*)fluxhist->Clone(); eventhist->Reset(); TH1D* xsechist = (TH1D*)eventhist->Clone(); // Create maps std::map modexsec; std::map modecount; std::vector genieids; std::vector targetids; std::vector interids; // Loop over all events for (int i = 0; i < nevt; i++) { tn->GetEntry(i); StopTalking(); EventRecord& event = *(genientpl->event); GHepParticle* neu = event.Probe(); StartTalking(); // Get XSec From Spline GHepRecord genie_record = static_cast(event); double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); // Parse Interaction String std::string mode = genie_record.Summary()->AsString(); std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); std::string targ = (modevec[0] + ";" + modevec[1]); std::string inter = mode; // Fill lists of Unique IDS if (std::find(targetids.begin(), targetids.end(), targ) == targetids.end()) { targetids.push_back(targ); } if (std::find(interids.begin(), interids.end(), inter) == interids.end()) { interids.push_back(inter); } // Create entries Mode Maps if (modexsec.find(mode) == modexsec.end()) { genieids.push_back(mode); modexsec[mode] = (TH1D*)xsechist->Clone(); modecount[mode] = (TH1D*)xsechist->Clone(); modexsec[mode]->GetYaxis()->SetTitle("d#sigma/dE_{#nu} #times 10^{-38} (events weighted by #sigma)"); modecount[mode]->GetYaxis()->SetTitle("Number of events in file"); } // Fill XSec Histograms modexsec[mode]->Fill(neu->E(), xsec); modecount[mode]->Fill(neu->E()); // Fill total event hist eventhist->Fill(neu->E()); // Clear Event genientpl->Clear(); size_t freq = nevt / 20; if (freq && !(i % freq)) { LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events (E: " << neu->E() << " GeV, xsec: " << xsec << " E-38 cm^2/nucleon)" << std::endl; } } LOG(FIT) << "Processed all events" << std::endl; TFile* outputfile; // If no output is specified just append to the file if (!gOutputFile.length()) { tn->GetEntry(0); outputfile = tn->GetFile(); outputfile->cd(); } else { outputfile = new TFile(gOutputFile.c_str(), "RECREATE"); outputfile->cd(); QLOG(FIT, "Cloning input vector to output file: " << gOutputFile); TTree* cloneTree = tn->CloneTree(); cloneTree->SetDirectory(outputfile); cloneTree->Write(); QLOG(FIT, "Cloning input nova_wgts to output file: " << gOutputFile); // *********************************** // *********************************** // FUDGE FOR NOVA MINERVA WORKSHOP // Also check for the nova_wgts tree from Jeremy TChain *nova_chain = new TChain("nova_wgts"); nova_chain->AddFile(input.c_str()); TTree* nova_tree = nova_chain->GetTree(); if (!nova_tree) { QLOG(FIT, "Could not find nova_wgts tree in " << gOutputFile); } else { QLOG(FIT, "Found nova_wgts tree in " << gOutputFile); } if (nova_tree) { nova_tree->SetDirectory(outputfile); nova_tree->Write(); } QLOG(FIT, "Done cloning tree."); } LOG(FIT) << "Getting splines in mono-energetic..." << std::endl; // Save each of the reconstructed splines to file std::map modeavg; TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines"); if (!inddir) inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines"); inddir->cd(); // Loop over GENIE ID's and get MEC count int MECcount = 0; bool MECcorrect = FitPar::Config().GetParB("CorrectGENIEMECNorm"); for (UInt_t i = 0; i < genieids.size(); i++) { if (genieids[i].find("MEC") != std::string::npos) { MECcount++; } } LOG(FIT) << "Found " << MECcount << " repeated MEC instances." << std::endl; for (UInt_t i = 0; i < genieids.size(); i++) { std::string mode = genieids[i]; modexsec[mode]->Write((mode + "_summed_xsec").c_str(), TObject::kOverwrite); modecount[mode]->Write((mode + "_summed_evt").c_str(), TObject::kOverwrite); // Form extra avg xsec map -> Reconstructed spline modeavg[mode] = (TH1D*)modexsec[mode]->Clone(); modeavg[mode]->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)"); modeavg[mode]->Divide(modecount[mode]); if (MECcorrect && (mode.find("MEC") != std::string::npos)) { modeavg[mode]->Scale(1.0 / double(MECcount)); } modeavg[mode]->Write((mode + "_rec_spline").c_str(), TObject::kOverwrite); } TDirectory* targdir = (TDirectory*)outputfile->Get("TargetGENIESplines"); if (!targdir) targdir = (TDirectory*)outputfile->mkdir("TargetGENIESplines"); targdir->cd(); LOG(FIT) << "Getting Target Splines" << std::endl; // For each target save a total spline std::map targetsplines; for (uint i = 0; i < targetids.size(); i++) { std::string targ = targetids[i]; LOG(FIT) << "Getting target " << i << ": " << targ << std::endl; targetsplines[targ] = (TH1D*)xsechist->Clone(); targetsplines[targ]->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)"); LOG(FIT) << "Created target spline for " << targ << std::endl; for (uint j = 0; j < genieids.size(); j++) { std::string mode = genieids[j]; if (mode.find(targ) != std::string::npos) { LOG(FIT) << " Mode " << mode << " contains " << targ << " target" << std::endl; targetsplines[targ]->Add(modeavg[mode]); LOG(FIT) << "Finished with Mode " << mode << " " << modeavg[mode]->Integral() << std::endl; } } LOG(FIT) << "Saving target spline:" << targ << std::endl; targetsplines[targ]->Write(("Total_" + targ).c_str(), TObject::kOverwrite); } LOG(FIT) << "Getting total splines" << std::endl; // Now we have each of the targets we need to create a total cross-section. int totalnucl = 0; // Get the targets specified by the user, separated by commas + // This has structure target1[fraction1], target2[fraction2] std::vector targprs = GeneralUtils::ParseToStr(target, ","); + + std::vector targ_list; + std::vector frac_list; + + // Chop up the target string which has format TARGET1[fraction1],TARGET2[fraction2] + + //std::cout << "Targets: " << std::endl; + // Loop over the vector of strings "TARGET1[fraction1]" "TARGET2[fraction2]" + for (std::vector::iterator it = targprs.begin(); it != targprs.end(); ++it) { + // Cut into "TARGET1" and "fraction1]" + std::vector targind = GeneralUtils::ParseToStr(*it, "["); + //std::cout << " " << *it << std::endl; + // Cut into "TARGET1" and "fraction1" + for (std::vector::iterator jt = targind.begin(); jt != targind.end(); ++jt) { + if ((*jt).find("]") != std::string::npos) { + (*jt) = (*jt).substr(0, (*jt).find("]")); + //*jt = "hello"; + frac_list.push_back(*jt); + // Won't find bracket for target + } else { + targ_list.push_back(*jt); + } + } + } + + targprs = targ_list; + + std::vector targ_fractions; + double minimum = 1.0; + for (std::vector::iterator it = frac_list.begin(); it != frac_list.end(); it++) { + //std::cout << " " << *it << std::endl; + double frac = std::atof((*it).c_str()); + targ_fractions.push_back(frac); + if (frac < minimum) minimum = frac; + } + + std::vector::iterator it = targ_fractions.begin(); + std::vector::iterator jt = targ_list.begin(); + double scaling = 0; + for (; it != targ_fractions.end(); it++, jt++) { + // First get the mass number from the targ_list + int nucl = atoi((*jt).c_str()); + nucl = (nucl%10000)/10; + // Gets the relative portions right + *it = (*it)/minimum; + // Scale relative the atomic mass + //(*it) *= (double(nucl)/(*it)); + double tempscaling = double(nucl)/(*it); + if (tempscaling > scaling) scaling=tempscaling; + } + it = targ_fractions.begin(); + for (; it != targ_fractions.end(); it++) { + // Round the scaling to nearest integer and multiply + *it *= int(scaling+0.5); + // Round to nearest integer + *it = int(*it+0.5); + totalnucl += *it; + } + + if (totalnucl == 0) { + THROW("Didn't find any nucleons in input file. Did you really specify the target ratios?\ne.g. TARGET1[fraction1],TARGET2[fraction2]" << std::endl); + } TH1D* totalxsec = (TH1D*)xsechist->Clone(); for (uint i = 0; i < targprs.size(); i++) { std::string targpdg = targprs[i]; // Check that we found the user requested target in GENIE bool FoundTarget = false; - for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) { std::string targstr = iter->first; TH1D* xsec = iter->second; + // Match the user targets to the targets found in GENIE if (targstr.find(targpdg) != std::string::npos) { FoundTarget = true; LOG(FIT) << "Adding target spline " << targstr << " Integral = " << xsec->Integral("width") << std::endl; totalxsec->Add(xsec); - int nucl = atoi(targpdg.c_str()); - totalnucl += int((nucl % 10000) / 10); + //int nucl = atoi(targpdg.c_str()); + //totalnucl += int((nucl % 10000) / 10); } } // Check that targets were all found if (!FoundTarget) { ERR(WRN) << "Didn't find target " << targpdg << " in the list of targets recorded by GENIE" << std::endl; ERR(WRN) << " The list of targets you requested is: " << std::endl; for (uint i = 0; i < targprs.size(); ++i) ERR(WRN) << " " << targprs[i] << std::endl; ERR(WRN) << " The list of targets found in GENIE is: " << std::endl; for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) ERR(WRN) << " " << iter->first<< std::endl; } } outputfile->cd(); totalxsec->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/nucleon)"); totalxsec->Write("nuisance_xsec", TObject::kOverwrite); eventhist = (TH1D*)fluxhist->Clone(); eventhist->Multiply(totalxsec); eventhist->GetYaxis()->SetTitle((std::string("Event rate (N = #sigma #times #Phi) #times 10^{-38} (cm^{2}/nucleon) #times ")+eventhist->GetYaxis()->GetTitle()).c_str()); LOG(FIT) << "Dividing by Total Nucl = " << totalnucl << std::endl; eventhist->Scale(1.0 / double(totalnucl)); eventhist->Write("nuisance_events", TObject::kOverwrite); fluxhist->Write("nuisance_flux", TObject::kOverwrite); LOG(FIT) << "Inclusive XSec Per Nucleon = " << eventhist->Integral("width") * 1E-38 / fluxhist->Integral("width") << std::endl; LOG(FIT) << "XSec Hist Integral = " << totalxsec->Integral("width") << std::endl; outputfile->Close(); return; } void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output) { LOG(FIT) << "Running GENIE Prepare with flux..." << std::endl; // Get Flux Hist std::vector fluxvect = GeneralUtils::ParseToStr(flux, ","); TH1* fluxhist = NULL; if (fluxvect.size() == 3) { double from = GeneralUtils::StrToDbl(fluxvect[0]); double to = GeneralUtils::StrToDbl(fluxvect[1]); double step = GeneralUtils::StrToDbl(fluxvect[2]); int nstep = ceil((to - from) / step); to = from + step * nstep; QLOG(FIT, "Generating flat flux histogram from " << from << " to " << to << " with bins " << step << " wide (NBins = " << nstep << ")."); fluxhist = new TH1D("spectrum", ";E_{#nu} (GeV);Count (A.U.)", nstep, from, to); for (Int_t bi_it = 1; bi_it < fluxhist->GetXaxis()->GetNbins(); ++bi_it) { fluxhist->SetBinContent(bi_it, 1.0 / double(step * nstep)); } fluxhist->SetDirectory(0); } else if (fluxvect.size() == 2) { TFile* fluxfile = new TFile(fluxvect[0].c_str(), "READ"); if (!fluxfile->IsZombie()) { fluxhist = dynamic_cast(fluxfile->Get(fluxvect[1].c_str())); if (!fluxhist) { ERR(FTL) << "Couldn't find histogram named: \"" << fluxvect[1] << "\" in file: \"" << fluxvect[0] << std::endl; throw; } fluxhist->SetDirectory(0); } } else if (fluxvect.size() == 1) { MonoEnergy = GeneralUtils::StrToDbl(fluxvect[0]); RunGENIEPrepareMono(input, target, output); return; } else { LOG(FTL) << "Bad flux specification: \"" << flux << "\"." << std::endl; throw; } // Setup TTree TChain* tn = new TChain("gtree"); if (input.find_first_of(',') != std::string::npos) { std::vector inputvect = GeneralUtils::ParseToStr(input, ","); for (size_t iv_it = 0; iv_it < inputvect.size(); ++iv_it) { tn->AddFile(inputvect[iv_it].c_str()); QLOG(FIT, "Added input file: " << inputvect[iv_it]); } } else { // The Add form can accept wildcards. tn->Add(input.c_str()); } if (tn->GetFile() == NULL) { tn->Print(); ERROR(FTL, "gtree not located in GENIE file: " << input); THROW("Check your inputs, they may need to be completely regenerated!"); throw; } int nevt = tn->GetEntries(); if (gNEvents != -999) { LOG(FIT) << "Overriding number of events by user from " << nevt << " to " << gNEvents << std::endl; nevt = gNEvents; } if (!nevt) { THROW("Couldn't load any events from input specification: \"" << input.c_str() << "\""); } else { QLOG(FIT, "Found " << nevt << " input entries in " << input); } NtpMCEventRecord* genientpl = NULL; tn->SetBranchAddress("gmcrec", &genientpl); // Make Event and xsec Hist TH1D* eventhist = (TH1D*)fluxhist->Clone(); eventhist->Reset(); TH1D* xsechist = (TH1D*)eventhist->Clone(); // Create maps std::map modexsec; std::map modecount; std::vector genieids; std::vector targetids; std::vector interids; // Loop over all events for (int i = 0; i < nevt; i++) { tn->GetEntry(i); // Hussssch GENIE StopTalking(); // Get the event EventRecord& event = *(genientpl->event); // Get the neutrino GHepParticle* neu = event.Probe(); StartTalking(); // Get XSec From Spline // Get the GHepRecord GHepRecord genie_record = static_cast(event); double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); // Parse Interaction String std::string mode = genie_record.Summary()->AsString(); std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); std::string targ = (modevec[0] + ";" + modevec[1]); std::string inter = mode; // Get target nucleus // Alternative ways of getting the summaries //GHepParticle *target = genie_record.TargetNucleus(); //int pdg = target->Pdg(); // Fill lists of Unique IDS (neutrino and target) if (std::find(targetids.begin(), targetids.end(), targ) == targetids.end()) { targetids.push_back(targ); } // The full interaction list if (std::find(interids.begin(), interids.end(), inter) == interids.end()) { interids.push_back(inter); } // Create entries Mode Maps if (modexsec.find(mode) == modexsec.end()) { genieids.push_back(mode); modexsec[mode] = (TH1D*)xsechist->Clone(); modecount[mode] = (TH1D*)xsechist->Clone(); modexsec[mode]->GetYaxis()->SetTitle("d#sigma/dE_{#nu} #times 10^{-38} (events weighted by #sigma)"); modecount[mode]->GetYaxis()->SetTitle("Number of events in file"); } // Fill XSec Histograms modexsec[mode]->Fill(neu->E(), xsec); modecount[mode]->Fill(neu->E()); // Fill total event hist eventhist->Fill(neu->E()); if (i % (nevt / 20) == 0) { LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events (E: " << neu->E() << " GeV, xsec: " << xsec << " E-38 cm^2/nucleon)" << std::endl; } // Clear Event genientpl->Clear(); } LOG(FIT) << "Processed all events" << std::endl; // Once event loop is done we can start saving stuff into the file TFile* outputfile; if (!gOutputFile.length()) { tn->GetEntry(0); outputfile = tn->GetFile(); outputfile->cd(); } else { outputfile = new TFile(gOutputFile.c_str(), "RECREATE"); outputfile->cd(); QLOG(FIT, "Cloning input vector to output file: " << gOutputFile); TTree* cloneTree = tn->CloneTree(); cloneTree->SetDirectory(outputfile); cloneTree->Write(); // ******************************** // CLUDGE KLUDGE KLUDGE FOR NOVA QLOG(FIT, "Cloning input nova_wgts to output file: " << gOutputFile); // Also check for the nova_wgts tree from Jeremy TChain *nova_chain = new TChain("nova_wgts"); nova_chain->AddFile(input.c_str()); TTree* nova_tree = nova_chain->CloneTree(); if (!nova_tree) { QLOG(FIT, "Could not find nova_wgts tree in " << input); } else { QLOG(FIT, "Found nova_wgts tree in " << input); nova_tree->SetDirectory(outputfile); nova_tree->Write(); } QLOG(FIT, "Done cloning tree."); } LOG(FIT) << "Getting splines..." << std::endl; // Save each of the reconstructed splines to file std::map modeavg; TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines"); if (!inddir) inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines"); inddir->cd(); // Loop over GENIE ID's and get MEC count int MECcount = 0; bool MECcorrect = FitPar::Config().GetParB("CorrectGENIEMECNorm"); for (UInt_t i = 0; i < genieids.size(); i++) { if (genieids[i].find("MEC") != std::string::npos) { MECcount++; } } LOG(FIT) << "Found " << MECcount << " repeated MEC instances." << std::endl; for (UInt_t i = 0; i < genieids.size(); i++) { std::string mode = genieids[i]; modexsec[mode]->Write((mode + "_summed_xsec").c_str(), TObject::kOverwrite); modecount[mode]->Write((mode + "_summed_evt").c_str(), TObject::kOverwrite); // Form extra avg xsec map -> Reconstructed spline modeavg[mode] = (TH1D*)modexsec[mode]->Clone(); modeavg[mode]->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)"); modeavg[mode]->Divide(modecount[mode]); if (MECcorrect && (mode.find("MEC") != std::string::npos)) { modeavg[mode]->Scale(1.0 / double(MECcount)); } modeavg[mode]->Write((mode + "_rec_spline").c_str(), TObject::kOverwrite); } TDirectory* targdir = (TDirectory*)outputfile->Get("TargetGENIESplines"); if (!targdir) targdir = (TDirectory*)outputfile->mkdir("TargetGENIESplines"); targdir->cd(); LOG(FIT) << "Getting Target Splines" << std::endl; // For each target save a total spline std::map targetsplines; for (uint i = 0; i < targetids.size(); i++) { std::string targ = targetids[i]; LOG(FIT) << "Getting target " << i << ": " << targ << std::endl; targetsplines[targ] = (TH1D*)xsechist->Clone(); targetsplines[targ]->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)"); LOG(FIT) << "Created target spline for " << targ << std::endl; for (uint j = 0; j < genieids.size(); j++) { std::string mode = genieids[j]; // Look at all matching modes/targets if (mode.find(targ) != std::string::npos) { LOG(FIT) << " Mode " << mode << " contains " << targ << " target" << std::endl; targetsplines[targ]->Add(modeavg[mode]); LOG(FIT) << "Finished with Mode " << mode << " " << modeavg[mode]->Integral() << std::endl; } } LOG(FIT) << "Saving target spline: " << targ << std::endl; targetsplines[targ]->Write(("Total_" + targ).c_str(), TObject::kOverwrite); } LOG(FIT) << "Getting total splines" << std::endl; // Now we have each of the targets we need to create a total cross-section. int totalnucl = 0; // This has structure target1[fraction1], target2[fraction2] std::vector targprs = GeneralUtils::ParseToStr(target, ","); std::vector targ_list; std::vector frac_list; // Chop up the target string which has format TARGET1[fraction1],TARGET2[fraction2] //std::cout << "Targets: " << std::endl; // Loop over the vector of strings "TARGET1[fraction1]" "TARGET2[fraction2]" for (std::vector::iterator it = targprs.begin(); it != targprs.end(); ++it) { // Cut into "TARGET1" and "fraction1]" std::vector targind = GeneralUtils::ParseToStr(*it, "["); //std::cout << " " << *it << std::endl; // Cut into "TARGET1" and "fraction1" for (std::vector::iterator jt = targind.begin(); jt != targind.end(); ++jt) { if ((*jt).find("]") != std::string::npos) { (*jt) = (*jt).substr(0, (*jt).find("]")); //*jt = "hello"; frac_list.push_back(*jt); // Won't find bracket for target } else { targ_list.push_back(*jt); } } } targprs = targ_list; std::vector targ_fractions; double minimum = 1.0; for (std::vector::iterator it = frac_list.begin(); it != frac_list.end(); it++) { //std::cout << " " << *it << std::endl; double frac = std::atof((*it).c_str()); targ_fractions.push_back(frac); if (frac < minimum) minimum = frac; } std::vector::iterator it = targ_fractions.begin(); std::vector::iterator jt = targ_list.begin(); double scaling = 0; for (; it != targ_fractions.end(); it++, jt++) { // First get the mass number from the targ_list int nucl = atoi((*jt).c_str()); nucl = (nucl%10000)/10; // Gets the relative portions right *it = (*it)/minimum; // Scale relative the atomic mass //(*it) *= (double(nucl)/(*it)); double tempscaling = double(nucl)/(*it); if (tempscaling > scaling) scaling=tempscaling; } it = targ_fractions.begin(); for (; it != targ_fractions.end(); it++) { // Round the scaling to nearest integer and multiply *it *= int(scaling+0.5); // Round to nearest integer *it = int(*it+0.5); totalnucl += *it; } if (totalnucl == 0) { THROW("Didn't find any nucleons in input file. Did you really specify the target ratios?\ne.g. TARGET1[fraction1],TARGET2[fraction2]" << std::endl); } - TH1D* totalxsec = (TH1D*)xsechist->Clone(); // Loop over the specified targets by the user for (uint i = 0; i < targprs.size(); i++) { std::string targpdg = targprs[i]; // Check that we found the user requested target in GENIE bool FoundTarget = false; - for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) { std::string targstr = iter->first; TH1D* xsec = iter->second; // Match the user targets to the targets found in GENIE if (targstr.find(targpdg) != std::string::npos) { FoundTarget = true; LOG(FIT) << "Adding target spline " << targstr << " Integral = " << xsec->Integral("width") << std::endl; totalxsec->Add(xsec); //int nucl = atoi(targpdg.c_str()); //totalnucl += int((nucl % 10000) / 10); } } // Looped over target splines // Check that targets were all found if (!FoundTarget) { ERR(WRN) << "Didn't find target " << targpdg << " in the list of targets recorded by GENIE" << std::endl; ERR(WRN) << " The list of targets you requested is: " << std::endl; for (uint i = 0; i < targprs.size(); ++i) ERR(WRN) << " " << targprs[i] << std::endl; ERR(WRN) << " The list of targets found in GENIE is: " << std::endl; for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) ERR(WRN) << " " << iter->first<< std::endl; } } outputfile->cd(); totalxsec->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/nucleon)"); totalxsec->Write("nuisance_xsec", TObject::kOverwrite); eventhist = (TH1D*)fluxhist->Clone(); eventhist->Multiply(totalxsec); eventhist->GetYaxis()->SetTitle((std::string("Event rate (N = #sigma #times #Phi) #times 10^{-38} (cm^{2}/nucleon) #times ")+eventhist->GetYaxis()->GetTitle()).c_str()); LOG(FIT) << "Dividing by Total Nucl = " << totalnucl << std::endl; eventhist->Scale(1.0 / double(totalnucl)); eventhist->Write("nuisance_events", TObject::kOverwrite); fluxhist->Write("nuisance_flux", TObject::kOverwrite); LOG(FIT) << "Inclusive XSec Per Nucleon = " << eventhist->Integral("width") * 1E-38 / fluxhist->Integral("width") << std::endl; LOG(FIT) << "XSec Hist Integral = " << totalxsec->Integral() << std::endl; outputfile->Close(); return; }; void PrintOptions() { std::cout << "PrepareGENIEEvents NUISANCE app. " << std::endl << "Takes GHep Outputs and prepares events for NUISANCE." << std::endl << std::endl << "PrepareGENIE [-h,-help,--h,--help] [-i " "inputfile1.root,inputfile2.root,inputfile3.root,...] " << "[-f flux_root_file.root,flux_hist_name] [-t " "target1[frac1],target2[frac2],...]" << "[-n number_of_events (experimental)]" << std::endl << std::endl; std::cout << "Prepare Mode [Default] : Takes a single GHep file, " "reconstructs the original GENIE splines, " << " and creates a duplicate file that also contains the flux, " "event rate, and xsec predictions that NUISANCE needs. " << std::endl; std::cout << "Following options are required for Prepare Mode:" << std::endl; std::cout << " [ -i inputfile.root ] : Reads in a single GHep input file " "that needs the xsec calculation ran on it. " << std::endl; std::cout << " [ -f flux_file.root,hist_name ] : Path to root file " "containing the flux histogram the GHep records were generated " "with." << " A simple method is to point this to the flux histogram genie " "generatrs '-f /path/to/events/input-flux.root,spectrum'. " << std::endl; std::cout << " [ -f elow,ehigh,estep ] : Energy range specification when no " "flux file was used." << std::endl; std::cout << " [ -t target ] : Target that GHepRecords were generated with. " "Comma seperated list. E.g. for CH2 " "target=1000060120,1000010010,1000010010" << std::endl; std::cout << " [ -o outputfile.root ] : File to write prepared input file to." << std::endl; std::cout << " [ -m Mono_E_nu_GeV ] : Run in mono-energetic mode with m GeV neutrino energy." << std::endl; std::cout << " [ -n number_of_evt ] : Run with a reduced number of events for debugging purposes" << std::endl; } void ParseOptions(int argc, char* argv[]) { bool flagopt = false; // If No Arguments print commands for (int i = 1; i < argc; ++i) { if (!std::strcmp(argv[i], "-h")) { flagopt = true; break; } if (i + 1 != argc) { // Cardfile if (!std::strcmp(argv[i], "-h")) { flagopt = true; break; } else if (!std::strcmp(argv[i], "-i")) { gInputFiles = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-o")) { gOutputFile = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-f")) { gFluxFile = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-t")) { gTarget = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-n")) { gNEvents = GeneralUtils::StrToInt(argv[i + 1]); ++i; } else if (!std::strcmp(argv[i], "-m")) { MonoEnergy = GeneralUtils::StrToDbl(argv[i + 1]); IsMonoE = true; ++i; } else { ERR(FTL) << "ERROR: unknown command line option given! - '" << argv[i] << " " << argv[i + 1] << "'" << std::endl; PrintOptions(); break; } } } if (gInputFiles == "" && !flagopt) { ERR(FTL) << "No input file(s) specified!" << std::endl; flagopt = true; } if (gFluxFile == "" && !flagopt && !IsMonoE) { ERR(FTL) << "No flux input specified for Prepare Mode" << std::endl; flagopt = true; } if (gTarget == "" && !flagopt) { ERR(FTL) << "No target specified for Prepare Mode" << std::endl; flagopt = true; } if (gTarget.find("[") == std::string::npos || gTarget.find("]") == std::string::npos) { ERR(FTL) << "Didn't specify target ratios in Prepare Mode" << std::endl; ERR(FTL) << "Are you sure you gave it as -t \"TARGET1[fraction1],TARGET2[fraction]\"?" << std::endl; flagopt = true; } if (argc < 1 || flagopt) { PrintOptions(); exit(-1); } return; } diff --git a/cmake/GENIESetup.cmake b/cmake/GENIESetup.cmake index 23f0841..fe49732 100644 --- a/cmake/GENIESetup.cmake +++ b/cmake/GENIESetup.cmake @@ -1,181 +1,186 @@ # 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 . ################################################################################ # TODO # check system for libxml2 # check whether we need the includes # check if we can use a subset of the GENIE libraries ################################################################################ # Check Dependencies ################################################################################ ################################# GENIE ###################################### if(GENIE STREQUAL "") cmessage(FATAL_ERROR "Variable GENIE is not defined. " "The location of a pre-built GENIE install must be defined either as" " $ cmake -DGENIE=/path/to/GENIE or as and environment vairable" " $ export GENIE=/path/to/GENIE") endif() if (BUILD_GEVGEN) cmessage(STATUS "Building custom gevgen") LIST(APPEND EXTRA_CXX_FLAGS -D__GEVGEN_ENABLED__) endif() +if(GENIE_EMPMEC_REWEIGHT) + cmessage(STATUS "Enable EMPMEC dials") + LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_EMP_MECRW_ENABLED) +endif() + # Extract GENIE VERSION if (GENIE_VERSION STREQUAL "AUTO") execute_process (COMMAND ${CMAKE_SOURCE_DIR}/cmake/getgenieversion.sh ${GENIE} OUTPUT_VARIABLE GENIE_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE) endif() execute_process (COMMAND genie-config --libs OUTPUT_VARIABLE GENIE_LD_FLAGS_STR OUTPUT_STRIP_TRAILING_WHITESPACE) execute_process (COMMAND genie-config --topsrcdir OUTPUT_VARIABLE GENIE_INCLUDES_DIR OUTPUT_STRIP_TRAILING_WHITESPACE) string(REGEX MATCH "-L\([^ ]+\) \(.*\)$" PARSE_GENIE_LIBS_MATCH ${GENIE_LD_FLAGS_STR}) cmessage(DEBUG "genie-config --libs: ${GENIE_LD_FLAGS_STR}") if(NOT PARSE_GENIE_LIBS_MATCH) cmessage(FATAL_ERROR "Expected to be able to parse the result of genie-config --libs to a lib directory and a list of libraries to include, but got: \"${GENIE_LD_FLAGS_STR}\"") endif() set(GENIE_LIB_DIR ${CMAKE_MATCH_1}) set(GENIE_LIBS_RAW ${CMAKE_MATCH_2}) string(REPLACE "-l" "" GENIE_LIBS_STRIPED "${GENIE_LIBS_RAW}") cmessage(STATUS "GENIE version : ${GENIE_VERSION}") cmessage(STATUS "GENIE libdir : ${GENIE_LIB_DIR}") cmessage(STATUS "GENIE libs : ${GENIE_LIBS_STRIPED}") string(REGEX MATCH "ReinSeghal" WASMATCHED ${GENIE_LIBS_STRIPED}) if(WASMATCHED AND GENIE_VERSION STREQUAL "210") set(GENIE_SEHGAL ${GENIE_LIBS_STRIPED}) STRING(REPLACE "ReinSeghal" "ReinSehgal" GENIE_LIBS_STRIPED ${GENIE_SEHGAL}) cmessage(DEBUG "Fixed inconsistency in library naming: ${GENIE_LIBS_STRIPED}") endif() string(REGEX MATCH "ReWeight" WASMATCHED ${GENIE_LIBS_STRIPED}) if(NOT WASMATCHED) set(GENIE_LIBS_STRIPED "GReWeight ${GENIE_LIBS_STRIPED}") cmessage(DEBUG "Force added ReWeight library: ${GENIE_LIBS_STRIPED}") endif() -string(REPLACE " " ";" GENIE_LIBS_LIST "-Wl,--start-group ${GENIE_LIBS_STRIPED} -Wl,--end-group") +string(REPLACE " " ";" GENIE_LIBS_LIST "-Wl,--no-as-needed -Wl,--start-group ${GENIE_LIBS_STRIPED} -Wl,--end-group") cmessage(DEBUG "genie-config --libs -- MATCH1: ${CMAKE_MATCH_1}") cmessage(DEBUG "genie-config --libs -- MATCH2: ${CMAKE_MATCH_2}") cmessage(DEBUG "genie-config --libs -- libs stripped: ${GENIE_LIBS_STRIPED}") cmessage(DEBUG "genie-config --libs -- libs list: ${GENIE_LIBS_LIST}") ################################ LHAPDF ###################################### if(LHAPDF_LIB STREQUAL "") cmessage(FATAL_ERROR "Variable LHAPDF_LIB is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_LIB=/path/to/LHAPDF_libraries or as and environment vairable $ export LHAPDF_LIB=/path/to/LHAPDF_libraries") endif() if(LHAPDF_INC STREQUAL "") cmessage(FATAL_ERROR "Variable LHAPDF_INC is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_INC=/path/to/LHAPDF_includes or as and environment vairable $ export LHAPDF_INC=/path/to/LHAPDF_includes") endif() if(LHAPATH STREQUAL "") cmessage(FATAL_ERROR "Variable LHAPATH is not defined. The location of a the LHAPATH directory must be defined either as $ cmake -DLHAPATH=/path/to/LHAPATH or as and environment variable $ export LHAPATH=/path/to/LHAPATH") endif() ################################ LIBXML ###################################### if(LIBXML2_LIB STREQUAL "") # Check for xml2-config find_program(LIBXMLCFGLIBS xml2-config) if(NOT LIBXMLCFGLIBS STREQUAL "LIBXMLLIBS-NOTFOUND") execute_process (COMMAND ${LIBXMLCFGLIBS} --libs OUTPUT_VARIABLE LIBXML2_LIB OUTPUT_STRIP_TRAILING_WHITESPACE) else() message(FATAL_ERROR "Variable LIBXML2_LIB is not defined and could not find xml2-config. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_LIB=/path/to/LIBXML2_libraries or as and environment vairable $ export LIBXML2_LIB=/path/to/LIBXML2_libraries") endif() endif() if(LIBXML2_INC STREQUAL "") # Check for xml2-config find_program(LIBXMLCFGINCS xml2-config) if(NOT LIBXMLCFGINCS STREQUAL "LIBXMLINCS-NOTFOUND") execute_process (COMMAND ${LIBXMLCFGINCS} --cflags OUTPUT_VARIABLE LIBXML2_INC OUTPUT_STRIP_TRAILING_WHITESPACE) else() message(FATAL_ERROR "Variable LIBXML2_INC is not defined and could not find xml2-config. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_INC=/path/to/LIBXML2_libraries or as and environment vairable $ export LIBXML2_INC=/path/to/LIBXML2_libraries") endif() endif() ############################### log4cpp ###################################### if(LOG4CPP_LIB STREQUAL "") find_program(LOG4CPPCFG log4cpp-config) if(NOT LOG4CPPCFG STREQUAL "LOG4CPPCFG-NOTFOUND") execute_process (COMMAND ${LOG4CPPCFG} --pkglibdir OUTPUT_VARIABLE LOG4CPP_LIB OUTPUT_STRIP_TRAILING_WHITESPACE) else() message(FATAL_ERROR "Variable LOG4CPP_LIB is not defined. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_LIB=/path/to/LOG4CPP_libraries or as and environment vairable $ export LOG4CPP_LIB=/path/to/LOG4CPP_libraries") endif() endif() if(LOG4CPP_INC STREQUAL "") find_program(LOG4CPPCFG log4cpp-config) if(NOT LOG4CPPCFG STREQUAL "LOG4CPPCFG-NOTFOUND") execute_process (COMMAND ${LOG4CPPCFG} --pkgincludedir OUTPUT_VARIABLE LOG4CPP_INC OUTPUT_STRIP_TRAILING_WHITESPACE) else() message(FATAL_ERROR "Variable LOG4CPP_INC is not defined. The location of a pre-built log4cpp install must be defined either as $ cmake -DGENIE_LOG4CPP_INC=/path/to/LOG4CPP_includes or as and environment vairable $ export LOG4CPP_INC=/path/to/LOG4CPP_includes") endif() endif() ################################################################################ LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_ENABLED__ -D__GENIE_VERSION__=${GENIE_VERSION}) LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GENIE_INCLUDES_DIR} ${GENIE_INCLUDES_DIR}/GHEP ${GENIE_INCLUDES_DIR}/Ntuple ${GENIE_INCLUDES_DIR}/ReWeight ${GENIE_INCLUDES_DIR}/Apps ${GENIE_INCLUDES_DIR}/FluxDrivers ${GENIE_INCLUDES_DIR}/EVGDrivers ${LHAPDF_INC} ${LIBXML2_INC} ${LOG4CPP_INC}) SAYVARS() LIST(APPEND EXTRA_LINK_DIRS ${GENIE_LIB_DIR} ${LHAPDF_LIB} ${LIBXML2_LIB} ${LOG4CPP_LIB}) LIST(APPEND EXTRA_LIBS ${GENIE_LIBS_LIST}) LIST(APPEND EXTRA_LIBS LHAPDF xml2 log4cpp) if(USE_PYTHIA8) set(NEED_PYTHIA8 TRUE) set(NEED_ROOTPYTHIA8 TRUE) else() set(NEED_PYTHIA6 TRUE) set(NEED_ROOTPYTHIA6 TRUE) endif() set(NEED_ROOTEVEGEN TRUE) SET(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. " FORCE) diff --git a/cmake/c++CompilerSetup.cmake b/cmake/c++CompilerSetup.cmake index 518dc24..bc6b473 100644 --- a/cmake/c++CompilerSetup.cmake +++ b/cmake/c++CompilerSetup.cmake @@ -1,113 +1,121 @@ # 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 . ################################################################################ if(USE_OMP) LIST(APPEND EXTRA_CXX_FLAGS -fopenmp) endif() if(USE_DYNSAMPLES) LIST(APPEND EXTRA_LIBS dl) LIST(APPEND EXTRA_CXX_FLAGS -D__USE_DYNSAMPLES__) endif() set(CXX_WARNINGS -Wall ) cmessage(DEBUG "EXTRA_CXX_FLAGS: ${EXTRA_CXX_FLAGS}") string(REPLACE ";" " " STR_EXTRA_CXX_FLAGS "${EXTRA_CXX_FLAGS}") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${STR_EXTRA_CXX_FLAGS} ${CXX_WARNINGS}") set(CMAKE_Fortran_FLAGS_RELEASE "-fPIC") set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0") if(USE_DYNSAMPLES) set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fPIC") set(CMAKE_Fortran_FLAGS_DEBUG "-fPIC") endif() set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fPIC -O3") if(CMAKE_BUILD_TYPE MATCHES DEBUG) set(CURRENT_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS_DEBUG}) elseif(CMAKE_BUILD_TYPE MATCHES RELEASE) set(CURRENT_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS_RELEASE}) else() cmessage(FATAL_ERROR "[ERROR]: Unknown CMAKE_BUILD_TYPE (\"${CMAKE_BUILD_TYPE}\"): Should be \"DEBUG\" or \"RELEASE\".") endif() SET(STR_EXTRA_LINK_DIRS) if(NOT EXTRA_LINK_DIRS STREQUAL "") string(REPLACE ";" " -L" STR_EXTRA_LINK_DIRS "-L${EXTRA_LINK_DIRS}") endif() SET(STR_EXTRA_LIBS) if(NOT EXTRA_LIBS STREQUAL "") SET(STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS) string(REPLACE ";" " -l" STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS "-l${EXTRA_LIBS}") string(REPLACE "-l-" "-" STR_EXTRA_LIBS ${STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS}) endif() SET(STR_EXTRA_SHAREDOBJS) if(NOT EXTRA_SHAREDOBJS STREQUAL "") string(REPLACE ";" " " STR_EXTRA_SHAREDOBJS "${EXTRA_SHAREDOBJS}") endif() SET(STR_EXTRA_LINK_FLAGS) if(NOT EXTRA_LINK_FLAGS STREQUAL "") string(REPLACE ";" " " STR_EXTRA_LINK_FLAGS "${EXTRA_LINK_FLAGS}") endif() cmessage(DEBUG "EXTRA_LINK_DIRS: ${STR_EXTRA_LINK_DIRS}") cmessage(DEBUG "EXTRA_LIBS: ${STR_EXTRA_LIBS}") cmessage(DEBUG "EXTRA_SHAREDOBJS: ${STR_EXTRA_SHAREDOBJS}") cmessage(DEBUG "EXTRA_LINK_FLAGS: ${STR_EXTRA_LINK_FLAGS}") if(NOT STR_EXTRA_LINK_DIRS STREQUAL "" AND NOT STR_EXTRA_LIBS STREQUAL "") SET(CMAKE_DEPENDLIB_FLAGS "${STR_EXTRA_LINK_DIRS} ${STR_EXTRA_LIBS}") endif() if(NOT EXTRA_SHAREDOBJS STREQUAL "") if(NOT STR_EXTRA_LINK_FLAGS STREQUAL "") SET(STR_EXTRA_LINK_FLAGS "${STR_EXTRA_SHAREDOBJS} ${STR_EXTRA_LINK_FLAGS}") else() SET(STR_EXTRA_LINK_FLAGS "${STR_EXTRA_SHAREDOBJS}") endif() endif() if(NOT EXTRA_LINK_FLAGS STREQUAL "") if(NOT CMAKE_LINK_FLAGS STREQUAL "") SET(CMAKE_LINK_FLAGS "${CMAKE_LINK_FLAGS} ${STR_EXTRA_LINK_FLAGS}") else() SET(CMAKE_LINK_FLAGS "${STR_EXTRA_LINK_FLAGS}") endif() endif() if(USE_OMP) cmessage(FATAL_ERROR "No OMP features currently enabled so this is a FATAL_ERROR to let you know that you don't gain anything with this declaration.") endif() + +string(STRIP ${CMAKE_CXX_FLAGS} CMAKE_CXX_FLAGS) +string(STRIP ${CMAKE_CXX_FLAGS_RELEASE} CMAKE_CXX_FLAGS_RELEASE) +string(STRIP ${CMAKE_CXX_FLAGS_DEBUG} CMAKE_CXX_FLAGS_DEBUG) + +string(STRIP ${CMAKE_LINK_FLAGS} CMAKE_LINK_FLAGS) +string(STRIP ${CMAKE_DEPENDLIB_FLAGS} CMAKE_DEPENDLIB_FLAGS) + if (VERBOSE) cmessage (STATUS "C++ Compiler : ${CXX_COMPILER_NAME}") cmessage (STATUS " flags : ${CMAKE_CXX_FLAGS}") cmessage (STATUS " Release flags : ${CMAKE_CXX_FLAGS_RELEASE}") cmessage (STATUS " Debug flags : ${CMAKE_CXX_FLAGS_DEBUG}") cmessage (STATUS " Link Flags : ${CMAKE_LINK_FLAGS}") cmessage (STATUS " Lib Flags : ${CMAKE_DEPENDLIB_FLAGS}") endif() diff --git a/cmake/cacheVariables.cmake b/cmake/cacheVariables.cmake index cd486db..c7155b3 100644 --- a/cmake/cacheVariables.cmake +++ b/cmake/cacheVariables.cmake @@ -1,216 +1,217 @@ # 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 . ################################################################################ function(CheckAndSetDefaultEnv VARNAME DEFAULT CACHETYPE DOCSTRING ENVNAME) #cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.") if(NOT DEFINED ${VARNAME}) if(DEFINED ENV{${ENVNAME}} AND NOT $ENV{${ENVNAME}} STREQUAL "") set(${VARNAME} $ENV{${ENVNAME}} CACHE ${CACHETYPE} ${DOCSTRING}) cmessage(DEBUG " Read ${VARNAME} from ENVVAR ${ENVNAME} as $ENV{${ENVNAME}}.") else() set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING}) endif() else() set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING}) unset(${VARNAME}) endif() cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.") endfunction() function(CheckAndSetDefaultCache VARNAME DEFAULT CACHETYPE DOCSTRING) # cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.") if(NOT DEFINED ${VARNAME}) set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING}) else() set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING}) unset(${VARNAME}) endif() cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.") endfunction() function(CheckAndSetDefault VARNAME DEFAULT) # cmessage(DEBUG "Trying to assign variable ${VARNAME}.") if(NOT DEFINED ${VARNAME}) set(${VARNAME} ${DEFAULT} PARENT_SCOPE) set(${VARNAME} ${DEFAULT}) endif() cmessage(CACHE "--Set variable: \"${VARNAME}\" to \"${${VARNAME}}\".") endfunction() CheckAndSetDefaultCache(VERBOSE TRUE BOOL "Whether to configure loudly.") set (CMAKE_SKIP_BUILD_RPATH TRUE) #Changes default install path to be a subdirectory of the build dir. #Can set build dir at configure time with -DCMAKE_INSTALL_PREFIX=/install/path if(CMAKE_INSTALL_PREFIX STREQUAL "" OR CMAKE_INSTALL_PREFIX STREQUAL "/usr/local") set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}") elseif(NOT DEFINED CMAKE_INSTALL_PREFIX) set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}") endif() if(CMAKE_BUILD_TYPE STREQUAL "") set(CMAKE_BUILD_TYPE DEBUG) elseif(NOT DEFINED CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE DEBUG) endif() CheckAndSetDefaultCache(EXTRA_SETUP_SCRIPT "" PATH "The path to an extra script to inject into the NUISANCE setup script. <>") CheckAndSetDefaultCache(USE_MINIMIZER TRUE INTERNAL "Whether we are using the ROOT minimization libraries. ") CheckAndSetDefaultCache(USE_ROOT6 FALSE INTERNAL "Whether we are using the ROOT 6. ") CheckAndSetDefaultCache(USE_HEPMC FALSE BOOL "Whether to enable HepMC input support. ") CheckAndSetDefaultEnv(HEPMC "" PATH "Path to HepMC source tree root directory. Overrides environment variable \$HEPMC <>" HEPMC) CheckAndSetDefaultCache(HEPMC_MOMUNIT "GEV" STRING "HepMC momentum units [MEV|GEV]. ") CheckAndSetDefaultCache(HEPMC_LENUNIT "CM" STRING "HepMC momentum units [MM|CM]. ") CheckAndSetDefaultCache(HEPMC_USED_EP FALSE INTERNAL "Whether we built HepMC or not. ") CheckAndSetDefaultCache(USE_NEUT FALSE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. ") CheckAndSetDefaultEnv(NEUT_VERSION FALSE STRING "NEUT version string, e.g. 5.4.0. <5.4.0>" NEUT_VERSION) CheckAndSetDefaultEnv(NEUT_ROOT "" PATH "Path to NEUT source tree root directory. Overrides environment variable \$NEUT_ROOT <>" NEUT_ROOT) CheckAndSetDefaultEnv(CERN "" PATH "Path to CERNLIB source tree root directory that NEUT was built against. Overrides environment variable \$CERN <>" CERN) CheckAndSetDefaultEnv(CERN_LEVEL "" STRING "CERNLIB Library version. Overrides environment variable \$CERN_LEVEL <>" CERN_LEVEL) CheckAndSetDefaultCache(USE_NuWro FALSE BOOL "Whether to enable NuWro support. ") CheckAndSetDefaultEnv(NUWRO "" PATH "Path to NuWro source tree root directory. Overrides environment variable \$NUWRO <>" NUWRO) CheckAndSetDefaultEnv(NUWRO_INC "" PATH "Path to NuWro installed includes directory, needs to contain \"params_all.h\". Overrides environment variable \$NUWRO_INC <>" NUWRO_INC) CheckAndSetDefaultCache(NUWRO_INPUT_FILE "" FILEPATH "Path to an input NuWro event vector, which can be used to build NuWro i/o libraries. <>") CheckAndSetDefaultCache(NUWRO_BUILT_FROM_FILE FALSE INTERNAL "Whether the NuWro libraries were built by NUISANCE. ") CheckAndSetDefaultCache(USE_NuWro_RW FALSE BOOL "Whether to try and build support for NuWro reweighting. ") CheckAndSetDefaultCache(USE_NuWro_SRW_Event FALSE BOOL "Whether to use cut down NuWro reweight event format. Requires NuWro reweight. ") CheckAndSetDefaultCache(USE_GENIE FALSE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. ") CheckAndSetDefaultCache(GENIE_VERSION "AUTO" STRING "GENIE Version ") CheckAndSetDefaultEnv(GENIE "" PATH "Path to GENIE source tree root directory. Overrides environment variable \$GENIE <>" GENIE) +CheckAndSetDefaultCache(GENIE_EMPMEC_REWEIGHT FALSE BOOL "Whether to use GENIE EMP MEC reweight (requires custom GENIE) ") CheckAndSetDefaultEnv(LHAPDF_LIB "" PATH "Path to pre-built LHAPDF libraries. Overrides environment variable \$LHAPDF_LIB. <>" LHAPDF_LIB) CheckAndSetDefaultEnv(LHAPDF_INC "" PATH "Path to installed LHAPDF headers. Overrides environment variable \$LHAPDF_INC. <>" LHAPDF_INC) CheckAndSetDefaultEnv(LHAPATH "" PATH "Path to LHA PDF inputs. Overrides environment variable \$LHAPATH. <>" LHAPATH) CheckAndSetDefaultEnv(LIBXML2_LIB "" PATH "Path to pre-built LIBXML2 libraries. Overrides environment variable \$LIBXML2_LIB. <>" LIBXML2_LIB) CheckAndSetDefaultEnv(LIBXML2_INC "" PATH "Path to installed LIBXML2 headers. Overrides environment variable \$LIBXML2_INC. <>" LIBXML2_INC) CheckAndSetDefaultEnv(LOG4CPP_LIB "" PATH "Path to pre-built LOG4CPP libraries. Overrides environment variable \$LOG4CPP_LIB. <>" LOG4CPP_LIB) CheckAndSetDefaultEnv(LOG4CPP_INC "" PATH "Path to installed LOG4CPP headers. Overrides environment variable \$LOG4CPP_INC. <>" LOG4CPP_INC) CheckAndSetDefaultCache(BUILD_GEVGEN FALSE BOOL "Whether to build nuisance_gevgen app.") CheckAndSetDefaultCache(USE_T2K FALSE BOOL "Whether to enable T2KReWeight support. Requires external libraries. ") CheckAndSetDefaultEnv(T2KREWEIGHT "" PATH "Path to installed T2KREWEIGHTReWeight. Overrides environment variable \$T2KREWEIGHT. <>" T2KREWEIGHT) CheckAndSetDefaultCache(USE_NIWG FALSE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. ") CheckAndSetDefaultEnv(NIWG_ROOT "" PATH "Path to installed NIWGReWeight. Overrides environment variable \$NIWG. <>" NIWG) CheckAndSetDefaultCache(USE_MINERvA_RW FALSE BOOL "Whether to enable MINERvA ReWeight support. ") CheckAndSetDefaultEnv(PYTHIA6 "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$PYTHIA6 <>" PYTHIA6) CheckAndSetDefaultEnv(PYTHIA8 "" PATH "Path to directory containing libPythia8.so. Overrides environment variable \$PYTHIA8 <>" PYTHIA8) CheckAndSetDefaultCache(USE_PYTHIA8 FALSE BOOL "Whether to enable PYTHIA8 event support. ") CheckAndSetDefaultCache(USE_GiBUU TRUE BOOL "Whether to enable GiBUU event support. ") CheckAndSetDefaultCache(BUILD_GiBUU FALSE BOOL "Whether to build supporting GiBUU event tools along with a patched version of GiBUU. ") CheckAndSetDefaultCache(USE_NUANCE TRUE BOOL "Whether to enable NUANCE event support. ") CheckAndSetDefaultCache(USE_PROB3PP FALSE BOOL "Whether to download and compile in Prob3++ support. ") CheckAndSetDefaultCache(NO_EXTERNAL_UPDATE FALSE BOOL "Whether to perform the update target for external dependencies. Note this may produce errors for CMake < 3.8 where a bug was fixed for the feature that this option invokes. ") CheckAndSetDefaultCache(USE_GPERFTOOLS FALSE BOOL "Whether to compile in google performance tools. ") CheckAndSetDefault(NEED_PYTHIA6 FALSE) CheckAndSetDefault(NEED_PYTHIA8 FALSE) CheckAndSetDefault(NEED_ROOTEVEGEN FALSE) CheckAndSetDefault(NEED_ROOTPYTHIA6 FALSE) CheckAndSetDefaultCache(USE_OMP FALSE BOOL "Whether to enable multicore features (there currently are none...). ") CheckAndSetDefaultCache(USE_DYNSAMPLES TRUE BOOL "Whether to enable the dynamic sample loader. ") CheckAndSetDefault(NO_EXPERIMENTS FALSE) cmessage(STATUS "NO_EXPERIMENTS: ${NO_EXPERIMENTS}") CheckAndSetDefaultCache(NO_ANL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ANL samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_ArgoNeuT ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ArgoNeuT samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_BEBC ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BEBC samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_BNL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BNL samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_FNAL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build FNAL samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_GGM ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build GGM samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_K2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build K2K samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_MINERvA ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MINERvA samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_MiniBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MiniBooNE samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_T2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build T2K samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_SciBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build SciBooNE samples. <-DNO_EXPERIMENTS=FALSE>") function(SAYVARS) LIST(APPEND VARS USE_HEPMC HEPMC HEPMC_MOMUNIT HEPMC_LENUNIT HEPMC_USED_EP USE_NEUT NEUT_ROOT CERN CERN_LEVEL USE_NuWro NUWRO NUWRO_INC NUWRO_INPUT_FILE NUWRO_BUILT_FROM_FILE USE_GENIE GENIE LHAPDF_LIB LHAPDF_INC LIBXML2_LIB LIBXML2_INC LOG4CPP_LIB GENIE_LOG4CPP_INC BUILD_GEVGEN USE_T2K USE_NIWG USE_GiBUU BUILD_GiBUU USE_NUANCE NO_EXTERNAL_UPDATE USE_GPERFTOOLS NO_ANL NO_ArgoNeuT NO_BEBC NO_BNL NO_FNAL NO_GGM NO_K2K NO_MINERvA NO_MiniBooNE NO_T2K NO_SciBooNE) foreach(v ${VARS}) if(DEFINED ${v}) cmessage(DEBUG "VARIABLE: \"${v}\" = \"${${v}}\"") endif() endforeach(v) endfunction() diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx index 6703430..34f7083 100755 --- a/src/FCN/JointFCN.cxx +++ b/src/FCN/JointFCN.cxx @@ -1,1122 +1,1130 @@ #include "JointFCN.h" -#include #include "FitUtils.h" +#include //*************************************************** -JointFCN::JointFCN(TFile* outfile) { +JointFCN::JointFCN(TFile *outfile) { //*************************************************** fOutputDir = gDirectory; - if (outfile) Config::Get().out = outfile; + 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) { +JointFCN::JointFCN(std::vector samplekeys, TFile *outfile) { //*************************************************** fOutputDir = gDirectory; - if (outfile) Config::Get().out = outfile; + 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; + MeasurementBase *exp = *iter; delete exp; } for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { - ParamPull* pull = *iter; + ParamPull *pull = *iter; delete pull; } // Sort Tree - if (fIterationTree) DestroyIterationTree(); - if (fDialVals) delete fDialVals; - if (fSampleLikes) delete fSampleLikes; + if (fIterationTree) + DestroyIterationTree(); + if (fDialVals) + delete fDialVals; + if (fSampleLikes) + delete fSampleLikes; }; //*************************************************** -void JointFCN::CreateIterationTree(std::string name, FitWeight* rw) { +void JointFCN::CreateIterationTree(std::string name, FitWeight *rw) { //*************************************************** LOG(FIT) << " Creating new iteration container! " << std::endl; DestroyIterationTree(); fIterationTreeName = name; // Add sample likelihoods and ndof for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { - MeasurementBase* exp = *iter; + MeasurementBase *exp = *iter; std::string name = exp->GetName(); std::string liketag = name + "_likelihood"; fNameValues.push_back(liketag); fCurrentValues.push_back(0.0); std::string ndoftag = name + "_ndof"; fNameValues.push_back(ndoftag); fCurrentValues.push_back(0.0); } // Add Pull terms for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { - ParamPull* pull = *iter; + ParamPull *pull = *iter; std::string name = pull->GetName(); std::string liketag = name + "_likelihood"; fNameValues.push_back(liketag); fCurrentValues.push_back(0.0); std::string ndoftag = name + "_ndof"; fNameValues.push_back(ndoftag); fCurrentValues.push_back(0.0); } // Add Likelihoods fNameValues.push_back("total_likelihood"); fCurrentValues.push_back(0.0); fNameValues.push_back("total_ndof"); fCurrentValues.push_back(0.0); // Setup Containers fSampleN = fSamples.size() + fPulls.size(); fSampleLikes = new double[fSampleN]; fSampleNDOF = new int[fSampleN]; // Add Dials std::vector dials = rw->GetDialNames(); for (size_t i = 0; i < dials.size(); i++) { fNameValues.push_back(dials[i]); fCurrentValues.push_back(0.0); } fNDials = dials.size(); fDialVals = new double[fNDials]; // Set IterationTree Flag fIterationTree = true; } //*************************************************** void JointFCN::DestroyIterationTree() { //*************************************************** fIterationCount.clear(); fCurrentValues.clear(); fNameValues.clear(); fIterationValues.clear(); } //*************************************************** void JointFCN::WriteIterationTree() { //*************************************************** LOG(FIT) << "Writing iteration tree" << std::endl; // Make a new TTree - TTree* itree = + TTree *itree = new TTree(fIterationTreeName.c_str(), fIterationTreeName.c_str()); - double* vals = new double[fNameValues.size()]; + double *vals = new double[fNameValues.size()]; int count = 0; itree->Branch("iteration", &count, "Iteration/I"); for (size_t i = 0; i < fNameValues.size(); i++) { itree->Branch(fNameValues[i].c_str(), &vals[i], (fNameValues[i] + "/D").c_str()); } // Fill Iterations for (size_t i = 0; i < fIterationValues.size(); i++) { std::vector itervals = fIterationValues[i]; // Fill iteration state count = fIterationCount[i]; for (size_t j = 0; j < itervals.size(); j++) { vals[j] = itervals[j]; } // Save to TTree itree->Fill(); } // Write to file itree->Write(); } //*************************************************** -void JointFCN::FillIterationTree(FitWeight* rw) { +void JointFCN::FillIterationTree(FitWeight *rw) { //*************************************************** // Loop over samples count int count = 0; for (int i = 0; i < fSampleN; i++) { fCurrentValues[count++] = fSampleLikes[i]; fCurrentValues[count++] = double(fSampleNDOF[i]); } // Fill Totals fCurrentValues[count++] = fLikelihood; fCurrentValues[count++] = double(fNDOF); // Loop Over Parameter Counts rw->GetAllDials(fDialVals, fNDials); for (int i = 0; i < fNDials; i++) { fCurrentValues[count++] = double(fDialVals[i]); } // Push Back Into Container fIterationCount.push_back(fCurIter); fIterationValues.push_back(fCurrentValues); } //*************************************************** -double JointFCN::DoEval(const double* x) { +double JointFCN::DoEval(const double *x) { //*************************************************** // WEIGHT ENGINE fDialChanged = FitBase::GetRW()->HasRWDialChanged(x); FitBase::GetRW()->UpdateWeightEngine(x); if (fDialChanged) { FitBase::GetRW()->Reconfigure(); FitBase::EvtManager().ResetWeightFlags(); } if (LOG_LEVEL(REC)) { FitBase::GetRW()->Print(); } // SORT SAMPLES ReconfigureSamples(); // GET TEST STAT fLikelihood = GetLikelihood(); fNDOF = GetNDOF(); // PRINT PROGRESS LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood << std::endl; // UPDATE TREE - if (fIterationTree) FillIterationTree(FitBase::GetRW()); + 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; + 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; + 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; + 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; + 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); + 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; + 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; + 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; +std::vector JointFCN::GetInputList() { + std::vector InputList; fIsAllSplines = true; MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { - MeasurementBase* exp = (*iterSam); + MeasurementBase *exp = (*iterSam); - std::vector subsamples = exp->GetSubSamples(); + std::vector subsamples = exp->GetSubSamples(); for (size_t i = 0; i < subsamples.size(); i++) { - InputHandlerBase* inp = subsamples[i]->GetInput(); + 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; +std::vector JointFCN::GetSubSampleList() { + std::vector SampleList; MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { - MeasurementBase* exp = (*iterSam); + MeasurementBase *exp = (*iterSam); - std::vector subsamples = exp->GetSubSamples(); + 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); + 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(); + std::vector::iterator inp_iter = fInputList.begin(); if (fIsAllSplines) { for (; inp_iter != fInputList.end(); inp_iter++) { - InputHandlerBase* curinput = (*inp_iter); + InputHandlerBase *curinput = (*inp_iter); // Tell reader in each BaseEvent it needs a Reconfigure next weight calc. - BaseFitEvt* curevent = curinput->FirstBaseEvent(); + 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); + InputHandlerBase *curinput = (*inp_iter); // Get event information - FitEvent* curevent = curinput->FirstNuisanceEvent(); + FitEvent *curevent = curinput->FirstNuisanceEvent(); curinput->CreateCache(); int i = 0; int nevents = curinput->GetNEvents(); int countwidth = nevents / 10; // Start event loop iterating until we get a NULL pointer. while (curevent) { // Get Event Weight // The reweighting weight curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent); // The Custom weight and reweight - curevent->Weight = curevent->RWWeight * curevent->InputWeight * curevent->CustomWeight; - //double rwweight = curevent->Weight; + curevent->Weight = + curevent->RWWeight * curevent->InputWeight * curevent->CustomWeight; + // double rwweight = curevent->Weight; // std::cout << "RWWeight = " << curevent->RWWeight << " " << // curevent->InputWeight << std::endl; // Logging // std::cout << CHECKLOG(1) << std::endl; if (LOGGING(REC)) { if (countwidth && (i % countwidth == 0)) { QLOG(REC, curinput->GetName() << " : Processed " << i << " events. [M, W] = [" << curevent->Mode << ", " << curevent->Weight << "]"); } } // Setup flag for if signal found in at least one sample bool foundsignal = false; // Create a new signal bitset for this event std::vector signalbitset(fSubSampleList.size()); // Create a new signal box vector for this event - std::vector signalboxes; + std::vector signalboxes; // Start measurement iterator size_t measitercount = 0; - std::vector::iterator meas_iter = + std::vector::iterator meas_iter = fSubSampleList.begin(); // Loop over all subsamples (sub in JointMeas) for (; meas_iter != fSubSampleList.end(); meas_iter++) { - MeasurementBase* curmeas = (*meas_iter); + 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); + 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); + 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; + int mem = ( // sizeof(fSignalEventBoxes) + + // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) + + sizeof(MeasurementVariableBox1D) * fillcount) * + 1E-6; LOG(REC) << " -> Saved " << fillcount << " signal boxes for faster access. (~" << mem << " MB)" << std::endl; if (fIsAllSplines and !fSignalEventSplines.empty()) { int splmem = sizeof(float) * fSignalEventSplines.size() * fSignalEventSplines[0].size() * 1E-6; LOG(REC) << " -> Saved " << fillcount << " " << fSignalEventSplines.size() << " spline sets into memory. (~" << splmem << " MB)" << std::endl; } } LOG(REC) << "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart << std::endl; // Check SignalReconfigures works for all samples if (savesignal) { double likefull = GetLikelihood(); ReconfigureFastUsingManager(); double likefast = GetLikelihood(); if (fabs(likefull - likefast) > 0.0001) { ERROR(FTL, "Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast); - ERROR(FTL, - "This means some samples you are using are not setup to use " - "SignalReconfigures=1"); + ERROR(FTL, "This means some samples you are using are not setup to use " + "SignalReconfigures=1"); ERROR(FTL, "Please turn OFF signal reconfigures."); throw; } else { LOG(FIT) << "Likelihoods for FULL and FAST match. Will use FAST next time." << std::endl; } } // End of reconfigure return; }; //*************************************************** void JointFCN::ReconfigureFastUsingManager() { //*************************************************** LOG(FIT) << " -> Doing FAST using manager" << std::endl; // Get Start time for profilling int timestart = time(NULL); // Reset all samples MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { - MeasurementBase* exp = (*iterSam); + MeasurementBase *exp = (*iterSam); exp->ResetAll(); } // Check for saved variables if not do a full reconfigure. if (fSignalEventFlags.empty()) { ERR(WRN) << "Signal Flags Empty! Using normal manager." << std::endl; ReconfigureUsingManager(); return; } bool fFillNuisanceEvent = FitPar::Config().GetParB("FullEventOnSignalReconfigure"); // Setup fast vector iterators. std::vector::iterator inpsig_iter = fSignalEventFlags.begin(); - std::vector >::iterator box_iter = + std::vector >::iterator box_iter = fSignalEventBoxes.begin(); std::vector >::iterator spline_iter = fSignalEventSplines.begin(); std::vector >::iterator samsig_iter = fSampleSignalFlags.begin(); int splinecount = 0; // Setup stuff for logging int fillcount = 0; int nevents = fSignalEventFlags.size(); int countwidth = nevents / 10; // If All Splines tell splines they need a reconfigure. - std::vector::iterator inp_iter = fInputList.begin(); + 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); + InputHandlerBase *curinput = (*inp_iter); // Tell each fSplineRead in BaseFitEvent to reconf next weight calc - BaseFitEvt* curevent = curinput->FirstBaseEvent(); + BaseFitEvt *curevent = curinput->FirstBaseEvent(); if (curevent->fSplineRead) curevent->fSplineRead->SetNeedsReconfigure(true); } } // Loop over all possible spline inputs - double* coreeventweights = new double[fSignalEventBoxes.size()]; + 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(); + 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) { curevent = curinput->GetNuisanceEvent(i); } else { curevent = curinput->GetBaseEvent(i); } } else { curevent->fSplineCoeff = &fSignalEventSplines[splinecount][0]; } curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent); - curevent->Weight = curevent->RWWeight * curevent->InputWeight*curevent->CustomWeight; + curevent->Weight = + curevent->RWWeight * curevent->InputWeight * curevent->CustomWeight; 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 = + std::vector::iterator subbox_iter = (*box_iter).begin(); // Loop over all sub measurements. - std::vector::iterator meas_iter = fSubSampleList.begin(); + std::vector::iterator meas_iter = fSubSampleList.begin(); for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) { - MeasurementBase* curmeas = (*meas_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); + 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; + MeasurementBase *exp = *iter; double like = exp->GetLikelihood(); double ndof = exp->GetNDOF(); std::string name = exp->GetName(); likes.push_back(like); ndofs.push_back(ndof); names.push_back(name); } - TH1D likehist = TH1D("likelihood_hist", "likelihood_hist", likes.size(), 0.0, - double(likes.size())); - TH1D ndofhist = - TH1D("ndof_hist", "ndof_hist", ndofs.size(), 0.0, double(ndofs.size())); - TH1D divhist = TH1D("likedivndof_hist", "likedivndof_hist", likes.size(), 0.0, - double(likes.size())); - for (int i = 0; i < likehist.GetNbinsX(); i++) { - likehist.SetBinContent(i + 1, likes[i]); - ndofhist.SetBinContent(i + 1, ndofs[i]); - if (ndofs[i] != 0.0) { - divhist.SetBinContent(i + 1, likes[i] / ndofs[i]); + if (likes.size()) { + TH1D likehist = TH1D("likelihood_hist", "likelihood_hist", likes.size(), + 0.0, double(likes.size())); + TH1D ndofhist = + TH1D("ndof_hist", "ndof_hist", ndofs.size(), 0.0, double(ndofs.size())); + TH1D divhist = TH1D("likedivndof_hist", "likedivndof_hist", likes.size(), + 0.0, double(likes.size())); + for (int i = 0; i < likehist.GetNbinsX(); i++) { + likehist.SetBinContent(i + 1, likes[i]); + ndofhist.SetBinContent(i + 1, ndofs[i]); + if (ndofs[i] != 0.0) { + divhist.SetBinContent(i + 1, likes[i] / ndofs[i]); + } + likehist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); + ndofhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); + divhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); } - likehist.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(); } - 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; + MeasurementBase *exp = *iter; exp->Write(); } // Save Pull Terms for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { - ParamPull* pull = *iter; + ParamPull *pull = *iter; pull->Write(); } if (FitPar::Config().GetParB("EventManager")) { // Get list of inputs - std::map fInputs = + std::map fInputs = FitBase::EvtManager().GetInputs(); - std::map::const_iterator iterInp; + std::map::const_iterator iterInp; for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) { - InputHandlerBase* input = (iterInp->second); + 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; + MeasurementBase *exp = *iter; exp->SetFakeDataValues(fakeinput); } return; } //*************************************************** void JointFCN::ThrowDataToy() { //*************************************************** for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { - MeasurementBase* exp = *iter; + MeasurementBase *exp = *iter; exp->ThrowDataToy(); } return; } //*************************************************** std::vector JointFCN::GetAllNames() { //*************************************************** // Vect of all likelihoods and total std::vector namevect; // Loop over samples first for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { - MeasurementBase* exp = *iter; + MeasurementBase *exp = *iter; // Get Likelihoods and push to vector namevect.push_back(exp->GetName()); } // Loop over pulls second for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { - ParamPull* pull = *iter; + 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(); iter++) { - MeasurementBase* exp = *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++) { - ParamPull* pull = *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(); iter++) { - MeasurementBase* exp = *iter; + MeasurementBase *exp = *iter; // Get Likelihoods and push to vector int singlendof = exp->GetNDOF(); ndofvect.push_back(singlendof); total_ndof += singlendof; } // Loop over pulls second for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { - ParamPull* pull = *iter; + 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/JointMeas1D.cxx b/src/FitBase/JointMeas1D.cxx index 6c3e699..6d99713 100644 --- a/src/FitBase/JointMeas1D.cxx +++ b/src/FitBase/JointMeas1D.cxx @@ -1,2303 +1,2303 @@ // 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 "JointMeas1D.h" //******************************************************************** JointMeas1D::JointMeas1D(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; 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"; fIsFix = false; fIsShape = false; fIsFree = false; fIsDiag = false; fIsFull = false; fAddNormPen = false; fIsMask = false; fIsChi2SVD = false; fIsRawEvents = false; fIsDifXSec = false; fIsEnu1D = false; // Inputs fInput = NULL; fRW = NULL; // Extra Histograms fMCHist_Modes = NULL; for (std::vector::const_iterator iter = fSubChain.begin(); iter != fSubChain.end(); iter++) { MeasurementBase* exp = *iter; if (exp) delete exp; } fSubChain.clear(); // Flags for Joint Measurements fIsRatio = false; fIsSummed = false; fSaveSubMeas = false; fIsJoint = true; } //******************************************************************** void JointMeas1D::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); return; } //******************************************************************** JointMeas1D::~JointMeas1D(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 (fCovar) delete fCovar; if (fInvert) delete fInvert; if (fDecomp) delete fDecomp; for (std::vector::const_iterator iter = fSubChain.begin(); iter != fSubChain.end(); iter++) { MeasurementBase* exp = *iter; if (exp) delete exp; } fSubChain.clear(); } //******************************************************************** SampleSettings JointMeas1D::LoadSampleSettings(nuiskey samplekey){ //******************************************************************** SampleSettings s = MeasurementBase::LoadSampleSettings(samplekey); // Parse Inputs fSubInFiles.clear(); std::vector entries = GeneralUtils::ParseToStr(s.GetS("input"), ";"); if (entries.size() < 2) { ERR(FTL) << "Joint measurement expected to recieve at least two semi-colon " "separated input files, but recieved: \"" << s.GetS("input") << "\"" << std::endl; throw; } std::vector first_file_descriptor = GeneralUtils::ParseToStr(entries.front(), ":"); if (first_file_descriptor.size() != 2) { ERR(FTL) << "Found Joint measurement where the input file had no type: \"" << s.GetS("input") << "\", expected \"INPUTTYPE:File.root;File2.root\"." << std::endl; throw; } std::string inpType = first_file_descriptor[0]; for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { if (GeneralUtils::ParseToStr(*iter, ":").size() != 2) { std::stringstream ss(""); ss << inpType << ":" << (*iter); fSubInFiles.push_back(ss.str()); } else { fSubInFiles.push_back(*iter); } } return s; } //******************************************************************** void JointMeas1D::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); } // Parse Inputs fSubInFiles.clear(); std::vector entries = GeneralUtils::ParseToStr(fSettings.GetS("input"), ";"); if (entries.size() < 2) { ERR(FTL) << "Joint measurement expected to recieve at least two semi-colon " "separated input files, but recieved: \"" << fSettings.GetS("input") << "\"" << std::endl; throw; } std::vector first_file_descriptor = GeneralUtils::ParseToStr(entries.front(), ":"); if (first_file_descriptor.size() != 2) { ERR(FTL) << "Found Joint measurement where the input file had no type: \"" << fSettings.GetS("input") << "\", expected \"INPUTTYPE:File.root;File2.root\"." << std::endl; throw; } std::string inpType = first_file_descriptor[0]; for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { if (GeneralUtils::ParseToStr(*iter, ":").size() != 2) { std::stringstream ss(""); ss << inpType << ":" << (*iter); fSubInFiles.push_back(ss.str()); } else { fSubInFiles.push_back(*iter); } } // 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; } } if (!fRW) fRW = FitBase::GetRW(); LOG(SAM) << "Finalised Sample Settings" << std::endl; } //******************************************************************** void JointMeas1D::SetDataFromTextFile(std::string datafile) { //******************************************************************** LOG(SAM) << "Reading data from text file: " << datafile << std::endl; fDataHist = PlotUtils::GetTH1DFromFile(datafile, fSettings.GetName() + "_data", fSettings.GetFullTitles()); } //******************************************************************** void JointMeas1D::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 JointMeas1D::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 JointMeas1D::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 JointMeas1D::SetCovarFromTextFile(std::string covfile, int dim) { //******************************************************************** LOG(SAM) << "Reading covariance from text file: " << covfile << std::endl; fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim); covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void JointMeas1D::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 JointMeas1D::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 JointMeas1D::SetCovarInvertFromTextFile(std::string covfile, int dim) { //******************************************************************** 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 JointMeas1D::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 JointMeas1D::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 JointMeas1D::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++) { // Note that there is a factor of 1E76 here. It is very silly indeed. // However, if you remove it, you also need to fix the factors of 1E38 added to the chi2 calculations! (*temp_cov)(i, j) = (*temp_cov)(i, j) * fDataHist->GetBinError(i + 1) * fDataHist->GetBinError(j + 1) * 1E76; } } (*fFullCovar) += (*temp_cov); delete temp_cov; } covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } //******************************************************************** void JointMeas1D::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 JointMeas1D::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 JointMeas1D::SetCholDecompFromTextFile(std::string covfile, int dim) { //******************************************************************** 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 JointMeas1D::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 JointMeas1D::ScaleData(double scale) { //******************************************************************** fDataHist->Scale(scale); } //******************************************************************** void JointMeas1D::ScaleCovar(double scale) { //******************************************************************** (*fFullCovar) *= scale; (*covar) *= 1.0 / scale; (*fDecomp) *= sqrt(scale); } //******************************************************************** void JointMeas1D::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) << "JointMeas1D::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 JointMeas1D::FinaliseMeasurement() { //******************************************************************** LOG(SAM) << "Finalising Measurement: " << fName << std::endl; // 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 out until scaling is used consistently... StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1E-38); // 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(), 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); } // 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 JointMeas1D::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("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 JointMeas1D::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 JointMeas1D::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 JointMeas1D::ResetAll() { //******************************************************************** fMCHist->Reset(); fMCFine->Reset(); fMCStat->Reset(); return; }; //******************************************************************** void JointMeas1D::FillHistograms() { //******************************************************************** if (Signal) { fMCHist->Fill(fXVar, Weight); fMCFine->Fill(fXVar, Weight); fMCStat->Fill(fXVar, 1.0); if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, fXVar, Weight); } return; }; //******************************************************************** void JointMeas1D::ScaleEvents() { //******************************************************************** LOG(FIT) << "Scaling JointMeas1D" << std::endl; // 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); // } // 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 JointMeas1D::ApplyNormScale(double norm) { //******************************************************************** fCurrentNorm = norm; fMCHist->Scale(1.0 / norm); fMCFine->Scale(1.0 / norm); return; }; /* Statistic Functions - Outsources to StatUtils */ //******************************************************************** int JointMeas1D::GetNDOF() { //******************************************************************** int ndof = fDataHist->GetNbinsX(); if (fMaskHist) ndof -= fMaskHist->Integral(); return ndof; } //******************************************************************** double JointMeas1D::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; 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 JointMeas1D::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 JointMeas1D::ResetFakeData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH1D*)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str()); } } //******************************************************************** void JointMeas1D::ResetData() { //******************************************************************** if (fIsFakeData) { if (fDataHist) delete fDataHist; fDataHist = (TH1D*)fDataOrig->Clone((fSettings.GetName() + "_data").c_str()); } fIsFakeData = false; } //******************************************************************** void JointMeas1D::ThrowCovariance() { //******************************************************************** // Take a fDecomposition and use it to throw the current dataset. // Requires fDataTrue also be set incase used repeatedly. if (fDataHist) delete fDataHist; fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); return; }; //******************************************************************** void JointMeas1D::ThrowDataToy(){ //******************************************************************** if (!fDataTrue) fDataTrue = (TH1D*) fDataHist->Clone(); if (fMCHist) delete fMCHist; fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar); } /* Access Functions */ //******************************************************************** TH1D* JointMeas1D::GetMCHistogram() { //******************************************************************** if (!fMCHist) return fMCHist; std::ostringstream chi2; - chi2 << std::setprecision(5) << this->GetLikelihood(); + chi2 << "#chi^{2}=" << 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* JointMeas1D::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 JointMeas1D::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 GetDataHistogram()->Write(); GetMCHistogram()->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()) GetEventHistogram()->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(); if (fSaveSubMeas) { for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; exp->Write(drawOpt); } } // Returning LOG(SAM) << "Written Histograms: " << fName << std::endl; return; } //******************************************************************** void JointMeas1D::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 JointMeas1D::WriteShapePlot() { //******************************************************************** 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); std::stringstream ss; - ss << shapeScale; + ss << "Scale=" << shapeScale; mcShape->SetTitle(ss.str().c_str()); mcShape->SetLineWidth(3); mcShape->SetLineStyle(7); mcShape->Write(); delete mcShape; } //******************************************************************** void JointMeas1D::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; } // 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 "JointMeas1D.h" /* Constructor/Deconstuctor */ /* Setup Functions */ //******************************************************************** void JointMeas1D::SetupMeasurement(std::string input, std::string type, FitWeight* rw, std::string fkdt) { //******************************************************************** // For joint samples, input files are given as a semi-colon seperated list. // Parse this list and save it for later, and set up the types etc. if (FitPar::Config().GetParB("EventManager")) { ERR(FTL) << "Event Manager does not yet work with JointMeas1D Samples" << std::endl; ERR(FTL) << "If you want good predictions for " << fName << " then run with it turned off! (-q EventManager=0)" << std::endl; } fSubInFiles.clear(); std::vector entries = GeneralUtils::ParseToStr(input, ";"); if (entries.size() < 2) { ERR(FTL) << "Joint measurement expected to recieve at least two semi-colon " "separated input files, but recieved: \"" << input << "\"" << std::endl; throw; } std::vector first_file_descriptor = GeneralUtils::ParseToStr(entries.front(), ":"); if (first_file_descriptor.size() != 2) { ERR(FTL) << "Found Joint measurement where the input file had no type: \"" << input << "\", expected \"INPUTTYPE:File.root;File2.root\"." << std::endl; throw; } std::string inpType = first_file_descriptor[0]; for (std::vector::iterator iter = entries.begin(); iter != entries.end(); iter++) { if (GeneralUtils::ParseToStr(*iter, ":").size() != 2) { std::stringstream ss(""); ss << inpType << ":" << (*iter); fSubInFiles.push_back(ss.str()); } else { fSubInFiles.push_back(*iter); } } // Set Engine and Fake Data fRW = rw; fFakeDataInput = fkdt; // Set Fit Options SetFitOptions(type); return; } /* XSec Functions */ //******************************************************************** double JointMeas1D::TotalIntegratedFlux(std::string intOpt, double low, double high) { //******************************************************************** double totalflux = 0.0; // Destroy the job for sub samples for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; double expflux = exp->TotalIntegratedFlux(intOpt, low, high); // Fill flux options if (fIsRatio) { totalflux = expflux; break; } if (fIsSummed) { totalflux += expflux; } } return totalflux; } /* Reconfigure Functions */ //******************************************************************** void JointMeas1D::Reconfigure() { //******************************************************************** for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; exp->Reconfigure(); } ConvertEventRates(); return; } //******************************************************************** void JointMeas1D::ConvertEventRates() { //******************************************************************** // Apply Event Scaling for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = static_cast(*expIter); exp->ScaleEvents(); } // Joint function called by top level class MakePlots(); // Do Final Normalisation ApplyNormScale(fRW->GetSampleNorm(this->fName)); } //******************************************************************** void JointMeas1D::MakePlots() { //******************************************************************** // Reset the 1D histograms but not the subClasses ResetAll(); // If Summed if (fIsSummed) { for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = static_cast(*expIter); this->fMCHist->Add(exp->GetMCList().at(0)); this->fMCFine->Add(exp->GetFineList().at(0)); } return; } // If Ratio if (fIsRatio) { int sample = 0; for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; if (sample == 0) { this->fMCHist->Add(exp->GetMCList().at(0)); this->fMCFine->Add(exp->GetFineList().at(0)); } else if (sample == 1) { this->fMCHist->Divide(exp->GetMCList().at(0)); this->fMCFine->Divide(exp->GetFineList().at(0)); } else { break; } sample++; } return; } return; } /* Access Functions */ //******************************************************************** std::vector JointMeas1D::GetMCList() { //******************************************************************** // Make Default Vector std::vector tempVect; tempVect.push_back(this->fMCHist); // Return vector from all sub samples for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; std::vector subTempVect = exp->GetMCList(); for (UInt_t i = 0; i < subTempVect.size(); i++) { tempVect.push_back(subTempVect.at(i)); } } return tempVect; } //******************************************************************** std::vector JointMeas1D::GetDataList() { //******************************************************************** // Make Default Vector std::vector tempVect; tempVect.push_back(this->fDataHist); // Return vector from all sub samples for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; std::vector subTempVect = exp->GetDataList(); for (UInt_t i = 0; i < subTempVect.size(); i++) { tempVect.push_back(subTempVect.at(i)); } } return tempVect; } //******************************************************************** std::vector JointMeas1D::GetFineList() { //******************************************************************** // Make Default Vector std::vector tempVect; tempVect.push_back(this->fMCFine); // Return vector from all sub samples for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; std::vector subTempVect = exp->GetFineList(); for (UInt_t i = 0; i < subTempVect.size(); i++) { tempVect.push_back(subTempVect.at(i)); } } return tempVect; } //******************************************************************** std::vector JointMeas1D::GetMaskList() { //******************************************************************** // Make Default Vector std::vector tempVect; tempVect.push_back(this->fMaskHist); // Return vector from all sub samples for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; std::vector subTempVect = exp->GetMaskList(); for (UInt_t i = 0; i < subTempVect.size(); i++) { tempVect.push_back(subTempVect.at(i)); } } return tempVect; } //******************************************************************** std::vector JointMeas1D::GetFluxList() { //******************************************************************** // Make Default Vector std::vector tempVect; tempVect.push_back(MeasurementBase::GetFluxHistogram()); // Return vector from all sub samples for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; std::vector subTempVect = exp->GetFluxList(); for (UInt_t i = 0; i < subTempVect.size(); i++) { tempVect.push_back(subTempVect.at(i)); } } return tempVect; } //******************************************************************** std::vector JointMeas1D::GetEventRateList() { //******************************************************************** // Make Default Vector std::vector tempVect; tempVect.push_back(MeasurementBase::GetEventHistogram()); // Return vector from all sub samples for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; std::vector subTempVect = exp->GetEventRateList(); for (UInt_t i = 0; i < subTempVect.size(); i++) { tempVect.push_back(subTempVect.at(i)); } } return tempVect; } //******************************************************************** std::vector JointMeas1D::GetXSecList() { //******************************************************************** // Make Default Vector std::vector tempVect; tempVect.push_back(MeasurementBase::GetXSecHistogram()); // Return vector from all sub samples for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; std::vector subTempVect = exp->GetXSecList(); for (UInt_t i = 0; i < subTempVect.size(); i++) { tempVect.push_back(subTempVect.at(i)); } } return tempVect; } //******************************************************************** TH1D* JointMeas1D::GetCombinedFlux() { //******************************************************************** TH1D* newflux = NULL; int sample = 0; for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; // Get flux from experiment std::vector fluxVect = exp->GetFluxList(); // Setup newflux if (sample == 0) { newflux = (TH1D*)fluxVect.at(0); newflux->Reset(); } // Add all fluxes for (UInt_t i = 0; i < fluxVect.size(); i++) { newflux->Add((TH1D*)fluxVect.at(i)); sample++; } } if (!newflux){ ERR(FTL) << "No combined flux setup in JointMeas1D" << std::endl; } return newflux; } //******************************************************************** TH1D* JointMeas1D::GetCombinedEventRate() { //******************************************************************** TH1D* newflux = NULL; int sample = 0; for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { MeasurementBase* exp = *expIter; // Get flux from experiment std::vector fluxVect = exp->GetFluxList(); // Setup newflux if (sample == 0) { newflux = (TH1D*)fluxVect.at(0); newflux->Reset(); } // Add all fluxes for (UInt_t i = 0; i < fluxVect.size(); i++) { newflux->Add(fluxVect.at(i)); sample++; } } if (!newflux){ ERR(FTL) << "No combined event rate setup in JointMeas1D" << std::endl; } return newflux; } //******************************************************************** std::vector JointMeas1D::GetSubSamples() { //******************************************************************** std::vector exps; for (std::vector::const_iterator expIter = fSubChain.begin(); expIter != fSubChain.end(); expIter++) { exps.push_back(*expIter); } return exps; } //// CRAP TO BE REMOVED //******************************************************************** void JointMeas1D::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 JointMeas1D::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 JointMeas1D::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 JointMeas1D::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 JointMeas1D::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 JointMeas1D::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 JointMeas1D::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; }; // std::vector JointMeas1D::GetMCList(void){ // std::vector temp; // return temp; // } // std::vector JointMeas1D::GetDataList(void){ // std::vector temp; // return temp; // } // std::vector JointMeas1D::GetMaskList(void){ // std::vector temp; // return temp; // } // std::vector JointMeas1D::GetFineList(void){ // std::vector temp; // return temp; // } diff --git a/src/InputHandler/InputHandler.cxx b/src/InputHandler/InputHandler.cxx index dd609c4..d606944 100644 --- a/src/InputHandler/InputHandler.cxx +++ b/src/InputHandler/InputHandler.cxx @@ -1,297 +1,298 @@ // 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 "InputHandler.h" #include "InputUtils.h" InputHandlerBase::InputHandlerBase() { fName = ""; fFluxHist = NULL; fEventHist = NULL; fNEvents = 0; fNUISANCEEvent = NULL; fBaseEvent = NULL; kRemoveUndefParticles = FitPar::Config().GetParB("RemoveUndefParticles"); kRemoveFSIParticles = FitPar::Config().GetParB("RemoveFSIParticles"); kRemoveNuclearParticles = FitPar::Config().GetParB("RemoveNuclearParticles"); fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); fTTreePerformance = NULL; }; InputHandlerBase::~InputHandlerBase() { if (fFluxHist) delete fFluxHist; if (fEventHist) delete fEventHist; // if (fXSecHist) delete fXSecHist; // if (fNUISANCEEvent) delete fNUISANCEEvent; jointfluxinputs.clear(); jointeventinputs.clear(); jointindexlow.clear(); jointindexhigh.clear(); jointindexallowed.clear(); jointindexscale.clear(); // if (fTTreePerformance) { // fTTreePerformance->SaveAs(("ttreeperfstats_" + fName + // ".root").c_str()); // } } void InputHandlerBase::Print(){}; TH1D* InputHandlerBase::GetXSecHistogram(void) { fXSecHist = (TH1D*)fFluxHist->Clone(); fXSecHist->Divide(fEventHist); return fXSecHist; }; double InputHandlerBase::PredictedEventRate(double low, double high, std::string intOpt) { Int_t minBin = fEventHist->GetXaxis()->FindFixBin(low); Int_t maxBin = fEventHist->GetXaxis()->FindFixBin(high); if ((fEventHist->IsBinOverflow(minBin) && (low != -9999.9))) { minBin = 1; } if ((fEventHist->IsBinOverflow(maxBin) && (high != -9999.9))) { maxBin = fEventHist->GetXaxis()->GetNbins() + 1; } // If we are within a single bin if (minBin == maxBin) { // Get the contained fraction of the single bin's width return ((high - low) / fEventHist->GetXaxis()->GetBinWidth(minBin)) * fEventHist->Integral(minBin, minBin, intOpt.c_str()); } double lowBinUpEdge = fEventHist->GetXaxis()->GetBinUpEdge(minBin); double highBinLowEdge = fEventHist->GetXaxis()->GetBinLowEdge(maxBin); double lowBinfracIntegral = ((lowBinUpEdge - low) / fEventHist->GetXaxis()->GetBinWidth(minBin)) * fEventHist->Integral(minBin, minBin, intOpt.c_str()); double highBinfracIntegral = ((high - highBinLowEdge) / fEventHist->GetXaxis()->GetBinWidth(maxBin)) * fEventHist->Integral(maxBin, maxBin, intOpt.c_str()); // If they are neighbouring bins if ((minBin + 1) == maxBin) { std::cout << "Get lowfrac + highfrac" << std::endl; // Get the contained fraction of the two bin's width return lowBinfracIntegral + highBinfracIntegral; } double ContainedIntegral = fEventHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str()); // If there are filled bins between them return lowBinfracIntegral + highBinfracIntegral + ContainedIntegral; }; double InputHandlerBase::TotalIntegratedFlux(double low, double high, std::string intOpt) { Int_t minBin = fFluxHist->GetXaxis()->FindFixBin(low); Int_t maxBin = fFluxHist->GetXaxis()->FindFixBin(high); if ((fFluxHist->IsBinOverflow(minBin) && (low != -9999.9))) { minBin = 1; } if ((fFluxHist->IsBinOverflow(maxBin) && (high != -9999.9))) { - maxBin = fFluxHist->GetXaxis()->GetNbins() + 1; + maxBin = fFluxHist->GetXaxis()->GetNbins(); + high = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin+1); } // If we are within a single bin if (minBin == maxBin) { // Get the contained fraction of the single bin's width return ((high - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) * fFluxHist->Integral(minBin, minBin, intOpt.c_str()); } double lowBinUpEdge = fFluxHist->GetXaxis()->GetBinUpEdge(minBin); double highBinLowEdge = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin); double lowBinfracIntegral = ((lowBinUpEdge - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) * fFluxHist->Integral(minBin, minBin, intOpt.c_str()); double highBinfracIntegral = ((high - highBinLowEdge) / fFluxHist->GetXaxis()->GetBinWidth(maxBin)) * fFluxHist->Integral(maxBin, maxBin, intOpt.c_str()); // If they are neighbouring bins if ((minBin + 1) == maxBin) { std::cout << "Get lowfrac + highfrac" << std::endl; // Get the contained fraction of the two bin's width return lowBinfracIntegral + highBinfracIntegral; } double ContainedIntegral = fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str()); // If there are filled bins between them return lowBinfracIntegral + highBinfracIntegral + ContainedIntegral; } std::vector InputHandlerBase::GetFluxList(void) { return std::vector(1, fFluxHist); }; std::vector InputHandlerBase::GetEventList(void) { return std::vector(1, fEventHist); }; std::vector InputHandlerBase::GetXSecList(void) { return std::vector(1, GetXSecHistogram()); }; FitEvent* InputHandlerBase::FirstNuisanceEvent() { fCurrentIndex = 0; return GetNuisanceEvent(fCurrentIndex); }; FitEvent* InputHandlerBase::NextNuisanceEvent() { fCurrentIndex++; if ((fMaxEvents != -1) && (fCurrentIndex > fMaxEvents)) { return NULL; } return GetNuisanceEvent(fCurrentIndex); }; BaseFitEvt* InputHandlerBase::FirstBaseEvent() { fCurrentIndex = 0; return GetBaseEvent(fCurrentIndex); }; BaseFitEvt* InputHandlerBase::NextBaseEvent() { fCurrentIndex++; if (jointinput and fMaxEvents != -1) { while (fCurrentIndex < jointindexlow[jointindexswitch] || fCurrentIndex >= jointindexhigh[jointindexswitch]) { jointindexswitch++; // Loop Around if (jointindexswitch == jointindexlow.size()) { jointindexswitch = 0; } } if (fCurrentIndex > jointindexlow[jointindexswitch] + jointindexallowed[jointindexswitch]) { fCurrentIndex = jointindexlow[jointindexswitch]; } } return GetBaseEvent(fCurrentIndex); }; void InputHandlerBase::RegisterJointInput(std::string input, int n, TH1D* f, TH1D* e) { if (jointfluxinputs.size() == 0) { jointindexswitch = 0; fNEvents = 0; } // Push into individual input vectors jointfluxinputs.push_back((TH1D*)f->Clone()); jointeventinputs.push_back((TH1D*)e->Clone()); jointindexlow.push_back(fNEvents); jointindexhigh.push_back(fNEvents + n); fNEvents += n; // Add to the total flux/event hist if (!fFluxHist) fFluxHist = (TH1D*)f->Clone(); else fFluxHist->Add(f); if (!fEventHist) fEventHist = (TH1D*)e->Clone(); else fEventHist->Add(e); } void InputHandlerBase::SetupJointInputs() { if (jointeventinputs.size() <= 1) { jointinput = false; } else if (jointeventinputs.size() > 1) { jointinput = true; jointindexswitch = 0; } fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); if (fMaxEvents != -1 and jointeventinputs.size() > 1) { THROW("Can only handle joint inputs when config MAXEVENTS = -1!"); } if (jointeventinputs.size() > 1) { ERROR(WRN, "GiBUU sample contains multiple inputs. This will only work for " "samples that expect multi-species inputs. If this sample does, you " "can ignore this warning."); } for (size_t i = 0; i < jointeventinputs.size(); i++) { double scale = double(fNEvents) / fEventHist->Integral("width"); scale *= jointeventinputs.at(i)->Integral("width"); scale /= double(jointindexhigh[i] - jointindexlow[i]); jointindexscale.push_back(scale); } fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str()); fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str()); // Setup Max Events if (fMaxEvents > 1 && fMaxEvents < fNEvents) { if (LOG_LEVEL(SAM)) { std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl; } fNEvents = fMaxEvents; } // Print out Status if (LOG_LEVEL(SAM)) { std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl << "\t\t|-> Event Integral : " << fEventHist->Integral("width") * 1.E-38 << " events/nucleon" << std::endl << "\t\t|-> Flux Integral : " << fFluxHist->Integral("width") << " /cm2" << std::endl << "\t\t|-> Event/Flux : " << fEventHist->Integral("width") * 1.E-38 / fFluxHist->Integral("width") << " cm2/nucleon" << std::endl; } } BaseFitEvt* InputHandlerBase::GetBaseEvent(const UInt_t entry) { return static_cast(GetNuisanceEvent(entry, true)); } double InputHandlerBase::GetInputWeight(int entry) { if (!jointinput) return 1.0; // Find Switch Scale while (entry < jointindexlow[jointindexswitch] || entry >= jointindexhigh[jointindexswitch]) { jointindexswitch++; // Loop Around if (jointindexswitch >= jointindexlow.size()) { jointindexswitch = 0; } } return jointindexscale[jointindexswitch]; }; diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx index 8e3d5a2..7df81d1 100644 --- a/src/InputHandler/NEUTInputHandler.cxx +++ b/src/InputHandler/NEUTInputHandler.cxx @@ -1,499 +1,503 @@ #ifdef __NEUT_ENABLED__ #include "NEUTInputHandler.h" #include "InputUtils.h" NEUTGeneratorInfo::~NEUTGeneratorInfo() { DeallocateParticleStack(); } void NEUTGeneratorInfo::AddBranchesToTree(TTree* tn) { tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I"); tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode, "NEUTParticleStatusCode[NEUTParticleN]/I"); tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode, "NEUTParticleAliveCode[NEUTParticleN]/I"); } void NEUTGeneratorInfo::SetBranchesFromTree(TTree* tn) { tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN); tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode); tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode); } void NEUTGeneratorInfo::AllocateParticleStack(int stacksize) { fNEUTParticleN = 0; fNEUTParticleStatusCode = new int[stacksize]; fNEUTParticleStatusCode = new int[stacksize]; } void NEUTGeneratorInfo::DeallocateParticleStack() { delete fNEUTParticleStatusCode; delete fNEUTParticleAliveCode; } void NEUTGeneratorInfo::FillGeneratorInfo(NeutVect* nevent) { Reset(); for (int i = 0; i < nevent->Npart(); i++) { fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus; fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive; fNEUTParticleN++; } } void NEUTGeneratorInfo::Reset() { for (int i = 0; i < fNEUTParticleN; i++) { fNEUTParticleStatusCode[i] = -1; fNEUTParticleAliveCode[i] = 9; } fNEUTParticleN = 0; } NEUTInputHandler::NEUTInputHandler(std::string const& handle, std::string const& rawinputs) { LOG(SAM) << "Creating NEUTInputHandler : " << handle << std::endl; // Run a joint input handling fName = handle; // Setup the TChain fNEUTTree = new TChain("neuttree"); fSaveExtra = FitPar::Config().GetParB("SaveExtraNEUT"); fCacheSize = FitPar::Config().GetParI("CacheSize"); fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); // Loop over all inputs and grab flux, eventhist, and nevents std::vector inputs = InputUtils::ParseInputFileList(rawinputs); for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) { // Open File for histogram access TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ"); if (!inp_file or inp_file->IsZombie()) { THROW("NEUT File IsZombie() at : '" << inputs[inp_it] << "'" << std::endl << "Check that your file paths are correct and the file exists!" << std::endl << "$ ls -lh " << inputs[inp_it]); } // Get Flux/Event hist TH1D* fluxhist = (TH1D*)inp_file->Get( (PlotUtils::GetObjectWithName(inp_file, "flux")).c_str()); TH1D* eventhist = (TH1D*)inp_file->Get( (PlotUtils::GetObjectWithName(inp_file, "evt")).c_str()); if (!fluxhist or !eventhist) { ERROR(FTL, "Input File Contents: " << inputs[inp_it]); inp_file->ls(); THROW( "NEUT FILE doesn't contain flux/xsec info. You may have to " "regenerate your MC!"); } // Get N Events TTree* neuttree = (TTree*)inp_file->Get("neuttree"); if (!neuttree) { ERROR(FTL, "neuttree not located in NEUT file: " << inputs[inp_it]); THROW("Check your inputs, they may need to be completely regenerated!"); throw; } int nevents = neuttree->GetEntries(); if (nevents <= 0) { THROW("Trying to a TTree with " << nevents << " to TChain from : " << inputs[inp_it]); } // Register input to form flux/event rate hists RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist); // Add To TChain fNEUTTree->AddFile(inputs[inp_it].c_str()); } // Registor all our file inputs SetupJointInputs(); // Assign to tree fEventType = kNEUT; fNeutVect = NULL; fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect); fNEUTTree->GetEntry(0); // Create Fit Event fNUISANCEEvent = new FitEvent(); fNUISANCEEvent->SetNeutVect(fNeutVect); if (fSaveExtra) { fNeutInfo = new NEUTGeneratorInfo(); fNUISANCEEvent->AddGeneratorInfo(fNeutInfo); } fNUISANCEEvent->HardReset(); }; NEUTInputHandler::~NEUTInputHandler(){ // if (fNEUTTree) delete fNEUTTree; // if (fNeutVect) delete fNeutVect; // if (fNeutInfo) delete fNeutInfo; }; void NEUTInputHandler::CreateCache() { if (fCacheSize > 0) { // fNEUTTree->SetCacheEntryRange(0, fNEvents); fNEUTTree->AddBranchToCache("vectorbranch", 1); fNEUTTree->SetCacheSize(fCacheSize); } } void NEUTInputHandler::RemoveCache() { // fNEUTTree->SetCacheEntryRange(0, fNEvents); fNEUTTree->AddBranchToCache("vectorbranch", 0); fNEUTTree->SetCacheSize(0); } FitEvent* NEUTInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) { // Catch too large entries if (entry >= (UInt_t)fNEvents) return NULL; // Read Entry from TTree to fill NEUT Vect in BaseFitEvt; fNEUTTree->GetEntry(entry); // Run NUISANCE Vector Filler if (!lightweight) { CalcNUISANCEKinematics(); } #ifdef __PROB3PP_ENABLED__ else { UInt_t npart = fNeutVect->Npart(); for (size_t i = 0; i < npart; i++) { NeutPart* part = fNUISANCEEvent->fNeutVect->PartInfo(i); if ((part->fIsAlive == false) && (part->fStatus == -1) && std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4, part->fPID)) { fNUISANCEEvent->probe_E = part->fP.T(); fNUISANCEEvent->probe_pdg = part->fPID; break; } else { continue; } } } #endif // Setup Input scaling for joint inputs fNUISANCEEvent->InputWeight = GetInputWeight(entry); // Return event pointer return fNUISANCEEvent; } // From NEUT neutclass/neutpart.h // Bool_t fIsAlive; // Particle should be tracked or not // ( in the detector simulator ) // // Int_t fStatus; // Status flag of this particle // -2: Non existing particle // -1: Initial state particle // 0: Normal // 1: Decayed to the other particle // 2: Escaped from the detector // 3: Absorped // 4: Charge exchanged // 5: Pauli blocked // 6: N/A // 7: Produced child particles // 8: Inelastically scattered // int NEUTInputHandler::GetNeutParticleStatus(NeutPart* part) { // State int state = kUndefinedState; // Remove Pauli blocked events, probably just single pion events if (part->fStatus == 5) { state = kFSIState; // fStatus == -1 means initial state } else if (part->fIsAlive == false && part->fStatus == -1) { state = kInitialState; // NEUT has a bit of a strange convention for fIsAlive and fStatus // combinations // for NC and neutrino particle isAlive true/false and status 2 means // final state particle // for other particles in NC status 2 means it's an FSI particle // for CC it means it was an FSI particle } else if (part->fStatus == 2) { // NC case is a little strange... The outgoing neutrino might be alive or // not alive. Remaining particles with status 2 are FSI particles that // reinteracted if (abs(fNeutVect->Mode) > 30 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) { state = kFinalState; // The usual CC case } else if (part->fIsAlive == true) { state = kFSIState; } } else if (part->fIsAlive == true && part->fStatus == 2 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) { state = kFinalState; } else if (part->fIsAlive == true && part->fStatus == 0) { state = kFinalState; } else if (!part->fIsAlive && (part->fStatus == 1 || part->fStatus == 3 || part->fStatus == 4 || part->fStatus == 7 || part->fStatus == 8)) { state = kFSIState; // There's one hyper weird case where fStatus = -3. This apparently corresponds to a nucleon being ejected via pion FSI when there is "data available" } else if (!part->fIsAlive && (part->fStatus == -3)) { state = kUndefinedState; // NC neutrino outgoing } else if (!part->fIsAlive && part->fStatus == 0 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) { state = kFinalState; // Warn if we still find alive particles without classifying them } else if (part->fIsAlive == true) { ERR(WRN) << "Undefined NEUT state " << " Alive: " << part->fIsAlive << " Status: " << part->fStatus << " PDG: " << part->fPID << std::endl; throw; // Warn if we find dead particles that we haven't classified } else { ERR(WRN) << "Undefined NEUT state " << " Alive: " << part->fIsAlive << " Status: " << part->fStatus << " PDG: " << part->fPID << std::endl; throw; } return state; } void NEUTInputHandler::CalcNUISANCEKinematics() { // Reset all variables fNUISANCEEvent->ResetEvent(); // Fill Globals fNUISANCEEvent->Mode = fNeutVect->Mode; fNUISANCEEvent->fEventNo = fNeutVect->EventNo; fNUISANCEEvent->fTargetA = fNeutVect->TargetA; fNUISANCEEvent->fTargetZ = fNeutVect->TargetZ; fNUISANCEEvent->fTargetH = fNeutVect->TargetH; fNUISANCEEvent->fBound = bool(fNeutVect->Ibound); if (fNUISANCEEvent->fBound) { fNUISANCEEvent->fTargetPDG = TargetUtils::GetTargetPDGFromZA( fNUISANCEEvent->fTargetZ, fNUISANCEEvent->fTargetA); } else { fNUISANCEEvent->fTargetPDG = 1000010010; } // Check Particle Stack UInt_t npart = fNeutVect->Npart(); UInt_t kmax = fNUISANCEEvent->kMaxParticles; if (npart > kmax) { ERR(FTL) << "NEUT has too many particles. Expanding stack." << std::endl; fNUISANCEEvent->ExpandParticleStack(npart); throw; } // Fill Particle Stack for (size_t i = 0; i < npart; i++) { // Get Current Count int curpart = fNUISANCEEvent->fNParticles; // Get NEUT Particle NeutPart* part = fNeutVect->PartInfo(i); // State int state = GetNeutParticleStatus(part); // Remove Undefined if (kRemoveUndefParticles && state == kUndefinedState) continue; // Remove FSI if (kRemoveFSIParticles && state == kFSIState) continue; // Remove Nuclear if (kRemoveNuclearParticles && (state == kNuclearInitial || state == kNuclearRemnant)) continue; // State fNUISANCEEvent->fParticleState[curpart] = state; // Mom fNUISANCEEvent->fParticleMom[curpart][0] = part->fP.X(); fNUISANCEEvent->fParticleMom[curpart][1] = part->fP.Y(); fNUISANCEEvent->fParticleMom[curpart][2] = part->fP.Z(); fNUISANCEEvent->fParticleMom[curpart][3] = part->fP.T(); // PDG fNUISANCEEvent->fParticlePDG[curpart] = part->fPID; // Add up particle count fNUISANCEEvent->fNParticles++; } // Save Extra Generator Info if (fSaveExtra) { fNeutInfo->FillGeneratorInfo(fNeutVect); } // Run Initial, FSI, Final, Other ordering. fNUISANCEEvent->OrderStack(); FitParticle* ISNeutralLepton = fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos); if (ISNeutralLepton) { fNUISANCEEvent->probe_E = ISNeutralLepton->E(); fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG(); } return; } void NEUTUtils::FillNeutCommons(NeutVect* nvect) { // WARNING: This has only been implemented for a neuttree and not GENIE // This should be kept in sync with T2KNIWGUtils::GetNIWGEvent(TTree) // NEUT version info. Can't get it to compile properly with this yet // neutversion_.corev = nvect->COREVer; // neutversion_.nucev = nvect->NUCEVer; // neutversion_.nuccv = nvect->NUCCVer; // Documentation: See nework.h nework_.modene = nvect->Mode; nework_.numne = nvect->Npart(); +#ifdef NEUT_COMMON_QEAV + nemdls_.mdlqeaf = nvect->QEAVForm; +#else nemdls_.mdlqeaf = nvect->QEVForm; +#endif nemdls_.mdlqe = nvect->QEModel; nemdls_.mdlspi = nvect->SPIModel; nemdls_.mdldis = nvect->DISModel; nemdls_.mdlcoh = nvect->COHModel; neutcoh_.necohepi = nvect->COHModel; nemdls_.xmaqe = nvect->QEMA; nemdls_.xmvqe = nvect->QEMV; nemdls_.kapp = nvect->KAPPA; // nemdls_.sccfv = SCCFVdef; // nemdls_.sccfa = SCCFAdef; // nemdls_.fpqe = FPQEdef; nemdls_.xmaspi = nvect->SPIMA; nemdls_.xmvspi = nvect->SPIMV; nemdls_.xmares = nvect->RESMA; nemdls_.xmvres = nvect->RESMV; neut1pi_.xmanffres = nvect->SPIMA; neut1pi_.xmvnffres = nvect->SPIMV; neut1pi_.xmarsres = nvect->RESMA; neut1pi_.xmvrsres = nvect->RESMV; neut1pi_.neiff = nvect->SPIForm; neut1pi_.nenrtype = nvect->SPINRType; neut1pi_.rneca5i = nvect->SPICA5I; neut1pi_.rnebgscl = nvect->SPIBGScale; nemdls_.xmacoh = nvect->COHMA; nemdls_.rad0nu = nvect->COHR0; // nemdls_.fa1coh = nvect->COHA1err; // nemdls_.fb1coh = nvect->COHb1err; // neutdis_.nepdf = NEPDFdef; // neutdis_.nebodek = NEBODEKdef; neutcard_.nefrmflg = nvect->FrmFlg; neutcard_.nepauflg = nvect->PauFlg; neutcard_.nenefo16 = nvect->NefO16; neutcard_.nemodflg = nvect->ModFlg; // neutcard_.nenefmodl = 1; // neutcard_.nenefmodh = 1; // neutcard_.nenefkinh = 1; // neutpiabs_.neabspiemit = 1; nenupr_.iformlen = nvect->FormLen; neutpiless_.ipilessdcy = nvect->IPilessDcy; neutpiless_.rpilessdcy = nvect->RPilessDcy; neutpiless_.ipilessdcy = nvect->IPilessDcy; neutpiless_.rpilessdcy = nvect->RPilessDcy; neffpr_.fefqe = nvect->NuceffFactorPIQE; neffpr_.fefqeh = nvect->NuceffFactorPIQEH; neffpr_.fefinel = nvect->NuceffFactorPIInel; neffpr_.fefabs = nvect->NuceffFactorPIAbs; neffpr_.fefcx = nvect->NuceffFactorPICX; neffpr_.fefcxh = nvect->NuceffFactorPICXH; neffpr_.fefcoh = nvect->NuceffFactorPICoh; neffpr_.fefqehf = nvect->NuceffFactorPIQEHKin; neffpr_.fefcxhf = nvect->NuceffFactorPICXKin; neffpr_.fefcohf = nvect->NuceffFactorPIQELKin; for (int i = 0; i < nework_.numne; i++) { nework_.ipne[i] = nvect->PartInfo(i)->fPID; nework_.pne[i][0] = (float)nvect->PartInfo(i)->fP.X() / 1000; // VC(NE)WORK in M(G)eV nework_.pne[i][1] = (float)nvect->PartInfo(i)->fP.Y() / 1000; // VC(NE)WORK in M(G)eV nework_.pne[i][2] = (float)nvect->PartInfo(i)->fP.Z() / 1000; // VC(NE)WORK in M(G)eV } // fsihist.h // neutroot fills a dummy object for events with no FSI to prevent memory leak // when // reading the TTree, so check for it here if ((int)nvect->NfsiVert() == 1) { // An event with FSI must have at least two vertices // if (nvect->NfsiPart()!=1 || nvect->Fsiprob!=-1) // ERR(WRN) << "T2KNeutUtils::fill_neut_commons(TTree) NfsiPart!=1 or // Fsiprob!=-1 when NfsiVert==1" << std::endl; fsihist_.nvert = 0; fsihist_.nvcvert = 0; fsihist_.fsiprob = 1; } else { // Real FSI event fsihist_.nvert = (int)nvect->NfsiVert(); for (int ivert = 0; ivert < fsihist_.nvert; ivert++) { fsihist_.iflgvert[ivert] = nvect->FsiVertInfo(ivert)->fVertID; fsihist_.posvert[ivert][0] = (float)nvect->FsiVertInfo(ivert)->fPos.X(); fsihist_.posvert[ivert][1] = (float)nvect->FsiVertInfo(ivert)->fPos.Y(); fsihist_.posvert[ivert][2] = (float)nvect->FsiVertInfo(ivert)->fPos.Z(); } fsihist_.nvcvert = nvect->NfsiPart(); for (int ip = 0; ip < fsihist_.nvcvert; ip++) { fsihist_.abspvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomLab; fsihist_.abstpvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomNuc; fsihist_.ipvert[ip] = nvect->FsiPartInfo(ip)->fPID; fsihist_.iverti[ip] = nvect->FsiPartInfo(ip)->fVertStart; fsihist_.ivertf[ip] = nvect->FsiPartInfo(ip)->fVertEnd; fsihist_.dirvert[ip][0] = (float)nvect->FsiPartInfo(ip)->fDir.X(); fsihist_.dirvert[ip][1] = (float)nvect->FsiPartInfo(ip)->fDir.Y(); fsihist_.dirvert[ip][2] = (float)nvect->FsiPartInfo(ip)->fDir.Z(); } fsihist_.fsiprob = nvect->Fsiprob; } neutcrscom_.crsx = nvect->Crsx; neutcrscom_.crsy = nvect->Crsy; neutcrscom_.crsz = nvect->Crsz; neutcrscom_.crsphi = nvect->Crsphi; neutcrscom_.crsq2 = nvect->Crsq2; neuttarget_.numbndn = nvect->TargetA - nvect->TargetZ; neuttarget_.numbndp = nvect->TargetZ; neuttarget_.numfrep = nvect->TargetH; neuttarget_.numatom = nvect->TargetA; posinnuc_.ibound = nvect->Ibound; // put empty nucleon FSI history (since it is not saved in the NeutVect // format) // Comment out as NEUT does not have the necessary proton FSI information yet // nucleonfsihist_.nfnvert = 0; // nucleonfsihist_.nfnstep = 0; } #endif diff --git a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx index b29e4d2..ab9d511 100644 --- a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx +++ b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx @@ -1,327 +1,335 @@ // 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 . *******************************************************************************/ // Implementation of 2018 MINERvA numu CC0pi STV analysis // arxiv 1805.05486.pdf // Clarence Wret // cwret@fnal.gov // Stephen Dolan // Stephen.Dolan@llr.in2p3.fr #include "MINERvA_SignalDef.h" #include "MINERvA_CC0pinp_STV_XSec_1D_nu.h" //******************************************************************** void MINERvA_CC0pinp_STV_XSec_1D_nu::SetupDataSettings() { //******************************************************************** // Set Distribution // See header file for enum and some descriptions std::string name = fSettings.GetS("name"); if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpmu_nu")) fDist= kMuonMom; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthmu_nu")) fDist= kMuonTh; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpprot_nu")) fDist= kPrMom; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthprot_nu")) fDist= kPrTh; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpnreco_nu")) fDist= kNeMom; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddalphat_nu")) fDist= kDalphaT; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddpt_nu")) fDist= kDpT; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddphit_nu")) fDist= kDphiT; // Location of data, correlation matrices and the titles std::string titles = "MINERvA_CC0pinp_STV_XSec_1D"; std::string foldername; std::string distdescript; // Data release is a single file std::string rootfile = "MINERvA_1805.05486.root"; fMin = -999; fMax = 999; switch (fDist) { case (kMuonMom): titles += "pmu"; foldername = "muonmomentum"; distdescript = "Muon momentum in lab frame"; + /* fMin = 2.0; fMax = 6.0; + */ + fMin = 1.5; + fMax = 10.0; break; case (kMuonTh): titles += "thmu"; foldername = "muontheta"; distdescript = "Muon angle relative neutrino in lab frame"; fMin = 0.0; fMax = 20.0; break; case (kPrMom): titles += "pprot"; foldername = "protonmomentum"; distdescript = "Proton momentum in lab frame"; - fMin = 0.5; + //fMin = 0.5; + fMin = 0.45; fMax = 1.2; break; case (kPrTh): titles += "thprot"; foldername = "protontheta"; distdescript = "Proton angle relative neutrino in lab frame"; fMin = 0.0; fMax = 70.0; break; case (kNeMom): titles += "pnreco"; foldername = "neutronmomentum"; distdescript = "Neutron momentum in lab frame"; fMin = 0.0; - fMax = 0.9; + //fMax = 0.9; + fMax = 2.0; break; case (kDalphaT): foldername = "dalphat"; titles += foldername; distdescript = "Delta Alpha_T"; fMin = 0.0; - fMax = 170; + //fMax = 170; + fMax = 180; break; case (kDpT): foldername = "dpt"; titles += foldername; distdescript = "Delta p_T"; fMin = 0.0; - fMax = 170; + fMax = 2.0; break; case (kDphiT): foldername = "dphit"; titles += foldername; distdescript = "Delta phi_T"; fMin = 0.0; - fMax = 60.0; + //fMax = 60.0; + fMax = 180.0; break; default: ERR(FTL) << "Did not find your specified distribution implemented, exiting" << std::endl; ERR(FTL) << "You gave " << fName << std::endl; throw; } titles += "_nu"; // All have the same name std::string dataname = foldername; // Sample overview --------------------------------------------------- std::string descrip = distdescript + \ "Target: CH \n" \ "Flux: MINERvA Forward Horn Current numu ONLY \n" \ "Signal: Any event with 1 muon, and 0pi0 in FS, no mesons, at least one proton with: \n" \ "1.5GeV < p_mu < 10 GeV\n"\ "theta_mu < 20 degrees\n"\ "0.45GeV < p_prot < 1.2 GeV\n"\ "theta_prot < 70 degrees\n" \ "arXiv 1805.05486"; fSettings.SetDescription(descrip); std::string filename = GeneralUtils::GetTopLevelDir()+"/data/MINERvA/CC0pi/CC0pi_STV/MINERvA_1805.05486.root"; // Specify the data fSettings.SetDataInput(filename+";"+dataname); // And the correlations fSettings.SetCovarInput(filename+";"+dataname); // Set titles fSettings.SetTitle(titles); }; //******************************************************************** MINERvA_CC0pinp_STV_XSec_1D_nu::MINERvA_CC0pinp_STV_XSec_1D_nu(nuiskey samplekey) { //******************************************************************** // A few different distributinos // Muon momentum, muon angle, proton momentum, proton angle, neutron momentum, dalphat, dpt, dphit // Setup common settings fSettings = LoadSampleSettings(samplekey); // Load up the data paths and sample descriptions SetupDataSettings(); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); // No Enu cut fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); // Finalise the settings FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / (double(fNEvents) * TotalIntegratedFlux("width")); // Set the data and covariance matrix SetDataFromRootFile( fSettings.GetDataInput() ); SetCovarianceFromRootFile(fSettings.GetCovarInput() ); fSettings.SetXTitle(fDataHist->GetXaxis()->GetTitle()); fSettings.SetYTitle(fDataHist->GetYaxis()->GetTitle()); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** // Data comes in a TList // Some bins need stripping out because of zero bin content. Why oh why void MINERvA_CC0pinp_STV_XSec_1D_nu::SetDataFromRootFile(std::string filename) { //******************************************************************** std::vector tempfile = GeneralUtils::ParseToStr(filename, ";"); TFile *File = new TFile(tempfile[0].c_str(), "READ"); // First object is the data TH1D* temp = (TH1D*)(((TList*)(File->Get(tempfile[1].c_str())))->At(0)); // Garh, some of the data points are zero in the TH1D (WHY?!) so messes with the covariance entries to data bins check // Skim through the data and check for zero bins std::vector CrossSection; std::vector Error; std::vector BinEdges; int lastbin = 0; startbin = 0; for (int i = 0; i < temp->GetXaxis()->GetNbins()+2; ++i) { - if (temp->GetBinContent(i+1) > 0 && temp->GetBinLowEdge(i+1) > fMin && temp->GetBinLowEdge(i+1) < fMax) { + if (temp->GetBinContent(i+1) > 0 && temp->GetBinLowEdge(i+1) >= fMin && temp->GetBinLowEdge(i+1) <= fMax) { if (startbin == 0) startbin = i; lastbin = i; CrossSection.push_back(temp->GetBinContent(i+1)); BinEdges.push_back(temp->GetXaxis()->GetBinLowEdge(i+1)); Error.push_back(temp->GetBinError(i+1)); } } BinEdges.push_back(temp->GetXaxis()->GetBinLowEdge(lastbin+2)); fDataHist = new TH1D((fSettings.GetName()+"_data").c_str(), (fSettings.GetFullTitles()).c_str(), BinEdges.size()-1, &BinEdges[0]); fDataHist->SetDirectory(0); for (unsigned int i = 0; i < BinEdges.size()-1; ++i) { fDataHist->SetBinContent(i+1, CrossSection[i]); fDataHist->SetBinError(i+1, Error[i]); } fDataHist->GetXaxis()->SetTitle(temp->GetXaxis()->GetTitle()); fDataHist->GetYaxis()->SetTitle(temp->GetYaxis()->GetTitle()); fDataHist->SetTitle(temp->GetTitle()); File->Close(); } //******************************************************************** // Covariance also needs stripping out // There's padding (two bins...) and overflow (last bin before the two empty bins) void MINERvA_CC0pinp_STV_XSec_1D_nu::SetCovarianceFromRootFile(std::string filename) { //******************************************************************** std::vector tempfile = GeneralUtils::ParseToStr(filename, ";"); TFile *File = new TFile(tempfile[0].c_str(), "READ"); // First object is the data, second is data with statistical error only, third is the covariance matrix TMatrixDSym *tempcov = (TMatrixDSym*)((TList*)File->Get(tempfile[1].c_str()))->At(2); // Count the number of zero entries int ngood = 0; int nstart = -1; int nend = -1; // Scan through the middle bin and look for entries int middle = tempcov->GetNrows()/2; int nbinsdata = fDataHist->GetXaxis()->GetNbins(); for (int j = 0; j < tempcov->GetNrows(); ++j) { if ((*tempcov)(middle,j) > 0 && ngood < nbinsdata) { ngood++; if (nstart == -1) nstart = j; if (j > nend) nend = j; } } fFullCovar = new TMatrixDSym(ngood); for (int i = 0; i < fFullCovar->GetNrows(); ++i) { for (int j = 0; j < fFullCovar->GetNrows(); ++j) { (*fFullCovar)(i,j) = (*tempcov)(i+nstart+startbin-1, j+nstart+startbin-1); } } (*fFullCovar) *= 1E38*1E38; File->Close(); // Fill other covars. covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } void MINERvA_CC0pinp_STV_XSec_1D_nu::FillEventVariables(FitEvent *event) { fXVar = -999.99; // Need a proton and a muon if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(13) == 0) { return; } TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 int HMFSProton = 0; double HighestMomentum = 0.0; // Get the stack of protons std::vector Protons = event->GetAllFSProton(); for (size_t i = 0; i < Protons.size(); ++i) { if (Protons[i]->p() > 450 && Protons[i]->p() < 1200 && Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && Protons[i]->p() > HighestMomentum) { HMFSProton = i; } } // Now get the proton TLorentzVector Pprot = Protons[HMFSProton]->fP; switch (fDist) { case (kMuonMom): fXVar = Pmu.Vect().Mag()/1000.0; break; case (kMuonTh): fXVar = Pmu.Vect().Angle(Pnu.Vect()) * (180.0/M_PI); break; case (kPrMom): fXVar = Pprot.Vect().Mag()/1000.0; break; case (kPrTh): fXVar = Pprot.Vect().Angle(Pnu.Vect()) * (180.0/M_PI); break; // Use Stephen's validated functions case (kNeMom): fXVar = FitUtils::Get_pn_reco_C(event, 14, true); break; case (kDalphaT): fXVar = FitUtils::Get_STV_dalphat(event, 14, true)*(180.0/M_PI); break; case (kDpT): fXVar = FitUtils::Get_STV_dpt(event, 14, true)/1000.0; break; case (kDphiT): fXVar = FitUtils::Get_STV_dphit(event, 14, true)*(180.0/M_PI); break; } return; }; //******************************************************************** bool MINERvA_CC0pinp_STV_XSec_1D_nu::isSignal(FitEvent *event) //******************************************************************** { return SignalDef::isCC0piNp_MINERvA_STV(event, EnuMin, EnuMax); } diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx index 3587b06..52b433f 100644 --- a/src/Reweight/GENIEWeightEngine.cxx +++ b/src/Reweight/GENIEWeightEngine.cxx @@ -1,251 +1,267 @@ #include "GENIEWeightEngine.h" +#ifdef __GENIE_EMP_MECRW_ENABLED +#include "ReWeight/GReWeightXSecEmpiricalMEC.h" +#endif + GENIEWeightEngine::GENIEWeightEngine(std::string name) { #ifdef __GENIE_ENABLED__ - // Setup the NEUT Reweight engien fCalcName = name; LOG(FIT) << "Setting up GENIE RW : " << fCalcName << std::endl; // Create RW Engine suppressing cout StopTalking(); fGenieRW = new genie::rew::GReWeight(); // Get List of Vetos (Just for debugging) - std::string rw_engine_list = FitPar::Config().GetParS("FitWeight_fGenieRW_veto"); + std::string rw_engine_list = + FitPar::Config().GetParS("FitWeight_fGenieRW_veto"); bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos; bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos; bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos; bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos; bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos; - bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos; + bool xsec_resdec = + rw_engine_list.find("hadro_res_decay") == std::string::npos; bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos; bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos; bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos; bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos; bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos; bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos; bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos; bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos; bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos; - bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos; + bool xsec_qeaxial = + rw_engine_list.find("xsec_ccqe_axial") == std::string::npos; +#ifdef __GENIE_EMP_MECRW_ENABLED + bool xsec_empMEC = rw_engine_list.find("xsec_empMEC") == std::string::npos; +#endif // Now actually add the RW Calcs if (xsec_ncel) fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL); - if (xsec_ccqe){ + if (xsec_ccqe) { fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE); // (dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe"))) // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") ); } - if (xsec_coh){ +#ifdef __GENIE_EMP_MECRW_ENABLED + if (xsec_empMEC) { + fGenieRW->AdoptWghtCalc("xsec_empMEC", + new genie::rew::GReWeightXSecEmpiricalMEC); + } +#endif + if (xsec_coh) { fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH()); // (dynamic_cast (fGenieRW->WghtCalc("xsec_coh"))) // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") ); } if (xsec_nnres) fGenieRW->AdoptWghtCalc("xsec_nonresbkg", - new genie::rew::GReWeightNonResonanceBkg); + new genie::rew::GReWeightNonResonanceBkg); if (xsec_nudis) fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod); if (xsec_resdec) fGenieRW->AdoptWghtCalc("hadro_res_decay", - new genie::rew::GReWeightResonanceDecay); + new genie::rew::GReWeightResonanceDecay); if (xsec_fzone) fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone); if (xsec_intra) fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke); if (xsec_agky) fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY); if (xsec_qevec) fGenieRW->AdoptWghtCalc("xsec_ccqe_vec", - new genie::rew::GReWeightNuXSecCCQEvec); + new genie::rew::GReWeightNuXSecCCQEvec); #if __GENIE_VERSION__ >= 212 if (xsec_qeaxial) fGenieRW->AdoptWghtCalc("xsec_ccqe_axial", - new genie::rew::GReWeightNuXSecCCQEaxial); + new genie::rew::GReWeightNuXSecCCQEaxial); #endif if (xsec_dis) fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS); if (xsec_nc) fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC); - if (xsec_ccres){ + if (xsec_ccres) { #if __GENIE_VERSION__ < 213 fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES); #else - fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES( FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default")); + fGenieRW->AdoptWghtCalc( + "xsec_ccres", + new genie::rew::GReWeightNuXSecCCRES( + FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default")); #endif } if (xsec_ncres) fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES); if (xsec_nucqe) fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM); - GReWeightNuXSecCCQE * rwccqe = - dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe")); - rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa); - - // Default to include shape and normalization changes for CCRES (can be changed downstream if desired) - GReWeightNuXSecCCRES * rwccres = dynamic_cast (fGenieRW->WghtCalc("xsec_ccres")); - std::string marestype = FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode"); - if (!marestype.compare("kModeNormAndMaMvShape")){ - LOG(FIT) << "Setting GENIE ReWeight CCRES to kModeNormAndMaMvShape" << std::endl; - rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); - } else if (!marestype.compare("kModeMaMv")) { - LOG(FIT) << "Setting GENIE ReWeight CCRES to kModeMaMv" << std::endl; - rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); - } else { - THROW("Unkown MARES Mode in GENIE Weight Engine : " << marestype ); + if (xsec_ccqe) { + GReWeightNuXSecCCQE *rwccqe = + dynamic_cast(fGenieRW->WghtCalc("xsec_ccqe")); + rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa); } - // Default to include shape and normalization changes for NCRES (can be changed downstream if desired) - GReWeightNuXSecNCRES * rwncres = dynamic_cast (fGenieRW->WghtCalc("xsec_ncres")); - rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv); + if (xsec_ccres) { + // Default to include shape and normalization changes for CCRES (can be + // changed downstream if desired) + GReWeightNuXSecCCRES *rwccres = + dynamic_cast(fGenieRW->WghtCalc("xsec_ccres")); + + std::string marestype = + FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode"); + if (!marestype.compare("kModeNormAndMaMvShape")) { + rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); + } else if (!marestype.compare("kModeMaMv")) { + rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); + } else { + THROW("Unkown MARES Mode in GENIE Weight Engine : " << marestype); + } + } + if (xsec_ncres) { + // Default to include shape and normalization changes for NCRES (can be + // changed downstream if desired) + GReWeightNuXSecNCRES *rwncres = + dynamic_cast(fGenieRW->WghtCalc("xsec_ncres")); + rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv); + } - // Default to include shape and normalization changes for DIS (can be changed downstream if desired) - GReWeightNuXSecDIS * rwdis = dynamic_cast (fGenieRW->WghtCalc("xsec_dis")); - rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u); + if (xsec_dis) { + // Default to include shape and normalization changes for DIS (can be + // changed downstream if desired) + GReWeightNuXSecDIS *rwdis = + dynamic_cast(fGenieRW->WghtCalc("xsec_dis")); + rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u); - // Set Abs Twk Config - fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); - if (fIsAbsTwk) { - LOG(FIT) << "GENIE ReWeight values are ABSOLUTE tweak dials" << std::endl; - } else { - LOG(FIT) << "GENIE ReWeight values are RELATIVE tweak dials" << std::endl; + // Set Abs Twk Config + fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); } // allow cout again StartTalking(); #else ERR(FTL) << "GENIE RW NOT ENABLED" << std::endl; #endif }; - void GENIEWeightEngine::IncludeDial(std::string name, double startval) { #ifdef __GENIE_ENABLED__ // Get First enum int nuisenum = Reweight::ConvDial(name, kGENIE); // Setup Maps - fEnumIndex[nuisenum];// = std::vector(0); - fNameIndex[name]; // = std::vector(0); + fEnumIndex[nuisenum]; // = std::vector(0); + fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get RW genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fGENIESysts.push_back(rwsyst); // Initialize dial - LOG(FIT) << "Registering " << singlename << " from " << name << std::endl; - fGenieRW->Systematics().Init( fGENIESysts[index] ); + std::cout << "Registering " << singlename << " from " << name << std::endl; + fGenieRW->Systematics().Init(fGENIESysts[index]); // If Absolute if (fIsAbsTwk) { - GSystUncertainty::Instance()->SetUncertainty( rwsyst, 1.0, 1.0 ); + GSystUncertainty::Instance()->SetUncertainty(rwsyst, 1.0, 1.0); } // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); - } // Set Value if given if (startval != -999.9) { SetDialValue(nuisenum, startval); } #endif }; - void GENIEWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef __GENIE_ENABLED__ std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; - fGenieRW->Systematics().Set( fGENIESysts[indices[i]], val); + fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val); } #endif } void GENIEWeightEngine::SetDialValue(std::string name, double val) { #ifdef __GENIE_ENABLED__ std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val); } #endif } void GENIEWeightEngine::Reconfigure(bool silent) { #ifdef __GENIE_ENABLED__ // Hush now... - LOG(MIN) << "Reconfiguring" << std::endl; - if (silent) StopTalking(); + if (silent) + StopTalking(); // Reconf fGenieRW->Reconfigure(); fGenieRW->Print(); // Shout again - if (silent) StartTalking(); + if (silent) + StartTalking(); #endif } - -double GENIEWeightEngine::CalcWeight(BaseFitEvt* evt) { +double GENIEWeightEngine::CalcWeight(BaseFitEvt *evt) { double rw_weight = 1.0; #ifdef __GENIE_ENABLED__ // Skip Non GENIE - if (evt->fType != kGENIE) return 1.0; + if (evt->fType != kGENIE) + return 1.0; // Make nom weight if (!evt) { THROW("evt not found : " << evt); } if (!(evt->genie_event)) { THROW("evt->genie_event not found!" << evt->genie_event); } if (!(evt->genie_event->event)) { - THROW("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event)); + THROW("evt->genie_event->event GHepRecord not found!" + << (evt->genie_event->event)); } if (!fGenieRW) { THROW("GENIE RW Not Found!" << fGenieRW); } rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event)); - // std::cout << "Returning GENIE Weight for electron scattering = " << rw_weight << std::endl; + // std::cout << "Returning GENIE Weight for electron scattering = " << + // rw_weight << std::endl; #endif // Return rw_weight return rw_weight; } - - - - - - - - - - diff --git a/src/Reweight/weightRPA.h b/src/Reweight/weightRPA.h index d0c564d..ea875df 100644 --- a/src/Reweight/weightRPA.h +++ b/src/Reweight/weightRPA.h @@ -1,418 +1,417 @@ -#ifndef weightRPA_h -#define weightRPA_h - -#include //ifstream -#include //cout - -#include -#include -#include -#include -#include "math.h" -#include "assert.h" -//For Compatibility with ROOT compiler -//uncomment the following: -//#define ROOT - -/*! - * Code example to extract the RPA effect central weight - * and its uncertainties from the prepared files. - * Heidi Schellman (Oregon State) and Rik Gran (Minnesota Duluth) - * for use in MINERvA experiment analysis - * must compile with the ROOT libraries - * g++ `root-config --glibs --cflags` -O3 weightRPAtest.cxx -o weightRPA - - - * The underlying model is from the IFIC Valencia group - * see (public) minerva docdb:12688 for the full physics discussion - - */ - - -// NOTE UNITS ARE GEV in the calculation -// make sure you convert MeV to GeV before calling these functions - - -// Class for getting RPA paramaters inputs from a given file -// Includes methods that return all five RPA weights at once -// (is the most cpu efficient way to get them) -// Or return just the central value -// (skipping the uncertainties code completely if only cv wanted) -// Or return each CV and uncertainty one at a time -// (is the least cpu efficient, repeats some calculations 5 times) - -class weightRPA { -public: - //Constructor: Read in params from a filename - weightRPA(const TString f) { read(f); } //Read in params from file - - TString filename; - TFile* fRPAratio; - TH2D *hRPArelratio; - TH2D *hRPAnonrelratio; - TH1D *hQ2relratio; - TH1D *hQ2nonrelratio; - TArrayD *TADrpapolyrel; - Double_t *rpapolyrel ; - TArrayD *TADrpapolynonrel; - Double_t *rpapolynonrel; - static const int CENTRAL=0; - static const int LOWQ2 = 1; - static const int HIGHQ2 = 2; - - // true to take from histogram, false use parameterization - static const bool Q2histORparam = true; - - // MINERvA holds kinematics in MeV, but all these functions require GeV - // So make sure you pass them in GeV. - inline double getWeight(const double mc_q0, const double mc_q3, double * weights); //in GeV - inline double getWeight(const double mc_q0, const double mc_q3); //in GeV - inline double getWeight(const double mc_q0, const double mc_q3, int type, int sign); //in GeV - inline double getWeight(const double Q2); //in GeV^2 - inline double getWeightLowQ2(const double mc_q0, const double mc_q3, const int sign); - inline double getWeightHighQ2(const double mc_q0, const double mc_q3, const int sign); - inline double getWeightQ2(const double mc_Q2, const bool relORnonrel=true); - //Initializer - inline void read(const TString f); - - // q0 and q3 in GeV, type = 1 for low Q2, 2 for high Q2, 0 for central - //double getWeightInternal(const double mc_q0, const double mc_q3,int type, int sign); - - private: - inline double getWeightInternal(const double mc_Q2); - - inline double getWeightInternal(const double mc_q0, const double mc_q3, const int type, const int sign); - inline double getWeightInternal(const double mc_q0, const double mc_q3, double *weights=0); - inline double getWeightQ2parameterization(const double mc_Q2, const bool relORnonrel); - inline double getWeightQ2fromhistogram(const double mc_Q2, const bool relORnonrel); - - -}; - -void weightRPA::read(const TString f) -//Read in the params doubles from a file -//argument: valid filename -{ - fRPAratio = TFile::Open(f,"READONLY"); - if (fRPAratio){ - hRPArelratio = (TH2D*)fRPAratio->Get("hrelratio"); - hRPAnonrelratio = (TH2D*)fRPAratio->Get("hnonrelratio"); - hQ2relratio = (TH1D*)fRPAratio->Get("hQ2relratio"); - hQ2nonrelratio = (TH1D*)fRPAratio->Get("hQ2nonrelratio"); - TADrpapolyrel = (TArrayD*)fRPAratio->Get("rpapolyrel"); - rpapolyrel = TADrpapolyrel->GetArray(); - TADrpapolynonrel = (TArrayD*)fRPAratio->Get("rpapolynonrel"); - rpapolynonrel = TADrpapolynonrel->GetArray(); - hRPArelratio->Print(); - std::cout << "have read in ratios from file " << f <= gevlimit) q0bin = rpamevlimit - 1; - if(mc_q3 >= gevlimit) q3bin = rpamevlimit - 1; - - // Nieves does not calculate anything below binding energy. - // I don't know what GENIE does, but lets be soft about this. - // Two things lurking here at once. - // One, we are taking out a 10 MeV offset between GENIE and Valencia. - // Second, we are protecting against being asked for a weight that is too low in q0. - // It actually shouldn't happen for real GENIE events, - // but this protection does something that doesn't suck, just in case. - // you would see the artifact in a plot for sure, but better than writing 1.0. - - // CWret after talking to Rik in Nov 2018 at MINERvA CM - // 2.12.8 sets this to 27 because change of Q definition: need to offset GENIE and Nieves Eb even more -#if __GENIE_VERSION__ >= 210 - Int_t q0offsetValenciaGENIE = 27; -#else - Int_t q0offsetValenciaGENIE = 10; -#endif - - - if(mc_q0 < 0.018) q0bin = 18+q0offsetValenciaGENIE; - Double_t thisrwtemp = hRPArelratio->GetBinContent(q3bin,q0bin-q0offsetValenciaGENIE); - - // now trap bogus entries. Not sure why they happen, but set to 1.0 not 0.0 - if(thisrwtemp <= 0.001)thisrwtemp = 1.0; - - // events in genie but not in valencia should get a weight - // related to a similar q0 from the bulk distribution. - if(mc_q0 < 0.15 && thisrwtemp > 0.9){ - thisrwtemp = hRPArelratio->GetBinContent(q3bin+150, q0bin-q0offsetValenciaGENIE); - } - - - - //Double_t *mypoly; - //mypoly = rpapoly; - - if(Q2gev >= 9.0){ - thisrwtemp = 1.0; - } - else if(Q2gev > 3.0) { - // hiding option, use the Q2 parameterization everywhere - // } else if(Q2gev > 3.0 || rwRPAQ2) { - // turn rwRPAQ2 all the way on to override the 2D histogram - // illustrates the old-style Q2 suppression many folks still use. - - thisrwtemp = getWeightQ2(Q2gev,true); - - // double powerQ2 = 1.0; - //thisrwtemp = 0.0; - //for(int ii=0; ii<10; ii++){ - // thisrwtemp += rpapoly[ii]*powerQ2; - // powerQ2 *= Q2gev; - //} - //std::cout << "test temp " << thisrwtemp << " " << rpamypoly[2] << std::endl; - } - - if(!(thisrwtemp >= 0.001 && thisrwtemp <= 2.0))thisrwtemp = 1.0; - - // hiding option, turn off the enhancement. - //if(rwoffSRC && thisrwtemp > 1.0)thisrwtemp = 1.0; - - if (0 == weights) return thisrwtemp; - // if this was called without passing an array, - // the user didn't want us to calculate the +/- 1-sigma bounds - // so the above line returned before even trying. - - weights[0] = thisrwtemp; - - //if (type == 0) return thisrwtemp; - - //if (type == 1) { - // Construct the error bands on the low Q2 suppression. - // Double_t thisrwSupP1 = 1.0; - // Double_t thisrwSupM1 = 1.0; - - if( thisrwtemp < 1.0){ - // make the suppression stronger or weaker to muon capture uncertainty - // rwRPAonesig is either +1 or -1, which is 0.25 (25%). - // possible to be re-written to produce 2 and 3 sigma. - - weights[1] = thisrwtemp + 1.0 * (0.25)*(1.0 - thisrwtemp); - weights[2] = thisrwtemp - 1.0 * (0.25)*(1.0 - thisrwtemp); - - - - } - else{ - weights[1] = thisrwtemp; - weights[2] = thisrwtemp; - } - - //std::cout << "check " << thisrwtemp << " " << weights[1] << " " << weights[2] << std::endl; - - - - - // Construct the rest of the error bands on the low Q2 suppression. - // this involves getting the weight from the non-relativistic ratio - - //if (type == 2){ - - Double_t thisrwEnhP1 = 1.0; - Double_t thisrwEnhM1 = 1.0; - - // make enhancement stronger or weaker to Federico Sanchez uncertainty - // this does NOT mean two sigma, its overloading the option. - Double_t thisrwextreme = hRPAnonrelratio->GetBinContent(q3bin,q0bin-q0offsetValenciaGENIE); - // now trap bogus entries. Not sure why they happen, but set to 1.0 not 0.0 - if(thisrwextreme <= 0.001)thisrwextreme = 1.0; - - if(mc_q0 < 0.15 && thisrwextreme > 0.9){ - thisrwextreme = hRPAnonrelratio->GetBinContent(q3bin+150, q0bin-q0offsetValenciaGENIE); - } - - //std::cout << "ext " << thisrwextreme << " " << thisrwtemp << std::endl; - - // get the same for the Q2 dependent thing, - // but from the nonrelativistic polynomial - - if(Q2gev >= 9.0){ - thisrwextreme = 1.0; - } - else if(Q2gev > 3.0 ) { - thisrwextreme = getWeightQ2(Q2gev,false); - //double powerQ2 = 1.0; - //thisrwextreme = 0.0; - //for(int ii=0; ii<10; ii++){ - // thisrwextreme += rpapolynonrel[ii]*powerQ2; - // powerQ2 *= Q2gev; - //} - //std::cout << "test extreme " << thisrwextreme << " " << mypolynonrel[2] << std::endl; - } - - if(!(thisrwextreme >= 0.001 && thisrwextreme <= 2.0))thisrwextreme = 1.0; - - //std::cout << "test extreme " << Q2gev << " " << thisrwextreme << " " << thisrwtemp << std::endl; - - Double_t RelToNonRel = 0.6; - - // based on some distance between central value and extreme - thisrwEnhP1 = thisrwtemp + RelToNonRel * (thisrwextreme-thisrwtemp); - Double_t thisrwEnhP1max = thisrwextreme; - - if(Q2gev < 0.9)thisrwEnhP1 += 1.5*(0.9 - Q2gev)*(thisrwEnhP1max - thisrwEnhP1); - // sanity check, don't let the upper error bound go above the nonrel limit. - if(thisrwEnhP1 > thisrwEnhP1max)thisrwEnhP1 = thisrwEnhP1max; - // don't let positive error bound be closer than 3% above the central value - // will happen at very high Q2 and very close to Q2 = 0 - if(thisrwEnhP1 < thisrwtemp + 0.03)thisrwEnhP1 = thisrwtemp + 0.03; - - thisrwEnhM1 = thisrwtemp - RelToNonRel * (thisrwextreme-thisrwtemp); - // don't let negative error bound be closer than 3% below the central value - if(thisrwEnhM1 > thisrwtemp - 0.03)thisrwEnhM1 = thisrwtemp - 0.03; - // even still, don't let the lower error bound go below 1.0 at high-ish Q2 - if(Q2gev > 1.0 && thisrwEnhM1 < 1.0)thisrwEnhM1 = 1.0; - - // whew. so now return the main weight - // and return the array of all five weights in some array - // thisrwtemp, thisrwSupP1, thisrwSupM1, thisrwEnhP1, thisrwEnhM1 - - //if (sign == 1) return thisrwEnhP1; - //if (sign == -1) return thisrwEnhM1; - - weights[3] = thisrwEnhP1; - weights[4] = thisrwEnhM1; - - // still return the central value - return thisrwtemp; - -} - - -double weightRPA::getWeight(const double mc_q0, const double mc_q3){ - - return getWeightInternal(mc_q0, mc_q3); - -} - -double weightRPA::getWeight(const double mc_q0, const double mc_q3, double *weights){ - - return getWeightInternal(mc_q0, mc_q3, weights); - -} - -double weightRPA::getWeight(const double mc_q0, const double mc_q3, int type, int sign){ - - return getWeightInternal(mc_q0, mc_q3, type, sign); - -} - -double weightRPA::getWeight(const double mc_Q2){ - - return getWeightQ2(mc_Q2); - -} - -double weightRPA::getWeightInternal(const double mc_q0, const double mc_q3, int type, int sign){ - - double weights[5] = {1., 1., 1., 1., 1.}; - double cv = getWeightInternal(mc_q0, mc_q3, weights); - - if(type==0)return cv; - else if(type==weightRPA::LOWQ2 && sign == 1)return weights[1]; - else if(type==weightRPA::LOWQ2 && sign == -1)return weights[2]; - else if(type==weightRPA::HIGHQ2 && sign == 1)return weights[3]; - else if(type==weightRPA::HIGHQ2 && sign == -1)return weights[4]; - //else { - // // should never happen? Bork ? - // return cv; //? - //} - - return cv; -} - -double weightRPA::getWeightQ2(const double mc_Q2, const bool relORnonrel){ - - if(mc_Q2 < 0.0)return 1.0; // this is Q2 actually, not sure hw - if(mc_Q2 > 9.0)return 1.0; - - // this function needs to know two options. - // does user want rel (cv) or nonrel - // does user want to use the histogram or parameterization - - if(Q2histORparam)return getWeightQ2fromhistogram(mc_Q2, relORnonrel); - else return getWeightQ2parameterization(mc_Q2, relORnonrel); - -} - -double weightRPA::getWeightQ2parameterization(const double mc_Q2, const bool relORnonrel){ - - if(mc_Q2 < 0.0)return 1.0; - if(mc_Q2 > 9.0)return 1.0; - - // this one returns just the polynomial Q2 version - // for special tests. Poor answer for baseline MINERvA QE events. - // No uncertainty assigned to this usecase at this time. - //double gevmev = 0.001; // minerva sends in MeV. - double Q2gev = mc_Q2; - double powerQ2 = 1.0; - double thisrwtemp = 0.0; - thisrwtemp = 0.0; - for(int ii=0; ii<10; ii++){ - if(relORnonrel)thisrwtemp += rpapolyrel[ii]*powerQ2; - else thisrwtemp += rpapolynonrel[ii]*powerQ2; - powerQ2 *= Q2gev; - } - return thisrwtemp; - - -} - -double weightRPA::getWeightQ2fromhistogram(const double mc_Q2, const bool relORnonrel){ - - if(mc_Q2 < 0.0)return 1.0; - if(mc_Q2 > 9.0) return 1.0; - - if(relORnonrel)return hQ2relratio->GetBinContent( hQ2relratio->FindBin(mc_Q2) ); - else return hQ2nonrelratio->GetBinContent( hQ2nonrelratio->FindBin(mc_Q2) ); - - // interpolation might be overkill for such a finely binned histogram, 0.01% - // but the extra cpu cycles maybe small. - // save it here for some future use. - //if(relORnonrel)return hQ2relratio->Interpolate(mc_Q2); - //else return hQ2nonrelratio->Interpolate(mc_Q2); - -} - - -double weightRPA::getWeightLowQ2(const double mc_q0, const double mc_q3, int sign){ - return getWeightInternal(mc_q0,mc_q3,weightRPA::LOWQ2,sign); -} - -double weightRPA::getWeightHighQ2(const double mc_q0, const double mc_q3, int sign){ - return getWeightInternal(mc_q0,mc_q3,weightRPA::HIGHQ2,sign); -} - - -#endif - +#ifndef weightRPA_h +#define weightRPA_h + +#include //ifstream +#include //cout + +#include +#include +#include +#include +#include "math.h" +#include "assert.h" +//For Compatibility with ROOT compiler +//uncomment the following: +//#define ROOT + +/*! + * Code example to extract the RPA effect central weight + * and its uncertainties from the prepared files. + * Heidi Schellman (Oregon State) and Rik Gran (Minnesota Duluth) + * for use in MINERvA experiment analysis + * must compile with the ROOT libraries + * g++ `root-config --glibs --cflags` -O3 weightRPAtest.cxx -o weightRPA + + + * The underlying model is from the IFIC Valencia group + * see (public) minerva docdb:12688 for the full physics discussion + + */ + + +// NOTE UNITS ARE GEV in the calculation +// make sure you convert MeV to GeV before calling these functions + + +// Class for getting RPA paramaters inputs from a given file +// Includes methods that return all five RPA weights at once +// (is the most cpu efficient way to get them) +// Or return just the central value +// (skipping the uncertainties code completely if only cv wanted) +// Or return each CV and uncertainty one at a time +// (is the least cpu efficient, repeats some calculations 5 times) + +class weightRPA { +public: + //Constructor: Read in params from a filename + weightRPA(const TString f) { read(f); } //Read in params from file + + TString filename; + TFile* fRPAratio; + TH2D *hRPArelratio; + TH2D *hRPAnonrelratio; + TH1D *hQ2relratio; + TH1D *hQ2nonrelratio; + TArrayD *TADrpapolyrel; + Double_t *rpapolyrel ; + TArrayD *TADrpapolynonrel; + Double_t *rpapolynonrel; + static const int CENTRAL=0; + static const int LOWQ2 = 1; + static const int HIGHQ2 = 2; + + // true to take from histogram, false use parameterization + static const bool Q2histORparam = true; + + // MINERvA holds kinematics in MeV, but all these functions require GeV + // So make sure you pass them in GeV. + inline double getWeight(const double mc_q0, const double mc_q3, double * weights); //in GeV + inline double getWeight(const double mc_q0, const double mc_q3); //in GeV + inline double getWeight(const double mc_q0, const double mc_q3, int type, int sign); //in GeV + inline double getWeight(const double Q2); //in GeV^2 + inline double getWeightLowQ2(const double mc_q0, const double mc_q3, const int sign); + inline double getWeightHighQ2(const double mc_q0, const double mc_q3, const int sign); + inline double getWeightQ2(const double mc_Q2, const bool relORnonrel=true); + //Initializer + inline void read(const TString f); + + // q0 and q3 in GeV, type = 1 for low Q2, 2 for high Q2, 0 for central + //double getWeightInternal(const double mc_q0, const double mc_q3,int type, int sign); + + private: + inline double getWeightInternal(const double mc_Q2); + + inline double getWeightInternal(const double mc_q0, const double mc_q3, const int type, const int sign); + inline double getWeightInternal(const double mc_q0, const double mc_q3, double *weights=0); + inline double getWeightQ2parameterization(const double mc_Q2, const bool relORnonrel); + inline double getWeightQ2fromhistogram(const double mc_Q2, const bool relORnonrel); + + +}; + +void weightRPA::read(const TString f) +//Read in the params doubles from a file +//argument: valid filename +{ + fRPAratio = TFile::Open(f,"READONLY"); + if (fRPAratio){ + hRPArelratio = (TH2D*)fRPAratio->Get("hrelratio"); + hRPAnonrelratio = (TH2D*)fRPAratio->Get("hnonrelratio"); + hQ2relratio = (TH1D*)fRPAratio->Get("hQ2relratio"); + hQ2nonrelratio = (TH1D*)fRPAratio->Get("hQ2nonrelratio"); + TADrpapolyrel = (TArrayD*)fRPAratio->Get("rpapolyrel"); + rpapolyrel = TADrpapolyrel->GetArray(); + TADrpapolynonrel = (TArrayD*)fRPAratio->Get("rpapolynonrel"); + rpapolynonrel = TADrpapolynonrel->GetArray(); + hRPArelratio->Print(); + std::cout << "have read in ratios from file " << f <= gevlimit) q0bin = rpamevlimit - 1; + if(mc_q3 >= gevlimit) q3bin = rpamevlimit - 1; + + // Nieves does not calculate anything below binding energy. + // I don't know what GENIE does, but lets be soft about this. + // Two things lurking here at once. + // One, we are taking out a 10 MeV offset between GENIE and Valencia. + // Second, we are protecting against being asked for a weight that is too low in q0. + // It actually shouldn't happen for real GENIE events, + // but this protection does something that doesn't suck, just in case. + // you would see the artifact in a plot for sure, but better than writing 1.0. + + // CWret after talking to Rik in Nov 2018 at MINERvA CM + // 2.12.8 sets this to 27 because change of Q definition: need to offset GENIE and Nieves Eb even more +#if __GENIE_VERSION__ >= 210 + Int_t q0offsetValenciaGENIE = 27; +#else + Int_t q0offsetValenciaGENIE = 10; +#endif + + if(mc_q0 < 0.018) q0bin = 18+q0offsetValenciaGENIE; + Double_t thisrwtemp = hRPArelratio->GetBinContent(q3bin,q0bin-q0offsetValenciaGENIE); + + // now trap bogus entries. Not sure why they happen, but set to 1.0 not 0.0 + if(thisrwtemp <= 0.001)thisrwtemp = 1.0; + + // events in genie but not in valencia should get a weight + // related to a similar q0 from the bulk distribution. + if(mc_q0 < 0.15 && thisrwtemp > 0.9){ + thisrwtemp = hRPArelratio->GetBinContent(q3bin+150, q0bin-q0offsetValenciaGENIE); + } + + + + //Double_t *mypoly; + //mypoly = rpapoly; + + if(Q2gev >= 9.0){ + thisrwtemp = 1.0; + } + else if(Q2gev > 3.0) { + // hiding option, use the Q2 parameterization everywhere + // } else if(Q2gev > 3.0 || rwRPAQ2) { + // turn rwRPAQ2 all the way on to override the 2D histogram + // illustrates the old-style Q2 suppression many folks still use. + + thisrwtemp = getWeightQ2(Q2gev,true); + + // double powerQ2 = 1.0; + //thisrwtemp = 0.0; + //for(int ii=0; ii<10; ii++){ + // thisrwtemp += rpapoly[ii]*powerQ2; + // powerQ2 *= Q2gev; + //} + //std::cout << "test temp " << thisrwtemp << " " << rpamypoly[2] << std::endl; + } + + if(!(thisrwtemp >= 0.001 && thisrwtemp <= 2.0))thisrwtemp = 1.0; + + // hiding option, turn off the enhancement. + //if(rwoffSRC && thisrwtemp > 1.0)thisrwtemp = 1.0; + + if (0 == weights) return thisrwtemp; + // if this was called without passing an array, + // the user didn't want us to calculate the +/- 1-sigma bounds + // so the above line returned before even trying. + + weights[0] = thisrwtemp; + + //if (type == 0) return thisrwtemp; + + //if (type == 1) { + // Construct the error bands on the low Q2 suppression. + // Double_t thisrwSupP1 = 1.0; + // Double_t thisrwSupM1 = 1.0; + + if( thisrwtemp < 1.0){ + // make the suppression stronger or weaker to muon capture uncertainty + // rwRPAonesig is either +1 or -1, which is 0.25 (25%). + // possible to be re-written to produce 2 and 3 sigma. + + weights[1] = thisrwtemp + 1.0 * (0.25)*(1.0 - thisrwtemp); + weights[2] = thisrwtemp - 1.0 * (0.25)*(1.0 - thisrwtemp); + + + + } + else{ + weights[1] = thisrwtemp; + weights[2] = thisrwtemp; + } + + //std::cout << "check " << thisrwtemp << " " << weights[1] << " " << weights[2] << std::endl; + + + + + // Construct the rest of the error bands on the low Q2 suppression. + // this involves getting the weight from the non-relativistic ratio + + //if (type == 2){ + + Double_t thisrwEnhP1 = 1.0; + Double_t thisrwEnhM1 = 1.0; + + // make enhancement stronger or weaker to Federico Sanchez uncertainty + // this does NOT mean two sigma, its overloading the option. + Double_t thisrwextreme = hRPAnonrelratio->GetBinContent(q3bin,q0bin-q0offsetValenciaGENIE); + // now trap bogus entries. Not sure why they happen, but set to 1.0 not 0.0 + if(thisrwextreme <= 0.001)thisrwextreme = 1.0; + + if(mc_q0 < 0.15 && thisrwextreme > 0.9){ + thisrwextreme = hRPAnonrelratio->GetBinContent(q3bin+150, q0bin-q0offsetValenciaGENIE); + } + + //std::cout << "ext " << thisrwextreme << " " << thisrwtemp << std::endl; + + // get the same for the Q2 dependent thing, + // but from the nonrelativistic polynomial + + if(Q2gev >= 9.0){ + thisrwextreme = 1.0; + } + else if(Q2gev > 3.0 ) { + thisrwextreme = getWeightQ2(Q2gev,false); + //double powerQ2 = 1.0; + //thisrwextreme = 0.0; + //for(int ii=0; ii<10; ii++){ + // thisrwextreme += rpapolynonrel[ii]*powerQ2; + // powerQ2 *= Q2gev; + //} + //std::cout << "test extreme " << thisrwextreme << " " << mypolynonrel[2] << std::endl; + } + + if(!(thisrwextreme >= 0.001 && thisrwextreme <= 2.0))thisrwextreme = 1.0; + + //std::cout << "test extreme " << Q2gev << " " << thisrwextreme << " " << thisrwtemp << std::endl; + + Double_t RelToNonRel = 0.6; + + // based on some distance between central value and extreme + thisrwEnhP1 = thisrwtemp + RelToNonRel * (thisrwextreme-thisrwtemp); + Double_t thisrwEnhP1max = thisrwextreme; + + if(Q2gev < 0.9)thisrwEnhP1 += 1.5*(0.9 - Q2gev)*(thisrwEnhP1max - thisrwEnhP1); + // sanity check, don't let the upper error bound go above the nonrel limit. + if(thisrwEnhP1 > thisrwEnhP1max)thisrwEnhP1 = thisrwEnhP1max; + // don't let positive error bound be closer than 3% above the central value + // will happen at very high Q2 and very close to Q2 = 0 + if(thisrwEnhP1 < thisrwtemp + 0.03)thisrwEnhP1 = thisrwtemp + 0.03; + + thisrwEnhM1 = thisrwtemp - RelToNonRel * (thisrwextreme-thisrwtemp); + // don't let negative error bound be closer than 3% below the central value + if(thisrwEnhM1 > thisrwtemp - 0.03)thisrwEnhM1 = thisrwtemp - 0.03; + // even still, don't let the lower error bound go below 1.0 at high-ish Q2 + if(Q2gev > 1.0 && thisrwEnhM1 < 1.0)thisrwEnhM1 = 1.0; + + // whew. so now return the main weight + // and return the array of all five weights in some array + // thisrwtemp, thisrwSupP1, thisrwSupM1, thisrwEnhP1, thisrwEnhM1 + + //if (sign == 1) return thisrwEnhP1; + //if (sign == -1) return thisrwEnhM1; + + weights[3] = thisrwEnhP1; + weights[4] = thisrwEnhM1; + + // still return the central value + return thisrwtemp; + +} + + +double weightRPA::getWeight(const double mc_q0, const double mc_q3){ + + return getWeightInternal(mc_q0, mc_q3); + +} + +double weightRPA::getWeight(const double mc_q0, const double mc_q3, double *weights){ + + return getWeightInternal(mc_q0, mc_q3, weights); + +} + +double weightRPA::getWeight(const double mc_q0, const double mc_q3, int type, int sign){ + + return getWeightInternal(mc_q0, mc_q3, type, sign); + +} + +double weightRPA::getWeight(const double mc_Q2){ + + return getWeightQ2(mc_Q2); + +} + +double weightRPA::getWeightInternal(const double mc_q0, const double mc_q3, int type, int sign){ + + double weights[5] = {1., 1., 1., 1., 1.}; + double cv = getWeightInternal(mc_q0, mc_q3, weights); + + if(type==0)return cv; + else if(type==weightRPA::LOWQ2 && sign == 1)return weights[1]; + else if(type==weightRPA::LOWQ2 && sign == -1)return weights[2]; + else if(type==weightRPA::HIGHQ2 && sign == 1)return weights[3]; + else if(type==weightRPA::HIGHQ2 && sign == -1)return weights[4]; + //else { + // // should never happen? Bork ? + // return cv; //? + //} + + return cv; +} + +double weightRPA::getWeightQ2(const double mc_Q2, const bool relORnonrel){ + + if(mc_Q2 < 0.0)return 1.0; // this is Q2 actually, not sure hw + if(mc_Q2 > 9.0)return 1.0; + + // this function needs to know two options. + // does user want rel (cv) or nonrel + // does user want to use the histogram or parameterization + + if(Q2histORparam)return getWeightQ2fromhistogram(mc_Q2, relORnonrel); + else return getWeightQ2parameterization(mc_Q2, relORnonrel); + +} + +double weightRPA::getWeightQ2parameterization(const double mc_Q2, const bool relORnonrel){ + + if(mc_Q2 < 0.0)return 1.0; + if(mc_Q2 > 9.0)return 1.0; + + // this one returns just the polynomial Q2 version + // for special tests. Poor answer for baseline MINERvA QE events. + // No uncertainty assigned to this usecase at this time. + //double gevmev = 0.001; // minerva sends in MeV. + double Q2gev = mc_Q2; + double powerQ2 = 1.0; + double thisrwtemp = 0.0; + thisrwtemp = 0.0; + for(int ii=0; ii<10; ii++){ + if(relORnonrel)thisrwtemp += rpapolyrel[ii]*powerQ2; + else thisrwtemp += rpapolynonrel[ii]*powerQ2; + powerQ2 *= Q2gev; + } + return thisrwtemp; + + +} + +double weightRPA::getWeightQ2fromhistogram(const double mc_Q2, const bool relORnonrel){ + + if(mc_Q2 < 0.0)return 1.0; + if(mc_Q2 > 9.0) return 1.0; + + if(relORnonrel)return hQ2relratio->GetBinContent( hQ2relratio->FindBin(mc_Q2) ); + else return hQ2nonrelratio->GetBinContent( hQ2nonrelratio->FindBin(mc_Q2) ); + + // interpolation might be overkill for such a finely binned histogram, 0.01% + // but the extra cpu cycles maybe small. + // save it here for some future use. + //if(relORnonrel)return hQ2relratio->Interpolate(mc_Q2); + //else return hQ2nonrelratio->Interpolate(mc_Q2); + +} + + +double weightRPA::getWeightLowQ2(const double mc_q0, const double mc_q3, int sign){ + return getWeightInternal(mc_q0,mc_q3,weightRPA::LOWQ2,sign); +} + +double weightRPA::getWeightHighQ2(const double mc_q0, const double mc_q3, int sign){ + return getWeightInternal(mc_q0,mc_q3,weightRPA::HIGHQ2,sign); +} + + +#endif +