calibrations/mbd/plot_calibgains.C

0001 // Make plots of the values of the gains
0002 #include "mbd/MbdCalib.h"
0003 //#include "get_runstr.h"
0004
0005 R__LOAD_LIBRARY(libmbd.so)
0006 R__LOAD_LIBRARY(libmbd_io.so)
0007
0008 void read_calibgains(const char *flist);
0009 void plot_relwidth_onerun(const int irun = 0);
0010
0011 TFile *savefile;
0012
0013 const int MAX_RUNS = 2500;
0014 int nruns = 0;
0015
0016 MbdCalib *bcal[MAX_RUNS];
0017
0018 Double_t bqmean[128][MAX_RUNS];
0019 Double_t bqmeanerr[128][MAX_RUNS];
0020 Double_t bqwidth[128][MAX_RUNS];
0021 Double_t bqwidtherr[128][MAX_RUNS];
0022 Double_t bqrelwidth[128][MAX_RUNS];
0023 Double_t bqrelwidtherr[128][MAX_RUNS];
0024 Double_t listofruns[MAX_RUNS];
0025 Double_t runindex[MAX_RUNS];
0026
0027 TGraphErrors *gainvals[128];
0028 TGraphErrors *g_relwidth[128];  // relative width
0029 TH1 *h_relwidth{nullptr};
0030
0031 const int update_qfit = 1;    // whether to write update gain files
0032 const int write_plots = 1;    // whether to save plots to png
0033
0034
0035 // flist is a list of run-seq or run number directories
0036 // where the calibrations are stored
0037 void read_calibgains(const char *flist)
0038 {
0039   cout << "read_calibgains()" << endl;
0040   ifstream inflist(flist);
0041
0042   ifstream cal_mip_file;
0043
0044   TString calrunseq;
0045   TString calfile;
0046   TString name;
0047   while ( inflist >> calrunseq )
0048   {
0049     calfile = "results/" + calrunseq + "/mbd_qfit.calib";
0050     cout << calfile << endl;
0051     cal_mip_file.open( calfile );
0052
0053     TString runtext = calrunseq;
0054     runtext.ReplaceAll("-0000","");
0055     listofruns[nruns] = runtext.Atof();
0056     cout << listofruns[nruns] << endl;
0057
0058     runindex[nruns] = nruns;
0059
0060     //bcal[nruns] = new MbdCalib();
0061     //bcal[nruns]->Download_Gains( calfile.Data() );
0062
0063     int    temp_pmt;
0064     double integ;
0065     double best_peak;
0066     double width;
0067     double integerr;
0068     double best_peakerr;
0069     double widtherr;
0070     double chi2ndf;
0071
0072     double corrected_peak;
0073
0074     for (int ipmt=0; ipmt<128; ipmt++)
0075     {
0076       //float gain = bcal[nruns]->get_qgain(ipmt);
0077       //float gainerr = bcal[nruns]->get_qgainerr(ipmt);
0078
0079       cal_mip_file >> temp_pmt >> integ >> best_peak >> width
0080         >> integerr >> best_peakerr >> widtherr >> chi2ndf;
0081
0082       if ( ipmt==0 ) cout << temp_pmt << "\t" << integ << "\t" << best_peak << "\t" << width
0083         << "\t" << integerr << "\t" << best_peakerr << "\t" << widtherr << "\t" << chi2ndf << endl;
0084
0085       bqmean[temp_pmt][nruns] = best_peak;
0086       bqmeanerr[temp_pmt][nruns] = best_peakerr;
0087       bqwidth[temp_pmt][nruns] = fabs(width);
0088       bqwidtherr[temp_pmt][nruns] = widtherr;
0089
0090       /*
0091       if ( listofruns[nruns] == 21520 )
0092       {
0093         cout << "XXX " << calfile << "\t" << ipmt << "\t" << best_peakerr << "\t"
0094           << best_peak << "\t" << bqmean[temp_pmt][nruns-1] << endl;
0095       }
0096       */
0097
0098       // check for bad fit
0099       //if ( integ < 100. || bqmean[temp_pmt][nruns]<0. || chi2ndf>4.0 )
0100       if ( integ < 100. || bqmean[temp_pmt][nruns]<0. || chi2ndf>10.0 )
0101       {
0102         cout << "BAD " << calfile << "\t" << ipmt << "\t" << best_peak
0103           << "\t" << integ << "\t" << bqmean[temp_pmt][nruns] << "\t" << chi2ndf << endl;
0104         bqmean[temp_pmt][nruns] = NAN;
0105       }
0106
0107       if ( bqmeanerr[temp_pmt][nruns]>10. )
0108       {
0109         cout << "BADERR " << calfile << "\t" << ipmt << "\t" << best_peakerr << "\t"
0110           << best_peak << "\t" << bqmean[temp_pmt][nruns-1] << endl;
0111         bqmeanerr[temp_pmt][nruns] = NAN;
0112       }
0113     }
0114
0115     cal_mip_file.close();
0116
0117     nruns++;
0118   }
0119   inflist.close();
0120 }
0121
0122 void plot_relwidth_onerun(const int irun = 0)
0123 {
0124   if ( h_relwidth==nullptr )
0125   {
0126     h_relwidth = new TH1F("h_relwidth","relative width",100,0.1,0.36);
0127   }
0128   cout << "== Relative Widths ==" << endl;
0129   for (int ipmt=0; ipmt<128; ipmt++)
0130   {
0131     double relwidth = bqwidth[ipmt][irun]/bqmean[ipmt][irun];
0132     double frac_merr = bqmeanerr[ipmt][irun]/bqmean[ipmt][irun];
0133     double frac_werr = bqwidtherr[ipmt][irun]/bqwidth[ipmt][irun];
0134     double relwidtherr = frac_merr*frac_merr + frac_werr*frac_werr;
0135     relwidtherr = relwidth*sqrt(relwidtherr);
0136
0137     cout << ipmt << "\t" << relwidth << "\t" << relwidtherr << endl;
0138
0139     h_relwidth->Fill( relwidth );
0140   }
0141   h_relwidth->Draw();
0142 }
0143
0144
0145 // Find any values that deviate greatly from prior run, and set to nan
0146 // Note that 1st run must be good
0147 //void CheckForLargeDeviation( const double max_deviation = 0.05 ) // 5% deviation
0148 void CheckForLargeDeviation( const double max_deviation = 0.10 ) // 10% deviation
0149 {
0150   for (int ipmt=0; ipmt<128; ipmt++)
0151   {
0152     for (int irun=1; irun<nruns; irun++)
0153     {
0154       if ( isnan(bqmean[ipmt][irun]) ) continue;  // skip ones already bad
0155
0156       // search for prev good val
0157       double prevgoodmean = NAN;
0158       for (int prevrun=irun-1; prevrun>=0; prevrun--)
0159       {
0160         if ( !isnan( bqmean[ipmt][prevrun] ) )
0161         {
0162           prevgoodmean = bqmean[ipmt][prevrun];
0163           break;
0164         }
0165       }
0166
0167       double deviation = fabs((bqmean[ipmt][irun] - prevgoodmean)/prevgoodmean);
0168       if ( deviation > max_deviation )
0169       {
0170         cout << "BADDEV " << listofruns[irun] << "\t" << ipmt << "\t" << bqmean[ipmt][irun]
0171           << "\t" << prevgoodmean << "\t" << deviation << endl;
0172         bqmean[ipmt][irun] = NAN;
0173       }
0174     }
0175   }
0176
0177 }
0178
0179 void plot_calibgains(const char *flist = "runs.list")
0180 {
0181   // Read in all the calibrations from flist
0182   read_calibgains(flist);
0183
0184   if ( write_plots )
0185   {
0186     savefile = new TFile("results/calibgains.root","RECREATE");
0187   }
0188
0189   // sanity check that gains aren't bad
0190   CheckForLargeDeviation();
0191
0192   TString name;
0193   TString title;
0194   for (int ipmt=0; ipmt<128; ipmt++)
0195   {
0196     // do corrections
0197     for (int irun=0; irun<nruns; irun++)
0198     {
0199
0200       if ( isnan(bqmean[ipmt][irun]) )
0201       {
0202         cout << "FOUND BAD " << listofruns[irun] << " " << ipmt << " " << bqmean[ipmt][irun] << endl;
0203         // search for next good val
0204         double nextgood = NAN;
0205         for (int nextrun=irun+1; nextrun<nruns; nextrun++)
0206         {
0207           if ( !isnan( bqmean[ipmt][nextrun] ) )
0208           {
0209             nextgood = bqmean[ipmt][nextrun];
0210             break;
0211           }
0212         }
0213
0214
0215         // search for prev good val
0216         double prevgood = NAN;
0217         for (int prevrun=irun-1; prevrun>=0; prevrun--)
0218         {
0219           if ( !isnan( bqmean[ipmt][prevrun] ) )
0220           {
0221             prevgood = bqmean[ipmt][prevrun];
0222             break;
0223           }
0224         }
0225
0226         if ( !isnan(nextgood) && !isnan(prevgood) )
0227         {
0228           bqmean[ipmt][irun] = (nextgood+prevgood)/2.0;
0229         }
0230         else if ( !isnan(nextgood) )
0231         {
0232           bqmean[ipmt][irun] = nextgood;
0233         }
0234         else if ( !isnan(prevgood) )
0235         {
0236           bqmean[ipmt][irun] = prevgood;
0237         }
0238         else
0239         {
0240           cout << "ERROR, no good run to interpolate from" << endl;
0241         }
0242
0243         cout << "NEW VALUE FOR BAD " << irun << " " << ipmt << " " << bqmean[ipmt][irun] << endl;
0244       }
0245
0246       if ( isnan(bqmeanerr[ipmt][irun]) )
0247       {
0248         // search for next good val
0249         double nextgood = NAN;
0250         for (int nextrun=irun+1; nextrun<nruns; nextrun++)
0251         {
0252           if ( !isnan( bqmeanerr[ipmt][nextrun] ) )
0253           {
0254             nextgood = bqmeanerr[ipmt][nextrun];
0255             break;
0256           }
0257         }
0258
0259
0260         // search for prev good val
0261         double prevgood = NAN;
0262         for (int prevrun=irun-1; prevrun>=0; prevrun--)
0263         {
0264           if ( !isnan( bqmeanerr[ipmt][prevrun] ) )
0265           {
0266             prevgood = bqmeanerr[ipmt][prevrun];
0267             break;
0268           }
0269         }
0270
0271         if ( irun==0 )
0272         {
0273           bqmeanerr[ipmt][irun] = nextgood;
0274         }
0275         else if ( irun==nruns-1 )
0276         {
0277           bqmeanerr[ipmt][irun] = prevgood;
0278         }
0279         else
0280         {
0281           bqmeanerr[ipmt][irun] = (nextgood+prevgood)/2.0;
0282         }
0283       }
0284     }
0285
0286     name = "gainvals"; name += ipmt;
0287     title = "gain, ch"; title += ipmt;
0288     //gainvals[ipmt] = new TGraphErrors(nruns,runindex,bqmean[ipmt],0,bqmeanerr[ipmt]);
0289     gainvals[ipmt] = new TGraphErrors(nruns,listofruns,bqmean[ipmt],0,bqmeanerr[ipmt]);
0290     gainvals[ipmt]->SetName( name );
0291     gainvals[ipmt]->SetTitle( title );
0292
0293     gainvals[ipmt]->Draw("ap");
0294     gPad->Modified();
0295     gPad->Update();
0296
0297     if ( write_plots )
0298     {
0299       name += ".png";
0300       gPad->Print( name );
0301       gainvals[ipmt]->Write();
0302     }
0303   }
0304
0305   // Get mean gain for each channel, and write out the values
0306   TString meangainfname = "results/"; meangainfname += flist;
0307   meangainfname.ReplaceAll(".list","_meangains.calib");
0308   ofstream meangainfile( meangainfname );
0309   TF1 *f_meangain[128] {nullptr};
0310   Double_t grp_mean[128];
0311   Double_t grp_meanerr[128];
0312   for (int ipmt=0; ipmt<128; ipmt++)
0313   {
0314     name = "f_meangain"; name += ipmt;
0315     f_meangain[ipmt] = new TF1(name,"pol0",0,1e9);
0316     gainvals[ipmt]->Fit(f_meangain[ipmt]);
0317     gPad->Update();
0318     gPad->Modified();
0319     grp_mean[ipmt] = f_meangain[ipmt]->GetParameter(0);
0320     grp_meanerr[ipmt] = f_meangain[ipmt]->GetParError(0);
0321
0322     meangainfile << ipmt << "\t" << grp_mean[ipmt] << "\t" << grp_meanerr[ipmt] << endl;
0323     cout << ipmt << "\t" << grp_mean[ipmt] << "\t" << grp_meanerr[ipmt] << endl;
0324   }
0325   meangainfile.close();
0326
0327   // generate new gains
0328   ifstream inflist;
0329   ifstream cal_mip_file;
0330   TString calrunseq;
0331   TString calfile;
0332
0333   if ( update_qfit )
0334   {
0335     inflist.open(flist);
0336     ofstream newcal_mip_file;
0337
0338     int    temp_pmt;
0339     double integ;
0340     double best_peak;
0341     double width;
0342     double integerr;
0343     double best_peakerr;
0344     double widtherr;
0345     double chi2ndf;
0346
0347     int irun = 0;
0348     while ( inflist >> calrunseq )
0349     {
0350       calfile = "results/" + calrunseq + "/mbd_qfit.calib";
0351       cout << calfile << endl;
0352       cal_mip_file.open( calfile );
0353
0354       // create new gain corr file
0355       calfile = "results/" + calrunseq + "/newgainvals_mbd_qfit.calib";
0356       newcal_mip_file.open( calfile );
0357
0358       for (int ipmt=0; ipmt<128; ipmt++)
0359       {
0360
0361         cal_mip_file >> temp_pmt >> integ >> best_peak >> width
0362           >> integerr >> best_peakerr >> widtherr >> chi2ndf;
0363
0364         newcal_mip_file << temp_pmt << "\t" << integ << "\t" << bqmean[ipmt][irun] << "\t" << width << "\t"
0365           << integerr << "\t" << bqmeanerr[ipmt][irun] << "\t" << widtherr << "\t"
0366           << chi2ndf << endl;
0367       }
0368
0369       cal_mip_file.close();
0370       newcal_mip_file.close();
0371       irun++;
0372     }
0373   }
0374
0375   if ( write_plots )
0376   {
0377     savefile->Write();
0378   }
0379 }
0380