DCMOffsetCalculator.C
Go to the documentation of this file.
1 #include <vector>
2 #include <algorithm>
3 #include <utility>
4 #include <map>
5 #include <iterator>
6 #include <iostream>
7 #include <fstream>
8 #include <stdio.h>
9 #include <string>
10 #include <cmath>
11 #include "TTree.h"
12 #include "TChain.h"
13 #include "TH1F.h"
14 #include "TH2F.h"
15 #include "TGraphErrors.h"
16 #include "TFile.h"
17 #include "TMath.h"
18 #include "TStyle.h"
19 #include "TColor.h"
20 #include <sstream>
21 #include "TCanvas.h"
22 #include "TMatrixD.h"
23 #include "TVectorD.h"
24 #include "TDecompSVD.h"
25 #include "TDecompLU.h"
26 #include <stdint.h>
27 
28 
29 using namespace std;
30 
31 struct DCMOffset{
32  int dcm; // dcm number
33  double offset; // dcm timing offset
34  double sigma; // uncertainty in offset
35  double stats; // number of track segments for dcm
36  uint32_t sTime; // start time offset is valid for
37  uint32_t eTime; // end time offset is valid for
38  };
39 
40 
42 
43  //****************************************************************
44  //Define constants needed for DCM calibration
45  const static int fNumDCM = 48; // Number of DCMs calibration will be performed on
46  const static int fMinDCM = 1; // Number of lowest DCM in calibration set
47  const static int fMaxDCM = 48; // Number of highest DCM in calibration set
48  const static int fDCMPerDB = 12; // Number of DCMs per diblock
49  const static int fMinCount = 20; // Minimum number of counts for a DCM pair
50  const static int fRefDCM = 6; // DCM number to use as refence, default is 6
51  const static double fRefOffset = 0; // absolute offset, in ns, to fix refence DCM to, default is 0.0
52  const static int fMinRun = 0; // minimum run number to use in study
53  const static int fMaxRun = 9999999; // maximum run number to use in study
54  const static int fMinSubrun = 20; // minimum number of subruns to keep run
55  bool fIsMC = false; // influences name of output csv file
56  //****************************************************************
57 
58 
59  //create a chain for the timing calibration
60  TChain *timingchain = new TChain("timingcal/TimingCalib");
61 
62  //load list of timing calibration files to take as input
63  string line;
64  ifstream inFile("/nova/ana/users/edniner/timingcalibration_Official2/file_list.txt");
65  if (inFile.is_open()){
66  while (inFile.good()){
67  getline(inFile,line);
68  if (line == "") continue;
69  //std::string test = line.substr(63,5);
70  //if (atoi(test.c_str()) < fMinRun) continue;
71  //if (atoi(test.c_str()) >= fMaxRun) continue;
72  std::cout <<"Adding file: " << line << std::endl;
73  timingchain->AddFile(line.c_str());
74  }
75  inFile.close();
76  }
77 
78  TTree* timingcal = timingchain;
79 
80  //ntuple variables
81  //event statistics
82  int run;
83  int subrun;
84  int event;
85  uint32_t evtTime;
86 
87  //hold dcm level hit information
88  vector<float> *dcmnum = 0;
89  vector<float> *dcmhits = 0;
90  vector<float> *dcmsimtime = 0;
91 
92  //turn off all branches, only use the ones we care about
93  timingcal->SetBranchStatus("*", 0);
94  //turn back on branches of interest
95  timingcal->SetBranchStatus("run",1);
96  timingcal->SetBranchStatus("subrun",1);
97  timingcal->SetBranchStatus("event",1);
98  timingcal->SetBranchStatus("evtTime",1);
99  timingcal->SetBranchStatus("dcmnum",1);
100  timingcal->SetBranchStatus("dcmhits",1);
101  timingcal->SetBranchStatus("dcmsimtime",1);
102 
103  timingcal->SetBranchAddress("run",&run);
104  timingcal->SetBranchAddress("subrun",&subrun);
105  timingcal->SetBranchAddress("event",&event);
106  timingcal->SetBranchAddress("evtTime",&evtTime);
107  timingcal->SetBranchAddress("dcmnum", &dcmnum);
108  timingcal->SetBranchAddress("dcmhits", &dcmhits);
109  timingcal->SetBranchAddress("dcmsimtime", &dcmsimtime);
110 
111 
112 
113 
114  //keep stats for entire group of runs
115  int countDCMOffsetArray[fNumDCM][fNumDCM] = {0};
116  double sumDCMOffsetArray[fNumDCM][fNumDCM] = {0.};
117  double sumSquareDCMOffsetArray[fNumDCM][fNumDCM] = {0.};
118  double meanDCMOffsetArray[fNumDCM][fNumDCM] = {0.};
119  double weightDCMOffsetArray[fNumDCM][fNumDCM] = {0.};
120  double absoluteDCMOffset[fNumDCM] = {0.};
121  double absoluteDCMErrors[fNumDCM] = {0.};
122  int countDCMUse[fNumDCM] = {0};
123  uint32_t minTime;
124  uint32_t maxTime;
125  //keep stats for individual run
126  int countDCMOffsetArrayRun[fNumDCM][fNumDCM] = {0};
127  double sumDCMOffsetArrayRun[fNumDCM][fNumDCM] = {0.};
128  double sumSquareDCMOffsetArrayRun[fNumDCM][fNumDCM] = {0.};
129  double meanDCMOffsetArrayRun[fNumDCM][fNumDCM] = {0.};
130  double weightDCMOffsetArrayRun[fNumDCM][fNumDCM] = {0.};
131  double absoluteDCMErrorsRun[fNumDCM] = {0.};
132  int countDCMUseRun[fNumDCM] = {0};
133  uint32_t minTimeRun;
134  uint32_t maxTimeRun;
135 
136  std::map<int, std::vector<DCMOffset> > offsetMap;
137 
138  minTimeRun = std::numeric_limits<uint32_t>::max()-1;
139  maxTimeRun = std::numeric_limits<uint32_t>::min();
142  int crun = -1;
143  int csrun = -1;
144  int srun = 0;
145 
146  for(unsigned int k=0; k<timingcal->GetEntries(); ++k){
147 
148  timingcal->GetEntry(k);
149  if (k % 100000 == 0) cout<<"processed: "<<k<<" entries"<<endl;
150 
151  //Files at this point should have passed all selection cuts,
152  //so no addtional cuts are made here.
153  if (csrun == -1) csrun = subrun;
154 
155  if (crun == -1) {
156  crun = run;
157  std::cout<<"Currently processing run: "<<run<<std::endl;
158  }
159 
160  if (crun != run){
161  std::cout<<"Completed run: "<<crun<<" moving on to run: "<<run<<std::endl;
162  crun = run;
163 
164  if (run < fMinRun) continue;
165  if (run >= fMaxRun) continue;
166 
167  if (srun < fMinSubrun){
168  srun = 0;
169  //reset run arrays
170  for (int i=0; i<fNumDCM; ++i){
171  absoluteDCMErrorsRun[i] = 0;
172  countDCMUseRun[i] = 0;
173  for (int j=0; j<fNumDCM; ++j){
174  countDCMOffsetArrayRun[i][j] = 0;
175  sumDCMOffsetArrayRun[i][j] = 0;
176  sumSquareDCMOffsetArrayRun[i][j] = 0;
177  meanDCMOffsetArrayRun[i][j] = 0;
178  weightDCMOffsetArrayRun[i][j] = 0;
179  }
180  }
181  continue;
182  }
183  srun = 0;
184 
185  //solve run level offsets
186  for(int i=0; i<fNumDCM-1; ++i){
187  for(int j=i+1; j<fNumDCM; ++j){
188  double sum = sumDCMOffsetArrayRun[i][j];
189  int count = countDCMOffsetArrayRun[i][j];
190  double sumSqr = sumSquareDCMOffsetArrayRun[i][j];
191  if (count >= fMinCount){
192  meanDCMOffsetArrayRun[i][j] = sum/count;
193  double sigma = sqrt((count*sumSqr - sum*sum)/(count*(count-1.0)));
194  weightDCMOffsetArrayRun[i][j] = 1.0/(sigma*sigma);
195  }
196  else {
197  meanDCMOffsetArrayRun[i][j] = 0;
198  weightDCMOffsetArrayRun[i][j] = 0;
199  }
200  //std::cout<<i<<" "<<j<<" "<<count<<" "<<sum<<" "<<sum/count<<std::endl;
201  }
202  }
203 
204  const int num= fNumDCM;
205  int refIndx = fRefDCM - fMinDCM;
206  TMatrixD A;
207  A.ResizeTo(fNumDCM, fNumDCM);
208  double y[num];
209  for (int x=0; x<num; ++x) y[x]=0.;
210 
211 
212  for(int i=0; i<fNumDCM; ++i){
213  double temp[num];
214  for (int x=0; x<num; ++x) temp[x]=0.;
215  for(int j=0; j<fNumDCM; ++j){
216  if (i==j) continue;
217  if (j<i){
218  temp[i] += weightDCMOffsetArrayRun[j][i]*2.0;
219  temp[j] += -weightDCMOffsetArrayRun[j][i]*2.0;
220  y[i] += meanDCMOffsetArrayRun[j][i]*weightDCMOffsetArrayRun[j][i]*2.0;
221  }
222  if (j>i){
223  temp[i] += weightDCMOffsetArrayRun[i][j]*2.0;
224  temp[j] += -weightDCMOffsetArrayRun[i][j]*2.0;
225  y[i] += -meanDCMOffsetArrayRun[i][j]*weightDCMOffsetArrayRun[i][j]*2.0;
226  }
227  }
228 
229  if (i == refIndx){
230  for (int j=0; j<fNumDCM; ++j) temp[j] = 0.0;
231  temp[refIndx]=1.0;
232  }
233 
234  y[i]*=-1.0;
235  TVectorD tempvec;
236  tempvec.Use(fNumDCM,temp);
237  TMatrixDRow(A,i) = tempvec;
238  }
239 
240  y[refIndx] = fRefOffset;
241  TVectorD yvec;
242  yvec.Use(fNumDCM,y);
243 
244  //calculate uncertainties
245  for(int i=0; i<fNumDCM; ++i){
246  for(int j=0; j<fNumDCM; ++j){
247  if (i==j) continue;
248  if (j<i) absoluteDCMErrorsRun[i] += 2.0*weightDCMOffsetArrayRun[j][i];
249  else if (j>i) absoluteDCMErrorsRun[i] += 2.0*weightDCMOffsetArrayRun[i][j];
250  }
251  absoluteDCMErrorsRun[i] = sqrt(2.0/absoluteDCMErrorsRun[i]);
252  }
253  absoluteDCMErrorsRun[refIndx] = 0.0;
254 
255  //solve for absolute offsets
256  TDecompLU lu(A);
257  bool inv = lu.Decompose();
258  if (inv) lu.Invert(A);
259  //std::cout<<A.IsValid()<<std::endl;
260  if (inv) yvec *= A;
261 
262 
263  for (int i=0; i<fNumDCM; ++i){
264  if (!inv) break;
265  DCMOffset dcmOffset;
266  dcmOffset.dcm = i+fMinDCM;
267  dcmOffset.offset = yvec[i];
268  dcmOffset.sigma = absoluteDCMErrorsRun[i];
269  dcmOffset.stats = countDCMUseRun[i];
270  dcmOffset.sTime = minTimeRun;
271  dcmOffset.eTime = maxTimeRun;
272  offsetMap[run].push_back(dcmOffset);
273  }
274 
275  //reset run arrays
276  for (int i=0; i<fNumDCM; ++i){
277  absoluteDCMErrorsRun[i] = 0;
278  countDCMUseRun[i] = 0;
279  for (int j=0; j<fNumDCM; ++j){
280  countDCMOffsetArrayRun[i][j] = 0;
281  sumDCMOffsetArrayRun[i][j] = 0;
282  sumSquareDCMOffsetArrayRun[i][j] = 0;
283  meanDCMOffsetArrayRun[i][j] = 0;
284  weightDCMOffsetArrayRun[i][j] = 0;
285  }
286  }
287 
288  minTimeRun = std::numeric_limits<uint32_t>::max()-1;
289  maxTimeRun = std::numeric_limits<uint32_t>::min();
290 
291  crun = run;
292  } //end condition for changing run numbers
293 
294  if (run < fMinRun) continue;
295  if (run >= fMaxRun) continue;
296 
297  if (csrun != subrun){csrun = subrun; srun++;}
298 
299  if (evtTime>maxTime) maxTime = evtTime;
300  if (evtTime<minTime) minTime = evtTime;
301  if (evtTime>maxTimeRun) maxTimeRun = evtTime;
302  if (evtTime<minTimeRun) minTimeRun = evtTime;
303 
304  vector<pair<unsigned int, double> > dcmtimes;
305  vector<pair<unsigned int, unsigned int> > dcmHits;
306  //loop over all dcms used in this track
307  //a cut has already been applied in the ana module
308  //to select the minumum number of hits allowed in a dcm
309  //the default cut is 10 hits
310  for (unsigned int j=0; j<dcmnum->size(); ++j){
311  unsigned int num = (*dcmnum)[j];
312  double stime = (*dcmsimtime)[j];
313  dcmtimes.push_back(make_pair(num,stime));
314  //std::cout<<"dcm "<<j<<" "<<num<<" "<<stime<<std::endl;
315  }
316  std::sort(dcmtimes.begin(), dcmtimes.end());
317  for (unsigned int dcm1=0; dcm1<dcmtimes.size(); ++dcm1){
318  if ((dcmtimes[dcm1].first < fMinDCM) || (dcmtimes[dcm1].first > fMaxDCM)) continue;
319  unsigned int ind1 = dcmtimes[dcm1].first-fMinDCM;
320  countDCMUse[ind1]++;
321  //countDCMUseRun[ind1]++;
322  for (unsigned int dcm2=dcm1+1; dcm2 <dcmtimes.size(); ++dcm2){
323  if ((dcmtimes[dcm2].first < fMinDCM) || (dcmtimes[dcm2].first > fMaxDCM)) continue;
324  double offset = dcmtimes[dcm1].second-dcmtimes[dcm2].second;
325  unsigned int ind2 = dcmtimes[dcm2].first-fMinDCM;
326  //std::cout<<ind1<<" "<<ind2<<" "<<offset<<std::endl;
327  countDCMOffsetArray[ind1][ind2]++;
328  sumDCMOffsetArray[ind1][ind2] += offset;
329  sumSquareDCMOffsetArray[ind1][ind2] += offset*offset;
330  countDCMOffsetArrayRun[ind1][ind2]++;
331  sumDCMOffsetArrayRun[ind1][ind2] += offset;
332  sumSquareDCMOffsetArrayRun[ind1][ind2] += offset*offset;
333  }
334  }
335 
336  }//end loop over events
337 
338  //now process the last run
339  for(int i=0; i<fNumDCM-1; ++i){
340  for(int j=i+1; j<fNumDCM; ++j){
341  double sum = sumDCMOffsetArrayRun[i][j];
342  double count = countDCMOffsetArrayRun[i][j];
343  double sumSqr = sumSquareDCMOffsetArrayRun[i][j];
344  if (count >= fMinCount){
345  meanDCMOffsetArrayRun[i][j] = sum/count;
346  double sigma = sqrt((count*sumSqr - sum*sum)/(count*(count-1.0)));
347  weightDCMOffsetArrayRun[i][j] = 1.0/(sigma*sigma);
348  }
349  else {
350  meanDCMOffsetArrayRun[i][j] = 0;
351  weightDCMOffsetArrayRun[i][j] = 0;
352  }
353  }
354  }
355 
356  const int num= fNumDCM;
357  int refIndx = fRefDCM - fMinDCM;
358  TMatrixD A;
359  A.ResizeTo(fNumDCM, fNumDCM);
360  double y[num];
361  for (int x=0; x<num; ++x) y[x]=0.;
362 
363 
364  for(int i=0; i<fNumDCM; ++i){
365  double temp[num];
366  for (int x=0; x<num; ++x) temp[x]=0.;
367  for(int j=0; j<fNumDCM; ++j){
368  if (i==j) continue;
369  if (j<i){
370  temp[i] += weightDCMOffsetArrayRun[j][i]*2.0;
371  temp[j] += -weightDCMOffsetArrayRun[j][i]*2.0;
372  y[i] += meanDCMOffsetArrayRun[j][i]*weightDCMOffsetArrayRun[j][i]*2.0;
373  }
374  if (j>i){
375  temp[i] += weightDCMOffsetArrayRun[i][j]*2.0;
376  temp[j] += -weightDCMOffsetArrayRun[i][j]*2.0;
377  y[i] += -meanDCMOffsetArrayRun[i][j]*weightDCMOffsetArrayRun[i][j]*2.0;
378  }
379  }
380 
381  if (i == refIndx){
382  for (int j=0; j<fNumDCM; ++j) temp[j] = 0.0;
383  temp[refIndx]=1.0;
384  }
385 
386  y[i]*=-1.0;
387  TVectorD tempvec;
388  tempvec.Use(fNumDCM,temp);
389  TMatrixDRow(A,i) = tempvec;
390  }
391 
392  y[refIndx] = fRefOffset;
393  TVectorD yvec;
394  yvec.Use(fNumDCM,y);
395 
396  //calculate uncertainties
397  for(int i=0; i<fNumDCM; ++i){
398  for(int j=0; j<fNumDCM; ++j){
399  if (i==j) continue;
400  if (j<i) absoluteDCMErrorsRun[i] += 2.0*weightDCMOffsetArrayRun[j][i];
401  else if (j>i) absoluteDCMErrorsRun[i] += 2.0*weightDCMOffsetArrayRun[i][j];
402  }
403  absoluteDCMErrorsRun[i] = sqrt(2.0/absoluteDCMErrorsRun[i]);
404  }
405  absoluteDCMErrorsRun[refIndx] = 0.0;
406 
407  //solve for absolute offsets
408  TDecompLU lu(A);
409  bool inv = lu.Decompose();
410  if (inv) lu.Invert(A);
411  //std::cout<<A.IsValid()<<std::endl;
412  if (inv) yvec *= A;
413 
414 
415  for (int i=0; i<fNumDCM; ++i){
416  if (!inv) break;
417  DCMOffset dcmOffset;
418  dcmOffset.dcm = i+fMinDCM;
419  dcmOffset.offset = yvec[i];
420  dcmOffset.sigma = absoluteDCMErrorsRun[i];
421  dcmOffset.stats = countDCMUseRun[i];
422  dcmOffset.sTime = minTimeRun;
423  dcmOffset.eTime = maxTimeRun;
424  offsetMap[run].push_back(dcmOffset);
425  }
426 
427 
428  //const int num= fNumDCM;
429  //int refIndx = fRefDCM - fMinDCM;
430  //double y[num];
431  for (int i=0; i<num; ++i) y[i] = 0;
432  TFile* file = new TFile("dcmoffsets_summary_plots.root","recreate");
433  //now process final stats
434  const int runs = offsetMap.size();
435  const int ndcm= fNumDCM;
436  TGraphErrors* grAbsDCM[ndcm];
437  TGraphErrors* grAbsRuns[runs];
438  TGraph* grStatsDCM[ndcm];
439  TGraph* grStatsRuns[runs];
440  TGraphErrors* grAbs;
441  TGraph* grStats;
442  TH2F* fHistoRelMap;
443 
444  double runDCM[ndcm][runs];
445  double runDCMOffset[ndcm][runs];
446  double runDCMOffsetE[ndcm][runs];
447  double runDCMStats[ndcm][runs];
448 
449  grAbs = new TGraphErrors;
450  gDirectory->Append(grAbs);
451  grStats = new TGraph;
452  gDirectory->Append(grStats);
453  fHistoRelMap = new TH2F("fHistoRelMap","Relative offsets (ns), all runs, DCM Y- DCM X;dcm X;dcm Y",
454  fNumDCM,(double)fMinDCM-0.5,(double)fMaxDCM+0.5,fNumDCM,(double)fMinDCM-0.5,(double)fMaxDCM+0.5);
455  for (int i=0; i< runs; ++i){
456  grAbsRuns[i] = new TGraphErrors;
457  gDirectory->Append(grAbsRuns[i]);
458  grStatsRuns[i] = new TGraph;
459  gDirectory->Append(grStatsRuns[i]);
460  }
461  for (int i=0; i< ndcm; ++i){
462  grAbsDCM[i] = new TGraphErrors;
463  gDirectory->Append(grAbsDCM[i]);
464  grStatsDCM[i] = new TGraph;
465  gDirectory->Append(grStatsDCM[i]);
466  }
467 
468 
469  for(int i=0; i<fNumDCM-1; ++i){
470  for(int j=i+1; j<fNumDCM; ++j){
471  double sum = sumDCMOffsetArray[i][j];
472  double count = countDCMOffsetArray[i][j];
473  double sumSqr = sumSquareDCMOffsetArray[i][j];
474  if (count >= fMinCount){
475  meanDCMOffsetArray[i][j] = sum/count;
476  double sigma = sqrt((count*sumSqr - sum*sum)/(count*(count-1.0)));
477  weightDCMOffsetArray[i][j] = 1.0/(sigma*sigma);
478  fHistoRelMap->SetBinContent(j+1,i+1,sum/count);
479  }
480  else {
481  meanDCMOffsetArray[i][j] = 0;
482  weightDCMOffsetArray[i][j] = 0;
483  }
484  //std::cout<<i<<" "<<j<<" "<<count<<" "<<sum<<" "<<sum/count<<std::endl;
485  }
486  }
487 
488 
489  TMatrixD A1;
490  A1.ResizeTo(fNumDCM, fNumDCM);
491  double y1[ndcm];
492  for (int x=0; x<ndcm; ++x) y1[x]=0.;
493 
494  for(int i=0; i<fNumDCM; ++i){
495  double temp[ndcm];
496  for (int x=0; x<ndcm; ++x) temp[x]=0.;
497  for(int j=0; j<fNumDCM; ++j){
498  if (i==j) continue;
499  if (j<i){
500  temp[i] += weightDCMOffsetArray[j][i]*2.0;
501  temp[j] += -weightDCMOffsetArray[j][i]*2.0;
502  y1[i] += meanDCMOffsetArray[j][i]*weightDCMOffsetArray[j][i]*2.0;
503  }
504  if (j>i){
505  temp[i] += weightDCMOffsetArray[i][j]*2.0;
506  temp[j] += -weightDCMOffsetArray[i][j]*2.0;
507  y1[i] += -meanDCMOffsetArray[i][j]*weightDCMOffsetArray[i][j]*2.0;
508  }
509  }
510 
511  if (i == refIndx){
512  for (int j=0; j<fNumDCM; ++j) temp[j] = 0.0;
513  temp[refIndx]=1.0;
514  }
515 
516  y1[i]*=-1.0;
517  TVectorD tempvec;
518  tempvec.Use(fNumDCM,temp);
519  TMatrixDRow(A1,i) = tempvec;
520  }
521 
522  y1[refIndx] = fRefOffset;
523  TVectorD yvec1;
524  yvec1.Use(fNumDCM,y1);
525 
526  //calculate uncertainties
527  for(int i=0; i<fNumDCM; ++i){
528  for(int j=0; j<fNumDCM; ++j){
529  if (i==j) continue;
530  if (j<i) absoluteDCMErrors[i] += 2.0*weightDCMOffsetArray[j][i];
531  else if (j>i) absoluteDCMErrors[i] += 2.0*weightDCMOffsetArray[i][j];
532  }
533  absoluteDCMErrors[i] = sqrt(2.0/absoluteDCMErrors[i]);
534  }
535  absoluteDCMErrors[refIndx] = 0.0;
536 
537  //solve for absolute offsets
538  A1.Invert();
539  yvec1 *= A1;
540 
541  //write results to csv file
542  const TString mcStr = fIsMC ? "_mc" : "";
543  FILE* fConstsCSV = fopen("calib_abs_dcm_delay_consts"+mcStr+".csv","w");
544 
545  fprintf(fConstsCSV, "#dcmoffset,sigma,stats,channel,tv\n");
546 
547  //Multipy by -1 since the correction facto is the opposite sign of the offset
548  //store channel as diblock*1000+dcm
549  //make CVS file with constants to put in DB
550  for (int i=0; i<fNumDCM; ++i){
551  int db= ((i+fMinDCM)/fDCMPerDB)+1;
552  int dcm = (i+fMinDCM)%fDCMPerDB;
553  if ((i+fMinDCM)%fDCMPerDB == 0){
554  db --;
555  dcm = fDCMPerDB;
556  }
557  absoluteDCMOffset[i] = yvec1[i];
558  int chan = db*1000 + dcm;
559  fprintf(fConstsCSV, "%g,%g,%d,%d,%d\n",
560  -1*absoluteDCMOffset[i], absoluteDCMErrors[i], countDCMUse[i], chan, minTime);
561  }
562 
563  fclose(fConstsCSV);
564 
565  int runcount= 0;
566  for (std::map<int, std::vector<DCMOffset> >::iterator itr = offsetMap.begin(); itr != offsetMap.end(); ++itr){
567  int rnum = itr->first;
568  std::vector<DCMOffset> offsets = itr->second;
569 
570  TString dummy;
571  std::ostringstream convert;
572  convert << rnum;
573  dummy = convert.str();
574  grAbsRuns[runcount]->SetName("grAbsOffset_run_"+dummy);
575  grAbsRuns[runcount]->SetTitle("Absolute DCM timing offsets (ns) for run: "+dummy);
576  grAbsRuns[runcount]->Set(offsets.size());
577  grStatsRuns[runcount]->SetName("grStatsOffset_run_"+dummy);
578  grStatsRuns[runcount]->SetTitle("Number of track segments used by dcm for run: "+dummy);
579  grStatsRuns[runcount]->Set(offsets.size());
580  for (unsigned int i=0; i<offsets.size(); ++i){
581  int tdcm = offsets[i].dcm;
582  double toffset = offsets[i].offset;
583  double tsigma = offsets[i].sigma;
584  int tstats = offsets[i].stats;
585  grStatsRuns[runcount]->SetPoint(i,tdcm,tstats);
586  grAbsRuns[runcount]->SetPoint(i,tdcm,toffset);
587  grAbsRuns[runcount]->SetPointError(i,0.0,tsigma);
588  runDCM[tdcm-fMinDCM][runcount]=rnum;
589  runDCMOffset[tdcm-fMinDCM][runcount]=toffset;
590  runDCMOffsetE[tdcm-fMinDCM][runcount]=tsigma;
591  runDCMStats[tdcm-fMinDCM][runcount]=tstats;
592  }
593  runcount++;
594  }
595 
596 
597  for(int i=0; i<ndcm; ++i){
598  int tdcm = i + fMinDCM;
599  TString dummy;
600  std::ostringstream convert;
601  convert << tdcm;
602  dummy = convert.str();
603  grAbsDCM[i]->SetName("grAbsOffset_dcm_"+dummy);
604  grAbsDCM[i]->SetTitle("Absolute DCM timing offsets (ns) across all run for dcm: "+dummy);
605  grAbsDCM[i]->Set(runs);
606  grStatsDCM[i]->SetName("grStatsOffset_dcm_"+dummy);
607  grStatsDCM[i]->SetTitle("Number of track segments used for all runs by dcm: "+dummy);
608  grStatsDCM[i]->Set(runs);
609  for(int j=0; j<runs; ++j){
610  grAbsDCM[i]->SetPoint(j,runDCM[i][j],runDCMOffset[i][j]);
611  grAbsDCM[i]->SetPointError(j,0.0,runDCMOffsetE[i][j]);
612  grStatsDCM[i]->SetPoint(j,runDCM[i][j],runDCMStats[i][j]);
613  }
614  }
615 
616  grAbs->SetName("grAbsOffset_total");
617  grAbs->SetTitle("Absolute DCM timing offsets (ns), all runs");
618  grAbs->Set(ndcm);
619  grStats->SetName("grStatsOffset_total");
620  grStats->SetTitle("Number of track segments used by dcm, all runs");
621  grStats->Set(ndcm);
622 
623  for(int i=0; i< fNumDCM; ++i){
624  grAbs->SetPoint(i,i+fMinDCM,absoluteDCMOffset[i]);
625  grAbs->SetPointError(i,0.0,absoluteDCMErrors[i]);
626  grStats->SetPoint(i,i+fMinDCM,countDCMUse[i]);
627  }
628 
629  file->Write();
630  file->Close();
631 
632 
633 }
void DCMOffsetCalculator()
Float_t y1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
void convert(std::string dir="cc_numu/C12")
Definition: convert.C:107
fclose(fg1)
ifstream inFile
Definition: AnaPlotMaker.h:34
const double j
Definition: BetheBloch.cxx:29
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
static const double A
Definition: Units.h:82
int num
Definition: f2_nu.C:119
TFile * file
Definition: cellShifts.C:17
Double_t sum
Definition: plot.C:31
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
list runs
Definition: shutoffs.py:47
fvar< T > inv(const fvar< T > &x)
Definition: inv.hpp:12