tpc/TpcDVCalib/fit.C

0001 #include "TFile.h"
0002 #include "TTree.h"
0003 #include "RooRealVar.h"
0004 #include "RooDataSet.h"
0005 #include "RooGaussian.h"
0006 #include "RooFitResult.h"
0007 #include "RooPlot.h"
0008 #include "TCanvas.h"
0009
0010 #include <stdio.h>
0011 #include <iostream>
0012
0013 #include <cdbobjects/CDBTTree.h>
0014
0015 // cppcheck-suppress unknownMacro
0016 R__LOAD_LIBRARY(libcdbobjects.so)
0017
0018 std::vector<float> vec_peak_pos_val;
0019 std::vector<float> vec_peak_pos_err;
0020
0021 int fit_dz(const char* filename = "data.root",
0022     const char* treename = "tree",
0023     const char* branchname = "dz",
0024     TString cut = "1",
0025     TString outfile = "fit_dz.pdf",
0026     TString title = "Run 51103, fit dz",
0027     int runnumber = 51103)
0028 {
0029     TChain* chain = new TChain(treename);
0030     chain->Add(filename);
0031     int nevent = chain->GetEntries();
0032
0033     if (nevent > 0) {
0034         std::cout << "Successfully added file to the TChain." << std::endl;
0035     } else {
0036         std::cout << "Failed to add file to the TChain." << std::endl;
0037         return -1;
0038     }
0039
0040     if (chain->GetEntries(cut) < 50)
0041     {
0042       std::cout<<"Cut: "<<cut<<" , nevent used in fit = "<<chain->GetEntries(cut)<<" . Too Low!!! Fit performance maybe not be good, skipping!!!"<<std::endl;
0043       return -1;
0044     }
0045
0046     double xmin = -20;
0047     double xmax = 20;
0048     int nbins = 50;
0049     RooRealVar x(branchname, branchname, nbins, xmin, xmax);
0050     double binwidth = (xmax-xmin)/nbins;
0051
0052     TFile* newfile = new TFile("tmp.root","recreate");
0053     TTree* newtree = (TTree*) chain->CopyTree(cut);
0054
0055     // Create a RooDataSet from the TTree
0056     RooDataSet data("data", "dataset of dz", newtree, x);
0057
0058     // Define the parameters for the Gaussian
0059     RooRealVar mean("mean", "mean of gaussian", data.mean(x), xmin, xmax);
0060     RooRealVar sigma("sigma", "width of gaussian", 0.1*data.sigma(x), 0.001, (xmax-xmin)/2);
0061
0062     // Create the Gaussian PDF
0063     RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma);
0064     RooRealVar ngauss("ngauss","number of gaussian",0.2*data.sumEntries(),0,data.sumEntries());
0065
0066     // Define the parameters for the linear background
0067     RooRealVar a0("a0", "constant", 0, -10, 10);
0068     RooRealVar a1("a1", "slope", 0, -1, 1);
0069     RooPolynomial poly("poly", "linear background", x, RooArgList(a1, a0));
0070     RooRealVar npoly("npoly","number of poly",0.3*data.sumEntries(),0,data.sumEntries());
0071
0072     // Combine the Gaussian and the polynomial into a single PDF
0073     RooAddPdf model("model", "gauss + poly", RooArgList(gauss, poly), RooArgList(ngauss, npoly));
0074
0075     // Fit the Gaussian to the data
0076     RooFitResult* fit_result = model.fitTo(data, RooFit::Save());
0077
0078     // Print the fit result
0079     fit_result->Print();
0080
0081     // Plot the data and the fit result
0082     RooPlot* frame = x.frame();
0083     data.plotOn(frame);
0084     model.plotOn(frame);
0085     model.plotOn(frame, RooFit::Components(gauss), RooFit::LineColor(kRed), RooFit::LineStyle(kDashed));
0086     model.plotOn(frame, RooFit::Components(poly), RooFit::LineColor(kGreen), RooFit::LineStyle(kDashed));
0087
0088     frame->SetTitle(title);
0089     frame->SetXTitle("#DeltaZ [cm]");
0090     frame->SetYTitle(Form("Events / (%1.0g cm)",binwidth));
0091     frame->GetXaxis()->SetTitleSize(0.05);
0092     frame->GetYaxis()->SetTitleSize(0.05);
0093
0094     double meanValue = mean.getVal();
0095     double meanError = mean.getError();
0096     double sigmaValue = sigma.getVal();
0097     double sigmaError = sigma.getError();
0098
0099     vec_peak_pos_val.push_back(meanValue);
0100     vec_peak_pos_err.push_back(meanError);
0101
0102 /*
0103     TPaveText *pt = new TPaveText(0.20, 0.73, 0.35, 0.88, "NDC");
0104     pt->SetFillColor(0);
0105     pt->SetTextAlign(12);
0106     pt->AddText(Form("#mu = %.2f #pm %.2f", meanValue, meanError));
0107     pt->AddText(Form("#sigma = %.2f #pm %.2f", sigmaValue, sigmaError));
0108
0109     // Draw the plot on a canvas
0110     TCanvas *canvas = new TCanvas("canvas", "fit result", 800, 600);
0111     canvas->SetLeftMargin(0.15);
0112     canvas->SetRightMargin(0.05);
0113     canvas->SetTopMargin(0.1);
0114     canvas->SetBottomMargin(0.15);
0115     frame->Draw();
0116     pt->Draw();
0117
0118     std::string dirName = "figure/" + std::to_string(runnumber);
0119     std::string command = "ls " + dirName + " 2>/dev/null";
0120     int result = system(command.c_str());
0121
0122     if (result != 0) {
0123         std::string mkdirCommand = "mkdir -p " + std::string(dirName);
0124         int mkdirResult = system(mkdirCommand.c_str());
0125
0126         if (mkdirResult == 0) {
0127             std::cout << "Directory '" << dirName << "' created successfully." << std::endl;
0128         } else {
0129             std::cerr << "Failed to create directory '" << dirName << "'." << std::endl;
0130         }
0131     } else {
0132         std::cout << "Directory '" << dirName << "' already exists." << std::endl;
0133     }
0134
0135     canvas->SaveAs(outfile);
0136     delete canvas;
0137 */
0138
0139     delete chain;
0140
0141     system("rm tmp.root");
0142     return 1;
0143 }
0144
0145 void CDBTTreeWrite(const std::string &fname = "test.root", float dv=0.007)
0146 {
0147   CDBTTree *cdbttree = new CDBTTree(fname);
0148   cdbttree->SetSingleFloatValue("tpc_drift_velocity",dv); // unit: cm/ns
0149   cdbttree->CommitSingle();
0150   cdbttree->Print();
0151   cdbttree->WriteCDBTTree();
0152   delete cdbttree;
0153   //gSystem->Exit(0);
0154 }
0155
0156 void fit(int runnumber)
0157 {
0158     gStyle->SetOptStat(0);
0159     std::vector<float> vec_trkr_z_val;
0160     std::vector<float> vec_trkr_z_err;
0161
0162     // fit dz in different trkr_z region
0163     vec_peak_pos_val.clear();
0164     vec_peak_pos_err.clear();
0165     vec_trkr_z_val.clear();
0166     vec_trkr_z_err.clear();
0167     float z_min = -140;
0168     float z_max = 140;
0169     float z_range = z_max - z_min;
0170     int nstep = 14;
0171     float z_stepsize = z_range / nstep;
0172     for (int i=0; i<nstep; i++)
0173     {
0174       float z_cut_min = z_min + i * z_stepsize;
0175       float z_cut_max = z_min + (i+1) * z_stepsize;
0176       int status = fit_dz(Form("Reconstructed/%d/TrackCalo_*_ana.root",runnumber), "tree", "dz", Form("fabs(dphi)<0.1 && track_z>%f && track_z<%f && fabs(dz)<20",z_cut_min,z_cut_max), Form("figure/%d/dz_fit_cuttrkrz_%d.pdf",runnumber,i), Form("Run %d, Track Z#in[%.0f,%.0f] cm",runnumber,z_cut_min,z_cut_max),runnumber);
0177       if (status==1)
0178       {
0179         vec_trkr_z_val.push_back((z_cut_min+z_cut_max)/2);
0180         vec_trkr_z_err.push_back((z_cut_max-z_cut_min)/2);
0181       }
0182     }
0183
0184     int npoint = vec_peak_pos_val.size();
0185
0186     if (npoint==0)
0187     {
0188       std::cout<<"Do not have anything to fit!!! DST is empty or Matching for this run is very bad"<<std::endl;
0189       return;
0190     }
0191
0192     int npoint_left=0;
0193     int npoint_right=0;
0194     for (int i=0; i<(vec_peak_pos_val.size()); i++)
0195     {
0196       if (vec_trkr_z_val.at(i)>0)
0197       {
0198         npoint_right++;
0199       }
0200       else if (vec_trkr_z_val.at(i)<0)
0201       {
0202         npoint_left++;
0203       }
0204     }
0205     if (npoint_left<3 || npoint_right<3)
0206     {
0207       std::cout<<"Do not have enough good point to fit!!! Matching for this run is very bad"<<std::endl;
0208       return;
0209     }
0210
0211     TCanvas *c1 = new TCanvas("c1", "Fitting Example", 800, 600);
0212     c1->SetLeftMargin(0.12);
0213     c1->SetRightMargin(0.05);
0214     TGraphErrors *graph = new TGraphErrors(vec_peak_pos_val.size(), vec_trkr_z_val.data(), vec_peak_pos_val.data(), vec_trkr_z_err.data(), vec_peak_pos_err.data());
0215     graph->SetTitle(Form("Run %d, TPC-EMCal matching;Track Z [cm];#DeltaZ=Z_{track}-Z_{calo} Peak [cm]",runnumber));
0216     graph->SetMarkerStyle(21);
0217     graph->SetMarkerColor(kBlue);
0218     graph->SetLineColor(kBlue);
0219     graph->GetYaxis()->SetRangeUser(-50,20);
0220
0221     double intercept_min = -20;
0222     double intercept_max= 20;
0223     double slope_min = -1;
0224     double slope_max = 1;
0225
0226     TF1 *fitFunc1 = new TF1("fitFunc1", "[0] + [1]*x", z_min, -10);
0227     fitFunc1->SetParameters(vec_peak_pos_val.at(npoint/2 - 2), (vec_peak_pos_val.at(npoint/2 - 2) - vec_peak_pos_val.at(0)) / (vec_trkr_z_val.at(npoint/2 - 2) - vec_trkr_z_val.at(0)));
0228     fitFunc1->SetParLimits(0, intercept_min, intercept_max);
0229     fitFunc1->SetParLimits(1, slope_min, slope_max);
0230     graph->Fit(fitFunc1, "R");
0231
0232     TF1 *fitFunc2 = new TF1("fitFunc2", "[0] + [1]*x", 10, z_max);
0233     fitFunc2->SetParameters(vec_peak_pos_val.at(npoint/2 + 1), (vec_peak_pos_val.at(npoint - 1) - vec_peak_pos_val.at(npoint/2 + 1)) / (vec_trkr_z_val.at(npoint - 1) - vec_trkr_z_val.at(npoint/2 + 1)));
0234     fitFunc2->SetParLimits(0, intercept_min, intercept_max);
0235     fitFunc2->SetParLimits(1, slope_min, slope_max);
0236     graph->Fit(fitFunc2, "R");
0237
0238     double chi2_1 = fitFunc1->GetChisquare();
0239     int ndf_1 = fitFunc1->GetNDF();
0240
0241     double chi2_2 = fitFunc2->GetChisquare();
0242     int ndf_2 = fitFunc2->GetNDF();
0243
0244     double p0_1 = fitFunc1->GetParameter(0);
0245     double p1_1 = fitFunc1->GetParameter(1);
0246     double p0_1_err = fitFunc1->GetParError(0);
0247     double p1_1_err = fitFunc1->GetParError(1);
0248     TString txt_fun1 = Form("y = %.3fx + (%.3f)",p1_1,p0_1);
0249     TString txt_quality1 = Form("#chi^2/NDF = %.2f/%d = %.2f",chi2_1,ndf_1,chi2_1/ndf_1);
0250     std::cout << "Fit left part: " << txt_fun1 << " , " << txt_quality1 << std::endl;
0251
0252     double p0_2 = fitFunc2->GetParameter(0);
0253     double p1_2 = fitFunc2->GetParameter(1);
0254     double p0_2_err = fitFunc2->GetParError(0);
0255     double p1_2_err = fitFunc2->GetParError(1);
0256     TString txt_fun2 = Form("y = %.3fx + (%.3f)",p1_2,p0_2);
0257     TString txt_quality2 = Form("#chi^2/NDF = %.2f/%d = %.2f",chi2_2,ndf_2,chi2_2/ndf_2);
0258     std::cout << "Fit right part: " << txt_fun2 << " , " << txt_quality2 << std::endl;
0259
0260     graph->Draw("AP");
0261
0262     fitFunc1->SetLineColor(kRed);
0263     fitFunc1->Draw("SAME");
0264     fitFunc2->SetLineColor(kBlue);
0265     fitFunc2->Draw("SAME");
0266
0267     double driftvelo_pre = 0.00710;
0268     double dz_separation_0_new = -(p1_1+p1_2)/2. * 2 * 105;
0269     double dz_separation_0_err = -(p1_1_err+p1_2_err)/2. * 2 * 105;
0270     double driftvelo_new = driftvelo_pre * (1+ dz_separation_0_new / 2. / 105.);
0271     double driftvelo_err = driftvelo_pre * (fabs(dz_separation_0_err) / 2. / 105.);
0272     TString printtxt1 = Form("DV used in reconstruction: %.5f cm/ns",driftvelo_pre);
0273     TString printtxt2 = Form("Calibrated DV from data: (%.5f#pm%.5f) cm/ns",driftvelo_new,driftvelo_err);
0274     cout<<printtxt1<<endl;
0275     cout<<printtxt2<<endl;
0276
0277     TPaveText *pt1 = new TPaveText(0.15, 0.275, 0.45, 0.40, "NDC");
0278     pt1->SetFillColor(0);
0279     pt1->SetTextAlign(12);
0280     pt1->AddText(txt_fun1);
0281     pt1->AddText(txt_quality1);
0282     pt1->Draw("same");
0283
0284     TPaveText *pt2 = new TPaveText(0.55, 0.275, 0.85, 0.40, "NDC");
0285     pt2->SetFillColor(0);
0286     pt2->SetTextAlign(12);
0287     pt2->AddText(txt_fun2);
0288     pt2->AddText(txt_quality2);
0289     pt2->Draw("same");
0290
0291     TPaveText *pt3 = new TPaveText(0.20, 0.15, 0.90, 0.275, "NDC");
0292     pt3->SetFillColor(0);
0293     pt3->SetTextAlign(12);
0294     pt3->AddText(printtxt1);
0295     pt3->AddText(printtxt2);
0296     pt3->Draw("same");
0297
0298     TPaveText *pt4 = new TPaveText(0.15, 0.40, 0.35, 0.50, "NDC");
0299     pt4->SetFillColor(0);
0300     pt4->SetTextAlign(12);
0301
0302     // checks if result is good to insert into CDB
0303     double fit_quality_thr = 10;
0304     bool fit_status = true;
0305     if (chi2_1/ndf_1>fit_quality_thr || chi2_2/ndf_2>fit_quality_thr)
0306     {
0307       cout<<"ERROR: fail fit quality cut"<<endl;
0308       fit_status = false;
0309     }
0310
0311     if (fabs(p0_1-intercept_min)<1e-3 || fabs(p0_1-intercept_max)<1e-3
0312      || fabs(p1_1-slope_min)<1e-3 || fabs(p1_1-slope_max)<1e-3
0313      || fabs(p0_2-intercept_min)<1e-3 || fabs(p0_2-intercept_max)<1e-3
0314      || fabs(p1_2-slope_min)<1e-3 || fabs(p1_2-slope_max)<1e-3)
0315     {
0316       cout<<"ERROR: parameters in the limit"<<endl;
0317       fit_status = false;
0318     }
0319
0320     if ( fabs(p1_1 - p1_2) > 0.1 )
0321     {
0322       cout<<"ERROR: left fit and right fit inconsistent"<<endl;
0323       fit_status = false;
0324     }
0325
0326     double dv_ref = driftvelo_pre;
0327     if (runnumber > 50800 && runnumber < 51200)
0328     {
0329       dv_ref = 0.00625; // for TPC gas mixture ratio changes during 08/09 - 08/12
0330     }
0331     if (driftvelo_new/dv_ref > 1.1 || driftvelo_new/dv_ref < 0.9)
0332     {
0333       fit_status = false;
0334     }
0335
0336     if (fit_status)
0337     {
0338       pt4->SetTextColor(3);
0339       pt4->AddText("GOOD");
0340       pt4->Draw("same");
0341       c1->Update();
0342       c1->SaveAs(Form("./figure/dz_fit_trkrz_Run%d.pdf",runnumber));
0343       CDBTTreeWrite(Form("cdbttree/tpc_drift_velocity_%d.root",runnumber), (float) driftvelo_new);
0344     }
0345     else
0346     {
0347       pt4->SetTextColor(2);
0348       pt4->AddText("BAD");
0349       pt4->Draw("same");
0350       c1->Update();
0351       c1->SaveAs(Form("./figure_failed/dz_fit_trkrz_Run%d.pdf",runnumber));
0352       CDBTTreeWrite(Form("cdbttree_failed/tpc_drift_velocity_%d.root",runnumber), (float) driftvelo_new);
0353     }
0354 }