mbd/SLEW/recal_mbd_slew.C

0001 //
0002 // Do a recalibration of the slew from the saved histograms
0003 //
0004 #include "get_runstr.h"
0005
0006 #include <mbd/MbdCalib.h>
0007 #include <mbd/MbdDefs.h>
0008 #include <mbd/MbdGeomV1.h>
0009
0010 #include <TCanvas.h>
0011 #include <TF1.h>
0012 #include <TFile.h>
0013 #include <TGraphErrors.h>
0014 #include <TH1F.h>
0015 #include <TH2F.h>
0016 #include <TString.h>
0017 #include <TSystem.h>
0018
0019 #include <fstream>
0020 #include <iostream>
0021
0022 R__LOAD_LIBRARY(libmbd_io.so)
0023
0024 const int NPOINTS = 16000;  // number of points in correction LUT
0025 const int MINADC = 0;       // subtracted adc
0026 const int MAXADC = 15999;
0027
0028 int verbose = 0;
0029
0030 TGraphErrors *g_slew[MbdDefs::MBD_N_PMT];
0031 TF1 *f_slewfit[MbdDefs::MBD_N_PMT];
0032
0033 // find the ridge of the TH2
0034 TGraphErrors *find_th2ridge(const TH2 *h2)
0035 {
0036   int nbinsx = h2->GetNbinsX();
0037   int nbinsy = h2->GetNbinsY();
0038   // double min_xrange = h2->GetXaxis()->GetBinLowEdge(1);
0039   // double max_xrange = h2->GetXaxis()->GetBinLowEdge(nbinsx+1);
0040   double min_yrange = h2->GetYaxis()->GetBinLowEdge(1);
0041   double max_yrange = h2->GetYaxis()->GetBinLowEdge(nbinsy + 1);
0042
0043   TString name;
0044   TString title;
0045   name = "aaa";
0046   title = "aaa";
0047   // TH1D prof(name,title,nbinsx,min_xrange,max_xrange);
0048   TGraphErrors *prof = new TGraphErrors();
0049   prof->SetName(name);
0050   prof->SetTitle(title);
0051
0052   TH1 *h_projx = h2->ProjectionX("projx");
0053   // std::unique_ptr<TH2D> h_projx( h2->ProjectionX("projx") );
0054
0055   TF1 gaussian("gaussian", "gaus", min_yrange, max_yrange);
0056   gaussian.SetLineColor(4);
0057
0058   TH1 *h_projy{nullptr};
0059   double adcmean = 0.;
0060   double adcnum = 0.;
0061
0062   for (int ibin = 1; ibin <= nbinsx; ibin++)
0063   {
0064     name = "hproj_";
0065     name += ibin;
0066     if (h_projy == nullptr)
0067     {
0068       h_projy = h2->ProjectionY(name, ibin, ibin);
0069       adcmean = h_projx->GetBinCenter(ibin);
0070       adcnum = 1.0;
0071     }
0072     else
0073     {
0074       TH1 *h_projyadd = h2->ProjectionY(name, ibin, ibin);
0075       h_projy->Add(h_projyadd);
0076       delete h_projyadd;
0077
0078       adcmean += h_projx->GetBinCenter(ibin);
0079       adcnum += 1.0;
0080     }
0081
0082     if (h_projy->Integral() > 2000 || ibin == nbinsx)
0083     {
0084       adcmean = adcmean / adcnum;
0085
0086       h_projy->Draw();
0087
0088       int maxbin = h_projy->GetMaximumBin();
0089       double xmax = h_projy->GetBinCenter(maxbin);
0090       double ymax = h_projy->GetBinContent(maxbin);
0091       gaussian.SetParameter(1, xmax);
0092       gaussian.SetParameter(0, ymax);
0093       gaussian.SetRange(xmax - 0.6, xmax + 0.6);
0094
0095       h_projy->Fit("gaussian", "RWW");
0096
0097       double mean = gaussian.GetParameter(1);
0098       double meanerr = gaussian.GetParError(1);
0099       if (meanerr < 1.0)
0100       {
0101         int n = prof->GetN();
0102         prof->SetPoint(n, adcmean, mean);
0103         prof->SetPointError(n, 0, meanerr);
0104       }
0105
0106       gPad->Modified();
0107       gPad->Update();
0108       /*
0109          std::string junk;
0110          std::cin >> junk;
0111          */
0112
0113       delete h_projy;
0114       h_projy = nullptr;
0115     }
0116   }
0117
0118   // interpolate last point out to ADC = 16000
0119   int n = prof->GetN();
0120   double x1;
0121   double x2;
0122   double y1;
0123   double y2;
0124   prof->GetPoint(n - 2, x1, y1);
0125   prof->GetPoint(n - 1, x2, y2);
0126
0127   delete h_projx;
0128
0129   prof->SetBit(TGraph::kIsSortedX);
0130   return prof;
0131 }
0132
0133 //
0134 // pass: should be the same as cal_mbd pass number
0135 //
0136 // need to fix run number issue...
0137 void recal_mbd_slew(const char *tfname = "calmbdslew_pass1-54321.root", const int pass = 1, const int /*nevt*/ = 0)
0138 {
0139   std::cout << "tfname " << tfname << std::endl;
0140   MbdGeom *mbdgeom = new MbdGeomV1();
0141
0142   // Read in TFile with h_q
0143   TFile *oldfile = new TFile(tfname, "READ");
0144
0145   TH2 *h2_slew[MbdDefs::MBD_N_FEECH] = {};
0146
0147   TString name;
0148   TString title;
0149   for (int ipmt = 0; ipmt < MbdDefs::MBD_N_PMT; ipmt++)
0150   {
0151     int feech = (ipmt / 8) * 16 + ipmt % 8;
0152     name = "h2_slew";
0153     name += ipmt;
0154     std::cout << name << "\t" << feech << std::endl;
0155     h2_slew[feech] = (TH2 *) oldfile->Get(name);
0156     if (h2_slew[feech] == nullptr)
0157     {
0158       std::cout << "ERROR, " << name << " not found in " << tfname << std::endl;
0159       return;
0160     }
0161   }
0162
0163   // Create new TFile
0164   TString dir = "results/";
0165   dir += get_runnumber(tfname);
0166   dir += "/";
0167   name = "mkdir -p ";
0168   name += dir;
0169   gSystem->Exec(name);
0170   name = dir;
0171   name += "recalmbdslew_pass2.";
0172   name += pass;
0173   name += ".root";
0174   std::cout << name << std::endl;
0175
0176   TFile *savefile = new TFile(name, "RECREATE");
0177
0178   TString pdfname = name;
0179   pdfname.ReplaceAll(".root", ".pdf");
0180
0181   // Load in calib constants
0182
0183   TCanvas *ac[100];
0184   int cvindex = 0;
0185
0186   // slew curves
0187   ac[cvindex] = new TCanvas("cal_slew", "slew", 425 * 1.5, 550 * 1.5);
0188   ac[cvindex]->Print(pdfname + "[");
0189
0190   for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++)
0191   {
0192     if (mbdgeom->get_type(ifeech) == 1)
0193     {
0194       continue;  // skip q-channels
0195     }
0196
0197     int pmtch = mbdgeom->get_pmt(ifeech);
0198
0199     /*
0200     // Use LUT instead of fit
0201     name = "f_slewfit"; name += pmtch;
0202     f_slewfit[pmtch] = new TF1(name,"[0]+([1]/x)+[2]*log(x)",20,16000);
0203     f_slewfit[pmtch]->SetParameters(4.0,1.,1.);
0204     f_slewfit[pmtch]->SetLineColor(2);
0205     */
0206
0207     // h2_slew[ifeech]->RebinX(20);
0208     // h2_slew[ifeech]->RebinY(10);
0209     name = "g_slew";
0210     name += pmtch;
0211     std::cout << name << std::endl;
0212     g_slew[pmtch] = find_th2ridge(h2_slew[ifeech]);
0213     g_slew[pmtch]->SetName(name);
0214     g_slew[pmtch]->SetMarkerStyle(20);
0215     g_slew[pmtch]->SetMarkerSize(0.25);
0216
0217     ac[cvindex]->cd();
0218     h2_slew[ifeech]->Draw("colz");
0219     g_slew[pmtch]->Draw("cp");
0220
0221     /*
0222        g_slew[pmtch]->Fit( f_slewfit[pmtch], "R" );
0223        f_slewfit[pmtch]->DrawCopy("same");
0224        double par1 = f_slewfit[pmtch]->GetParameter(1);
0225        double par2 = f_slewfit[pmtch]->GetParameter(2);
0226        f_slewfit[pmtch]->SetParameter(1,0.);
0227        f_slewfit[pmtch]->SetLineColor(6);
0228        f_slewfit[pmtch]->DrawCopy("same");
0229        f_slewfit[pmtch]->SetParameter(1,par1);
0230        f_slewfit[pmtch]->SetParameter(2,0.);
0231        f_slewfit[pmtch]->SetLineColor(4);
0232        f_slewfit[pmtch]->DrawCopy("same");
0233     */
0234
0235     gPad->Modified();
0236     gPad->Update();
0237
0238     if (verbose > 10)
0239     {
0240       std::string junk;
0241       std::cout << "? ";
0242       std::cin >> junk;
0243     }
0244
0245     // name = dir + "/h2_slewfit"; name += pmtch;
0246     // name += "_pass"; name += pass; name += ".png";
0247     name = "h2_slewfit";
0248     name += pmtch;
0249     name += "_pass";
0250     name += pass;
0251     std::cout << name << std::endl;
0252     ac[cvindex]->Print(pdfname, name);
0253   }
0254   ac[cvindex]->Print(pdfname + "]");
0255
0256   ++cvindex;
0257
0258   // Write out slew curves to temp calib file
0259   TString scorr_fname = dir;
0260   scorr_fname += "/pass";
0261   scorr_fname += pass;
0262   scorr_fname += "_mbd_slewcorr.calib";
0263   std::cout << scorr_fname << std::endl;
0264   std::ofstream scorr_file(scorr_fname);
0265   for (int ifeech = 0; ifeech < MbdDefs::MBD_N_FEECH; ifeech++)
0266   {
0267     if (mbdgeom->get_type(ifeech) == 1)
0268     {
0269       continue;  // skip q-channels
0270     }
0271     int pmtch = mbdgeom->get_pmt(ifeech);
0272
0273     scorr_file << ifeech << "\t" << NPOINTS << "\t" << MINADC << "\t" << MAXADC << std::endl;
0274     int step = (MAXADC - MINADC) / (NPOINTS - 1);
0275     // std::cout << "STEP " << step << std::endl;
0276     for (int iadc = MINADC; iadc <= MAXADC; iadc += step)
0277     {
0278       float slewcorr = g_slew[pmtch]->Eval(iadc);
0279       scorr_file << slewcorr << " ";
0280       if (iadc % 10 == 9)
0281       {
0282         scorr_file << std::endl;
0283       }
0284     }
0285   }
0286   scorr_file.close();
0287
0288   if (pass > 0)
0289   {
0290     // write out the slew curves
0291     for (auto &ipmt : g_slew)
0292     {
0293       ipmt->Write();
0294     }
0295   }
0296   savefile->Write();
0297   // savefile->Close();
0298 }