mbd/HV/calcnew_gainsbyhvgroup.C

0001 // Make plots of the values of the gains
0002 #include <mbd/MbdCalib.h>
0003 #include <mbd/MbdGeomV2.h>
0004
0005 #include <TCanvas.h>
0006 #include <TFile.h>
0007 #include <TGraph.h>
0008 #include <TGraphErrors.h>
0009 #include <TH1.h>
0010 #include <TPad.h>
0011
0012 #include <algorithm>
0013
0014 R__LOAD_LIBRARY(libmbd.so)
0015 R__LOAD_LIBRARY(libmbd_io.so)
0016
0017 int nruns = 0;
0018 MbdCalib *mcal{nullptr};
0019 MbdGeomV2 *mbdgeom{nullptr};
0020
0021 void read_calibgains(const char *fname)
0022 {
0023   mcal = new MbdCalib();
0024   mcal->Download_Gains(fname);
0025 }
0026
0027 TGraphErrors *g_gains{nullptr};
0028 TGraphErrors *g_gainsbyhv[16]{nullptr};
0029
0030 double gmean(TGraph *g)
0031 {
0032   int n = g->GetN();
0033   Double_t *y = g->GetY();
0034   double sum = 0.;
0035   for (int i = 0; i < n; i++)
0036   {
0037     sum += y[i];
0038   }
0039   double avg = sum / n;
0040
0041   return avg;
0042 }
0043
0044 // void calcnew_gainsbyhvgroup(const char *fname = "results/42045/mbd_qfit.calib")
0045 // void calcnew_gainsbyhvgroup(const char *fname = "results/53635/mbd_qfit.calib")
0046 void calcnew_gainsbyhvgroup(const char *fname = "results/54935/mbd_qfit.calib")
0047 {
0048   mbdgeom = new MbdGeomV2();
0049
0050   read_calibgains(fname);
0051
0052   TString name;
0053   // g_gains = new TGraph();
0054   for (int ihv = 0; ihv < 16; ihv++)
0055   {
0056     name = "g_gainsbyhv";
0057     name += ihv;
0058     g_gainsbyhv[ihv] = new TGraphErrors();
0059     g_gainsbyhv[ihv]->SetName(name);
0060     g_gainsbyhv[ihv]->SetLineColor(ihv % 8 + 1);
0061     g_gainsbyhv[ihv]->SetMarkerColor(ihv % 8 + 1);
0062     /*
0063     if ( ihv+1 == 10 )
0064     {
0065       g_gainsbyhv[ihv]->SetLineColor( 50 );
0066       g_gainsbyhv[ihv]->SetMarkerColor( 50 );
0067     }
0068     */
0069     g_gainsbyhv[ihv]->SetMarkerStyle(20);
0070     g_gainsbyhv[ihv]->GetHistogram()->SetXTitle("hv groups");
0071     g_gainsbyhv[ihv]->GetHistogram()->SetYTitle("gain [ADC/mip]");
0072   }
0073
0074   int npmts = 0;
0075   std::array<double, 16> gmin{};
0076   std::array<double, 16> gmax{};
0077   std::array<int, 16> minpmt{};
0078   std::array<int, 16> maxpmt{};
0079   gmin.fill(1e9);  // the gain settings that were used for run analyzed
0080   gmax.fill(0);
0081   minpmt.fill(-1);
0082   maxpmt.fill(-1);
0083   std::multimap hvmap = mbdgeom->get_hvmap();
0084   for (auto hv : hvmap)
0085   {
0086     int hvmod = hv.first;
0087     int pmt = hv.second;
0088     double gain = mcal->get_qgain(pmt);
0089     std::cout << hvmod << "\t" << pmt << "\t" << gain << std::endl;
0090
0091     int n = g_gainsbyhv[hvmod]->GetN();
0092     g_gainsbyhv[hvmod]->SetPoint(n, (double) npmts, gain);
0093
0094     if (gain > gmax[hvmod])
0095     {
0096       gmax[hvmod] = gain;
0097       maxpmt[hvmod] = pmt;
0098     }
0099     if (gain < gmin[hvmod])
0100     {
0101       gmin[hvmod] = gain;
0102       minpmt[hvmod] = pmt;
0103     }
0104
0105     npmts++;
0106   }
0107
0108   double &overall_gmin = *std::min_element(gmin.begin(), gmin.end());
0109   double &overall_gmax = *std::max_element(gmax.begin(), gmax.end());
0110
0111   // Draw the gains by HVMOD
0112   //  TCanvas *ac = new TCanvas("c_gainsbyhvmod","Gains by HV Mod",1250,600);
0113   g_gainsbyhv[0]->Draw("ap");
0114   g_gainsbyhv[0]->GetXaxis()->SetLimits(0, 128);
0115   g_gainsbyhv[0]->SetMaximum(overall_gmax * 1.1);
0116   g_gainsbyhv[0]->SetMinimum(overall_gmin * 0.9);
0117   gPad->Modified();
0118   gPad->Update();
0119   for (int ihv = 1; ihv < 16; ihv++)
0120   {
0121     g_gainsbyhv[ihv]->Draw("p");
0122   }
0123
0124   // return;
0125
0126   // calculate means, and the scaling to normalize all of them to the same
0127   std::array<double, 16> means{};
0128   means.fill(0.);
0129
0130   std::cout << "means in each hv group" << std::endl;
0131   std::array<double, 16> new_gmin{};
0132   std::array<double, 16> new_gmax{};
0133   std::array<double, 16> newscale{};
0134   double overall_newgmin{1e9};
0135   double overall_newgmax{0};
0136   int minhvmod = -1;  // hv mod with min gain
0137   int maxhvmod = -1;  // hv mod with max gain
0138   for (int ihv = 0; ihv < 16; ihv++)
0139   {
0140     means[ihv] = gmean(g_gainsbyhv[ihv]);
0141     newscale[ihv] = means[0] / means[ihv];
0142     new_gmin[ihv] = gmin[ihv] * newscale[ihv];
0143     new_gmax[ihv] = gmax[ihv] * newscale[ihv];
0144     std::cout << ihv << "\t" << means[ihv] << "\t" << newscale[ihv] << "\t" << new_gmin[ihv] << std::endl;
0145
0146     if (new_gmax[ihv] > overall_newgmax)
0147     {
0148       overall_newgmax = new_gmax[ihv];
0149       maxhvmod = ihv;
0150     }
0151     if (new_gmin[ihv] < overall_newgmin)
0152     {
0153       overall_newgmin = new_gmin[ihv];
0154       minhvmod = ihv;
0155     }
0156   }
0157
0158   std::cout << "MIN " << minpmt[minhvmod] << "\t" << minhvmod << "\t" << new_gmin[minhvmod] << std::endl;
0159   std::cout << "MAX " << maxpmt[maxhvmod] << "\t" << maxhvmod << "\t" << new_gmax[maxhvmod] << std::endl;
0160   std::cout << "SPREAD " << new_gmax[maxhvmod] / new_gmin[minhvmod] << std::endl;
0161
0162   // double minscale = 150./overall_newgmin;   // 150 is the min ADC we want in p+p
0163   // double minscale = 65./overall_newgmin;   // 65 is the min ADC we want in Au+Au
0164   double minscale = 80. / overall_newgmin;  // 80 is the min ADC we want in Au+Au
0165
0166   std::cout << "new means in each hv group" << std::endl;
0167   std::cout << "hvmod\tnewscale\tnew means [ADC/mip]" << std::endl;
0168   for (int ihv = 0; ihv < 16; ihv++)
0169   {
0170     newscale[ihv] *= minscale;
0171     std::cout << ihv << "\t" << newscale[ihv] << "\t" << means[ihv] * newscale[ihv] << std::endl;
0172   }
0173
0174   // Plot the gain curves by hvmod
0175   // TCanvas *bc = new TCanvas("c_gaincurvesbyhvmod","Gain curves by hvmod",1200,1200,.0001,.0001);
0176   TCanvas *bc = new TCanvas("c_gaincurvesbyhvmod", "Gain curves by hvmod", 1200, 1200);
0177   bc->Divide(4, 4);
0178   TFile *gcurvefile = new TFile("SCAN/mbdlaser_202405151654.root");
0179   TGraphErrors *g_norm_laserscan[128]{nullptr};
0180   TGraph *ginv_norm_laserscan[128]{nullptr};  // inverse curves
0181   for (int ipmt = 0; ipmt < 128; ipmt++)
0182   {
0183     name = "g_norm_laserscan";
0184     name += ipmt;
0185     g_norm_laserscan[ipmt] = (TGraphErrors *) gcurvefile->Get(name);
0186     int n = g_norm_laserscan[ipmt]->GetN();
0187     double *x = g_norm_laserscan[ipmt]->GetX();
0188     double *y = g_norm_laserscan[ipmt]->GetY();
0189     ginv_norm_laserscan[ipmt] = new TGraph(n, y, x);
0190     name = "ginv_norm_laserscan";
0191     name += ipmt;
0192     ginv_norm_laserscan[ipmt]->SetName(name);
0193   }
0194
0195   std::array<double, 16> ref_hvsetting{
0196       0.758629, 0.762746, 0.804495, 0.7953, 0.766233, 0.818163, 0.797741, 0.775627,   // south
0197       0.819641, 0.798468, 0.860734, 0.798129, 0.874018, 0.793739, 0.873797, 0.852621  // north
0198   };
0199   // ref_hvsetting.fill(0.81);      // the gain settings that were used for run analyzed
0200   std::array<double, 16> sumsetting{};
0201   sumsetting.fill(0.);
0202   std::array<int, 16> npmts_in_hvmod{};
0203   npmts_in_hvmod.fill(0.);
0204   std::array<double, 16> newsetting{};
0205   newsetting.fill(0.);
0206
0207   for (auto hv : hvmap)
0208   {
0209     int hvmod = hv.first;
0210     int pmt = hv.second;
0211
0212     bc->cd(hvmod + 1);
0213     if (npmts_in_hvmod[hvmod] == 0)
0214     {
0215       g_norm_laserscan[pmt]->Draw("acp");
0216     }
0217     else
0218     {
0219       g_norm_laserscan[pmt]->Draw("cp");
0220     }
0221     ++npmts_in_hvmod[hvmod];
0222
0223     // std::cout << "ref_hvset " << hvmod << "\t" << ref_hvsetting[hvmod] << std::endl;
0224     double ref_gain = g_norm_laserscan[pmt]->Eval(ref_hvsetting[hvmod]);
0225     double target_gain = newscale[hvmod] * ref_gain;
0226
0227     sumsetting[hvmod] += ginv_norm_laserscan[pmt]->Eval(target_gain);
0228     // std::cout << "xxx " << pmt << "\t" << ref_gain << "\t" << target_gain << "\t" << ginv_norm_laserscan[pmt]->Eval( target_gain ) << std::endl;
0229   }
0230
0231   std::cout << "New Settings" << std::endl;
0232   for (int ihv = 0; ihv < 16; ihv++)
0233   {
0234     newsetting[ihv] = sumsetting[ihv] / npmts_in_hvmod[ihv];
0235     std::cout << ihv << "\t" << newsetting[ihv] << std::endl;
0236   }
0237
0238   // to plug into hvscale.py
0239   std::cout << "hvscale = [ ";
0240   // first the north
0241   for (int ihv = 8; ihv < 16; ihv++)
0242   {
0243     std::cout << newsetting[ihv] << ", ";
0244   }
0245   for (int ihv = 0; ihv < 8; ihv++)
0246   {
0247     std::cout << newsetting[ihv] << ", ";
0248   }
0249   std::cout << " ]" << std::endl;
0250 }