mbd/HV/analaserscan.C

0001 //
0002 // Plot laser amplitudes, make gain curves, extract HV to run gain at
0003 //
0004 #include "read_dstmbd.C"
0005
0006 #include <TH1.h>
0007 #include <TF1.h>
0008 #include <TCanvas.h>
0009 #include <TTree.h>
0010 #include <TFile.h>
0011 #include <TGraphErrors.h>
0012
0013 #include <iostream>
0014 #include <fstream>
0015 #include <vector>
0016
0017
0018 const int MAXRUNS = 100;    // 11 is the typical
0019 const int MAXCH = 256;
0020 const int N_PMT = 128;
0021 const int MAXBOARDS = MAXCH/16;     // FEM boards
0022 int NCH;                 // total number of channels (including charge channels)
0023 int NBOARDS;
0024
0025 TH1 *h_laseramp[MAXRUNS][MAXCH];    //[run][ch]
0026 TF1 *fgaus[MAXRUNS][MAXCH];
0027 int nrun = 0;
0028 int run_number[MAXRUNS];
0029 float hvscale[MAXRUNS];             // Scale from final Au+Au
0030 float hvset[MAXRUNS];               // HV setting
0031 float laseramp[MAXCH][MAXRUNS];
0032 float laseramperr[MAXCH][MAXRUNS];
0033 float norm_laseramp[MAXCH][MAXRUNS];
0034 float norm_laseramperr[MAXCH][MAXRUNS];
0035 TGraphErrors *g_laserscan[MAXCH];
0036 TGraphErrors *g_norm_laserscan[MAXCH];
0037
0038 TF1 *f_gain[N_PMT];
0039
0040 const int verbose = 1;
0041
0042 TString name;
0043 TString title;
0044
0045 //
0046 //void anafile(const char *tfname = "prdf_478_times.root", const int nrun = 0)
0047 void anafile(const char *tfname = "DST_UNCALMBD-00042674-0000.root", const int nrun_local = 0)
0048 {
0049   std::cout << "tfname " << tfname << std::endl;
0050
0051   // Book Histograms, etc
0052   for (int ich=0; ich<N_PMT; ich++)
0053   {
0054     name = "h_laseramp"; name += ich; name += "_"; name += nrun_local;
0055     title = name;
0056     h_laseramp[nrun_local][ich] = new TH1F(name,title,1620,-100,16100);
0057     //h_laseramp[nrun_local][ich] = new TH1F(name,title,160,0,16000);
0058   }
0059
0060   /*
0061   Int_t   evt;
0062   Float_t f_tt[MAXCH];  // time, t-ch
0063   Float_t f_tq[MAXCH];  // time, q-ch
0064   Float_t f_q[MAXCH];   // charge
0065   */
0066
0067   std::cout << "tfname " << tfname << std::endl;
0068
0069   read_dstmbd(tfname);
0070   /*
0071   TFile *tfile = new TFile(tfname,"READ");
0072   TTree *tree = (TTree*)tfile->Get("t");
0073   tree->SetBranchAddress("evt",&evt);
0074   for (int ich=0; ich<N_PMT; ich++)
0075   {
0076     name = "tt"; name += ich;
0077     tree->SetBranchAddress(name,&f_tt[ich]);
0078     name = "tq"; name += ich;
0079     tree->SetBranchAddress(name,&f_tq[ich]);
0080     name = "q"; name += ich;
0081     tree->SetBranchAddress(name,&f_q[ich]);
0082   }
0083   */
0084
0085   int nentries = tree->GetEntries();
0086   // skip 1st 70 because those are possibly bad
0087   for (int ientry=70; ientry<nentries; ientry++)
0088   {
0089     tree->GetEntry(ientry);
0090
0091     for (int ich=0; ich<N_PMT; ich++)
0092     {
0093       //std::cout << evt << "\t" << t[1] << "\t" << f_q[1] << "\t" << t[14] << "\t" << f_q[14] << std::endl;
0094       if ( f_q[ich]>0 )
0095       {
0096         h_laseramp[nrun_local][ich]->Fill( f_q[ich] );
0097       }
0098     }
0099   }
0100
0101   // Fit gaussians to get the average amplitude
0102   for (int ich=0; ich<N_PMT; ich++)
0103   {
0104     //std::cout << "Fitting ch " << ich << std::endl;
0105     name = "fgaus"; name += ich; name += "_"; name += nrun_local;
0106     fgaus[nrun_local][ich] = new TF1(name,"gaus",-100,16000);
0107     fgaus[nrun_local][ich]->SetLineColor(2);
0108
0109     //h_laseramp[nrun_local][ich]->Rebin(2);
0110
0111     float ampl = 0.;
0112     float amplerr = 0.;
0113     if ( h_laseramp[nrun_local][ich]->GetEntries() > 10 )
0114     {
0115       double seed_mean = h_laseramp[nrun_local][ich]->GetMean();
0116       double seed_rms = h_laseramp[nrun_local][ich]->GetRMS();
0117       fgaus[nrun_local][ich]->SetParameters(1000,seed_mean,seed_rms);
0118       fgaus[nrun_local][ich]->SetRange( seed_mean-5*seed_rms, seed_mean+5*seed_rms );
0119       //h_laseramp[nrun_local][ich]->Fit( fgaus[nrun_local][ich], "QR" );
0120       h_laseramp[nrun_local][ich]->Fit( fgaus[nrun_local][ich], "LLRQ" );
0121     }
0122
0123     ampl = fgaus[nrun_local][ich]->GetParameter(1);
0124     amplerr = fgaus[nrun_local][ich]->GetParError(1);
0125     laseramp[ich][nrun_local] = ampl;
0126     laseramperr[ich][nrun_local] = amplerr;
0127   }
0128
0129   if ( verbose && nrun_local==0 )
0130   {
0131     const int NCANVAS = 4;
0132     TCanvas *c_laseramp[NCANVAS];
0133     for (int icv = 0; icv<NCANVAS; icv++)
0134     {
0135       name = "c_laser"; name += icv;
0136       c_laseramp[icv] = new TCanvas(name,name,1600,800);
0137       c_laseramp[icv]->Divide(8,4);
0138     }
0139
0140     for (int ipmt=0; ipmt<N_PMT; ipmt++)
0141     {
0142       int quad = ipmt/32;    // quadrant
0143       c_laseramp[quad]->cd( ipmt%32 + 1 );
0144       //h_laseramp[nrun_local][ipmt]->Rebin(100);
0145       h_laseramp[nrun_local][ipmt]->Draw();
0146       /*
0147       gPad->SetLogy(1);
0148       gPad->Modified();
0149       gPad->Update();
0150       */
0151     }
0152
0153   }
0154
0155 }
0156
0157 float prev_ampl[N_PMT];
0158
0159 void get_prev_gains(const char *prev_gains_file = "gains_at_9_zerofield.out")
0160 {
0161     std::ifstream infile( prev_gains_file );
0162     int ch;
0163     float hvscale_local;
0164     float ratio;
0165     while ( infile >> ch >> hvscale_local >> ratio )
0166     {
0167         if ( ch>=0 && ch<128 )
0168         {
0169           infile >> prev_ampl[ch];
0170           std::cout << ch << "\t" << std::endl;
0171         }
0172         else
0173         {
0174           std::cout << "Bad ch" << std::endl;
0175         }
0176     }
0177 }
0178
0179 void analaserscan(const char *fname = "hvscan.62880")
0180 {
0181   NCH = MAXCH;                 // total number of channels (including charge channels)
0182   NBOARDS = MAXBOARDS;
0183
0184   //get_prev_gains();
0185
0186   // Get the number of actual channels to process
0187   /*
0188   std::ifstream configfile("digsig.cfg");
0189   if ( configfile.is_open() )
0190   {
0191     string junk;
0192     configfile >> junk >> NCH;
0193     NBOARDS = NCH/16;
0194
0195     std::cout << "Found config file digsig.cfg" << std::endl;
0196     std::cout << "Setting NCH = " << NCH << std::endl;
0197   }
0198   */
0199
0200   TString rootfname;
0201   std::ifstream scanfile(fname);
0202   int nruns = 0;
0203   while ( scanfile >> run_number[nruns] >> hvscale[nruns] )
0204   {
0205     std::cout << run_number[nruns] << "\t" << hvscale[nruns] << std::endl;
0206
0207     // get name of root file
0208     //rootfname.Form( "calib_mbd-%08d-0000_mbd.root", run_number[nruns] );
0209     rootfname.Form( "DST_UNCALMBD-%08d-0000.root", run_number[nruns] ); // NOLINT(hicpp-vararg)
0210     std::cout << "Processing " << rootfname << std::endl;
0211
0212     anafile( rootfname, nruns );
0213
0214     nruns++;
0215   }
0216   std::cout << "Processed " << nruns << " runs" << std::endl;
0217
0218   // Make the canvases for the gain curves
0219   const int NCANVAS = 4;
0220   TCanvas *c_laserscan[NCANVAS];
0221   for (int icv = 0; icv<NCANVAS; icv++)
0222   {
0223     name = "c_laserscan"; name += icv;
0224     c_laserscan[icv] = new TCanvas(name,name,1600,800);
0225     c_laserscan[icv]->Divide(8,4);
0226   }
0227
0228   // Make and draw the gain curves, and
0229   // write out the HV scale needed to get XX gain
0230   name = "mbdlaser_"; name += fname; name += ".root";
0231   name.ReplaceAll("hvscan.","");
0232   TFile *savefile = new TFile(name,"RECREATE");
0233   std::ofstream gainfile("gains.out");
0234   std::ofstream gain9file("gains_at_9.out");
0235 //  float gain_wanted = 0.55;
0236   for (int ipmt=0; ipmt<N_PMT; ipmt++)
0237   {
0238     int quad = ipmt/32;  // quadrant
0239
0240     // Create the TGraph with the summary results
0241     name = "g_laserscan"; name += ipmt;
0242     title = "ch "; title += ipmt;
0243     g_laserscan[ipmt] = new TGraphErrors(nruns,hvscale,laseramp[ipmt],nullptr,laseramperr[ipmt]);
0244     g_laserscan[ipmt]->SetName(name);
0245     g_laserscan[ipmt]->SetTitle(title);
0246     g_laserscan[ipmt]->SetMarkerStyle(20);
0247     g_laserscan[ipmt]->SetMarkerSize(0.5);
0248
0249     //std::cout << "ch " << pmtch << std::endl;
0250     c_laserscan[quad]->cd( ipmt%32 + 1 );
0251     g_laserscan[ipmt]->Draw("acp");
0252     g_laserscan[ipmt]->GetHistogram()->SetTitleSize(4);
0253
0254     name = "f_gain"; name += ipmt;
0255     f_gain[ipmt] = new TF1(name,"[0]*exp([1]*x)-1",0,1.0);
0256     f_gain[ipmt]->SetParameters(1,1);
0257     //g_laserscan[ipmt]->Fit( f_gain[ipmt] );
0258
0259     // now step down until we find the gain we want
0260 //    double maxgain = g_laserscan[ipmt]->Eval(1.0);
0261
0262     /*
0263     double ref_gain = g_laserscan[ipmt]->Eval(0.81)/maxgain;
0264     gain9file << ipmt << "\t" << 0.81 << "\t" << ref_gain << "\t" << g_laserscan[ipmt]->Eval(0.81) << std::endl;
0265
0266     for (double iv=0.5; iv<1.0; iv+=0.001)
0267     {
0268         double ampl = g_laserscan[ipmt]->Eval(iv);
0269         //double ratio = g_laserscan[ipmt]->Eval(iv)/maxgain;
0270         //if ( ratio>gain_wanted )
0271         //{
0272         //    gainfile << ipmt << "\t" << iv << "\t" << ratio << std::endl;
0273         //    break;
0274         //}
0275
0276         if ( ampl > prev_ampl[ipmt] )
0277         {
0278             gainfile << ipmt << "\t" << iv << "\t" << ampl/prev_ampl[ipmt] << std::endl;
0279             break;
0280         }
0281     }
0282     */
0283   }
0284
0285   gainfile.close();
0286
0287   // Now create the normalized gain curves
0288   for (int ipmt=0; ipmt<N_PMT; ipmt++)
0289   {
0290     double ref_ampl = g_laserscan[ipmt]->Eval(1.0);
0291     for (int irun=0; irun<nruns; irun++)
0292     {
0293       norm_laseramp[ipmt][irun] = laseramp[ipmt][irun]/ref_ampl;
0294       norm_laseramperr[ipmt][irun] = laseramperr[ipmt][irun]/ref_ampl;
0295     }
0296     name = "g_norm_laserscan"; name += ipmt;
0297     title = "ch "; title += ipmt; title += ", norm";
0298     g_norm_laserscan[ipmt] = new TGraphErrors(nruns,hvscale,norm_laseramp[ipmt],nullptr,norm_laseramperr[ipmt]);
0299     g_norm_laserscan[ipmt]->SetName(name);
0300     g_norm_laserscan[ipmt]->SetTitle(title);
0301     g_norm_laserscan[ipmt]->SetMarkerStyle(20);
0302     g_norm_laserscan[ipmt]->SetMarkerSize(0.5);
0303   }
0304
0305   for (int ipmt=0; ipmt<N_PMT; ipmt++)
0306   {
0307     g_laserscan[ipmt]->Write();
0308     g_norm_laserscan[ipmt]->Write();
0309   }
0310
0311   savefile->Write();
0312
0313   // dump full gainfile
0314   std::ofstream gain2file("gains_full.csv");
0315
0316   // Header
0317   gain2file << "iv,";
0318   for (int ipmt=0; ipmt<N_PMT; ipmt++)
0319   {
0320       gain2file << ipmt << ",";
0321   }
0322   gain2file << std::endl;
0323
0324   // Gains
0325   for (double iv=1.0; iv>0.5; iv -= 0.05)
0326   {
0327     gain2file << iv << ",";
0328     for (int ipmt=0; ipmt<N_PMT; ipmt++)
0329     {
0330         gain2file << g_norm_laserscan[ipmt]->Eval(iv) << ",";
0331     }
0332     gain2file << std::endl;
0333   }
0334   gain2file.close();
0335
0336 }
0337