calibrations/mbd/calib_t0mean.C

0001 #include "make_cdbtree.C"
0002
0003 #include <TCanvas.h>
0004 #include <TFile.h>
0005 #include <TSpectrum.h>
0006 #include <TH2.h>
0007 #include <TH1.h>
0008 #include <TF1.h>
0009 #include <TStyle.h>
0010
0011 #include <filesystem>
0012 #include <fstream>
0013
0014 TSpectrum *tspec{nullptr};
0015 TF1 *FitHitTime{nullptr};
0016
0017 void FitGausToPeak(TH1 *h1, const Double_t sigma)
0018 {
0019 //  int npeak = tspec->Search(h1, sigma, "goff",0.2);  // finds the highest peak, draws marker
0020
0021   double *peakpos = tspec->GetPositionX();
0022   float centerpeak = peakpos[0];
0023   float rangemin = centerpeak - sigma*1.2;
0024   float rangemax = centerpeak + sigma*1.2;
0025
0026   /*
0027   std::cout << "NPEAKS " << npeak << std::endl;
0028   if (npeak < 3)
0029   {
0030     h1->Fit(FitHitTime,"QN0L");
0031     rangemin = FitHitTime->GetParameter(1) - 0.6*FitHitTime->GetParameter(2);
0032     rangemax = FitHitTime->GetParameter(1) + 0.6*FitHitTime->GetParameter(2);
0033   }
0034   else
0035   {
0036     float centerpeak = peakpos[0];
0037     float sidepeak[2];
0038     if (peakpos[2] > peakpos[1])
0039     {
0040       sidepeak[0] = peakpos[1];
0041       sidepeak[1] = peakpos[2];
0042     }
0043     else
0044     {
0045       sidepeak[1] = peakpos[1];
0046       sidepeak[0] = peakpos[2];
0047     }
0048     rangemin = centerpeak - (centerpeak - sidepeak[0]) / 2.;
0049     rangemax = centerpeak + (sidepeak[1] - centerpeak) / 2.;
0050   }
0051   */
0052
0053   FitHitTime->SetRange(rangemin, rangemax);
0054   h1->Fit("FitHitTime", "QRL");
0055
0056 }
0057
0058 void calib_t0mean(const char *fname = "results/48700/calmbdpass2.3_q_48700.root")
0059 {
0060   double gmean[2];
0061   double gmeanerr[2];
0062   double gsigma[2];
0063   double gsigmaerr[2];
0064   double gchi2ndf[2];
0065   double hmean[2];
0066   double hmeanerr[2];
0067   double hrms[2];
0068   double hrmserr[2];
0069
0070   gStyle->SetStatStyle(0);
0071
0072   TFile *tfile = new TFile(fname,"READ");
0073   TH2 *h2_tt = (TH2F*)tfile->Get("h2_tt");
0074
0075   TH1 *h_ttarm[2]{nullptr};
0076   h_ttarm[0] = h2_tt->ProjectionX("h_ttarm0",1,64);
0077   h_ttarm[1] = h2_tt->ProjectionX("h_ttarm1",65,128);
0078
0079   // Get Histogram results
0080   for (int iarm=0; iarm<2; iarm++)
0081   {
0082     //h_ttarm[iarm]->ResetStats();
0083
0084     hmean[iarm] = h_ttarm[iarm]->GetMean();
0085     hmeanerr[iarm] = h_ttarm[iarm]->GetMeanError();
0086     hrms[iarm] = h_ttarm[iarm]->GetRMS();
0087     hrms[iarm] = h_ttarm[iarm]->GetRMSError();
0088   }
0089
0090   // Now do gaus fit to central peak
0091   tspec = new TSpectrum(5);  // 5 peaks is enough - we have 4
0092   FitHitTime = new TF1("FitHitTime","gaus",-5,5);
0093
0094   TCanvas *ac = new TCanvas("ac","hittime",1200,600);
0095   ac->Divide(2,1);
0096   for (int iarm=0; iarm<2; iarm++)
0097   {
0098     ac->cd(iarm+1);
0099     Double_t sigma = 5.0; // start by assuming peaks are 5 ns sigma
0100     gchi2ndf[iarm] = 1e9;
0101     while ( gchi2ndf[iarm] > 2.0 )
0102     {
0103       FitGausToPeak( h_ttarm[iarm], sigma );
0104       gchi2ndf[iarm] = FitHitTime->GetChisquare()/FitHitTime->GetNDF();
0105       if ( gchi2ndf[iarm] > 2.0 )
0106       {
0107         sigma = sigma - 0.5;
0108         //std::cout << "ERROR, run iarm sigma gmean chi2ndf " << fname << "\t" << iarm << "\t" << sigma << "\t" << gchi2ndf[iarm] << std::endl;
0109         if ( sigma <= 0.5 )
0110         {
0111           break;
0112         }
0113       }
0114     }
0115     h_ttarm[iarm]->Draw();
0116
0117     gmean[iarm] = FitHitTime->GetParameter(1);
0118     gmeanerr[iarm] = FitHitTime->GetParError(1);
0119     gsigma[iarm] = FitHitTime->GetParameter(2);
0120     gsigmaerr[iarm] = FitHitTime->GetParError(2);
0121
0122   }
0123
0124   TString pdfname = fname;
0125   pdfname.ReplaceAll("_q-","_t0mean-");
0126   pdfname.ReplaceAll(".root",".pdf");
0127   ac->Print( pdfname );
0128
0129   // Get Mean of two peaks
0130   float t0corr_mean = (gmean[0]+gmean[1])/2.0;
0131   float t0corr_meanerr = sqrt(gmeanerr[0]*gmeanerr[0]+gmeanerr[1]*gmeanerr[1])/2.0;
0132
0133   // Save results to calib files
0134   std::filesystem::path dir(fname);
0135   //std::cout << dir.parent_path() << std::endl;
0136   std::string calibfname = dir.parent_path();
0137   calibfname += "/mbd_t0corr.calib";
0138   std::cout << "Saving calibs to " << calibfname << std::endl;
0139
0140   std::ofstream calibsavefile( calibfname );
0141   calibsavefile << t0corr_mean << "\t" << t0corr_meanerr << std::endl;
0142   calibsavefile << gmean[0] << "\t" << gmeanerr[0] << "\t" << gsigma[0] << "\t" << gsigmaerr[0]
0143         << "\t" << gmean[1] << "\t" << gmeanerr[1] << "\t" << gsigma[1] << "\t" << gsigmaerr[1] << std::endl;
0144   calibsavefile << hmean[0] << "\t" << hmeanerr[0] << "\t" << hrms[0] << "\t" << hrmserr[0]
0145         << "\t" << hmean[1] << "\t" << hmeanerr[1] << "\t" << hrms[1] << "\t" << hrmserr[1] << std::endl;
0146   calibsavefile.close();
0147   make_cdbtree( calibfname.c_str() );
0148 }