diff --git a/test/TestSplineDTAdivision.cc b/test/TestSplineDTAdivision.cc index 08b2ea5..8004e4f 100644 --- a/test/TestSplineDTAdivision.cc +++ b/test/TestSplineDTAdivision.cc @@ -1,211 +1,207 @@ #include #include #include #include #include #include "TFile.h" #include "TH1.h" #include "TF1.h" #include "TCanvas.h" #include "TVirtualPad.h" #include "TLegend.h" #include "TString.h" #include "TObjArray.h" #include "TRegexp.h" #include "TObjString.h" #include "TStyle.h" #include "TFitResult.h" #include "TGraph.h" #include "Lau1DCubicSpline.hh" void writeVecToFile(const std::vector& vec, std::ofstream& file) { for(auto it = vec.begin(); it < vec.end(); ++it) { file << *it; if(it != vec.end()-1){file << ", ";} } file << "\n\n"; return; } std::vector convertStringToVector(TString string, TString delim = "_") { std::vector vec = {}; TObjArray* ObjArray = string.Tokenize(delim); for(Int_t i = 0; i <= ObjArray->GetLast(); ++i) { vec.push_back( static_cast(ObjArray->At(i))->String() ); } delete ObjArray; return vec; } TLegend* makeLegend() { TLegend* leg = new TLegend; TH1D* l1 = new TH1D(); l1->SetLineColor(8); leg->AddEntry(l1, "#chi^{2} successful", "l"); TH1D* l2 = new TH1D(); l2->SetLineColor(9); leg->AddEntry(l2, "WL successful", "l"); TH1D* l3 = new TH1D(); l3->SetLineColor(2); leg->AddEntry(l3, "Neither successful", "l"); TGraph* g = new TGraph; g->SetMarkerColor(41); g->SetMarkerStyle(3); leg->AddEntry(g, "Knots", "p"); leg->SetLineColorAlpha(0,0); return leg; } int main() { std::ofstream outputFile; outputFile.open("MCsplines.dat"); TFile* DTA_file = TFile::Open("root://eoslhcb.cern.ch//eos/lhcb/wg/b2oc/B2D0pipi/Tuples22/DTAhists.root"); - //TFile* DTA_file = TFile::Open("/storage/epp2/phrdkf/DTAhists_fewerBins.root"); outputFile << "X vals : \n"; -// const std::vector dtvals = {0.2 ,0.25, 0.5,0.75 ,1., 2., 3., 4., 7., /*10.,*/ 15.}; const std::vector dtvals = {0.2 , 0.5,0.75 ,1., 2., 3., 4., 6., 15.}; writeVecToFile(dtvals, outputFile); const std::map > initialVals = { - {"Bd2DKpi", {0.978086,1.03688 ,0.814762,0.973555,0.901241,0.985174,1.03396 ,1.08291,1.00181,/*1.1311,*/0.1}}, - {"Bs2DKpi", {1.0256 ,0.906071,0.785794,0.767683,0.749571,0.88353 ,0.974613,1.04696,1.14629,/*1.15958,*/0.983078}}, - {"Lb2Dppi", {1.30118 ,0.295433,0.491633,0.619358,0.70069,1.02602 ,1.09221 ,1.15841,1.05802,/*1.06077,*/0.237281}}, - {"signalMC",{1.69052e-04, 8.99483e-05, 1.89572e-03, 4.41640e-03, 7.07254e-03, 1.42052e-02, 1.70844e-02, 1.79040e-02, 1.76894e-02, /*1.62783e-02,*/ 1.36895e-02}} + {"Bd2DKpi", {0.978086,1.03688 ,0.814762,0.973555,0.901241,0.985174,1.03396 ,1.08291,1.00181,0.1}}, + {"Bs2DKpi", {1.0256 ,0.906071,0.785794,0.767683,0.749571,0.88353 ,0.974613,1.04696,1.14629,0.983078}}, + {"Lb2Dppi", {1.30118 ,0.295433,0.491633,0.619358,0.70069,1.02602 ,1.09221 ,1.15841,1.05802,0.237281}}, + {"signalMC",{1.69052e-04, 8.99483e-05, 1.89572e-03, 4.41640e-03, 7.07254e-03, 1.42052e-02, 1.70844e-02, 1.79040e-02, 1.76894e-02, 1.36895e-02}} }; const std::map BmodeIndices = { {"Bd2DKpi", 3}, {"Bs2DKpi", 1}, {"Lb2Dppi", 2}, {"signalMC",0} }; const std::map DmodeIndices = { {"Kpi" , 0}, {"KK" , 1}, {"pipi", 2} }; std::array canvases = {new TCanvas(), new TCanvas()}; canvases[0]->Divide(4,3); canvases[1]->Divide(4,3); gStyle->SetOptStat(0); bool run1DTAappend = false, run2DTAappend = false; for( const auto&& obj : *(DTA_file->GetListOfKeys()) ) { TString histName = obj->GetName(); if(histName.Contains("Sideband")){continue;} std::cout << "histName : " << histName << std::endl; std::vector tokenizedName = convertStringToVector( histName ); const int run = tokenizedName.at(2)[3] - '0'; //Cool ascii hack converts char -> int outputFile << histName; std::vector effvals = initialVals.at( tokenizedName[0] ); effvals.assign(dtvals.size(), 1.); TH1D* thid = dynamic_cast( DTA_file->Get(histName) ); thid->SetDirectory(0); if(not thid) { std::cerr << "ERROR : Could not find histogram : " << histName << std::endl; std::cerr << "skipping..." << std::endl; continue; } - //Lau1DCubicSpline::SplineType splineType = Lau1DCubicSpline::SplineType::AkimaSpline; Lau1DCubicSpline::SplineType splineType = Lau1DCubicSpline::SplineType::StandardSpline; if(not histName.BeginsWith("signalMC_Kpi")) { TString denomName = Form( "signalMC%s", histName("_.+$").Data() ); denomName.ReplaceAll("_KK_", "_Kpi_"); denomName.ReplaceAll("_pipi_","_Kpi_"); TH1D* denom = dynamic_cast( DTA_file->Get(denomName) ); if(not denom){std::cerr << "ERROR : Could not find denominator!" << std::endl;} outputFile << " / " << denomName; thid->Divide( denom ); delete denom; } else { splineType = Lau1DCubicSpline::SplineType::AkimaSpline; } std::map fixedParams = {}; if( histName.Contains("signalMC") and not histName.Contains("Kpi") ) { fixedParams[0] = 1.; } - //fixedParams = {}; Lau1DCubicSpline* dtEffSpline = nullptr; for(size_t i = 0; i < effvals.size(); ++i) {effvals[i] = thid->GetBinContent(thid->FindFixBin(dtvals[i]));} dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,splineType); Int_t lineColor = 8; TFitResultPtr fit = dtEffSpline->fitToTH1(thid, LauOutputLevel::Quiet, false, fixedParams); if( fit->CovMatrixStatus() != 3 ) { lineColor = 9; std::cout << "chi2 fit failed : trying WL ..." << std::endl; fit = dtEffSpline->fitToTH1(thid, LauOutputLevel::Quiet, true, fixedParams); if( fit->CovMatrixStatus() != 3 ) { std::cout << " ... This also failed" << std::endl; lineColor = 2; } else{std::cout << " ... Success!" << std::endl;} } TVirtualPad* pad = canvases.at(run-1)->cd( BmodeIndices.at(tokenizedName.at(0)) + 4*DmodeIndices.at(tokenizedName.at(1)) + 1 ); pad->SetLogx(); thid->SetLineColor(1); thid->Draw("SCAT E"); TF1* tf1 = dtEffSpline->makeTF1(); tf1->SetLineColor( lineColor ); tf1->Draw("SAME"); TGraph* graph = new TGraph( dtEffSpline->getnKnots(), dtEffSpline->getXValues().data(), dtEffSpline->getYValues().data() ); graph->SetMarkerColor(41); graph->Draw("SAME*"); outputFile << "\n"; writeVecToFile( dtEffSpline->getYValues() , outputFile); switch(run) { case 1: dtEffSpline->writeToJson("run1DTAsplines.json", obj->GetName(), run1DTAappend); run1DTAappend = true; break; case 2: dtEffSpline->writeToJson("run2DTAsplines.json", obj->GetName(), run2DTAappend); run2DTAappend = true; break; default: std::cerr << "ERROR: Default case hit!" << std::endl; break; } } TLegend* legend = makeLegend(); canvases[0]->cd(1); legend->Draw(); canvases[1]->cd(1); legend->Draw(); canvases[0]->SaveAs("Run1DTAfits.pdf"); canvases[1]->SaveAs("Run2DTAfits.pdf"); outputFile.close(); DTA_file->Close(); return 0; }