calibrations/mbd/calib_t0mean.C

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