QM25_productions/K_S0/K_S0_Fitter.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 <RooHist.h>
0014 #include <RooPlot.h>
0015 #include <RooExponential.h>
0016 #include <RooAddPdf.h>
0017 #include <RooStats/SPlot.h>
0018 #include <TH1F.h>
0019 #include <TFile.h>
0020 #include <TTree.h>
0021 #include <TCut.h>
0022 #include <TChain.h>
0023 #include <TMath.h>
0024 #include <TPad.h>
0025 #include <TCanvas.h>
0026 #include <TBranch.h>
0027
0028 using namespace RooFit;
0029 using namespace std;
0030
0031 //https://stackoverflow.com/questions/997946/how-to-get-current-time-and-date-in-c
0032 std::string getDate()
0033 {
0034     std::time_t t = std::time(0);   // get time now
0035     std::tm* now = std::localtime(&t);
0036
0037     std::stringstream date;
0038     date << (now->tm_mon + 1) << '/'
0039          <<  now->tm_mday << '/'
0040          << (now->tm_year + 1900);
0041     return date.str();
0042 }
0043
0044 template <typename T>
0045 string to_string_with_precision(const T a_value, const int n = 2)
0046 {
0047     ostringstream out;
0048     out.precision(n);
0049     out << fixed << a_value;
0050     return out.str();
0051 }
0052
0053 void K_S0_Fitter()
0054 {
0055   int maxDigits = 3;
0056   TGaxis::SetMaxDigits(maxDigits);
0057
0058   string path, inputFile, plotTitle, saveName, xAxisTitle, branch, dataType;
0059   double minMass = 0.400;
0060   double maxMass = 0.600;
0061
0062   path = "./";
0063   inputFile = "../pipi_reco/outputKFParticle_pipi_reco_00053877_00000_00000.root";
0064   dataType = "pp";
0065
0066   branch = "K_S0_mass";
0067   plotTitle = "";
0068   saveName = path + "testFit_Kshort";
0069   xAxisTitle = "m(#pi^{+}#pi^{-}) [GeV]";
0070
0071   stringstream cutStream;
0072   cutStream << minMass << " <= " << branch << " && " << branch << " <= " << maxMass;
0073   TCut masscut = cutStream.str().c_str();
0074
0075   /*
0076    * Get files and data sets
0077    */
0078
0079   TFile* dataFile = new TFile(inputFile.c_str());
0080   TTree* dataTree = (TTree*)dataFile->Get("DecayTree");
0081
0082   string datasWeight =  inputFile.substr(0, inputFile.size()-5) + "_weighted.root";
0083   TFile* sWeightedDataFile  = new TFile(datasWeight.c_str(), "RECREATE");
0084   TTree* dataSWTree = dataTree->CopyTree(masscut);
0085   TTree* sWeightedDataTree = dataSWTree->CloneTree(-1);
0086
0087   RooRealVar mass(branch.c_str(), "mass", minMass, maxMass);
0088   RooDataSet dataSet(branch.c_str(), "data", mass, Import(*sWeightedDataTree));
0089
0090   /*
0091    * Signal Model
0092    */
0093
0094   RooRealVar  mean("K_S0_mean", "mean", 0.5, 0.46, 0.54);
0095
0096   //Base Gaussian model
0097   RooRealVar  sigma("sigma", "sigma", 0.01, 0.001, 0.02);
0098   RooGaussian signal("signal", "signal", mass, mean, sigma);
0099
0100   RooRealVar  fSig("fSig", "fSig", 0.*dataSet.numEntries(), 0.1, 1*dataSet.numEntries());
0101
0102   /*
0103    * Background Model
0104    */
0105
0106   RooRealVar expConst("expConst", "expConst", -1., -1.e4, 0.);         //creating the decay constant
0107   RooExponential background("background", "background", mass, expConst);          //creating the exponential background signal
0108
0109   /*
0110    * Fitting to the data
0111    */
0112
0113   RooArgList fitModelList(signal, background), fitFracList(fSig);
0114   RooAddPdf model("model", "model", fitModelList, fitFracList);
0115   model.fitTo(dataSet);
0116
0117   RooPlot* frame = mass.frame(Title(plotTitle.c_str()));       //creating the frame
0118
0119   RooBinning bins(minMass, maxMass);
0120   int nBins = 25;
0121   bins.addUniform(nBins, minMass, maxMass);
0122
0123   dataSet.plotOn(frame, Binning(bins), XErrorSize(0), DataError(RooAbsData::SumW2));          //plotting the raw unbinned data on the fram
0124   model.plotOn(frame, Components(background), LineColor(kGray), DrawOption("F"), FillColor(kGray));
0125   model.plotOn(frame, Components(RooArgSet(signal, background)), LineColor(kAzure+8), DrawOption("F"), FillColor(kAzure+8), MoveToBack());
0126   model.plotOn(frame, LineColor(kBlack));        //plotting the fit onto the frame
0127   dataSet.plotOn(frame, DrawOption("PE1"), Binning(bins),XErrorSize(0), DataError(RooAbsData::SumW2));          //plotting the raw unbinned data on the fram
0128
0129   // Create a new frame to draw the pull distribution and add the distribution to the frame
0130   RooHist* pull = frame->pullHist();
0131   RooPlot* frame2 = mass.frame(Title(""));
0132   frame2->addPlotable(pull,"PE1");
0133
0134   TCanvas* c = new TCanvas("massFitCanvas", "massFitCanvas",800, 800);
0135
0136   TPad mainPad("mainPad", "mainPad", 0., 0.3, 1., 1.);
0137   mainPad.SetBottomMargin(0);
0138   mainPad.Draw();
0139
0140   TPad pullPad("pullPad", "pullPad", 0., 0.0, 1., 0.3);
0141   pullPad.SetBottomMargin(0.5);
0142   pullPad.SetTopMargin(0);
0143
0144   pullPad.Draw();
0145
0146   mainPad.cd();
0147
0148   frame->SetMarkerStyle(kCircle);
0149   frame->SetMarkerSize(0.02);
0150   frame->SetLineWidth(1);
0151   frame->GetXaxis()->SetTitleSize(0);
0152   frame->GetXaxis()->SetLabelSize(0);
0153   frame->GetYaxis()->SetTitleFont(42);
0154   frame->GetYaxis()->SetLabelFont(42);
0155   frame->GetYaxis()->SetRangeUser(1, 50);
0156   float binWidth = 1000.*(maxMass - minMass) / nBins;
0157   string yAxisTitle = "Candidates / (" + to_string_with_precision(binWidth, 0) + " MeV)";
0158   frame->GetYaxis()->SetTitle(yAxisTitle.c_str());
0159   frame->Draw();
0160
0161   //Print signal yield
0162   string fittedSigYield = to_string_with_precision(dataSet.numEntries() * fSig.getValV(), 0);
0163   string fittedSigError = to_string_with_precision(dataSet.numEntries() * fSig.getError(), 0);
0164   //int nEvents = dataTree->GetEntries();
0165   int nEvents = sWeightedDataTree->GetEntries();
0166   cout << "*\n*\n* nEvents = " << nEvents << "\n* nSig = " << fittedSigYield << " +/- " << fittedSigError << endl;
0167   double chiSq = frame->chiSquare("model", "data", 4);
0168   cout << "* Fit chisq/ndf = " << chiSq << "\n*\n*" << endl;
0169
0170   TPaveText *pt;
0171   pt = new TPaveText(0.59,0.40,0.89,0.56, "NDC");
0172   pt->SetFillColor(0);
0173   pt->SetFillStyle(0);
0174   pt->SetTextFont(42);
0175   TText *pt_LaTex = pt->AddText("#it{#bf{sPHENIX}} Internal");
0176   pt_LaTex = pt->AddText("#it{p}+#it{p} #sqrt{s} = 200 GeV");
0177   pt->SetBorderSize(0);
0178   pt->DrawClone();
0179
0180   //date timestamp
0181   TPaveText *ptDate;
0182   ptDate = new TPaveText(0.67,0.92,1.05,1.00, "NDC");
0183   ptDate->SetFillColor(0);
0184   ptDate->SetFillStyle(0);
0185   ptDate->SetTextFont(42);
0186   std::string compilation_date = getDate();
0187   TText *pt_LaTexDate = ptDate->AddText(compilation_date.c_str());
0188   ptDate->SetBorderSize(0);
0189   ptDate->Draw();
0190   gPad->Modified();
0191
0192   gPad->SetTopMargin(0.08);
0193
0194   //frame->Print();
0195   TLegend *legend = new TLegend(0.57,0.60,0.85,0.85);
0196   legend->AddEntry(frame->findObject("h_K_S0_mass"),"Data","PE2");
0197   legend->AddEntry(frame->findObject("model_Norm[K_S0_mass]"),"Fit","L");
0198   legend->AddEntry(frame->findObject("model_Norm[K_S0_mass]_Comp[signal,background]"),"K_{S}^{0}#rightarrow#pi^{+}#pi^{-}","f");
0199   legend->AddEntry(frame->findObject("model_Norm[K_S0_mass]_Comp[background]"),"Comb. Bkg.","f");
0200   legend->SetFillColor(0);
0201   legend->SetFillStyle(0);
0202   legend->SetBorderSize(0);
0203   legend->SetTextSize(0.05);
0204   legend->Draw();
0205
0206   gPad->Modified();
0207   pullPad.cd();
0208
0209   frame2->SetMarkerStyle(kCircle);
0210   frame2->SetMarkerSize(0.02);
0211   frame2->SetTitle("");
0212   frame2->GetXaxis()->SetNdivisions(505);
0213   frame2->GetXaxis()->SetTitle(xAxisTitle.c_str());
0214   frame2->GetXaxis()->SetTitleOffset(0.9);
0215   frame2->GetXaxis()->SetTitleFont(42);
0216   frame2->GetXaxis()->SetTitleSize(0.12);
0217   frame2->GetXaxis()->SetLabelFont(42);
0218   frame2->GetXaxis()->SetLabelSize(0.12);
0219   frame2->GetYaxis()->SetTitle("Pull");
0220   frame2->GetYaxis()->SetTitleOffset(0.55);
0221   frame2->GetYaxis()->SetTitleFont(42);
0222   frame2->GetYaxis()->SetTitleSize(0.12);
0223   frame2->GetYaxis()->SetLabelFont(42);
0224   frame2->GetYaxis()->SetLabelSize(0.12);
0225   frame2->GetYaxis()->SetRangeUser(-6, 6);
0226   frame2->GetYaxis()->SetNdivisions(5);
0227   TF1 *plusThreeLine = new TF1("plusThreeLine",  "pol1", minMass, maxMass);
0228   plusThreeLine->SetParameters(3, 0);
0229   plusThreeLine->SetLineColor(1);
0230   plusThreeLine->SetLineStyle(2);
0231   plusThreeLine->SetLineWidth(2);
0232   TF1 *zeroLine = new TF1("zeroLine",  "pol1", minMass, maxMass);
0233   zeroLine->SetParameters(0, 0);
0234   zeroLine->SetLineColor(1);
0235   zeroLine->SetLineStyle(2);
0236   zeroLine->SetLineWidth(2);
0237   TF1 *minusThreeLine = new TF1("minusThreeLine", "pol1", minMass, maxMass);
0238   minusThreeLine->SetParameters(-3, 0);
0239   minusThreeLine->SetLineColor(1);
0240   minusThreeLine->SetLineStyle(2);
0241   minusThreeLine->SetLineWidth(2);
0242   frame2->Draw();
0243   plusThreeLine->Draw("same");
0244   zeroLine->Draw("same");
0245   minusThreeLine->Draw("same");
0246
0247   vector<string> extensions = {".C", ".pdf", ".png"};
0248   for (auto &extension : extensions)
0249   {
0250     string outputFile = saveName + "_" + dataType + extension;
0251     c->SaveAs(outputFile.c_str());
0252   }
0253 }