mbd/HV/analaserscan.C

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