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