QM25_productions/D0/D0_fitter_QM25_noPulls.C

0001 #include <string>
0002 #include <iostream>
0003 #include <sstream>
0004 #include <cmath>
0005 #ifndef __CINT__
0006 #include <RooGlobalFunc.h>
0007 #endif
0008 #include <RooRealVar.h>
0009 #include <RooDataSet.h>
0010 #include <RooDataHist.h>
0011 #include <RooGaussian.h>
0012 #include <RooBifurGauss.h>
0013 //#include <RooDoubleCB.h>
0014 #include <RooHist.h>
0015 #include <RooPlot.h>
0016 #include <RooExponential.h>
0017 #include <RooAddPdf.h>
0018 #include <RooStats/SPlot.h>
0019 #include <TH1F.h>
0020 #include <TFile.h>
0021 #include <TTree.h>
0022 #include <TCut.h>
0023 #include <TChain.h>
0024 #include <TMath.h>
0025 #include <TPad.h>
0026 #include <TCanvas.h>
0027 #include <TBranch.h>
0028
0029 using namespace RooFit;
0030 using namespace std;
0031
0032 //https://stackoverflow.com/questions/997946/how-to-get-current-time-and-date-in-c
0033 std::string getDate()
0034 {
0035     std::time_t t = std::time(0);   // get time now
0036     std::tm* now = std::localtime(&t);
0037
0038     std::stringstream date;
0039     date << (now->tm_mon + 1) << '/'
0040          <<  now->tm_mday << '/'
0041          << (now->tm_year + 1900);
0042     return date.str();
0043 }
0044
0045 template <typename T>
0046 string to_string_with_precision(const T a_value, const int n = 2)
0047 {
0048     ostringstream out;
0049     out.precision(n);
0050     out << fixed << a_value;
0051     return out.str();
0052 }
0053
0054 void D0_fitter_QM25_noPulls()
0055 {
0056   int maxDigits = 3;
0057   TGaxis::SetMaxDigits(maxDigits);
0058
0059   string path, inputFile, plotTitle, saveName, xAxisTitle, branch, dataType;
0060   double minMass = 1.7;
0061   double maxMass = 2.0;
0062
0063   path = "/sphenix/user/cdean/public/QM25/";
0064
0065   inputFile = path + "output_D0_oldKFP_ana475p017_optimized.root";
0066   dataType = "pp";
0067
0068   branch = "D0_mass";
0069   plotTitle = "";
0070   saveName = path + "fits/D02Kpi_ana475p01_noPull";
0071   xAxisTitle = "m(K^{#mp}#pi^{#pm}) [GeV]";
0072
0073   stringstream cutStream;
0074   cutStream << minMass << " <= " << branch << " && " << branch << " <= " << maxMass;
0075   TCut masscut = cutStream.str().c_str();
0076
0077   /*
0078    * Get files and data sets
0079    */
0080
0081   TFile* dataFile = new TFile(inputFile.c_str());
0082   TTree* dataTree = (TTree*)dataFile->Get("DecayTree");
0083
0084   string datasWeight =  inputFile.substr(0, inputFile.size()-5) + "_weighted.root";
0085   TFile* sWeightedDataFile  = new TFile(datasWeight.c_str(), "RECREATE");
0086   TTree* dataCopyTree = dataTree->CopyTree(masscut);
0087   TTree* dataCloneTree = dataCopyTree->CloneTree(-1);
0088
0089   RooRealVar mass(branch.c_str(), "mass", minMass, maxMass);
0090   RooDataSet dataSet(branch.c_str(), "data", mass, Import(*dataCloneTree));
0091
0092   /*
0093    * Signal Model
0094    */
0095
0096   RooRealVar  mean("D0_mean", "mean", 1.864, 1.855, 1.875);
0097
0098   //Base Gaussian model
0099   RooRealVar  sigma("sigma", "sigma", 0.01, 0.000, 0.02);
0100   RooGaussian signal("signal", "signal", mass, mean, sigma);
0101
0102   RooRealVar  fSig("fSig", "fSig", 0.*dataSet.numEntries(), 0., 0.1*dataSet.numEntries());
0103   RooRealVar  fBkg("fBkg", "fBkg", 0.*dataSet.numEntries(), 0., 1*dataSet.numEntries());
0104
0105   /*
0106    * Background Model
0107    */
0108
0109   RooRealVar expConstOne("expConstOne", "expConstOne", -10, -1e2, 0.);
0110   RooExponential background("background", "background", mass, expConstOne);
0111
0112   /*
0113    * Fitting to the data
0114    */
0115
0116   RooArgList fitModelList(signal, background), fitFracList(fSig);
0117   RooAddPdf model("model", "model", fitModelList, fitFracList);
0118   model.fitTo(dataSet);
0119
0120   RooPlot* frame = mass.frame(Title(plotTitle.c_str()));       //creating the frame
0121
0122   RooBinning bins(minMass, maxMass);
0123   int nBins = 25;
0124   bins.addUniform(nBins, minMass, maxMass);
0125
0126   dataSet.plotOn(frame, Binning(bins), XErrorSize(0), DataError(RooAbsData::SumW2));          //plotting the raw unbinned data on the fram
0127   model.plotOn(frame, Components(background), LineColor(kGray), DrawOption("F"), FillColor(kGray));
0128   model.plotOn(frame, Components(RooArgSet(signal, background)), LineColor(kAzure+8), DrawOption("F"), FillColor(kAzure+8), MoveToBack());
0129   model.plotOn(frame, LineColor(kBlack));        //plotting the fit onto the frame
0130   dataSet.plotOn(frame, DrawOption("PE1"), Binning(bins),XErrorSize(0), DataError(RooAbsData::SumW2));          //plotting the raw unbinned data on the fram
0131
0132   // Create a new frame to draw the pull distribution and add the distribution to the frame
0133   RooHist* pull = frame->pullHist();
0134   RooPlot* frame2 = mass.frame(Title(""));
0135   frame2->addPlotable(pull,"PE1");
0136
0137   TCanvas* c = new TCanvas("massFitCanvas", "massFitCanvas",800, 800);
0138
0139   TPad mainPad("mainPad", "mainPad", 0., 0., 1., 1.);
0140   mainPad.Draw();
0141
0142   mainPad.cd();
0143
0144   frame->SetMarkerStyle(kCircle);
0145   frame->SetMarkerSize(0.02);
0146   frame->SetLineWidth(1);
0147   frame->GetXaxis()->SetTitle(xAxisTitle.c_str());
0148   frame->GetXaxis()->SetNdivisions(505);
0149   frame->GetYaxis()->SetNdivisions(505);
0150   frame->GetYaxis()->SetTitleFont(42);
0151   frame->GetYaxis()->SetLabelFont(42);
0152   frame->GetYaxis()->SetRangeUser(181, 580);
0153   float binWidth = 1000.*(maxMass - minMass) / nBins;
0154   string yAxisTitle = "Candidates / (" + to_string_with_precision(binWidth, 0) + " MeV)";
0155   frame->GetYaxis()->SetTitle(yAxisTitle.c_str());
0156   frame->Draw();
0157
0158   //Print signal yield
0159   string fitSigYield = to_string_with_precision(dataSet.numEntries() * fSig.getValV(), 0);
0160   string fitSigError = to_string_with_precision(dataSet.numEntries() * fSig.getError(), 0);
0161   string fitYield = "Yield = " + fitSigYield + "#pm " + fitSigError;
0162   double significance = fSig.getValV()/fSig.getError();
0163   string fitSignificance = "Significance = " + to_string_with_precision(significance, 1);
0164   int nEvents = dataCloneTree->GetEntries();
0165   cout << "*\n*\n* nEvents = " << nEvents << "\n* nSig = " << fitSigYield << " +/- " << fitSigError << endl;
0166   double chiSq = frame->chiSquare("model", "data", 4);
0167   cout << "* Fit chisq/ndf = " << chiSq << "\n*\n*" << endl;
0168
0169   //Mean and width
0170   string fitMeanVal = to_string_with_precision(1e3 * mean.getValV(), 1);
0171   string fitMeanErr = to_string_with_precision(1e3 * mean.getError(), 1);
0172   string fitMean = "#mu = " + fitMeanVal + "#pm " + fitMeanErr + " MeV";
0173   string fitWidthVal = to_string_with_precision(1e3 * sigma.getValV(), 1);
0174   string fitWidthErr = to_string_with_precision(1e3 * sigma.getError(), 1);
0175   string fitWidth = "#sigma = " + fitWidthVal + "#pm " + fitWidthErr + " MeV";
0176
0177   TPaveText *pt;
0178   pt = new TPaveText(0.19,0.20,0.52,0.45, "NDC"); //For 4 GeV
0179   pt->SetFillColor(0);
0180   pt->SetFillStyle(0);
0181   pt->SetTextFont(42);
0182   TText *pt_LaTex = pt->AddText("#it{#bf{sPHENIX}} Preliminary");
0183   pt_LaTex = pt->AddText("#it{p}+#it{p} #sqrt{s} = 200 GeV");
0184   pt_LaTex = pt->AddText("Early Calibration");
0185   pt_LaTex = pt->AddText(fitMean.c_str());
0186   pt_LaTex = pt->AddText(fitWidth.c_str());
0187   pt_LaTex = pt->AddText(fitYield.c_str());
0188   //pt_LaTex = pt->AddText(fitSignificance.c_str());
0189   pt->SetBorderSize(0);
0190   pt->DrawClone();
0191
0192   //date timestamp
0193   TPaveText *ptDate;
0194   ptDate = new TPaveText(0.67,0.92,1.0,0.97, "NDC");
0195   ptDate->SetFillColor(0);
0196   ptDate->SetFillStyle(0);
0197   ptDate->SetTextFont(42);
0198   std::string compilation_date = getDate();
0199   TText *pt_LaTexDate = ptDate->AddText(compilation_date.c_str());
0200   ptDate->SetBorderSize(0);
0201   ptDate->Draw();
0202   gPad->Modified();
0203
0204   gPad->SetTopMargin(0.08);
0205
0206   //frame->Print();
0207   TLegend *legend = new TLegend(0.53,0.68,0.83,0.88);
0208   legend->AddEntry(frame->findObject("h_D0_mass"),"Data","PE2");
0209   legend->AddEntry(frame->findObject("model_Norm[D0_mass]"),"Fit","L");
0210   legend->AddEntry(frame->findObject("model_Norm[D0_mass]_Comp[signal,background]"),"D^{0}#rightarrow K^{-}#pi^{+} + c.c.","f");
0211   legend->AddEntry(frame->findObject("model_Norm[D0_mass]_Comp[background]"),"Comb. Bkg.","f");
0212   legend->SetFillColor(0);
0213   legend->SetFillStyle(0);
0214   legend->SetBorderSize(0);
0215   legend->SetTextSize(0.05);
0216   legend->Draw();
0217
0218   gPad->Modified();
0219
0220   vector<string> extensions = {".C", ".pdf", ".png"};
0221   for (auto &extension : extensions)
0222   {
0223     string outputFile = saveName + "_" + dataType + extension;
0224     c->SaveAs(outputFile.c_str());
0225   }
0226 }