QM25_productions/D0/D0_fitter_QM25.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()
0055 {
0056   //int changeType[] = {53534, 53631, 53687, 53716, 53738, 53739, 53741, 53742, 53743, 53744, 53756, 53783, 53871, 53876, 53877, 53879};
0057   string changeType[] = {"optimized"};
0058
0059   int maxDigits = 3;
0060   TGaxis::SetMaxDigits(maxDigits);
0061
0062   string path, inputFile, plotTitle, saveName, xAxisTitle, branch, dataType;
0063   double minMass = 1.7;
0064   double maxMass = 2.0;
0065
0066   path = "/sphenix/user/cdean/public/QM25/";
0067
0068   for (auto change : changeType)
0069   {
0070     inputFile = path + "output_D0_oldKFP_ana475p017_optimized.root";
0071     dataType = "pp";
0072
0073     branch = "D0_mass";
0074     plotTitle = "";
0075     saveName = path + "fits/D02Kpi_ana475p01_optimized";
0076     xAxisTitle = "m(K^{#mp}#pi^{#pm}) [GeV]";
0077
0078     stringstream cutStream;
0079     cutStream << minMass << " <= " << branch << " && " << branch << " <= " << maxMass;
0080     TCut masscut = cutStream.str().c_str();
0081
0082     /*
0083      * Get files and data sets
0084      */
0085
0086     TFile* dataFile = new TFile(inputFile.c_str());
0087     TTree* dataTree = (TTree*)dataFile->Get("DecayTree");
0088
0089     string datasWeight =  inputFile.substr(0, inputFile.size()-5) + "_weighted.root";
0090     TFile* sWeightedDataFile  = new TFile(datasWeight.c_str(), "RECREATE");
0091     TTree* dataCopyTree = dataTree->CopyTree(masscut);
0092     TTree* dataCloneTree = dataCopyTree->CloneTree(-1);
0093
0094     RooRealVar mass(branch.c_str(), "mass", minMass, maxMass);
0095     RooDataSet dataSet(branch.c_str(), "data", mass, Import(*dataCloneTree));
0096
0097     /*
0098      * Signal Model
0099      */
0100
0101     RooRealVar  mean("D0_mean", "mean", 1.864, 1.855, 1.875);
0102
0103     //Base Gaussian model
0104     RooRealVar  sigma("sigma", "sigma", 0.01, 0.000, 0.02);
0105     RooGaussian signal("signal", "signal", mass, mean, sigma);
0106
0107     RooRealVar  fSig("fSig", "fSig", 0.*dataSet.numEntries(), 0., 0.1*dataSet.numEntries());
0108     RooRealVar  fBkg("fBkg", "fBkg", 0.*dataSet.numEntries(), 0., 1*dataSet.numEntries());
0109
0110     /*
0111      * Background Model
0112      */
0113
0114     RooRealVar expConstOne("expConstOne", "expConstOne", -10, -1e2, 0.);
0115     RooExponential background("background", "background", mass, expConstOne);
0116
0117     /*
0118      * Fitting to the data
0119      */
0120
0121     RooArgList fitModelList(signal, background), fitFracList(fSig);
0122     RooAddPdf model("model", "model", fitModelList, fitFracList);
0123     model.fitTo(dataSet);
0124
0125     RooPlot* frame = mass.frame(Title(plotTitle.c_str()));       //creating the frame
0126
0127     RooBinning bins(minMass, maxMass);
0128     int nBins = 25;
0129     bins.addUniform(nBins, minMass, maxMass);
0130
0131     dataSet.plotOn(frame, Binning(bins), XErrorSize(0), DataError(RooAbsData::SumW2));          //plotting the raw unbinned data on the fram
0132     model.plotOn(frame, Components(background), LineColor(kGray), DrawOption("F"), FillColor(kGray));
0133     model.plotOn(frame, Components(RooArgSet(signal, background)), LineColor(kAzure+8), DrawOption("F"), FillColor(kAzure+8), MoveToBack());
0134     model.plotOn(frame, LineColor(kBlack));        //plotting the fit onto the frame
0135     dataSet.plotOn(frame, DrawOption("PE1"), Binning(bins),XErrorSize(0), DataError(RooAbsData::SumW2));          //plotting the raw unbinned data on the fram
0136
0137     // Create a new frame to draw the pull distribution and add the distribution to the frame
0138     RooHist* pull = frame->pullHist();
0139     RooPlot* frame2 = mass.frame(Title(""));
0140     frame2->addPlotable(pull,"PE1");
0141
0142     TCanvas* c = new TCanvas("massFitCanvas", "massFitCanvas",800, 800);
0143
0144     TPad mainPad("mainPad", "mainPad", 0., 0.3, 1., 1.);
0145     mainPad.SetBottomMargin(0);
0146     mainPad.Draw();
0147
0148     TPad pullPad("pullPad", "pullPad", 0., 0.0, 1., 0.3);
0149     pullPad.SetBottomMargin(0.5);
0150     pullPad.SetTopMargin(0);
0151
0152     pullPad.Draw();
0153
0154     mainPad.cd();
0155
0156     frame->SetMarkerStyle(kCircle);
0157     frame->SetMarkerSize(0.02);
0158     frame->SetLineWidth(1);
0159     frame->GetXaxis()->SetTitleSize(0);
0160     frame->GetXaxis()->SetLabelSize(0);
0161     frame->GetYaxis()->SetTitleFont(42);
0162     frame->GetYaxis()->SetLabelFont(42);
0163     frame->GetYaxis()->SetRangeUser(181, 580);
0164     float binWidth = 1000.*(maxMass - minMass) / nBins;
0165     string yAxisTitle = "Candidates / (" + to_string_with_precision(binWidth, 0) + " MeV)";
0166     frame->GetYaxis()->SetTitle(yAxisTitle.c_str());
0167     frame->Draw();
0168
0169     //Print signal yield
0170     string fitSigYield = to_string_with_precision(dataSet.numEntries() * fSig.getValV(), 0);
0171     string fitSigError = to_string_with_precision(dataSet.numEntries() * fSig.getError(), 0);
0172     string fitYield = "Yield = " + fitSigYield + "#pm " + fitSigError;
0173     double significance = fSig.getValV()/fSig.getError();
0174     string fitSignificance = "Significance = " + to_string_with_precision(significance, 1);
0175     int nEvents = dataCloneTree->GetEntries();
0176     cout << "*\n*\n* nEvents = " << nEvents << "\n* nSig = " << fitSigYield << " +/- " << fitSigError << endl;
0177     double chiSq = frame->chiSquare("model", "data", 4);
0178     cout << "* Fit chisq/ndf = " << chiSq << "\n*\n*" << endl;
0179
0180     //Mean and width
0181     string fitMeanVal = to_string_with_precision(1e3 * mean.getValV(), 1);
0182     string fitMeanErr = to_string_with_precision(1e3 * mean.getError(), 1);
0183     string fitMean = "#mu = " + fitMeanVal + "#pm " + fitMeanErr + " MeV";
0184     string fitWidthVal = to_string_with_precision(1e3 * sigma.getValV(), 1);
0185     string fitWidthErr = to_string_with_precision(1e3 * sigma.getError(), 1);
0186     string fitWidth = "#sigma = " + fitWidthVal + "#pm " + fitWidthErr + " MeV";
0187
0188     TPaveText *pt;
0189     pt = new TPaveText(0.14,0.05,0.52,0.35, "NDC"); //For 4 GeV
0190     pt->SetFillColor(0);
0191     pt->SetFillStyle(0);
0192     pt->SetTextFont(42);
0193     TText *pt_LaTex = pt->AddText("#it{#bf{sPHENIX}} Preliminary");
0194     pt_LaTex = pt->AddText("#it{p}+#it{p} #sqrt{s} = 200 GeV");
0195     pt_LaTex = pt->AddText("Early Calibration");
0196     pt_LaTex = pt->AddText(fitMean.c_str());
0197     pt_LaTex = pt->AddText(fitWidth.c_str());
0198     pt_LaTex = pt->AddText(fitYield.c_str());
0199     //pt_LaTex = pt->AddText(fitSignificance.c_str());
0200     pt->SetBorderSize(0);
0201     pt->DrawClone();
0202
0203     //date timestamp
0204     TPaveText *ptDate;
0205     ptDate = new TPaveText(0.67,0.92,1.05,1.00, "NDC");
0206     ptDate->SetFillColor(0);
0207     ptDate->SetFillStyle(0);
0208     ptDate->SetTextFont(42);
0209     std::string compilation_date = getDate();
0210     TText *pt_LaTexDate = ptDate->AddText(compilation_date.c_str());
0211     ptDate->SetBorderSize(0);
0212     ptDate->Draw();
0213     gPad->Modified();
0214
0215     gPad->SetTopMargin(0.08);
0216
0217     //frame->Print();
0218     TLegend *legend = new TLegend(0.62,0.60,0.90,0.85);
0219     legend->AddEntry(frame->findObject("h_D0_mass"),"Data","PE2");
0220     legend->AddEntry(frame->findObject("model_Norm[D0_mass]"),"Fit","L");
0221     legend->AddEntry(frame->findObject("model_Norm[D0_mass]_Comp[signal,background]"),"D^{0}#rightarrow K^{-}#pi^{+} + c.c.","f");
0222     legend->AddEntry(frame->findObject("model_Norm[D0_mass]_Comp[background]"),"Comb. Bkg.","f");
0223     legend->SetFillColor(0);
0224     legend->SetFillStyle(0);
0225     legend->SetBorderSize(0);
0226     legend->SetTextSize(0.05);
0227     legend->Draw();
0228
0229     gPad->Modified();
0230     pullPad.cd();
0231
0232     frame2->SetMarkerStyle(kCircle);
0233     frame2->SetMarkerSize(0.02);
0234     frame2->SetTitle("");
0235     frame2->GetXaxis()->SetNdivisions(505);
0236     frame2->GetXaxis()->SetTitle(xAxisTitle.c_str());
0237     frame2->GetXaxis()->SetTitleOffset(0.9);
0238     frame2->GetXaxis()->SetTitleFont(42);
0239     frame2->GetXaxis()->SetTitleSize(0.12);
0240     frame2->GetXaxis()->SetLabelFont(42);
0241     frame2->GetXaxis()->SetLabelSize(0.12);
0242     frame2->GetYaxis()->SetTitle("Pull");
0243     frame2->GetYaxis()->SetTitleOffset(0.55);
0244     frame2->GetYaxis()->SetTitleFont(42);
0245     frame2->GetYaxis()->SetTitleSize(0.12);
0246     frame2->GetYaxis()->SetLabelFont(42);
0247     frame2->GetYaxis()->SetLabelSize(0.12);
0248     frame2->GetYaxis()->SetRangeUser(-6, 6);
0249     frame2->GetYaxis()->SetNdivisions(5);
0250     TF1 *plusThreeLine = new TF1("plusThreeLine",  "pol1", minMass, maxMass);
0251     plusThreeLine->SetParameters(3, 0);
0252     plusThreeLine->SetLineColor(1);
0253     plusThreeLine->SetLineStyle(2);
0254     plusThreeLine->SetLineWidth(2);
0255     TF1 *zeroLine = new TF1("zeroLine",  "pol1", minMass, maxMass);
0256     zeroLine->SetParameters(0, 0);
0257     zeroLine->SetLineColor(1);
0258     zeroLine->SetLineStyle(2);
0259     zeroLine->SetLineWidth(2);
0260     TF1 *minusThreeLine = new TF1("minusThreeLine", "pol1", minMass, maxMass);
0261     minusThreeLine->SetParameters(-3, 0);
0262     minusThreeLine->SetLineColor(1);
0263     minusThreeLine->SetLineStyle(2);
0264     minusThreeLine->SetLineWidth(2);
0265     frame2->Draw();
0266     plusThreeLine->Draw("same");
0267     zeroLine->Draw("same");
0268     minusThreeLine->Draw("same");
0269
0270     vector<string> extensions = {".C", ".pdf", ".png"};
0271     for (auto &extension : extensions)
0272     {
0273       string outputFile = saveName + "_" + dataType + extension;
0274       c->SaveAs(outputFile.c_str());
0275     }
0276   }
0277 }