BBCStudy/macros/plot_checks.C

0001 #include "get_runstr.h"
0002 #include <iostream>
0003 #include <fstream>
0004 #include <TGraphErrors.h>
0005 #include <TString.h>
0006
0007 int do_zdc = 0;
0008 int do_emcal = 0;
0009 int do_ohcal = 0;
0010 int do_ihcal = 0;
0011 int calc_gain_corr = 1;
0012
0013 const int MAXRUNS = 20000;
0014
0015 int nruns = 0;
0016
0017 void fit_across_runs(TGraphErrors *gainvals[], const char *flist);
0018 void plot_bz();
0019
0020 Double_t bqmeansum[MAXRUNS];            // bq mean in sum of two arms
0021 Double_t bqmeansumerr[MAXRUNS];
0022 Double_t bqmeantot[2][MAXRUNS];         // bq mean in each arm
0023 Double_t bqmeantoterr[2][MAXRUNS];
0024 Double_t bqmean[128][MAXRUNS];
0025 Double_t bqmeanerr[128][MAXRUNS];
0026 Double_t runindex[MAXRUNS];
0027 TGraphErrors *g_meansum {nullptr};
0028 TGraphErrors *g_meantot[2] {nullptr};
0029 TGraphErrors *g_mean[128] {nullptr};
0030
0031 TFile *tfile[MAXRUNS];
0032 TH1 *h_bbcqsum[MAXRUNS];
0033 TH1 *h_bbcqtot[2][MAXRUNS];
0034 TH1 *h_bbcq[MAXRUNS][128];
0035 TH2 *h2_bbcq[128];
0036
0037 TH1 *h_zdce[MAXRUNS];
0038 TH1 *h_emcale[MAXRUNS];
0039 TH1 *h_ohcale[MAXRUNS];
0040 TH1 *h_ihcale[MAXRUNS];
0041
0042 TH1 *h_bz[MAXRUNS];
0043
0044 TCanvas *cv[100];
0045 int icv = 0;
0046
0047 // Get mean gain for each channel across runs, and write out the values
0048 void fit_across_runs(TGraphErrors *gainvals[], const char *flist)
0049 {
0050   TString meangainfname = "results/"; meangainfname.Append( flist );
0051   meangainfname.ReplaceAll(".list","_qmean.calib");
0052   ofstream meangainfile( meangainfname );
0053   TF1 *f_meangain[128] {nullptr};
0054   Double_t grp_mean[128];
0055   Double_t grp_meanerr[128];
0056   for (int ipmt=0; ipmt<128; ipmt++)
0057   {
0058     TString name = "f_meangain"; name += ipmt;
0059     f_meangain[ipmt] = new TF1(name,"pol0",0,1e9);
0060     gainvals[ipmt]->Fit(f_meangain[ipmt]);
0061     gPad->Update();
0062     gPad->Modified();
0063     grp_mean[ipmt] = f_meangain[ipmt]->GetParameter(0);
0064     grp_meanerr[ipmt] = f_meangain[ipmt]->GetParError(0);
0065
0066     meangainfile << ipmt << "\t" << grp_mean[ipmt] << "\t" << grp_meanerr[ipmt] << endl;
0067     cout << ipmt << "\t" << grp_mean[ipmt] << "\t" << grp_meanerr[ipmt] << endl;
0068   }
0069   meangainfile.close();
0070 }
0071
0072
0073 TH2 *h2_bzvsirun{nullptr};
0074 void plot_bz()
0075 {
0076   int nbinsx = h_bz[0]->GetNbinsX();
0077   double xmin = h_bz[0]->GetBinLowEdge(1);
0078   double xmax = h_bz[0]->GetBinLowEdge(nbinsx+1);
0079   h2_bzvsirun = new TH2F("h2_bzvsirun","MBD_z vs run idx",nruns,0,nruns,nbinsx,xmin,xmax);
0080
0081   for (int irun=0; irun<nruns; irun++)
0082   {
0083     if ( h_bz[irun]->Integral() == 0 )
0084     {
0085       cout << irun << endl;
0086     }
0087     for (int ibin=1; ibin<=nbinsx; ibin++)
0088     {
0089       float val = h_bz[irun]->GetBinContent(ibin);
0090       h2_bzvsirun->SetBinContent(irun+1,ibin,val);
0091     }
0092   }
0093
0094   cv[icv] = new TCanvas("bz","bz",1000,600);
0095   h2_bzvsirun->Draw("colz");
0096 }
0097
0098 void plot_checks(const char *flistname = "c.list")
0099 {
0100
0101   ifstream flist(flistname);
0102
0103   TString dir = "results/plot_checks_"; dir += flistname;
0104   dir.ReplaceAll(".list","/");
0105   TString name = "mkdir -p " + dir;
0106   cout << name << endl;
0107   gSystem->Exec( name );
0108
0109   TString dstfname;
0110   while ( flist >> dstfname )
0111   {
0112     tfile[nruns] = new TFile(dstfname,"READ");
0113
0114     h_bz[nruns] = (TH1*)tfile[nruns]->Get("h_bz");
0115
0116     h_bbcqsum[nruns] = (TH1*)tfile[nruns]->Get("h_bbcqsum");
0117     h_bbcqtot[0][nruns] = (TH1*)tfile[nruns]->Get("h_bbcqtot0");
0118     h_bbcqtot[1][nruns] = (TH1*)tfile[nruns]->Get("h_bbcqtot1");
0119     bqmeansum[nruns] = h_bbcqsum[nruns]->GetMean();
0120     bqmeansumerr[nruns] = h_bbcqsum[nruns]->GetMeanError();
0121     bqmeantot[0][nruns] = h_bbcqtot[0][nruns]->GetMean();
0122     bqmeantoterr[0][nruns] = h_bbcqtot[0][nruns]->GetMeanError();
0123     bqmeantot[1][nruns] = h_bbcqtot[0][nruns]->GetMean();
0124     bqmeantoterr[1][nruns] = h_bbcqtot[0][nruns]->GetMeanError();
0125
0126     h_zdce[nruns] = (TH1*)tfile[nruns]->Get("h_zdce");
0127     h_emcale[nruns] = (TH1*)tfile[nruns]->Get("h_emcale");
0128     h_ohcale[nruns] = (TH1*)tfile[nruns]->Get("h_ohcale");
0129     h_ihcale[nruns] = (TH1*)tfile[nruns]->Get("h_ihcale");
0130
0131     for (int ipmt=0; ipmt<128; ipmt++)
0132     {
0133       name = "h_q"; name += ipmt;
0134       h_bbcq[nruns][ipmt] = (TH1*)tfile[nruns]->Get(name);
0135       bqmean[ipmt][nruns] = h_bbcq[nruns][ipmt]->GetMean();
0136       bqmeanerr[ipmt][nruns] = h_bbcq[nruns][ipmt]->GetMeanError();
0137     }
0138
0139     runindex[nruns] = nruns;
0140     nruns++;
0141   }
0142   cout << "Processed " << nruns << " runs." << endl;
0143
0144   plot_bz();
0145   return;
0146
0147   TCanvas *bc = new TCanvas("bc","mean fits",800,600);
0148   TString title;
0149   TString meangainfname = dir + flistname;
0150   meangainfname.ReplaceAll(".list","_qmean.calib");
0151   ofstream meanfile( meangainfname );
0152   TF1 *meanfit = new TF1("meanfit","pol0",0,nruns);
0153   double runmean[128];
0154   double runmeanerr[128];
0155   double chi2ndf[128];
0156   for (int ipmt=0; ipmt<128; ipmt++)
0157   {
0158     g_mean[ipmt] = new TGraphErrors(nruns,runindex,bqmean[ipmt],0,bqmeanerr[ipmt]);
0159     name = "q_mean"; name += ipmt;
0160     title = name;
0161     g_mean[ipmt]->SetName( name );
0162     g_mean[ipmt]->SetTitle( title );
0163
0164     g_mean[ipmt]->Draw("ap");
0165     meanfit->SetParameter(0,100);
0166     g_mean[ipmt]->Fit(meanfit,"R");
0167     gPad->Modified();
0168     gPad->Update();
0169
0170     runmean[ipmt] = meanfit->GetParameter(0);
0171     runmeanerr[ipmt] = meanfit->GetParError(0);
0172     Double_t chi2 = meanfit->GetChisquare();
0173     Double_t ndf = meanfit->GetNDF();
0174     chi2ndf[ipmt] = chi2/ndf;
0175
0176     meanfile << ipmt << "\t" << runmean[ipmt] << "\t" << runmeanerr[ipmt] << "\t" << chi2ndf[ipmt] << endl;
0177     cout << ipmt << "\t" << runmean[ipmt] << "\t" << runmeanerr[ipmt] << "\t" << chi2ndf[ipmt] << endl;
0178   }
0179   meanfile.close();
0180
0181   // Make vs run index plots
0182   TH2 *h2_bbcqsum_vs_run = new TH2F("h2_bbcqsum_vs_run","bbc qsum vs run",nruns,0,nruns,3000,0,3000);
0183   TH2 *h2_bbcqs_vs_run = new TH2F("h2_bbcqs_vs_run","bbc south qsum vs run",nruns,0,nruns,1400,0,1400);
0184   TH2 *h2_bbcqn_vs_run = new TH2F("h2_bbcqn_vs_run","bbc north qsum vs run",nruns,0,nruns,1400,0,1400);
0185   for (int ipmt=0; ipmt<128; ipmt++)
0186   {
0187     name = "h2_bbcq"; name += ipmt;
0188     h2_bbcq[ipmt] = new TH2F(name,name,nruns,0,nruns,1200,0,60);
0189   }
0190   TH2 *h2_zdce_vs_run {nullptr};
0191   TH2 *h2_emcale_vs_run {nullptr};
0192   TH2 *h2_ohcale_vs_run {nullptr};
0193   TH2 *h2_ihcale_vs_run {nullptr};
0194   if ( do_zdc )
0195   {
0196     int nbinsx = h_zdce[0]->GetNbinsX();
0197     Double_t xmin = h_zdce[0]->GetBinLowEdge(1);
0198     Double_t xmax = h_zdce[0]->GetBinLowEdge(nbinsx+1);
0199     h2_zdce_vs_run = new TH2F("h2_zdce_vs_run","ZDC E vs run",nruns,0,nruns,nbinsx,xmin,xmax);
0200   }
0201   if ( do_emcal )
0202   {
0203     int nbinsx = h_emcale[0]->GetNbinsX();
0204     Double_t xmin = h_emcale[0]->GetBinLowEdge(1);
0205     Double_t xmax = h_emcale[0]->GetBinLowEdge(nbinsx+1);
0206     h2_emcale_vs_run = new TH2F("h2_emcale_vs_run","EMC E vs run",nruns,0,nruns,nbinsx,xmin,xmax);
0207   }
0208   if ( do_ohcal )
0209   {
0210     int nbinsx = h_ohcale[0]->GetNbinsX();
0211     Double_t xmin = h_ohcale[0]->GetBinLowEdge(1);
0212     Double_t xmax = h_ohcale[0]->GetBinLowEdge(nbinsx+1);
0213     h2_ohcale_vs_run = new TH2F("h2_ohcale_vs_run","OHCAL E vs run",nruns,0,nruns,nbinsx,xmin,xmax);
0214   }
0215   if ( do_ihcal )
0216   {
0217     int nbinsx = h_ihcale[0]->GetNbinsX();
0218     Double_t xmin = h_ihcale[0]->GetBinLowEdge(1);
0219     Double_t xmax = h_ihcale[0]->GetBinLowEdge(nbinsx+1);
0220     h2_ihcale_vs_run = new TH2F("h2_ihcale_vs_run","iHCAL E vs run",nruns,0,nruns,nbinsx,xmin,xmax);
0221   }
0222
0223
0224   for (int irun=0; irun<nruns; irun++)
0225   {
0226     int nbinsx = h_bbcqsum[irun]->GetNbinsX();
0227     for (int ibin=1; ibin<=nbinsx; ibin++)
0228     {
0229       Float_t val = h_bbcqsum[irun]->GetBinContent(ibin);
0230       Float_t qsum = h_bbcqsum[irun]->GetBinCenter(ibin);
0231       h2_bbcqsum_vs_run->Fill( irun, qsum, val );
0232
0233       val = h_bbcqtot[0][irun]->GetBinContent(ibin);
0234       qsum = h_bbcqtot[0][irun]->GetBinCenter(ibin);
0235       h2_bbcqs_vs_run->Fill( irun, qsum, val );
0236
0237       val = h_bbcqtot[1][irun]->GetBinContent(ibin);
0238       qsum = h_bbcqtot[1][irun]->GetBinCenter(ibin);
0239       h2_bbcqn_vs_run->Fill( irun, qsum, val );
0240     }
0241
0242     for (int ipmt=0; ipmt<128; ipmt++)
0243     {
0244       int nbinsx = h_bbcq[irun][ipmt]->GetNbinsX();
0245       for (int ibin=1; ibin<=nbinsx; ibin++)
0246       {
0247         Double_t val = h_bbcq[irun][ipmt]->GetBinContent(ibin);
0248         Double_t q = h_bbcq[irun][ipmt]->GetBinCenter(ibin);
0249         h2_bbcq[ipmt]->Fill( (Double_t)irun, q, val );
0250       }
0251     }
0252
0253     if ( do_zdc )
0254     {
0255       int nbinsx = h_zdce[irun]->GetNbinsX();
0256       for (int ibin=1; ibin<=nbinsx; ibin++)
0257       {
0258         Float_t val = h_zdce[irun]->GetBinContent(ibin);
0259         Float_t e = h_zdce[irun]->GetBinCenter(ibin);
0260         h2_zdce_vs_run->Fill( irun, e, val );
0261       }
0262     }
0263
0264     if ( do_emcal )
0265     {
0266       int nbinsx = h_emcale[irun]->GetNbinsX();
0267       for (int ibin=1; ibin<=nbinsx; ibin++)
0268       {
0269         Float_t val = h_emcale[irun]->GetBinContent(ibin);
0270         Float_t emce = h_emcale[irun]->GetBinCenter(ibin);
0271         h2_emcale_vs_run->Fill( irun, emce, val );
0272       }
0273     }
0274
0275     if ( do_ohcal )
0276     {
0277       int nbinsx = h_ohcale[irun]->GetNbinsX();
0278       for (int ibin=1; ibin<=nbinsx; ibin++)
0279       {
0280         Float_t val = h_ohcale[irun]->GetBinContent(ibin);
0281         Float_t e = h_ohcale[irun]->GetBinCenter(ibin);
0282         h2_ohcale_vs_run->Fill( irun, e, val );
0283       }
0284     }
0285
0286     if ( do_ihcal )
0287     {
0288       int nbinsx = h_ihcale[irun]->GetNbinsX();
0289       for (int ibin=1; ibin<=nbinsx; ibin++)
0290       {
0291         Float_t val = h_ihcale[irun]->GetBinContent(ibin);
0292         Float_t e = h_ihcale[irun]->GetBinCenter(ibin);
0293         h2_ihcale_vs_run->Fill( irun, e, val );
0294       }
0295     }
0296
0297
0298   }
0299
0300   TCanvas *ac[100];
0301   int icv = 0;
0302
0303   ac[icv] = new TCanvas("mbdq_vs_run2d","MBD Q vs Run",800,600);
0304   h2_bbcqsum_vs_run->SetMinimum(1e-5);
0305   h2_bbcqsum_vs_run->DrawCopy("colz");
0306   gPad->SetLogz(1);
0307
0308   ac[++icv] = new TCanvas("mbdqs_vs_run2d","MBD Q.S vs Run",800,600);
0309   h2_bbcqs_vs_run->SetMinimum(1e-6);
0310   h2_bbcqs_vs_run->DrawCopy("colz");
0311   gPad->SetLogz(1);
0312
0313   ac[++icv] = new TCanvas("mbdqn_vs_run2d","MBD Q.N vs Run",800,600);
0314   h2_bbcqn_vs_run->SetMinimum(1e-6);
0315   h2_bbcqn_vs_run->DrawCopy("colz");
0316   gPad->SetLogz(1);
0317
0318   ac[++icv] = new TCanvas("mbdq_vs_run","MBD Qsum vs Run",800,600);
0319   g_meansum = new TGraphErrors(nruns,runindex,bqmeansum,0,bqmeansumerr);
0320   name = "q_meansum";
0321   title = name;
0322   g_meansum->SetName( name );
0323   g_meansum->SetTitle( title );
0324   g_meansum->SetMarkerStyle(20);
0325   g_meansum->Draw("ap");
0326
0327   ac[++icv] = new TCanvas("mbdqs_vs_run","MBD Q.S vs Run",800,600);
0328   g_meantot[0] = new TGraphErrors(nruns,runindex,bqmeantot[0],0,bqmeantoterr[0]);
0329   name = "q_meantot0";
0330   title = name;
0331   g_meantot[0]->SetName( name );
0332   g_meantot[0]->SetTitle( title );
0333   g_meantot[0]->SetMarkerStyle(20);
0334   g_meantot[0]->Draw("ap");
0335
0336   ac[++icv] = new TCanvas("mbdqn_vs_run","MBD Q.N vs Run",800,600);
0337   g_meantot[1] = new TGraphErrors(nruns,runindex,bqmeantot[1],0,bqmeantoterr[1]);
0338   name = "q_meantot1";
0339   title = name;
0340   g_meantot[1]->SetName( name );
0341   g_meantot[1]->SetTitle( title );
0342   g_meantot[1]->SetMarkerStyle(20);
0343   g_meantot[1]->Draw("ap");
0344
0345   ac[++icv] = new TCanvas("mbdq_vs_run_ch","MBD Q vs Run",1200,800);
0346   for (int ipmt=0; ipmt<128; ipmt++)
0347   {
0348     h2_bbcq[ipmt]->Draw("colz");
0349     gPad->SetLogz(1);
0350
0351     name = dir; name += h2_bbcq[ipmt]->GetName(); name += ".png";
0352     //cout << name << endl;
0353     ac[icv]->Print( name );
0354   }
0355
0356   /*
0357   ac[++icv] = new TCanvas("mbdq_vs_run_ch_low","MBD Q vs Run, Low",1200,800);
0358   for (int ipmt=0; ipmt<128; ipmt++)
0359   {
0360     h2_bbcq[ipmt]->GetYaxis()->SetRangeUser(0,4);
0361     h2_bbcq[ipmt]->Draw("colz");
0362
0363     name = dir; name += h2_bbcq[ipmt]->GetName(); name += "_low.png";
0364     //cout << name << endl;
0365     ac[icv]->Print( name );
0366   }
0367   */
0368
0369   int run1 = 0;
0370   int run2 = 1;
0371   /*
0372      h_bbcqsum[run1]->Rebin(8);
0373      h_bbcqsum[run2]->Rebin(8);
0374      h_bbcqsum[run1]->Draw("hist");
0375      h_bbcqsum[run2]->SetLineColor(2);
0376      h_bbcqsum[run2]->Draw("histsame");
0377      gPad->SetLogy(1);
0378   */
0379
0380   if ( do_zdc )
0381   {
0382     ac[++icv] = new TCanvas("zdce_vs_run","ZDC E vs Run",800,600);
0383     h2_zdce_vs_run->SetMinimum(1e-6);
0384     h2_zdce_vs_run->Draw("colz");
0385     gPad->SetLogz(1);
0386   }
0387
0388   if ( do_emcal )
0389   {
0390     ac[++icv] = new TCanvas("emce_vs_run","EMC E vs Run",800,600);
0391     h2_emcale_vs_run->SetMinimum(1e-6);
0392     h2_emcale_vs_run->Draw("colz");
0393     gPad->SetLogz(1);
0394
0395     h_emcale[run1]->Rebin(4);
0396     h_emcale[run2]->Rebin(4);
0397     h_emcale[run2]->Scale(0.6);
0398     h_emcale[run1]->Draw("hist");
0399     h_emcale[run2]->SetLineColor(2);
0400     h_emcale[run2]->Draw("histsame");
0401     gPad->SetLogy(1);
0402   }
0403
0404   if ( do_ohcal )
0405   {
0406     ac[++icv] = new TCanvas("ohcal_vs_run","oHCAL E vs Run",800,600);
0407     h2_ohcale_vs_run->SetMinimum(1e-6);
0408     h2_ohcale_vs_run->Draw("colz");
0409     gPad->SetLogz(1);
0410
0411     h_ohcale[run1]->Rebin(4);
0412     h_ohcale[run2]->Rebin(4);
0413     h_ohcale[run2]->Scale(0.6);
0414     h_ohcale[run1]->Draw("hist");
0415     h_ohcale[run2]->SetLineColor(2);
0416     h_ohcale[run2]->Draw("histsame");
0417     gPad->SetLogy(1);
0418   }
0419
0420   if ( do_ihcal )
0421   {
0422     ac[++icv] = new TCanvas("ihcal_vs_run","iHCAL E vs Run",800,600);
0423     h2_ihcale_vs_run->SetMinimum(1e-6);
0424     h2_ihcale_vs_run->Draw("colz");
0425     gPad->SetLogz(1);
0426
0427     h_ihcale[run1]->Rebin(4);
0428     h_ihcale[run2]->Rebin(4);
0429     h_ihcale[run2]->Scale(0.6);
0430     h_ihcale[run1]->Draw("hist");
0431     h_ihcale[run2]->SetLineColor(2);
0432     h_ihcale[run2]->Draw("histsame");
0433     gPad->SetLogy(1);
0434   }
0435
0436
0437   // find gain corrections
0438   if ( calc_gain_corr )
0439   {
0440     //int refrun = 0;  // index of reference run
0441     flist.clear();
0442     flist.seekg(0);
0443     int irun = 0;
0444     ofstream gaincorr_calfile;
0445     while ( flist >> dstfname )
0446     {
0447       TString dir = "results/";
0448       dir += get_runstr(dstfname);
0449       dir += "/";
0450       name = "mkdir -p " + dir;
0451       gSystem->Exec( name );
0452       name = dir + "mbd_gaincorrmean.calib";
0453       gaincorr_calfile.open( name );
0454
0455       cout << name << endl;
0456
0457       for (int ipmt=0; ipmt<128; ipmt++)
0458       {
0459         double corr = bqmean[ipmt][irun] / runmean[ipmt];
0460         //double corr = bqmean[ipmt][irun] / bqmean[refun][ipmt];
0461         if ( fabs(corr-1.0)>0.01 )
0462         {
0463           cout << irun << "\t" << ipmt << "\t" << corr << endl;
0464         }
0465         gaincorr_calfile << ipmt << "\t" << corr << endl;
0466       }
0467
0468       gaincorr_calfile.close();
0469
0470       irun++;
0471     }
0472   }
0473
0474 }