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