mbd/TIMING/make_timecorr.C

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