MakeValidCalibPDF.C
Go to the documentation of this file.
1 
2 #include "CalibAnaPlot.h"
3 #include "CalibAnaTypes.h"
4 #include "TexBuilder.h"
5 
6 #include <iomanip>
7 #include <fstream>
8 #include <sstream>
9 
10 #include "TF1.h"
11 
12 const int wid1 = 17;
13 const int wid2 = 13;
14 const char sep = ' ';
15 
16 class ValidCalibPDF : public CalibAnaPlot{
17 
18 public:
19 
20  ValidCalibPDF( std::string outpath, std::string tex_cfg );
21  // vvv CalibAnaPlot implimentations
22  //void ScheduleAllDataPlots();
23  void SchedulePlots();
24  void ScheduleSamples();
25  void ScheduleHitVars();
26  void FillHitPlots();
27  // ^^^ CalibAnaPlot implimentations
28 
29  void MakeTex(std::string config_path);
30  void DiblockTex(std::string texdir, std::vector<std::string> configs);
31  void DriftTex(std::string texdir,
32  unsigned short bright_bin = std::numeric_limits<unsigned short>::max());
34  unsigned short bright_bin = std::numeric_limits<unsigned short>::max());
35  void CalibAllSamples();
36  void CalibSummaryTable();
37 
38  void SetApplyAbsCutToAll(bool set) { fApplyAbsCutToAll=set; }
39  void SetPerformAbsCal (bool set) { fPerformAbsCal=set; }
40  void SetShowRecoMeV (bool set) { fShowRecoMeV=set; }
41  void SetEpochList (std::string tex_cfg);
42 
44 
45  void GetCSVRow( float &meu, float &meuerr, float &gev,
46  std::string calibtag,
47  std::string det, std::string dataormc,
48  std::string runrange);
49 
50 private:
51  void InitRunsMap();
52  std::map< std::string, std::string > fRunsMap;
53  std::string RunRange(std::string sample, std::string oldcalibtag="");
54 
55  std::vector< std::string > fEpochList;
57 
58  std::map<std::string, AbsCal > fAbsCals;
59  void AbsoluteCalibration( std::string sample );
63  bool PassesAbsCalCut();
64 
65  std::vector< std::string > fMCSamples;
66 
67  std::string ToString(float num,unsigned short digits=0);
69  std::string ToString(unsigned int num);
70 
72 
73 };
74 
75 
76 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
78  : CalibAnaPlot(outpath)
79  , fTexConfig(tex_cfg)
80  , fPerformAbsCal(false)
81  , fApplyAbsCutToAll(false)
82  , fShowRecoMeV(false)
83 {
84  InitRunsMap();
85  SetEpochList(tex_cfg);
87 }
88 
89 
90 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
92  std::string tex_config,
93  std::string sample,
94  bool justdraw=false,
95  unsigned int stride=0,
96  unsigned int limit=0)
97 {
98  ValidCalibPDF bigpdf( outpath, tex_config );
99  bigpdf.Initialize();
100  if(stride) bigpdf.SetTreeLoopStride(stride);
101  if(limit) bigpdf.SetTreeLoopLimit(limit);
102 
103  // Unique tag for tex is first line, and a handy way to set cuts
104  std::ifstream cfg(tex_config);
105  if (!cfg.is_open()){ std::cout << "could not find config file " << tex_config << std::endl; return; }
106  std::string tex_tag;
107  std::getline (cfg,tex_tag);
108  bigpdf.SetApplyAbsCutToAll( tex_tag.find("abscal")!=std::string::npos );
109  bigpdf.SetShowRecoMeV( tex_tag.find("withreco")!=std::string::npos );
110 
111  // Enforce that Abs Cal be done all in the same pass
112  // to avoid unknown skews between histogram files
113  //if(sample=="all"&&!justdraw) bigpdf.SetPerformAbsCal(true);
114  if(sample=="all") bigpdf.SetPerformAbsCal(true);
115 
116  if(justdraw){
117  // If histogram files are already made, just make plots and tex...
118  bigpdf.MakeTex(tex_config);
119  } else {
120  // ...Otherwise, fill histograms
121 
122  if(sample=="all"){
123  for( auto& s : bigpdf.GetSamples() ){
124  bigpdf.ProcessTrees(s.first);
125  }
126  bigpdf.MakeTex(tex_config);
127  // may as well make the tex
128  } else {
129  bigpdf.ProcessTrees(sample);
130  // may not have all samples to make tex, so don't
131  }
132 
133  }
134 
135 }
136 
137 
138 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
140  // Figure out which data samples need to be processed from sweep through config file
141  std::ifstream cfg(tex_cfg);
142  if (!cfg.is_open()){ std::cout << "could not find config file " << tex_cfg << std::endl; return; }
144  std::vector< std::string > possible_epochs = {
145  "1a", "2a","2b","2c","2d","2e", "3a","3b","3c","3d","3e",
146  "4a", "4b.1", "4b.2", "4b", "4c", // split 4b up manually in ND
147  "5a","5b", "6a","6b", "7a","7b","7c","7d","7e",
148  "8a","8b", "9a","9b" };
149  std::string all_epochs;
150  while ( std::getline (cfg,line) ){
151  if(line.front()=='#') continue;
152  for(auto& ep : possible_epochs){
153  if( line.find("_data_ep"+ep)!=std::string::npos &&
154  all_epochs.find("_"+ep)==std::string::npos ){
155  fEpochList.push_back(ep);
156  std::cout << "Scheduling Epoch " << ep << std::endl;
157  all_epochs += "_"+ep;
158  }
159  }
160  }
161 }
162 
163 
164 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
166 { // table straight out for catching mistakes early
167  std::cout << std::left << std::setw(wid2) << std::setfill(sep) << "Meu(X)";
168  std::cout << std::left << std::setw(wid2) << std::setfill(sep) << "Meu(Y)";
169  std::cout << std::left << std::setw(wid2) << std::setfill(sep) << "Meu(PEc/cm)";
170  std::cout << std::left << std::setw(wid2) << std::setfill(sep) << "Meu_Err";
171  std::cout << std::left << std::setw(wid2) << std::setfill(sep) << "dEdx(X)";
172  std::cout << std::left << std::setw(wid2) << std::setfill(sep) << "dEdx(Y)";
173  std::cout << std::left << std::setw(wid2) << std::setfill(sep) << "dEdx(MeV/cm)";
174  std::cout << std::left << std::setw(wid2) << std::setfill(sep) << "dEdx_Err";
175  std::cout << std::left << std::setw(wid2+1) << std::setfill(sep) << "Scale(MeV/PEc)";
176  //std::cout << std::left << std::setw(wid2) << std::setfill(sep) << "Scale Err";
177  std::cout << std::left << std::setw(wid1) << std::setfill(sep) << "#_Reco_x";
178  std::cout << std::left << std::setw(wid1) << std::setfill(sep) << "#_Reco_y";
179  std::cout << std::endl;
180 
181  // Keep officially formated csvs along the way
182  MakeDir(fOutPath+"/calib_csvs");
183 
184  // Do absolute calibration, need MC first for the trueE/dx
185  for( auto& mcs : fMCSamples )
186  AbsoluteCalibration(mcs);
187  for( auto& s : fSamples )
188  if(s.first.find("mc")==std::string::npos)
189  AbsoluteCalibration(s.first);
190 }
191 
192 
193 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
195 {
196  std::cout << "Do Absolute Calibration for " << sample << std::endl;
197  bool isMC = sample.find("mc")!=std::string::npos;
198 
199  if(!isMC){
200  bool allmc = true;
201  for(auto& mcs : fMCSamples)
202  if( fAbsCals.count(mcs)!=1 ) allmc = false;
203 
204  if(!allmc){
205  std::cout << "Calibrate All MC First." << std::endl;
206  std::abort();
207  }
208  }
209 
210  AbsCal ac;
211  ac.sample = sample;
212 
213  TH1F* nhits_meu_x;
214  TH1F* nhits_meu_y;
215  TH1F* nhits_mev_x;
216  TH1F* nhits_mev_y;
217  TFile* sample_file = new TFile((fOutPath+"/"+sample+".root").c_str());
218  TFile* mc_file = new TFile((fOutPath+"/"+MCSample(sample)+".root").c_str());
219  std::cout << "printing useful info: " << fOutPath+"/"+sample+".root" << " BREAK " << fOutPath+"/"+MCSample(sample)+".root" << std::endl;
220  if(!sample_file) std::cout << "Did not find " << fOutPath+"/"+sample+".root" << std::endl;
221  if(!mc_file) std::cout << "Did not find " << fOutPath+"/"+MCSample(sample)+".root" << std::endl;
222  sample_file ->GetObject("nhits_meu_x",nhits_meu_x);
223  sample_file ->GetObject("nhits_meu_y",nhits_meu_y);
224  mc_file ->GetObject("nhits_mev_x",nhits_mev_x);
225  mc_file ->GetObject("nhits_mev_y",nhits_mev_y);
226  if(!nhits_meu_x) std::cout << "Did not find nhits_meu_x" << std::endl;
227  if(!nhits_meu_y) std::cout << "Did not find nhits_meu_y" << std::endl;
228  if(!nhits_mev_x) std::cout << "Did not find nhits_mev_x" << std::endl;
229  if(!nhits_mev_y) std::cout << "Did not find nhits_mev_y" << std::endl;
230 
231  nhits_meu_x->GetXaxis()->SetRangeUser(0,100);
232  nhits_meu_y->GetXaxis()->SetRangeUser(0,100);
233  ac.meu_x = nhits_meu_x->GetMean(1);
234  ac.meu_y = nhits_meu_y->GetMean(1);
235  double meuRMS_x = nhits_meu_x->GetRMS(1);
236  double meuRMS_y = nhits_meu_y->GetRMS(1);
237  ac.meuErr_x = meuRMS_x/std::sqrt(nhits_meu_x->GetEntries());
238  ac.meuErr_y = meuRMS_y/std::sqrt(nhits_meu_y->GetEntries());
239 
240 
241  nhits_mev_x->GetXaxis()->SetRangeUser(0,6);
242  nhits_mev_y->GetXaxis()->SetRangeUser(0,6);
243  ac.mev_x = nhits_mev_x->GetMean(1);
244  ac.mev_y = nhits_mev_y->GetMean(1);
245  double mevRMS_x = nhits_mev_x->GetRMS(1);
246  double mevRMS_y = nhits_mev_y->GetRMS(1);
247  ac.mevErr_x = mevRMS_x/std::sqrt(nhits_mev_x->GetEntries());
248  ac.mevErr_y = mevRMS_y/std::sqrt(nhits_mev_y->GetEntries());
249 
250 
251  ac.meu = (ac.meu_x + ac.meu_y)/2;
252  ac.meuErr = 1 / std::sqrt(std::pow(ac.meuErr_x,-2) + std::pow(ac.meuErr_y,-2));
253  ac.mev = (ac.mev_x + ac.mev_y)/2;
254  ac.mevErr = 1 / std::sqrt(std::pow(ac.mevErr_x,-2) + std::pow(ac.mevErr_y,-2));
255 
256  ac.numHits_x = nhits_meu_x->GetEntries();
257  ac.numHits_y = nhits_meu_y->GetEntries();
258 
259  // Done calculateing the calibration numbers
260  fAbsCals[sample] = ac;
261 
262  std::stringstream returnString;
263  returnString << std::left << std::setw(wid2) << std::setfill(sep) << ac.meu_x;
264  returnString << std::left << std::setw(wid2) << std::setfill(sep) << ac.meu_y;
265  returnString << std::left << std::setw(wid2) << std::setfill(sep) << ac.meu;
266  returnString << std::left << std::setw(wid2) << std::setfill(sep) << ac.meuErr;
267  returnString << std::left << std::setw(wid2) << std::setfill(sep) << ac.mev_x;
268  returnString << std::left << std::setw(wid2) << std::setfill(sep) << ac.mev_y;
269  returnString << std::left << std::setw(wid2) << std::setfill(sep) << ac.mev;
270  returnString << std::left << std::setw(wid2) << std::setfill(sep) << ac.mevErr;
271  returnString << std::left << std::setw(wid2+1) << std::setfill(sep) << 100*ac.mev/ac.meu;
272  //returnString << std::left << std::setw(wid2) << std::setfill(sep) << 0; // TODO: calculate error on En. Scale!
273  returnString << std::left << std::setw(wid1) << std::setfill(sep) << ac.numHits_x;
274  returnString << std::left << std::setw(wid1) << std::setfill(sep) << ac.numHits_y;
275  returnString << std::endl;
276  std::cout << returnString.str() << std::endl;
277 
278  // other CSV values regardless of method
279  double gev = 0.001*ac.mev;
280  int tv = 1;
281  std::string newtag = "v14";
282 
283 
284  Table tab;
285  tab.name = ac.sample+"_abscal_csv";
286  tab.csvdir = fOutPath+"/calib_csvs";
287  // Calibrator format:
288  tab.csvname = Form("calib_abs_consts.%s.%s.%s.%s.csv",
289  fDet->name().c_str(),
290  (isMC ? "mc" : "data"),
291  newtag.c_str(),
292  RunRange(sample).c_str());
293 
294  tab.title = tab.csvname;
295  tab.header = {"meu","meu_err","gev","channel","tv"};
296  for(unsigned int channel=1; channel<=fDet->NDiblocks(); channel++){
297  std::vector< std::string > row;
298 
299  row.push_back( std::to_string( ac.meu ) );
300  row.push_back( std::to_string( ac.meuErr ) );
301  row.push_back( std::to_string( gev ) );
302  row.push_back( std::to_string( channel ) );
303  row.push_back( std::to_string( tv ) );
304  // or format with Form() like before?
305  //csv << Form("%f,%f,%f,%i,%i", ac.meu, ac.meuErr, gev, channel, tv) << std::endl;
306 
307  tab.rows.push_back(row);
308  }
309  fTex->AddTable(tab, false,!(sample.find("2016")!=std::string::npos));
310 
311 }
312 
313 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
314 std::string ValidCalibPDF::ToString(float num,unsigned short digits)
315 {
316  std::stringstream ret;
317  if(digits!=0) ret << std::setprecision(digits) << num;
318  else ret << num;
319  return ret.str();
320 }
322 {
323  std::stringstream ret; ret << num; return ret.str();
324 }
326 {
327  std::stringstream ret; ret << num; return ret.str();
328 }
329 
330 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
332 {
333  Table tab;
334  tab.name = "calib_summary_table";
335  tab.csvdir = fOutPath+"/calib_csvs";
336  tab.csvname = "calib_summary_table.csv";
337  tab.title = "Abs. Cal. Summary Table";
338  tab.header = {"Sample",
339  "NHits$_x$","MEU$_x$",
340  "NHits$_y$","MEU$_y$",
341  "MEU","MEU Err",
342  "TrueE/dx","tE/dx Err"};
343 
344  unsigned short dig = 4;
345  for( auto& ac: fAbsCals ){
346  std::vector< std::string > row;
347 
348  row.push_back( ac.second.sample );
349  row.push_back( ToString( ac.second.numHits_x, dig) );
350  row.push_back( ToString( ac.second.meu_x, dig) );
351  row.push_back( ToString( ac.second.numHits_y, dig) );
352  row.push_back( ToString( ac.second.meu_y, dig) );
353  row.push_back( ToString( ac.second.meu, dig) );
354  row.push_back( ToString( ac.second.meuErr, dig) );
355  row.push_back( ToString( ac.second.mev, dig) );
356  row.push_back( ToString( ac.second.mevErr, dig) );
357 
358  tab.rows.push_back(row);
359  }
360  fTex->AddTable(tab, true);
361 
362  // New table
363  /*
364  tab.name = "version_compare_table";
365  tab.csvdir = fOutPath;
366  tab.csvname = "version_compare_table";
367  tab.title = "Compare to UPS Versions";
368  tab.header = {"Sample",
369  "MEU","MEU Err",
370  "TrueE/dx","tE/dx Err"};
371  std::vector< std::string > unitsrow = {
372  "","PEc/cm","\\%","PEc/cm","\\%","MeV/cm","\\%"};
373  tab.rows.push_back(units_row);
374 
375 
376  std::string calibtag = "v13";
377  for( auto& ac: fAbsCals ){
378 
379  std::string dataormc = (ac.sample.find("mc")!=std::string::npos) ? "mc" : "data";
380 
381  float meu, meuerr, gev;
382  GetCSVRow( meu,meuerr,gev,
383  calibtag,
384  fDet->name(),
385  dataormc,
386  RunRange(sample,calibtag) );
387  }
388  */
389 
390 }
391 
392 
393 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
394 void ValidCalibPDF::GetCSVRow( float &meu, float &meuerr, float &gev,
395  std::string calibtag,
396  std::string det, std::string dataormc,
397  std::string runrange)
398 {
399 
400  std::string csvname = Form("calib_abs_consts.%s.%s.%s.%s.csv",
401  det.c_str(), dataormc.c_str(),
402  calibtag.c_str(),
403  runrange.c_str() );
404  ifstream csvfile( csvname );
405 
407  unsigned short n=0;
408  while( getline( csvfile, line ) ){
409  if( n==0 ) continue;
410  else if( n>1 ) break;
411  n++;
412 
413  std::stringstream ss( line );
414 
415  unsigned short cell=0;
416  while(ss.good()){
418  if (!getline( ss, num, ',' )) break;
419  if( cell==0 ) meu = std::stof(num);
420  else if( cell==1 ) meuerr = std::stof(num);
421  else if( cell==2 ) gev = std::stof(num);
422  else return;
423  }
424 
425  }
426 
427  std::cout << "Warning, did not find line in CSV " << csvname << std::endl;
428 }
429 
430 
431 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
433 {
434  std::string samplekey;
435 
436  if( oldcalibtag=="" ){ samplekey=sample; }
437  else if ( oldcalibtag=="v13" ){
438  // Map new calib samples to old calibrator version
439  if(sample.find("nd_mc")!=std::string::npos) samplekey = oldcalibtag+"_nd_mc";
440  if(sample.find("fd_mc_lowgain")!=std::string::npos) samplekey = oldcalibtag+"_fd_mc_lowgain";
441  if(sample.find("fd_mc_highgain")!=std::string::npos) samplekey = oldcalibtag+"_fd_mc_highgain";
442  if(sample.find("nd_data_ep2")!=std::string::npos) samplekey = oldcalibtag+"_nd_data_p2";
443  if(sample.find("nd_data_ep3")!=std::string::npos) samplekey = oldcalibtag+"_nd_data_p3";
444  if(sample.find("nd_data_ep4")!=std::string::npos) samplekey = oldcalibtag+"_nd_data_p4";
445  if(sample.find("nd_data_ep5")!=std::string::npos) samplekey = oldcalibtag+"_nd_data_p56";
446  if(sample.find("nd_data_ep6")!=std::string::npos) samplekey = oldcalibtag+"_nd_data_p56";
447  if(sample.find("nd_data_ep7")!=std::string::npos) samplekey = oldcalibtag+"_nd_data_p7";
448  if(sample.find("nd_data_ep8")!=std::string::npos) samplekey = oldcalibtag+"_nd_data_p8";
449  // TODO: Add FD Data
450  } else {
451  std::cout << "Calibrator tag " << oldcalibtag << " has not been registered in fRunsMap, abort" << std::endl;
452  std::abort();
453  }
454 
455  if( fRunsMap.count(samplekey) != 1 ){
456  std::cout << "Sample " << samplekey << " not registered in fRunsMap, abort" << std::endl;
457  std::abort();
458  }
459 
460  return fRunsMap[samplekey];
461 }
462 
463 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
465 {
466  // Encoding the authoritative period/epoch naming page:
467  // https://cdcvs.fnal.gov/redmine/projects/novaart/wiki/Period_and_Epoch_Naming#Period-and-Epoch-Naming
468 
469  fRunsMap["tb_mc"] = "r-r";
470  fRunsMap["tb_data_p1"] = "r-r";
471  fRunsMap["nd_mc"] = "r-r";
472  fRunsMap["fd_mc_lowgain"] = "r-r20752";
473  fRunsMap["fd_mc_p2"] = "r-r20752";
474  fRunsMap["fd_mc_p2_2016"] = "r-r20752";
475  fRunsMap["fd_mc_highgain"] = "r20753-r";
476  fRunsMap["fd_mc_p3"] = "r20753-r";
477  fRunsMap["fd_mc_p3_2016"] = "r20753-r";
478  fRunsMap["nd_data_ep2a"] = "r10496-r10824";
479  fRunsMap["nd_data_ep2b"] = "r10825-r10932";
480  fRunsMap["nd_data_ep2c"] = "r10933-r10973";
481  fRunsMap["nd_data_ep2d"] = "r10974-r11018";
482  fRunsMap["nd_data_ep2e"] = "r11019-r11228";
483  fRunsMap["nd_data_ep3b"] = "r11229-r11301"; // 3b calibrating 3a+3b
484  fRunsMap["nd_data_ep3c"] = "r11302-r11383";
485  fRunsMap["nd_data_ep3d"] = "r11384-r11537";
486  fRunsMap["nd_data_ep3e"] = "r11538-r11631";
487  fRunsMap["nd_data_ep4a"] = "r11632-r11667";
488  fRunsMap["nd_data_ep4b.1"] = "r11668-r11725";
489  fRunsMap["nd_data_ep4b.2"] = "r11726-r11920";
490  fRunsMap["nd_data_ep4c"] = "r11921-r11925";
491  fRunsMap["nd_data_ep5a"] = "r11926-r12028";
492  fRunsMap["nd_data_ep5b"] = "r12029-r12086";
493  fRunsMap["nd_data_ep6a"] = "r12087-r12289";
494  fRunsMap["nd_data_ep6b"] = "r12290-r12516";
495  fRunsMap["nd_data_ep7a"] = "r12517-r12579";
496  fRunsMap["nd_data_ep7b"] = "r12580-r12667";
497  fRunsMap["nd_data_ep7c"] = "r12668-r12690";
498  fRunsMap["nd_data_ep7d"] = "r12691-r12955"; // 7d calibrating 7d+7e
499  fRunsMap["nd_data_ep8b"] = "r12956-r"; // 8b calibrating 8a+8b+9a-on
500 
501  // vor v13 csvs:
502  fRunsMap["v13_nd_mc"] = "r-r";
503  fRunsMap["v13_fd_mc_lowgain"] = "r-r20752";
504  fRunsMap["v13_fd_mc_highgain"] = "r20753-r";
505  fRunsMap["v13_nd_data_p2"] = "r-r11228";
506  fRunsMap["v13_nd_data_p3"] = "r11229-r11628";
507  fRunsMap["v13_nd_data_p4"] = "r11632-r11925";
508  fRunsMap["v13_nd_data_p56"] = "r11926-r12516";
509  fRunsMap["v13_nd_data_p7"] = "r12517-r12955";
510  fRunsMap["v13_nd_data_p8"] = "r12956-r";
511 
512 
513  ///// FD Data
514 
515  fRunsMap["fd_data_ep1a"] = "r-r17139";
516  fRunsMap["fd_data_ep2a"] = "r17891-r19096";
517  fRunsMap["fd_data_ep2b"] = "r19097-r19586";
518  //fRunsMap["fd_data_ep2c"] = "r19587-r19746";
519  //fRunsMap["fd_data_ep2d"] = "r19747-r19958"; // No-Horn-Current horn scans
520  //fRunsMap["fd_data_ep2e"] = "r19959-r20752"; // 2015 Summer shutdown
521  fRunsMap["fd_data_ep2c"] = "r19587-r20752"; // 2c to calibrating 2c+2d+2e
522 
523  // 3a: r20753-r20922 // no beam
524  // 3b: r20923-r21230 // ~3 weeks FHC
525  fRunsMap["fd_data_ep3b"] = "r20753-r21230"; // 3b calibrating 3a+3b
526  fRunsMap["fd_data_ep3c"] = "r21231-r22019"; // 2 months FHC
527  fRunsMap["fd_data_ep3d"] = "r22020-r22900"; // ~2.5 months FHC
528  fRunsMap["fd_data_ep3e"] = "r22901-r23419";
529 
530  // 4a: r23420-r23670 // one month RHC
531  //fRunsMap["fd_data_ep4b"] = "r23671-r24587"; // summer shutdown
532  fRunsMap["fd_data_ep4a"] = "r23420-r24587"; // 4a calibrating 4a+4b
533  //fRunsMap["fd_data_ep4c"] = "r24588-r24613"; // a few days of RHC
534 
535  // 5a: r24614-r25035 // 2 months FHC
536  fRunsMap["fd_data_ep5a"] = "r24588-r25035"; // 5a calibrating 4c+5a (same beam run, horn current swapped)
537  fRunsMap["fd_data_ep5b"] = "r25036-r25412"; // 5 weeks FHC
538 
539  //fRunsMap["fd_data_ep6a"] = "r25413-r26685"; // ~4.5 months RHC
540  //fRunsMap["fd_data_ep6b"] = "r26686-r28036"; // summer shutdown
541  fRunsMap["fd_data_ep6a"] = "r25413-r28036"; // 6a calibrating 6a+6b
542 
543  fRunsMap["fd_data_ep7a"] = "r28037-r28549";
544  fRunsMap["fd_data_ep7b"] = "r28550-r29225";
545  fRunsMap["fd_data_ep7c"] = "r29226-r29427";
546  //fRunsMap["fd_data_ep7d"] = "r29428-r30226";
547  //fRunsMap["fd_data_ep7e"] = "r30227-r31167";
548  fRunsMap["fd_data_ep7d"] = "r29428-r31167"; // 7d calibrating 7d+7e
549 
550  fRunsMap["fd_data_ep8a"] = "r31168-r31249";
551  fRunsMap["fd_data_ep8b"] = "r31250-r32397";
552  //fRunsMap["fd_data_ep9a"] = "r32398-r33541";
553  // 9b: r33542-r... // beam shutdown
554  fRunsMap["fd_data_ep9a"] = "r32398-r"; // 9a calibrating 9a+9b...
555 }
556 
557 
558 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
560 {
561  SchedulePlots(); // know the full list of plots, probably cleared after filling
562 
563  std::ifstream cfg(tex_config);
564  if (!cfg.is_open()){ std::cout << "could not find config file " << tex_config << std::endl; return; }
565 
566  // First 4 lines of config make the title page
567  std::string tex_tag, title, authors, abstract;
568  std::getline (cfg,tex_tag);
569  std::getline (cfg,title);
570  std::getline (cfg,authors);
571  std::getline (cfg,abstract);
572  tex_tag = fOutPath+"/"+tex_tag+".tex"; // tex_tag should be tag that is also in fOutPath
573  std::cout << "Title: " << title
574  << "\nAuthor: " << authors
575  << "\nAbstract: " << abstract
576  << "\nTex Path: " << tex_tag << std::endl;
577  fTex = new TexBuilder( tex_tag, title, authors, abstract );
578 
579  // Do the Absolute Calibration
580  // save the Tables to add to the Tex at any point
581  CalibAllSamples(); // Makes official csv
582 
583  // Summary table of all calibrated samples
584  fTex->NewPage();
585  fTex->Section("Absolute Calibration");
586  fTex->AddText("Standard absolute calibration cuts: track window, flat-response W, positive pe, pecorr, and pathlenght reco");
588  fTex->NewPage();
589 
590  // Each following line configures a section.
591  // Delineated with commas, so avoid using commas in section_title
592  // section_label,setion_title,ratio_label,denom_sample,sample1,sample2,sample3... etc
594  while ( std::getline (cfg,line) ){
595  if(line=="") continue;
596  if(line.front()=='#') continue;
597 
598  std::cout << "cfg line: " << line << std::endl;
599 
600  // Break line down around double-commas into vector
601  std::stringstream ss(line);
602  vector<string> configs;
603  while( ss.good() ){ std::string c; std::getline( ss, c, ',' );
604  configs.push_back( c );
605  }
606  if(configs.size()<5){
607  std::cout << "Need at least 5 entries in section config line: \n"
608  << " section_type,setion_title,ratio_label,denom_sample,sample1 " << std::endl;
609  std::abort();
610  }
611 
612  // Hard code intended order of configuration line
613  std::string section_label = configs[0];
614  // section_label: picks the Tex function to apply and
615  // sets the directory to keep the plots
616  fTex->Section(configs[1]);
617  std::string ratio_label = configs[2];
618  if(ratio_label=="noratio") ratio_label = ""; // tell CalibAnaPlot::Draw to not draw ratio
619 
620 
621  if(section_label.find("dbsec")!=std::string::npos){
622  DiblockTex(section_label,configs);
623  continue; // nonstandard section shich doesn't compare epoch samples
624  }
625 
626  for( auto& p : fPlots ){
627  // Just make every plot for reference, though only
628  // a small subset will likely go in a given section.
629 
630  Plot denom_plot( fOutPath+"/"+configs[3]+".root",
631  p.first,
632  LegendLabel(configs[3]) );
633  std::vector< Plot > plotstack;
634  for( size_t cfg_i=4; cfg_i<configs.size(); ++cfg_i )
635  plotstack.push_back( Plot( fOutPath+"/"+configs[cfg_i]+".root",
636  p.first,
637  LegendLabel(configs[cfg_i]) ) );
638 
639  Draw(section_label, plotstack, denom_plot, ratio_label );
640  } // all fPlots
641 
642  if(section_label.find("driftsec")!=std::string::npos) DriftTex(section_label);
643  if(section_label.find("essentialsec")!=std::string::npos) EssentialComparisonTex(section_label);
644  }
645 
646  // Official CSV tables
647  fTex->Section("CSV Tables: \\$CALIBCSVS\\_DIR/csv/");
648  fTex->AddText("Note that absolute calibration cuts have been applied to make these tables.");
649  unsigned short i(0);
650  for( auto& samp : fSamples ){
651  if(i%2==0 && i!=0) fTex->NewPage();
652  std::string samplename = samp.first;
653  std::replace(samplename.begin(), samplename.end(), '_', ' ');
654  fTex->SubSection(samplename);
655  fTex->DrawTable(samp.first+"_abscal_csv");
656  i++;
657  }
658 
659  fTex->Close();
660 }
661 
662 
663 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
664 void ValidCalibPDF::DiblockTex(std::string texdir, std::vector<std::string> cfgs)
665 {
666  // config line goes like:
667  // dbsec_*,setion_title,ratio_label,...diblock list...
668 
669  // Make sure we have an up-to-date all_data file
670  HaddAllDataFile();
671 
672  // only fit high-gain (threshold effects in low-gain aren't great to fit against)
673  std::map< std::string, uint32_t > epoch_end_map = fDet->EndOfEpochMap();
674  uint32_t fit_start = epoch_end_map["2e"];
675  uint32_t fit_end = epoch_end_map["9a"];
676  Table tab_pecm_xy;
677 
678  // We want each diblock for all time to be it's own entry in the plot,
679  // so we will use the <det>_all_data.root sample made under the hood
680  // and instead iterate the plotname for each plot to be stacked
681  std::vector< std::string > plotnames = {
682  "pecm_time", "pecorrcm_time", "pe_time", "pecorr_time", "cm_time"};
683  if(fShowRecoMeV) plotnames.push_back("recomevcm_time");
684 
685  std::vector<std::string> views = {"x", "y"};
686  for(auto& pl : plotnames){
687  std::vector< std::vector< TF1* > > fitstack;
688  fitstack.resize(2);
689  for(auto& vstr : views){
690  std::string pname = pl + "_" + vstr;
691  size_t v = (vstr=="x") ? 0 : 1;
692 
693  std::vector< Plot > plotstack;
694  for( size_t i=3; i<cfgs.size(); ++i ){
695  std::string plotname = pname+"_db"+cfgs[i];
696  plotstack.push_back( Plot( fOutPath+"/"+fDet->name()+"_all_data.root",
697  plotname,
698  "Diblock "+cfgs[i] ) );
699  fitstack[v].push_back(new TF1((plotname).c_str(),
700  "pol1",
701  fit_start,fit_end));
702  plotstack.back().h->Fit(fitstack[v].back(),"RC0Q"); // use range in function
703  plotstack.back().h->Rebin(10);
704 
705  }
706  // Make denominator the average for all diblocks by not appending _db to plotname
707  Plot denom_plot( fOutPath+"/"+fDet->name()+"_all_data.root",
708  pname,
709  "All Diblocks" );
710  denom_plot.h->Rebin(10);
711 
712  // Make the plot
713  Draw(cfgs[0], plotstack, denom_plot, cfgs[2] );
714  }
715 
716  // Make a table of drift slopes
717  // TODO: Add goodness of fit. Looks like it isn't always linear...
718  if ( pl=="pecm_time" ){
719  tab_pecm_xy.name = "diblock_drift";
720  tab_pecm_xy.title = "Drift by Diblock";
721  tab_pecm_xy.header = {"Diblock","PE/cm Slope (X)","PE/cm Slope (Y)"};
722  for(size_t j=0; j<fitstack.size(); ++j){
723  std::vector< std::string > row;
724  row.push_back( cfgs[j+3] );
725  row.push_back( std::to_string( fitstack[0][j]->GetParameter(1) ) );
726  row.push_back( std::to_string( fitstack[1][j]->GetParameter(1) ) );
727  tab_pecm_xy.rows.push_back(row);
728  }
729  }
730 
731  }
732 
733 
734  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
735  fTex->AddFigureRow({texdir+"/pecm_time_x", texdir+"/pecm_time_y"}, {"", ""}, {"label1", "label2"});
736  // fTex->AddTable(tab_pecm_xy, true, false);
737  fTex->AddFigureRow({texdir+"/pecorrcm_time_x", texdir+"/pecorrcm_time_y"}, {"", ""}, {"label1", "label2"});
738  fTex->NewPage();
739  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
740  fTex->AddFigureRow({texdir+"/pe_time_x", texdir+"/pe_time_y"}, {"", ""}, {"label1", "label2"});
741  fTex->AddFigureRow({texdir+"/pecorr_time_x", texdir+"/pecorr_time_y"}, {"", ""}, {"label1", "label2"});
742  fTex->NewPage();
743  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
744  fTex->AddFigureRow({texdir+"/cm_time_x", texdir+"/cm_time_y"}, {"", ""}, {"label1", "label2"});
745  if(fShowRecoMeV)
746  fTex->AddFigureRow({texdir+"/recomevcm_time_x", texdir+"/recomevcm_time_y"}, {"", ""}, {"label1", "label2"});
747  fTex->NewPage();
748 }
749 
750 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
752  unsigned short bright_bin)
753 {
754  std::string bb = "";
756  std::to_string(bright_bin);
757 
758  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
759  fTex->AddFigureRow({texdir+"/pecm_time_x"+bb, texdir+"/pecm_time_y"+bb}, {"", ""}, {"label1", "label2"});
760  fTex->AddFigureRow({texdir+"/pecorrcm_time_x"+bb, texdir+"/pecorrcm_time_y"+bb}, {"", ""}, {"label1", "label2"});
761  if(fShowRecoMeV)
762  fTex->AddFigureRow({texdir+"/recomevcm_time_x"+bb, texdir+"/recomevcm_time_y"+bb}, {"", ""}, {"label1", "label2"});
763  fTex->NewPage();
764  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
765  fTex->AddFigureRow({texdir+"/pe_time_x"+bb, texdir+"/pe_time_y"+bb }, {"", ""}, {"label1", "label2"});
766  fTex->AddFigureRow({texdir+"/pecorr_time_x"+bb, texdir+"/pecorr_time_y"+bb}, {"", ""}, {"label1", "label2"});
767  fTex->AddFigureRow({texdir+"/cm_time_x"+bb, texdir+"/cm_time_y"+bb }, {"", ""}, {"label1", "label2"});
768  fTex->NewPage();
769 }
770 
771 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
773  unsigned short bright_bin)
774 {
775  std::string bb = "";
777  std::to_string(bright_bin);
778 
779  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
780  fTex->AddFigureRow({texdir+"/pecm_w_x"+bb, texdir+"/pecm_w_y"+bb}, {"", ""}, {"label1", "label2"});
781  fTex->AddFigureRow({texdir+"/pecorrcm_w_x"+bb, texdir+"/pecorrcm_w_y"+bb}, {"", ""}, {"label1", "label2"});
782  if(fShowRecoMeV)
783  fTex->AddFigureRow({texdir+"/recomevcm_w_x"+bb, texdir+"/recomevcm_w_y"+bb}, {"", ""}, {"label1", "label2"});
784 
785  fTex->NewPage();
786  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
787  fTex->AddFigureRow({texdir+"/pe_w_x"+bb, texdir+"/pe_w_y"+bb }, {"", ""}, {"label1", "label2"});
788  fTex->AddFigureRow({texdir+"/pecorr_w_x"+bb, texdir+"/pecorr_w_y"+bb}, {"", ""}, {"label1", "label2"});
789  fTex->AddFigureRow({texdir+"/cm_w_x"+bb, texdir+"/cm_w_y"+bb }, {"", ""}, {"label1", "label2"});
790 
791  fTex->NewPage();
792 
793  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
794  fTex->AddFigureRow({texdir+"/pecm_cell_x"+bb, texdir+"/pecm_cell_y"+bb}, {"", ""}, {"label1", "label2"});
795  fTex->AddFigureRow({texdir+"/pecorrcm_cell_x"+bb, texdir+"/pecorrcm_cell_y"+bb}, {"", ""}, {"label1", "label2"});
796 
797  fTex->NewPage();
798  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
799  fTex->AddFigureRow({texdir+"/pe_cell_x"+bb, texdir+"/pe_cell_y"+bb }, {"", ""}, {"label1", "label2"});
800  fTex->AddFigureRow({texdir+"/pecorr_cell_x"+bb, texdir+"/pecorr_cell_y"+bb}, {"", ""}, {"label1", "label2"});
801  fTex->AddFigureRow({texdir+"/cm_cell_x"+bb, texdir+"/cm_cell_y"+bb }, {"", ""}, {"label1", "label2"});
802 
803  fTex->NewPage();
804 
805  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
806  fTex->AddFigureRow({texdir+"/pecm_plane_x"+bb, texdir+"/pecm_plane_y"+bb}, {"", ""}, {"label1", "label2"});
807  fTex->AddFigureRow({texdir+"/pecorrcm_plane_x"+bb, texdir+"/pecorrcm_plane_y"+bb}, {"", ""}, {"label1", "label2"});
808 
809  fTex->NewPage();
810  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
811  fTex->AddFigureRow({texdir+"/pe_plane_x"+bb, texdir+"/pe_plane_y"+bb }, {"", ""}, {"label1", "label2"});
812  fTex->AddFigureRow({texdir+"/pecorr_plane_x"+bb, texdir+"/pecorr_plane_y"+bb}, {"", ""}, {"label1", "label2"});
813  fTex->AddFigureRow({texdir+"/cm_plane_x"+bb, texdir+"/cm_plane_y"+bb }, {"", ""}, {"label1", "label2"});
814 
815  fTex->NewPage();
816 
817  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
818  fTex->AddFigureRow({texdir+"/pecm_diblock_x"+bb, texdir+"/pecm_diblock_y"+bb}, {"", ""}, {"label1", "label2"});
819  fTex->AddFigureRow({texdir+"/pecorrcm_diblock_x"+bb, texdir+"/pecorrcm_diblock_y"+bb}, {"", ""}, {"label1", "label2"});
820 
821  fTex->NewPage();
822  fTex->AddFigureRow({texdir+"/legend"}, {""}, {""});
823  fTex->AddFigureRow({texdir+"/pe_diblock_x"+bb, texdir+"/pe_diblock_y"+bb }, {"", ""}, {"label1", "label2"});
824  fTex->AddFigureRow({texdir+"/pecorr_diblock_x"+bb, texdir+"/pecorr_diblock_y"+bb}, {"", ""}, {"label1", "label2"});
825  fTex->AddFigureRow({texdir+"/cm_diblock_x"+bb, texdir+"/cm_diblock_y"+bb }, {"", ""}, {"label1", "label2"});
826 
827  fTex->NewPage();
828 }
829 
830 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
831 /*
832 void ValidCalibPDF::FillEventPlots()
833 {
834  // The filled branches are protected members of
835  // CalibAnaPlot, so can be used here
836 
837  evt_nhits_tricell
838 
839 }
840 */
841 
842 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
844 {
845  /*
846  std::cout << fDet->IsND() << ", " << fDet->IsFD() << ", "
847  << w << ", "
848  << path << ", "
849  << pecorr << ", "
850  << pe << ", "
851  << cmFromEnd << ", " << std::endl;
852  */
853  if(fDet->IsTB())
854  return ( w > -80 && w < 80 &&
855  path > 0. &&
856  pecorr > 0. && pecorr/path < 100. &&
857  pe > 0. &&
858  cmFromEnd > 100 && cmFromEnd < 200 );
859 
860  if(fDet->IsND())
861  return ( w > -100 && w < 100 &&
862  path > 0. &&
863  pecorr > 0. && pecorr/path < 100. &&
864  pe > 0. &&
865  cmFromEnd > 100 && cmFromEnd < 200 ); // --> in ND v14
866 
867  if(fDet->IsFD())
868  return ( w > 200 && w < 600 &&
869  path > 0. &&
870  pecorr > 0. && pecorr/path < 100. &&
871  pe > 0. &&
872  cmFromEnd > 100 && cmFromEnd < 200 );
873  return false;
874 }
875 
876 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
878 {
879  // The filled branches are protected members of
880  // CalibAnaPlot, so can be used here
881 
882  // Temporary hack to calibrate in a finer run division
883  if( fCurrentSample=="nd_data_ep4b.1" && run > 11725) return;
884  if( fCurrentSample=="nd_data_ep4b.2" && run < 11726) return;
885 
886  // It occurs to me the "no cuts" plots we make from the stopping cosmics
887  // are heavily biased in the lower part of the detector by the dE/dx
888  // increase near the end of the track.
889  // Look with the track window cut, instead, so the plots can be
890  // interpreted as coming from a dE/dx uniform accross the detector
891 
892  if(!fDet->IsTB()){
893  if( cmFromEnd < 100 || cmFromEnd > 200 ) return;
894  }
895 
896  bool abscal_hit = PassesAbsCalCut();
897  if( fApplyAbsCutToAll && !abscal_hit ) return;
898 
900  //std::cout << brightbin << std::endl;
901  std::string vstr;
902  if(view==0){ vstr="x"; }
903  if(view==1){ vstr="y"; }
904 
905  // hack out fb bin for now
906  //std::vector< std::string > all_and_perBin = {"", "_b"+b_i};
907  std::vector< std::string > all_and_perBin = {""};
908  //if(b_i=="0" || b_i=="11") all_and_perBin.push_back("_b"+b_i);
909  for( auto& bstr : all_and_perBin ){
910  std::string end = vstr + bstr;
912 
913  float pecm = pe/path;
914  fPlots["pecm_w_"+end] ->Fill(w,pecm);
915  fPlots["pecm_cell_"+end] ->Fill(cell,pecm);
916  fPlots["pecm_plane_"+end] ->Fill(plane,pecm);
917  fPlots["pecm_time_"+end] ->Fill(evt_time,pecm);
918  fPlots["pecm_time_"+end+"_db"+db] ->Fill(evt_time,pecm);
919  fPlots["pecm_diblock_"+end] ->Fill(diblock,pecm);
920  f2dPlots["pecm_cell_plane_"+end] ->Fill(plane,cell,pecm);
921  float pecorrcm = pecorr/path;
922  fPlots["pecorrcm_w_"+end] ->Fill(w,pecorrcm);
923  fPlots["pecorrcm_cell_"+end] ->Fill(cell,pecorrcm);
924  fPlots["pecorrcm_plane_"+end] ->Fill(plane,pecorrcm);
925  fPlots["pecorrcm_time_"+end] ->Fill(evt_time,pecorrcm);
926  fPlots["pecorrcm_diblock_"+end] ->Fill(diblock,pecorrcm);
927  f2dPlots["pecorrcm_cell_plane_"+end]->Fill(plane,cell,pecorrcm);
928  // pe
929  fPlots["pe_w_"+end] ->Fill(w,pe);
930  fPlots["pe_cell_"+end] ->Fill(cell,pe);
931  fPlots["pe_plane_"+end] ->Fill(plane,pe);
932  fPlots["pe_time_"+end] ->Fill(evt_time,pe);
933  fPlots["pe_diblock_"+end] ->Fill(diblock,pe);
934  f2dPlots["pe_cell_plane_"+end] ->Fill(plane,cell,pe);
935  // pecorr
936  fPlots["pecorr_w_"+end] ->Fill(w,pecorr);
937  fPlots["pecorr_cell_"+end] ->Fill(cell,pecorr);
938  fPlots["pecorr_plane_"+end] ->Fill(plane,pecorr);
939  fPlots["pecorr_time_"+end] ->Fill(evt_time,pecorr);
940  fPlots["pecorr_diblock_"+end] ->Fill(diblock,pecorr);
941  f2dPlots["pecorr_cell_plane_"+end]->Fill(plane,cell,pecorr);
942  // path
943  fPlots["cm_w_"+end] ->Fill(w,path);
944  fPlots["cm_cell_"+end] ->Fill(cell,path);
945  fPlots["cm_plane_"+end] ->Fill(plane,path);
946  fPlots["cm_time_"+end] ->Fill(evt_time,path);
947  fPlots["cm_diblock_"+end] ->Fill(diblock,path);
948  f2dPlots["cm_cell_plane_"+end] ->Fill(plane,cell,path);
949  // recoMeV
950  if(fShowRecoMeV){
951  float recomevcm = recoMeV/path;
952  fPlots["recomevcm_w_"+end] ->Fill(w,recomevcm);
953  fPlots["recomevcm_cell_"+end] ->Fill(cell,recomevcm);
954  fPlots["recomevcm_plane_"+end] ->Fill(plane,recomevcm);
955  fPlots["recomevcm_time_"+end] ->Fill(evt_time,recomevcm);
956  fPlots["recomevcm_diblock_"+end] ->Fill(diblock,recomevcm);
957  f2dPlots["recomevcm_cell_plane_"+end] ->Fill(plane,cell,recomevcm);
958  fPlots["recomevcm_time_"+end+"_db"+db] ->Fill(evt_time,recomevcm);
959  }
960 
961  // Diblock breakdown over time
962  fPlots["pecm_time_"+end+"_db"+db] ->Fill(evt_time,pecm);
963  fPlots["pecorrcm_time_"+end+"_db"+db] ->Fill(evt_time,pecorrcm);
964  fPlots["pe_time_"+end+"_db"+db] ->Fill(evt_time,pe);
965  fPlots["pecorr_time_"+end+"_db"+db] ->Fill(evt_time,pecorr);
966  fPlots["cm_time_"+end+"_db"+db] ->Fill(evt_time,path);
967 
968  // number of hits passing whatever cuts applied
969  fPlots["nhits_w_"+end] ->Fill(w);
970  fPlots["nhits_cell_"+end] ->Fill(cell);
971  fPlots["nhits_plane_"+end] ->Fill(plane);
972  fPlots["nhits_time_"+end] ->Fill(evt_time);
973  fPlots["nhits_diblock_"+end] ->Fill(diblock);
974  f2dPlots["nhits_cell_plane_"+end] ->Fill(plane,cell);
975 
976 
977  if(abscal_hit){
978  fPlots["nhits_meu_"+end]->Fill(pecorrcm);
979  fPlots["nhits_mev_"+end]->Fill(trueMeV/path);
980  }
981  }
982 
983 }
984 
985 
986 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
988 {
989  fHitVars="w, view, pe, pecorr, path, trueMeV, recoMeV, cell, plane, evt_time, cmFromEnd, run, diblock,";
990 }
991 
992 
993 //^~^~^~^~^~^~^~^~^~^~^~^~^<~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
995 {
996  std::map<unsigned short,bool> BrightnessBin;
997  BrightnessBin[12]=true; // 12 is all bins combined
998  //BrightnessBin[0]=true;
999  //BrightnessBin[6]=true;
1000  //BrightnessBin[11]=true;
1001  int nBins = 120;
1002 
1003  int ncells = fDet->NCells();
1004  int nplanes = fDet->NPlanes();
1005 
1006  std::vector< PlotAxis > xaxes = {
1007  {"w", "W (cm)", nBins, -fDet->HalfW(), fDet->HalfW() },
1008  {"cell", "Cell", (int)ncells, 0, (float)ncells },
1009  {"plane", "Plane", (int)nplanes/2, 0, (float)nplanes },
1010  {"time", "Unix Time (s)", 1200, 14e8, 16e8 },
1011  {"diblock", "Diblock", (int)fDet->NDiblocks(), 0, (float)fDet->NDiblocks() }
1012  };
1013 
1014  std::vector< ContentAxis > content_axes = {
1015  {"pecm", "PE/cm"},
1016  {"pecorrcm", "PECorr/cm"},
1017  {"pe", "PE"},
1018  {"pecorr", "PECorr"},
1019  {"cm", "Pathlength (cm)"},
1020  {"nhits", "N Hits"}
1021  };
1022  if(fShowRecoMeV) content_axes.push_back({"recomevcm", "Reco dE/dx (MeV/cm)"});
1023 
1024  std::vector< std::pair< PlotAxis, PlotAxis > > th2_axes = {
1025  { {"plane", "Plane", (int)nplanes/2, 0, (float)nplanes }, // cell vs plane
1026  {"cell", "Cell", (int)ncells, 0, (float)ncells } },
1027  };
1028 
1029 
1030  std::string vstr;
1031  std::vector< std::string > views = {"x", "y"};
1032  for( auto& v : views){
1033  if(v=="x") vstr = "X";
1034  if(v=="y") vstr = "Y";
1035 
1036  for(unsigned short i=0; i<=12; i++){
1037  if(BrightnessBin.count(i)!=1) continue;
1038  // 12th one is all bins
1039  std::string bstr = "";
1040  if(i<12) bstr = "_b"+std::to_string(i);
1041 
1042  for(auto& xax : xaxes){
1043  for(auto& yax : content_axes){
1044 
1045  std::string name = yax.label +"_"+ xax.label +"_"+ v + bstr;
1046 
1047  if(yax.label=="nhits")
1048  fPlots[name] = new TH1F( name.c_str(),
1049  (vstr+" View;" + xax.title +";"+ yax.title).c_str(),
1050  xax.nbins, xax.min, xax.max );
1051  else
1052  fPlots[name] = new TProfile( name.c_str(),
1053  (vstr+" View;" + xax.title +";"+ yax.title).c_str(),
1054  xax.nbins, xax.min, xax.max );
1055 
1056  // study drift over time per diblock
1057  if(xax.label=="time" && yax.label!="nhits")
1058  for(size_t db=1; db<=14; ++db)
1059  fPlots[name+"_db"+std::to_string(db)] =
1060  new TProfile( (name+"_db"+std::to_string(db)).c_str(),
1061  (vstr+" View;" + xax.title +";"+ yax.title).c_str(),
1062  xax.nbins, xax.min, xax.max );
1063 
1064  }
1065  }
1066 
1067 
1068  // The special absolute calibration PECorr/cm distributions
1069  std::string name = "nhits_meu_"+ v + bstr;
1070  fPlots[name] = new TH1F( name.c_str(),
1071  (vstr+" View Abs. Cal. Hits;PECorr/cm;N Hits").c_str(),
1072  1000, 0, 100. );
1073  name = "nhits_mev_"+ v + bstr;
1074  fPlots[name] = new TH1F( name.c_str(),
1075  (vstr+" View Abs. Cal. Hits;True dE/dx (MeV/cm);N Hits").c_str(),
1076  1000, 0, 6. );
1077 
1078 
1079  // 2d plots
1080  for( auto& th2ax : th2_axes ){
1081  for(auto& zax : content_axes){
1082  auto xax = th2ax.first;
1083  auto yax = th2ax.second;
1084  std::string name = zax.label +"_"+ yax.label +"_"+ xax.label +"_"+ v + bstr;
1085 
1086  if(zax.label=="nhits")
1087  f2dPlots[name] = new TH2F( name.c_str(),
1088  (vstr+" View;" + xax.title +";"+ yax.title).c_str(),
1089  xax.nbins, xax.min, xax.max,
1090  yax.nbins, yax.min, yax.max );
1091  else
1092  f2dPlots[name] = new TProfile2D( name.c_str(),
1093  (vstr+" View;" + xax.title +";"+ yax.title).c_str(),
1094  xax.nbins, xax.min, xax.max,
1095  yax.nbins, yax.min, yax.max );
1096  }
1097  } // all 2d plots
1098  } // all brightness bins
1099  } // each view
1100 
1101 }
1102 
1103 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
1105 {
1106 
1107  if(fDet->IsND()){
1108 
1109  for( auto& ep : fEpochList ){
1110  if(ep=="4b.1" || ep=="4b.2")
1111  fSamples["nd_data_ep"+ep]="/nova/ana/users/talion/miniprod5/pcliststop_miniprod5.b/nd_data_epoch4b/calibana/merged.root";
1112  else
1113  fSamples["nd_data_ep"+ep]="/nova/ana/users/talion/miniprod5/pcliststop_miniprod5.b/nd_data_epoch"+ep+"/calibana/merged.root";
1114  }
1115 
1116  fSamples["nd_mc"]="/nova/ana/users/abooth/miniprod5.b/nd_mc/merged.root";
1117  //fSamples["mc_2016"]="/nova/ana/users/abooth/R17-03-09-prod3genie.c/merged.root";
1118  }
1119 
1120  else if (fDet->IsFD()){
1121 
1122  for( auto& ep : fEpochList )
1123  fSamples["fd_data_ep"+ep]="/nova/ana/users/talion/miniprod5_calibana/fd_data_epoch"+ep+"/merged.root";
1124 
1125  // low gain:
1126  fSamples["fd_mc_p2"]="/nova/ana/users/talion/miniprod5_calibana/fd_mc_lowgain/merged.p2.root";
1127  // high gain:
1128  fSamples["fd_mc_p3"]="/nova/ana/users/talion/miniprod5_calibana/fd_mc_highgain/merged.p3.root";
1129  //fSamples["fd_mc_p5"]="/nova/ana/users/talion/miniprod5_calibana/fd_mc_highgain/merged.p3.root";
1130 
1131  //fSamples["fd_mc_p2_2016"]="/nova/ana/users/talion/miniprod5_calibana/fd_mc_lowgain_2016/merged.root";
1132  //fSamples["fd_mc_p3_2016"]="/nova/ana/users/talion/miniprod5_calibana/fd_mc_highgain_2016/merged.root";
1133 
1134  }
1135 
1136  else if (fDet->IsTB()){
1137 
1138  fSamples["tb_mc"]="/nova/ana/users/kmulder/28062019_testbeam_abs_calib/tb_mc_period_1/calibana_tb_mc_period_1.root";
1139  fSamples["tb_data_p1"]="/nova/ana/users/kmulder/28062019_testbeam_abs_calib/tb_data_period_1/calibana_tb_data_period_1.root";
1140 
1141  }
1142 
1143 
1144  for(auto &s: fSamples)
1145  if(s.first.find("mc")!=std::string::npos)
1146  fMCSamples.push_back(s.first);
1147 }
1148 
1149 
1150 //^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
1152 {
1153  // Configure which MC sample corresponds to what
1154  if(sample.find("nd_")!=std::string::npos) return "nd_mc";
1155  // FD data periods 1-2 are low gain (p2) MC
1156  if( sample.find("fd_data_ep1")!=std::string::npos ||
1157  sample.find("fd_data_ep2")!=std::string::npos ) return "fd_mc_p2";
1158  // The rest of the FD data periods are high gain (p3) MC
1159  if(sample.find("fd_data_ep")!=std::string::npos) return "fd_mc_p3";
1160  return sample;
1161 }
int numHits_x
Definition: CalibAnaTypes.h:20
std::vector< std::string > header
Definition: TexBuilder.h:7
const XML_Char * name
Definition: expat.h:151
void AddTable(Table table, bool drawnow=false, bool makecsv=true)
Definition: TexBuilder.h:237
std::map< std::string, TH1 * > fPlots
Definition: CalibAnaPlot.h:57
std::map< std::string, std::string > GetSamples()
Definition: CalibAnaPlot.h:28
int nBins
Definition: plotROC.py:16
unsigned int NCells()
Definition: Detector.h:12
void DrawTable(std::string tabname)
Definition: TexBuilder.h:276
std::string ToString(float num, unsigned short digits=0)
void NewPage()
Definition: TexBuilder.h:165
struct Plot Plot
Definition: DrawUtils.h:15
std::vector< std::string > fEpochList
double meu_y
Definition: CalibAnaTypes.h:6
std::map< std::string, TH2 * > f2dPlots
Definition: CalibAnaPlot.h:58
void SetPerformAbsCal(bool set)
void PrintEpochLengths()
Definition: CalibAnaPlot.h:517
std::string fOutPath
Definition: CalibAnaPlot.h:62
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
TexBuilder * fTex
bool IsTB()
Definition: Detector.h:19
Float_t ss
Definition: plot.C:24
void GetCSVRow(float &meu, float &meuerr, float &gev, std::string calibtag, std::string det, std::string dataormc, std::string runrange)
std::map< std::string, uint32_t > EndOfEpochMap()
Definition: Detector.h:24
bool IsFD()
Definition: Detector.h:21
std::string fTexConfig
double mevErr
Definition: CalibAnaTypes.h:18
void MakeDir(std::string path)
Definition: CalibAnaPlot.h:484
void Close()
Definition: TexBuilder.h:85
void MakeTex(std::string config_path)
uint32_t evt_time
Definition: CalibAnaPlot.h:80
void AddText(std::string text)
Definition: TexBuilder.h:92
std::string LegendLabel(std::string key)
Definition: CalibAnaPlot.h:239
double meuErr
Definition: CalibAnaTypes.h:17
void Initialize()
Definition: CalibAnaPlot.h:138
float HalfW()
Definition: Detector.h:15
void SetTreeLoopLimit(unsigned int limit)
Definition: CalibAnaPlot.h:31
std::string name
Definition: TexBuilder.h:4
TH1 * h
Definition: DrawUtils.h:19
const XML_Char * s
Definition: expat.h:262
std::string fHitVars
Definition: CalibAnaPlot.h:67
void AbsoluteCalibration(std::string sample)
std::vector< std::vector< std::string > > rows
Definition: TexBuilder.h:8
ValidCalibPDF(std::string outpath, std::string tex_cfg)
void SetApplyAbsCutToAll(bool set)
====================================================================== ///
Definition: CutFlow_Data.C:28
void SetEpochList(std::string tex_cfg)
std::map< std::string, AbsCal > fAbsCals
double mevErr_x
Definition: CalibAnaTypes.h:12
std::string sample
Definition: CalibAnaTypes.h:4
double meu_x
Definition: CalibAnaTypes.h:5
std::string csvdir
Definition: TexBuilder.h:11
int numHits_y
Definition: CalibAnaTypes.h:21
const double j
Definition: BetheBloch.cxx:29
double mev_y
Definition: CalibAnaTypes.h:8
double mev
Definition: CalibAnaTypes.h:16
void ProcessTrees(std::string sample)
Definition: CalibAnaPlot.h:329
string pname
Definition: eplot.py:33
std::string MCSample(std::string sample)
std::string name()
Definition: Detector.h:9
void MakeValidCalibPDF(std::string outpath, std::string tex_config, std::string sample, bool justdraw=false, unsigned int stride=0, unsigned int limit=0)
Definition: run.py:1
const char sep
void SetTreeLoopStride(unsigned int stride)
Definition: CalibAnaPlot.h:32
OStream cout
Definition: OStream.cxx:6
int nplanes
Definition: geom.C:145
const int wid1
void EssentialComparisonTex(std::string texdir, unsigned short bright_bin=std::numeric_limits< unsigned short >::max())
void HaddAllDataFile()
Definition: CalibAnaPlot.h:459
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
unsigned int NPlanes()
Definition: Detector.h:13
bool IsND()
Definition: Detector.h:20
void Draw(std::string pdfdir, std::vector< Plot > plots, Plot denom_plot, std::string ratio_title="")
Definition: CalibAnaPlot.h:168
void SetShowRecoMeV(bool set)
void DriftTex(std::string texdir, unsigned short bright_bin=std::numeric_limits< unsigned short >::max())
int num
Definition: f2_nu.C:119
std::string title
Definition: TexBuilder.h:5
std::map< std::string, std::string > fRunsMap
std::string csvname
Definition: TexBuilder.h:10
std::vector< std::string > fMCSamples
double meuErr_y
Definition: CalibAnaTypes.h:11
int ncells
Definition: geom.C:124
double mev_x
Definition: CalibAnaTypes.h:7
void Section(std::string secname)
Definition: TexBuilder.h:149
unsigned int NDiblocks()
Definition: Detector.h:14
double meuErr_x
Definition: CalibAnaTypes.h:10
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void SubSection(std::string subsecname)
Definition: TexBuilder.h:157
std::map< std::string, std::string > fSamples
Definition: CalibAnaPlot.h:59
void AddFigureRow(std::vector< std::string > figure, std::vector< std::string > caption, std::vector< std::string > label)
Definition: TexBuilder.h:98
std::string RunRange(std::string sample, std::string oldcalibtag="")
double mevErr_y
Definition: CalibAnaTypes.h:13
const int wid2
double meu
Definition: CalibAnaTypes.h:15
void DiblockTex(std::string texdir, std::vector< std::string > configs)
std::string fCurrentSample
Definition: CalibAnaPlot.h:60
Detector * fDet
Definition: CalibAnaPlot.h:68
enum BeamMode string