mbd/SLEW/recal_mbd_slew.C

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