diff --git a/MINERvA_NOvA_WS/FlatAnalysis.cc b/MINERvA_NOvA_WS/FlatAnalysis.cc new file mode 100644 index 0000000..e428144 --- /dev/null +++ b/MINERvA_NOvA_WS/FlatAnalysis.cc @@ -0,0 +1,518 @@ +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 new file mode 100644 index 0000000..9f71121 --- /dev/null +++ b/MINERvA_NOvA_WS/exp_list.xml @@ -0,0 +1,149 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/MINERvA_NOvA_WS/minerva_tune.xml b/MINERvA_NOvA_WS/minerva_tune.xml new file mode 100644 index 0000000..aecb234 --- /dev/null +++ b/MINERvA_NOvA_WS/minerva_tune.xml @@ -0,0 +1,37 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/MINERvA_NOvA_WS/plotscript.cc b/MINERvA_NOvA_WS/plotscript.cc new file mode 100644 index 0000000..acdeb05 --- /dev/null +++ b/MINERvA_NOvA_WS/plotscript.cc @@ -0,0 +1,245 @@ +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.md b/README.md new file mode 100644 index 0000000..fd753ae --- /dev/null +++ b/README.md @@ -0,0 +1,98 @@ +# 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. + +---- +## 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. + +---- +## 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) +