mbd/TIMING/make_timecorr.C

0001 //
0002 // make the timecorr calib files
0003 // input is the mbdtimecalib.C output
0004 //
0005 #include "get_runstr.h"
0006
0007 #include <mbd/MbdCalib.h>
0008 #include <mbd/MbdDefs.h>
0009 #include <mbd/MbdGeomV1.h>
0010
0011 #include <TF1.h>
0012 #include <TFile.h>
0013 #include <TGraphErrors.h>
0014 #include <TPad.h>
0015 #include <TSystem.h>
0016
0017 #include <iostream>
0018 // #include <GausProc.h>
0019 #include <fstream>
0020 #include <vector>
0021
0022 R__LOAD_LIBRARY(libmbd_io.so)
0023
0024 const int NPOINTS = 1000;  // number of points in correction LUT
0025 Int_t feech[MbdDefs::MBD_N_FEECH];
0026 Int_t npts[MbdDefs::MBD_N_FEECH];
0027 Double_t mintdc[MbdDefs::MBD_N_FEECH];  // subtracted tdc
0028 Double_t maxtdc[MbdDefs::MBD_N_FEECH];
0029 Double_t mintime[MbdDefs::MBD_N_FEECH];  // min time in range
0030 Double_t maxtime[MbdDefs::MBD_N_FEECH];  // max time in range
0031 Double_t timecorr[MbdDefs::MBD_N_FEECH][NPOINTS];
0032
0033 /*
0034 void gpr_interpolate_timecorr( const int ifeech, TGraphErrors *g )
0035 {
0036   Int_t n = g->GetN();
0037   Double_t *x = g->GetX();
0038   Double_t *y = g->GetY();
0039   Double_t *ey = g->GetEY();
0040   g->Print("ALL");
0041   g->SetMarkerStyle(20);
0042   g->SetMarkerSize(0);
0043   g->Draw("ap");
0044
0045   // parameters for calibrations output
0046   feech[ifeech] = ifeech;
0047   npts[ifeech] = NPOINTS;
0048   mintdc[ifeech] = 0.;
0049   maxtdc[ifeech] = 16000.;
0050
0051   // set up GPR to do interpolation
0052   // convert tp data to std::vectors
0053   std::vector<Double_t> vx( x, x + n );
0054   std::vector<Double_t> vy( y, y + n );
0055   std::vector<Double_t> vey( ey, ey + n );
0056   //std::fill(vey.begin(), vey.end(), 0.003);     // temporary
0057
0058   for (int i=100; i<110; i++)
0059   {
0060     std::cout << "vvv " << i << "\t" << n << "\t" << vx[i] << "\t" << vy[i] << "\t" << vey[i] << std::endl;
0061   }
0062
0063   unsigned int nPredictions = NPOINTS;
0064   Double_t xmin = mintdc[ifeech];
0065   Double_t xmax = maxtdc[ifeech];
0066   const char *outfile = "null.root";
0067   gSystem->Exec(Form("rm -f %s",outfile));
0068
0069   GausProc gp(vx, vy, vey, xmin, xmax, nPredictions, outfile);
0070   int verbosity = 1;
0071   gp.SetVerbosity(verbosity);
0072   gp.SetKernel(GausProc::RBF);
0073
0074   const int npar = 2;
0075   const int nparfixed = 0;
0076   // 2,10 = par[0],par[1], which are the start points for the minimizer
0077   GPOptimizer gpopt(&gp,2,10);
0078   gpopt.GPoptimize(npar,nparfixed);
0079
0080   std::cout << "gpopt pars " << gpopt.getPar(0) << " " << gpopt.getPar(1) << std::endl;
0081
0082   GausProc gp_result(vx, vy, vey, xmin, xmax, nPredictions, outfile);
0083   gp_result.SetPar(0,gpopt.getPar(0));
0084   gp_result.SetPar(1,gpopt.getPar(1));
0085   gp_result.process();
0086   gp_result.Write(-1);
0087
0088   // Get the resulting histogram
0089   TFile fin(outfile,"READ");
0090   TH1F *ho = (TH1F*)fin.Get("ho");
0091
0092   int nbinsx = ho->GetNbinsX();
0093   feech[ifeech] = ifeech;
0094   npts[ifeech] = NPOINTS;
0095   mintdc[ifeech] = 0.;
0096   maxtdc[ifeech] = 16000.;
0097   std::cout << feech[ifeech] << "\t" << npts[ifeech] << "\t" << mintdc[ifeech] << "\t" << maxtdc[ifeech] << std::endl;
0098   std::cout << nbinsx << std::endl;
0099   for (int ibin=1; ibin<=5; ibin++)
0100   {
0101     std::cout << ibin << " ";
0102     std::cout << ho->GetBinCenter(ibin) << " ";
0103     std::cout << ho->GetBinContent(ibin) << " ";
0104     //if ( ibin%10 == 0 ) std::cout << std::endl;
0105     std::cout << std::endl;
0106   }
0107   for (int ibin=nbinsx-5; ibin<=nbinsx; ibin++)
0108   {
0109     std::cout << ibin << " ";
0110     std::cout << ho->GetBinCenter(ibin) << " ";
0111     std::cout << ho->GetBinContent(ibin) << " ";
0112     //if ( ibin%10 == 0 ) std::cout << std::endl;
0113     std::cout << std::endl;
0114   }
0115
0116   fin.Close();
0117
0118 }
0119 */
0120
0121 // use poly fit to interpolate
0122 // now using spline fit
0123 void interpolate_timecorr(const int ifeech, TGraphErrors *g, TF1 *f)
0124 {
0125   Int_t n = g->GetN();
0126   Double_t *x = g->GetX();
0127   Double_t *y = g->GetY();
0128
0129   g->Fit(f, "R");
0130
0131   g->Draw("ap");
0132   gPad->Modified();
0133   gPad->Update();
0134
0135   float mintimelocal = TMath::MinElement(n, y);
0136   float maxtimelocal = TMath::MaxElement(n, y);
0137   float mintdclocal = TMath::MinElement(n, x);
0138   float maxtdclocal = TMath::MaxElement(n, x);
0139   std::cout << ifeech << "\t" << mintimelocal << "\t" << maxtimelocal << "\t" << maxtimelocal - mintimelocal << std::endl;
0140   std::cout << "tdc\t" << mintdclocal << "\t" << maxtdclocal << std::endl;
0141
0142   // parameters for calibrations output
0143   feech[ifeech] = ifeech;
0144   npts[ifeech] = NPOINTS;
0145   mintdc[ifeech] = mintdclocal;
0146   maxtdc[ifeech] = maxtdclocal;
0147   mintime[ifeech] = mintimelocal;
0148   maxtime[ifeech] = maxtimelocal;
0149
0150   // std::cout << feech[ifeech] << "\t" << npts[ifeech] << "\t" << mintdc[ifeech] << "\t" << maxtdc[ifeech] << std::endl;
0151
0152   // double dstep = (maxtdc[ifeech] - mintdc[ifeech])/(NPOINTS-1);
0153   // double tdc = mintdc[ifeech];
0154   double dstep = 14985. / (NPOINTS - 1);
0155   double tdc = 0;
0156   for (int istep = 0; istep < NPOINTS; istep++)
0157   {
0158     // timecorr[ifeech][istep] = f->Eval(tdc) - mintimelocal;
0159     timecorr[ifeech][istep] = g->Eval(tdc) - mintimelocal;
0160
0161     if (istep == 0 || istep == NPOINTS - 1)
0162     {
0163       std::cout << "ch " << ifeech << "\t" << istep << "\t" << tdc << "\t" << timecorr[ifeech][istep] << std::endl;
0164     }
0165
0166     tdc += dstep;
0167   }
0168 }
0169
0170 void make_timecorr(const char *rootfname = "calib_seb18-00029705-0000_mbdtimecalib.root")
0171 {
0172   MbdGeom *mbdgeom = new MbdGeomV1();
0173   // MbdCalib *mcal = new MbdCalib();
0174
0175   //== Create output directory (should already exist)
0176   TString dir = "results/";
0177   dir += get_runstr(rootfname);
0178   dir.ReplaceAll("_mbdtimecalib", "");
0179   // dir.ReplaceAll("_mbdtrigtimecalib","");
0180   dir += "/";
0181   // TString name = "mkdir -p " + dir;
0182   TString name = "echo " + dir;
0183   gSystem->Exec(name);
0184   std::cout << name << std::endl;
0185
0186   // TF1 *fpoly = new TF1("fpoly","pol3",0,15000);
0187   // fpoly->SetParameters(28.,-0.0018,0.,0.);
0188   TF1 *fpoly = new TF1("fpoly", "pol4", 0, 15000);
0189   fpoly->SetParameters(28., -0.0018, -1e-7, 1.3e-11, -5e-16);
0190   fpoly->SetLineColor(4);
0191   TGraphErrors *g_delaytime[MbdDefs::MBD_N_FEECH];  // delay curves, corrected for time
0192   TFile *tfile = new TFile(rootfname, "READ");
0193   for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++)
0194   {
0195     if (mbdgeom->get_type(ifeech) == 1)
0196     {
0197       continue;  // skip q-channels
0198     }
0199
0200     name = "g_delaytime";
0201     name += ifeech;
0202     g_delaytime[ifeech] = (TGraphErrors *) tfile->Get(name);
0203     std::cout << name << std::endl;
0204
0205     // do a prefit since for some reason the first fit fails
0206     if (ifeech == 0)
0207     {
0208       g_delaytime[ifeech]->Fit(fpoly, "R");
0209     }
0210
0211     // if ( ifeech!=118)
0212     //{
0213     interpolate_timecorr(ifeech, g_delaytime[ifeech], fpoly);
0214     //}
0215
0216     // if ( ifeech>4 ) break;
0217   }
0218
0219   // write out mbd_tpscan.txt file
0220   TString tscan_fname = dir;
0221   tscan_fname += "/mbd_tpscan.txt";
0222   std::cout << tscan_fname << std::endl;
0223   std::ofstream tscan_file(tscan_fname);
0224   for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++)
0225   {
0226     if (mbdgeom->get_type(ifeech) == 1)
0227     {
0228       continue;  // skip q-channels
0229     }
0230
0231     tscan_file << feech[ifeech] << "\t" << mintdc[ifeech] << "\t" << maxtdc[ifeech] << "\t"
0232                << maxtdc[ifeech] - mintdc[ifeech] << "\t"
0233                << mintime[ifeech] << "\t" << maxtime[ifeech] << "\t" << maxtime[ifeech] - mintime[ifeech] << std::endl;
0234   }
0235   tscan_file.close();
0236
0237   // write out mbd_timecorr.calib file
0238   TString tcorr_fname = dir;
0239   tcorr_fname += "/mbd_timecorr.calib";
0240   // TString tcorr_fname = dir; tcorr_fname += "/mbd_trigtimecorr.calib";
0241   std::cout << "tcorr_fname " << tcorr_fname << std::endl;
0242   std::ofstream tcorr_file(tcorr_fname);
0243   for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++)
0244   {
0245     if (mbdgeom->get_type(ifeech) == 1)
0246     {
0247       continue;  // skip q-channels
0248     }
0249
0250     // tcorr_file << feech[ifeech] << "\t" << npts[ifeech] << "\t" << mintdc[ifeech] << "\t" << maxtdc[ifeech] << std::endl;
0251     tcorr_file << feech[ifeech] << "\t" << npts[ifeech] << "\t" << 0 << "\t" << 14985 << std::endl;
0252     for (int ipt = 0; ipt < NPOINTS; ipt++)
0253     {
0254       tcorr_file << timecorr[ifeech][ipt] << " ";
0255       if (ipt % 10 == 9)
0256       {
0257         tcorr_file << std::endl;
0258       }
0259     }
0260   }
0261   tcorr_file.close();
0262 }