mbd/HV/rearrange_pmts.C

0001 // Make plots of the values of the gains
0002 // Use this to figure out which PMTs to swap around to normalize gains
0003 #include <mbd/MbdCalib.h>
0004 #include <mbd/MbdGeomV2.h>
0005
0006
0007 #include <TCanvas.h>
0008 #include <TFile.h>
0009 #include <TGraph.h>
0010 #include <TGraphErrors.h>
0011 #include <TH1.h>
0012 #include <TPad.h>
0013 #include <TString.h>
0014
0015 R__LOAD_LIBRARY(libmbd.so)
0016 R__LOAD_LIBRARY(libmbd_io.so)
0017
0018 int nruns = 0;
0019 MbdCalib *mcal{nullptr};
0020 MbdGeomV2 *mbdgeom{nullptr};
0021
0022 void read_calibgains(const char *fname)
0023 {
0024   mcal = new MbdCalib();
0025   mcal->Download_Gains( fname );
0026 }
0027
0028
0029 TGraphErrors *g_gainsbyhv[16]{nullptr};
0030
0031 double gmean(TGraph *g)
0032 {
0033   int n = g->GetN();
0034   Double_t *y = g->GetY();
0035   double sum = 0.;
0036   for (int i=0; i<n; i++)
0037   {
0038     sum += y[i];
0039   }
0040   double avg = sum / n;
0041
0042   return avg;
0043 }
0044
0045
0046 void rearrange_pmts(const char *fname = "results/54935/mbd_qfit.calib")
0047 {
0048   mbdgeom = new MbdGeomV2();
0049
0050   read_calibgains(fname);
0051
0052   TString name;
0053   for (int ihv=0; ihv<16; ihv++)
0054   {
0055     name = "g_gainsbyhv"; name += ihv;
0056     g_gainsbyhv[ihv] = new TGraphErrors();
0057     g_gainsbyhv[ihv]->SetName( name );
0058     g_gainsbyhv[ihv]->SetLineColor( ihv%8 + 1 );
0059     g_gainsbyhv[ihv]->SetMarkerColor( ihv%8 + 1 );
0060     /*
0061     if ( ihv+1 == 10 )
0062     {
0063       g_gainsbyhv[ihv]->SetLineColor( 50 );
0064       g_gainsbyhv[ihv]->SetMarkerColor( 50 );
0065     }
0066     */
0067     g_gainsbyhv[ihv]->SetMarkerStyle( 20 );
0068     g_gainsbyhv[ihv]->GetHistogram()->SetXTitle( "hv groups" );
0069     g_gainsbyhv[ihv]->GetHistogram()->SetYTitle( "gain [ADC/mip]" );
0070   }
0071
0072
0073   int npmts = 0;
0074   std::array<double,16> gmin{};
0075   std::array<double,16> gmax{};
0076   std::array<int,16> minpmt{};
0077   std::array<int,16> maxpmt{};
0078   gmin.fill(1e9);      // the gain settings that were used for run analyzed
0079   gmax.fill(0);
0080   minpmt.fill(-1);
0081   maxpmt.fill(-1);
0082   std::multimap hvmap = mbdgeom->get_hvmap();
0083   std::cout << "\thvmod\tpmt\tgain\tr" << std::endl;
0084   for ( auto hv : hvmap )
0085   {
0086     int hvmod = hv.first;
0087     int pmt = hv.second;
0088     double gain = mcal->get_qgain( pmt );
0089
0090     // print out the outliers
0091     if ( gain>250 || gain<100 )
0092     {
0093       float r = mbdgeom->get_r( pmt );
0094       std::cout << "outlier\t"<< hvmod << "\t" << pmt << "\t" << gain << "\t" << r << std::endl;
0095     }
0096
0097     int n = g_gainsbyhv[hvmod]->GetN();
0098     g_gainsbyhv[hvmod]->SetPoint(n,(double)npmts,gain);
0099
0100     if ( gain>gmax[hvmod] )
0101     {
0102       gmax[hvmod] = gain;
0103       maxpmt[hvmod] = pmt;
0104     }
0105     if ( gain<gmin[hvmod] )
0106     {
0107       gmin[hvmod] = gain;
0108       minpmt[hvmod] = pmt;
0109     }
0110
0111     npmts++;
0112   }
0113
0114   double& overall_gmin = *std::min_element(gmin.begin(), gmin.end());
0115   double& overall_gmax = *std::max_element(gmax.begin(), gmax.end());
0116
0117   // Draw the gains by HVMOD
0118 //  TCanvas *ac = new TCanvas("c_gainsbyhvmod","Gains by HV Mod",1250,600);
0119   g_gainsbyhv[0]->Draw("ap");
0120   g_gainsbyhv[0]->GetXaxis()->SetLimits(0,128);
0121   g_gainsbyhv[0]->SetMaximum( overall_gmax*1.1 );
0122   g_gainsbyhv[0]->SetMinimum( overall_gmin*0.9 );
0123   gPad->Modified();
0124   gPad->Update();
0125   for (int ihv=1; ihv<16; ihv++)
0126   {
0127     g_gainsbyhv[ihv]->Draw("p");
0128   }
0129
0130   //return;
0131
0132   /*
0133   // calculate means, and the scaling to normalize all of them to the same
0134   std::array<double,16> means;
0135   means.fill(0.);
0136
0137   std::cout << "means in each hv group" << std::endl;
0138   std::array<double,16> new_gmin;
0139   std::array<double,16> new_gmax;
0140   std::array<double,16> newscale;
0141   double overall_newgmin{1e9};
0142   double overall_newgmax{0};
0143   int   minhvmod = -1;  // hv mod with min gain
0144   int   maxhvmod = -1;  // hv mod with max gain
0145   for (int ihv=0; ihv<16; ihv++)
0146   {
0147     means[ihv] = gmean( g_gainsbyhv[ihv] );
0148     newscale[ihv] = means[0]/means[ihv];
0149     new_gmin[ihv] = gmin[ihv]*newscale[ihv];
0150     new_gmax[ihv] = gmax[ihv]*newscale[ihv];
0151     std::cout << ihv << "\t" << means[ihv] << "\t" << newscale[ihv] << "\t" << new_gmin[ihv] << std::endl;
0152
0153     if ( new_gmax[ihv]>overall_newgmax )
0154     {
0155       overall_newgmax = new_gmax[ihv];
0156       maxhvmod = ihv;
0157     }
0158     if ( new_gmin[ihv]<overall_newgmin )
0159     {
0160       overall_newgmin = new_gmin[ihv];
0161       minhvmod = ihv;
0162     }
0163   }
0164   */
0165
0166   // Starting voltages (from last known values for PHENIX BBC)
0167   std::array<double,16> orig_hv {
0168     -1950.0, -2150.0, -2000.0, -2308.6, -2300.0, -1668.8, -2000.0, -2200.0,      // south
0169     -1850.0, -2400.0, -1950.0, -2170.0, -1612.0, -2350.0, -2050.0, -1811.5       // north
0170   };
0171
0172   // Settings (relative) for HV in Run2 Au+Au
0173   std::array<double,16> ref_hvsetting {
0174       0.695889, 0.69948, 0.745683, 0.716363, 0.705564, 0.762345, 0.757273, 0.712111,   // south
0175       0.738386, 0.726954, 0.779238, 0.725638, 0.797808, 0.70735, 0.786134, 0.763983    // north
0176   };
0177
0178   std::array<double,16> hv_setting{};
0179
0180   std::cout << "HV SETTING" << std::endl;
0181   for (int imod=0; imod<orig_hv.size(); imod++)
0182   {
0183     hv_setting[imod] = orig_hv[imod]*ref_hvsetting[imod];
0184     std::cout << hv_setting[imod] << "\t";
0185     if ( imod%8==7 ) { std::cout << std::endl;
0186 }
0187   }
0188
0189   // Plot the gain curves by hvmod
0190   //TCanvas *bc = new TCanvas("c_gaincurvesbyhvmod","Gain curves by hvmod",1200,1200,.0001,.0001);
0191 //  TCanvas *bc = new TCanvas("c_gaincurvesbyhvmod","Gain curves by hvmod",1200,1200);
0192   //bc->Divide(4,4);
0193   TFile *gcurvefile = new TFile("SCAN/mbdlaser_202405151654.root");
0194   TGraphErrors *g_norm_laserscan[128]{nullptr};
0195   TGraph *ginv_norm_laserscan[128]{nullptr};  // inverse curves
0196   for (int ipmt=0; ipmt<128; ipmt++)
0197   {
0198     name = "g_norm_laserscan"; name += ipmt;
0199     g_norm_laserscan[ipmt] = (TGraphErrors*)gcurvefile->Get( name );
0200   }
0201
0202   std::array<double,16> sumsetting{};
0203   sumsetting.fill(0.);
0204   std::array<int,16> npmts_in_hvmod{};
0205   npmts_in_hvmod.fill(0.);
0206   std::array<double,16> newsetting{};
0207   newsetting.fill(0.);
0208
0209   int flag = 0;
0210   for ( auto hv : hvmap )
0211   {
0212     int hvmod = hv.first;
0213     int pmt = hv.second;
0214
0215     // convert x axis to HV and then create inverted graph
0216     int n = g_norm_laserscan[pmt]->GetN();
0217     double *x = g_norm_laserscan[pmt]->GetX();
0218     double *y = g_norm_laserscan[pmt]->GetY();
0219     for (int ip=0; ip<n; ip++)
0220     {
0221       x[ip] *= -orig_hv[hvmod];
0222     }
0223
0224     ginv_norm_laserscan[pmt] = new TGraph(n,y,x);
0225     name = "ginv_norm_laserscan"; name += pmt;
0226     ginv_norm_laserscan[pmt]->SetName( name );
0227
0228     //bc->cd( hvmod+1 );
0229     if ( flag == 0 )
0230     {
0231       g_norm_laserscan[pmt]->Draw("acp");
0232       g_norm_laserscan[pmt]->GetXaxis()->SetLimits(800.,2450.);
0233       gPad->Modified();
0234       gPad->Update();
0235       flag = 1;
0236     }
0237     else
0238     {
0239       g_norm_laserscan[pmt]->Draw("cp");
0240     }
0241
0242     // See what happens when we change HV
0243     double gain = mcal->get_qgain( pmt );
0244     if ( pmt==17 )
0245     {
0246       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1503. );
0247       float new_hv = g_norm_laserscan[pmt]->Eval( 1356. );
0248       float new_gain = (new_hv/orig_hv_loc)*gain;
0249       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0250     }
0251     if ( pmt==26 )
0252     {
0253       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1491. );
0254       float new_hv = g_norm_laserscan[pmt]->Eval( 1356. );
0255       float new_gain = (new_hv/orig_hv_loc)*gain;
0256       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0257     }
0258     if ( pmt==36 )
0259     {
0260       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1622. );
0261       float new_hv = g_norm_laserscan[pmt]->Eval( 1653. );
0262       float new_gain = (new_hv/orig_hv_loc)*gain;
0263       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0264     }
0265     if ( pmt==51 )
0266     {
0267       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1514. );
0268       //float new_hv = g_norm_laserscan[pmt]->Eval( 1566. );
0269       float new_hv = g_norm_laserscan[pmt]->Eval( 1622. );
0270       float new_gain = (new_hv/orig_hv_loc)*gain;
0271       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0272     }
0273     if ( pmt==70 )
0274     {
0275       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1366. );
0276       float new_hv = g_norm_laserscan[pmt]->Eval( 1519. );
0277       float new_gain = (new_hv/orig_hv_loc)*gain;
0278       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0279     }
0280     if ( pmt==79 )
0281     {
0282       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1574. );
0283       float new_hv = g_norm_laserscan[pmt]->Eval( 1611. );
0284       float new_gain = (new_hv/orig_hv_loc)*gain;
0285       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0286     }
0287     if ( pmt==85 )
0288     {
0289       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1574. );
0290       float new_hv = g_norm_laserscan[pmt]->Eval( 1662. );
0291       float new_gain = (new_hv/orig_hv_loc)*gain;
0292       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0293     }
0294     if ( pmt==87 )
0295     {
0296       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1774. );
0297       float new_hv = g_norm_laserscan[pmt]->Eval( 1662. );
0298       float new_gain = (new_hv/orig_hv_loc)*gain;
0299       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0300     }
0301     if ( pmt==89 )
0302     {
0303       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1519. );
0304       float new_hv = g_norm_laserscan[pmt]->Eval( 1366. );
0305       float new_gain = (new_hv/orig_hv_loc)*gain;
0306       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0307     }
0308     if ( pmt==91 )
0309     {
0310       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1774. );
0311       float new_hv = g_norm_laserscan[pmt]->Eval( 1662. );
0312       float new_gain = (new_hv/orig_hv_loc)*gain;
0313       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0314     }
0315     if ( pmt==107 )
0316     {
0317       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1611. );
0318       float new_hv = g_norm_laserscan[pmt]->Eval( 1662. );
0319       float new_gain = (new_hv/orig_hv_loc)*gain;
0320       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0321     }
0322     if ( pmt==113 )
0323     {
0324       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1383. );
0325       float new_hv = g_norm_laserscan[pmt]->Eval( 1286. );
0326       float new_gain = (new_hv/orig_hv_loc)*gain;
0327       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0328     }
0329     if ( pmt==116 )
0330     {
0331       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1611. );
0332       float new_hv = g_norm_laserscan[pmt]->Eval( 1574. );
0333       float new_gain = (new_hv/orig_hv_loc)*gain;
0334       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0335     }
0336     if ( pmt==121 )
0337     {
0338       float orig_hv_loc = g_norm_laserscan[pmt]->Eval( 1611. );
0339       //float new_hv = g_norm_laserscan[pmt]->Eval( 1519. );
0340       float new_hv = g_norm_laserscan[pmt]->Eval( 1574. );
0341       float new_gain = (new_hv/orig_hv_loc)*gain;
0342       std::cout << "PMT " << pmt << "\t" << orig_hv_loc << "\t" << new_hv << "\t" << new_hv/orig_hv_loc << "\t" << new_gain << std::endl;
0343     }
0344
0345   }
0346
0347 }