DCMTimingOffset_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: DCMTimingOffset
3 // Module Type: producer
4 // File: DCMTimingOffset_module.cc
5 //
6 // Generated at Mon Feb 11 09:01:32 2013 by Evan Niner using artmod
7 // from art v1_02_06.
8 ////////////////////////////////////////////////////////////////////////
9 
10 //NOvASoft includes
11 #include "Geometry/Geometry.h"
12 #include "RecoBase/CellHit.h"
13 #include "RecoBase/Cluster.h"
14 #include "RecoBase/Track.h"
15 #include "Utilities/func/MathUtil.h"
20 #include "NovaDAQConventions/DAQConventions.h"
22 #include "DAQChannelMap/DAQChannelMap.h"
23 
24 //framework includes
34 #include "Utilities/AssociationUtil.h"
35 #include "cetlib_except/exception.h"
37 
38 #include <stdint.h>
39 #include <sstream>
40 #include <map>
41 #include <iterator>
42 #include <iostream>
43 #include <fstream>
44 #include <stdio.h>
45 #include <string>
46 
47 //ROOT includes
48 #include "TTree.h"
49 #include "TH2F.h"
50 #include "TH1F.h"
51 #include "TGraph.h"
52 #include "TGraphErrors.h"
53 #include "TMatrixD.h"
54 #include "TVectorD.h"
55 #include "TDecompSVD.h"
56 #include "TDecompLU.h"
57 #include "TDirectory.h"
58 
59 namespace calib {
60 
61  class DCMTimingOffset : public art::EDFilter {
62  public:
63  explicit DCMTimingOffset(fhicl::ParameterSet const & p);
64  virtual ~DCMTimingOffset();
65 
66  //void produce(art::Event & e) override;
67  virtual bool filter(art::Event& evt) {return false;}
68 
69  void reconfigure(const fhicl::ParameterSet& pset);
70  void beginJob();
71  void endJob();
72  virtual bool beginRun(art::Run& run);
73  virtual bool endSubRun(art::SubRun& srun);
74  virtual bool endRun(art::Run & r);
75  //virtual void respondToCloseOutputFiles(art::FileBlock const& fb);
76 
77 
78  private:
79 
80  std::string fDCMLabel; ///< Where to find dcm stats
81  int fNumDCM; ///< Number of DCMs calibration will be performed on
82  int fMinDCM; ///< Number of lowest DCM in calibration set
83  int fMaxDCM; ///< Number of highest DCM in calibration set
84  int fMinCount; ///< Minimum number of counts for a DCM pair
85  bool fInitialized = false; ///<Flag if initialization has happened
86  bool fOverrideRef; ///<override default of fixing calibration of dcm 6 with an offset of 0
87  int fRefDCMTop = 6; ///< DCM number to use as reference for top chain, default is 6
88  int fRefDCMSide = 7; ///< DCM number to use as reference for side chain, default is 7
89  float fRefOffsetTop = 0.0; ///<absolute offset, in ns, to fix refence DCM to for top chain, default is 0.0
90  float fRefOffsetSide = 0.0; ///<absolute offset, in ns, to fix refence DCM to for side chain, default is 0.0
91  bool fIsMC;
94 
96 
97 
98  //keep stats for entire group of runs
99  std::vector<std::vector<int> > countDCMOffsetArray;
100  std::vector<std::vector<float> > sumSquareDCMOffsetArray;
101  std::vector<std::vector<float> > sumDCMOffsetArray;
102  std::vector<std::vector<float> > meanDCMOffsetArray;
103  std::vector<std::vector<float> > weightDCMOffsetArray;
104  std::vector<float> absoluteDCMOffset;
105  std::vector<float> absoluteDCMErrors;
106  std::vector<int> countDCMUse;
107  uint32_t minTime;
108  uint32_t maxTime;
109 
110  //keep stats for an individual run
111  std::vector<std::vector<int> > countDCMOffsetArrayRun;
112  std::vector<std::vector<float> > sumSquareDCMOffsetArrayRun;
113  std::vector<std::vector<float> > sumDCMOffsetArrayRun;
114  std::vector<std::vector<float> > meanDCMOffsetArrayRun;
115  std::vector<std::vector<float> > weightDCMOffsetArrayRun;
116  std::vector<float> absoluteDCMOffsetRun;
117  std::vector<float> absoluteDCMErrorsRun;
118  std::vector<int> countDCMUseRun;
119  uint32_t minTimeRun;
120  uint32_t maxTimeRun;
121 
123 
124  std::map<int, std::vector<caldp::DCMOffset> > offsetMap;
125 
126  };
127 }
128 
129 namespace calib {
130 
131  //define speed of light in cm/ns
132  static const double sCmPerNsec = 29.9792458;
133 
135  // :
136  // Initialize member data here.
137  {
138  this->reconfigure(pset);
139 
140  produces< std::vector<caldp::DCMOffset>, art::InRun >();
141  produces< caldp::DCMSummary, art::InRun >();
142 
143  }
144 
145 //---------------------------------------------------------------------------
147  {
148  fDCMLabel = pset.get< std::string >("DCMLabel");
149  fMinSubrun = pset.get< int>("MinSubrun");
150  fSaveRunCalibration = pset.get< bool >("SaveRunCalibration");
151  fOverrideRef = false;
152  fMinDCM = 0;
153  fMaxDCM = 0;
154  fMinCount = 0;
155  fInitialized = false;
156  fPSetND = pset.get<fhicl::ParameterSet>("nd");
157  fPSetFD = pset.get<fhicl::ParameterSet>("fd");
158  fPSetTB = pset.get<fhicl::ParameterSet>("tb");
159 
160  }
161 
162  //---------------------------------------------------------------------------
164  {
165  // Clean up dynamic memory and other resources here.
166  }
167 
168  //---------------------------------------------------------------------------
170  {
171 
172 
173 
174 
175  //We want to be able to handle any time possible.
178  }
179 
180  //---------------------------------------------------------------------------
182  {
183  //for first run determine detector, initialize arrays, do not repeat
184  if(!fInitialized){
185  fInitialized = true;
187  fIsMC = (rh->GetDataType() == nova::dbi::kMCOnly);
188 
191 
192  //get fcl parameters based on detector
194 
195  fhicl::ParameterSet pset;
196  switch(geom->DetId()){
198  pset = fPSetND;
199  break;
201  pset = fPSetFD;
202  break;
204  pset = fPSetTB;
205  break;
206  default:
207  assert(0 && "Unknown detector");
208  }
209 
210  fOverrideRef = pset.get< bool >("OverrideRef");
211  fMinDCM = pset.get< int >("MinDCM");
212  fMaxDCM = pset.get< int >("MaxDCM");
213  fMinCount = pset.get< int >("MinCount");
214  if (fOverrideRef){
215  fRefDCMTop = pset.get< int >("RefDCMTop");
216  fRefDCMSide = pset.get< int >("RefDCMSide");
217  fRefOffsetTop = pset.get< float >("RefOffsetTop");
218  fRefOffsetSide = pset.get< float >("RefOffsetSide");
219  }
220 
221  //assert parameters make sense
222  assert(fMaxDCM >= fMinDCM);
223  assert((fRefDCMTop >= fMinDCM) && (fRefDCMTop <= fMaxDCM));
224  assert((fRefDCMSide >= fMinDCM) && (fRefDCMSide <= fMaxDCM));
225  fNumDCM = fMaxDCM - fMinDCM + 1;
226 
227 
230  sumDCMOffsetArray.resize(fNumDCM);
231  meanDCMOffsetArray.resize(fNumDCM);
233  absoluteDCMOffset.resize(fNumDCM);
234  absoluteDCMErrors.resize(fNumDCM);
235  countDCMUse.resize(fNumDCM);
236 
237  for (int i=0; i<fNumDCM; ++i){
238  countDCMOffsetArray[i].resize(fNumDCM);
239  sumSquareDCMOffsetArray[i].resize(fNumDCM);
240  sumDCMOffsetArray[i].resize(fNumDCM);
241  meanDCMOffsetArray[i].resize(fNumDCM);
242  weightDCMOffsetArray[i].resize(fNumDCM);
243  }
244  }
245 
246  //clear arrays for the run to be safe
247  countDCMOffsetArrayRun.clear();
249  sumDCMOffsetArrayRun.clear();
250  meanDCMOffsetArrayRun.clear();
251  weightDCMOffsetArrayRun.clear();
252  absoluteDCMOffsetRun.clear();
253  absoluteDCMErrorsRun.clear();
254  countDCMUseRun.clear();
255 
256  //resize run arrays
264  countDCMUseRun.resize(fNumDCM);
265 
266  for (int i=0; i<fNumDCM; ++i){
267  countDCMOffsetArrayRun[i].resize(fNumDCM);
268  sumSquareDCMOffsetArrayRun[i].resize(fNumDCM);
269  sumDCMOffsetArrayRun[i].resize(fNumDCM);
270  meanDCMOffsetArrayRun[i].resize(fNumDCM);
271  weightDCMOffsetArrayRun[i].resize(fNumDCM);
272  }
273 
274  subrunCount = 0;
275 
276  return true;
277  }
278 
279  //---------------------------------------------------------------------------
281  {
282  //get out summary dataproduct
283  //get the cosmic tracks associated with each event
285  srun.getByLabel(fDCMLabel,dcmsum);
286 
287  subrunCount++;
288 
289  //keep running tally of relative dcm timing offsets for the run and overall
290  for (int j=0; j<fNumDCM; ++j){
291  int jabs = j + fMinDCM - 1;
292  countDCMUseRun[j] += dcmsum->DCMUseIndx(jabs);
293  countDCMUse[j] += dcmsum->DCMUseIndx(jabs);
294  for (int k=j+1; k<fNumDCM; ++k){
295  int kabs = k +fMinDCM - 1;
296  sumDCMOffsetArrayRun[j][k]+= dcmsum->DCMOffsetSumIndx(jabs,kabs);
297  sumDCMOffsetArray[j][k]+= dcmsum->DCMOffsetSumIndx(jabs,kabs);
298  sumSquareDCMOffsetArrayRun[j][k] += dcmsum->DCMOffsetSumSquareIndx(jabs,kabs);
299  sumSquareDCMOffsetArray[j][k] += dcmsum->DCMOffsetSumSquareIndx(jabs,kabs);
300  countDCMOffsetArrayRun[j][k] += dcmsum->DCMOffsetCountIndx(jabs,kabs);
301  countDCMOffsetArray[j][k] += dcmsum->DCMOffsetCountIndx(jabs,kabs);
302  }
303  }
304 
305  //adjust validity time
306  uint32_t min = dcmsum->StartTime();
307  uint32_t max = dcmsum->EndTime();
308  if (max>maxTime) maxTime = max;
309  if (min<minTime) minTime = min;
310  if (max>maxTimeRun) maxTimeRun = max;
311  if (min<minTimeRun) minTimeRun = min;
312 
313 
314  return false;
315  }
316 
317 
318 
319  //---------------------------------------------------------------------------
321  {
322  //perform timing calibration on a run basis if desired
323  if (!fSaveRunCalibration){
324  countDCMOffsetArrayRun.clear();
326  sumDCMOffsetArrayRun.clear();
327  meanDCMOffsetArrayRun.clear();
328  weightDCMOffsetArrayRun.clear();
329  absoluteDCMOffsetRun.clear();
330  absoluteDCMErrorsRun.clear();
331  countDCMUseRun.clear();
332  return false;
333  }
334 
335  if (subrunCount < fMinSubrun){
336  countDCMOffsetArrayRun.clear();
338  sumDCMOffsetArrayRun.clear();
339  meanDCMOffsetArrayRun.clear();
340  weightDCMOffsetArrayRun.clear();
341  absoluteDCMOffsetRun.clear();
342  absoluteDCMErrorsRun.clear();
343  countDCMUseRun.clear();
344  return false;
345  }
346  //Ptr to dcm timing offsets for the run
347  std::unique_ptr<std::vector<caldp::DCMOffset > >runOffsets(new std::vector<caldp::DCMOffset >);
348 
349  caldp::DCMSummary runsum;
350  runsum.SetStartTime(minTimeRun);
351  runsum.SetEndTime(maxTimeRun);
352 
353  //compute mean and standard deviation for each relative DCM pair in the matrix
354  for(int i=0; i<fNumDCM-1; ++i){
355  int iabs = i + fMinDCM - 1;
356  runsum.SetDCMUseIndx(iabs, countDCMUseRun[i]);
357  for(int j=i+1; j<fNumDCM; ++j){
358  int jabs = j + fMinDCM - 1;
359  float sum = sumDCMOffsetArrayRun[i][j];
361  float sumSqr = sumSquareDCMOffsetArrayRun[i][j];
362  runsum.SetDCMOffsetSumIndx(iabs,jabs,sum);
363  runsum.SetDCMOffsetSumSquareIndx(iabs,jabs,sumSqr);
364  runsum.SetDCMOffsetCountIndx(iabs,jabs,count);
365  if (count >= fMinCount){
366  meanDCMOffsetArrayRun[i][j] = sum/count;
367  float sigma = sqrt((count*sumSqr - sum*sum)/(count*(count-1.0)));
368  weightDCMOffsetArrayRun[i][j] = 1.0/(sigma*sigma);
369  }
370  else {
371  meanDCMOffsetArrayRun[i][j] = 0;
373  }
374  }
375  }
376 
377  const int num= fNumDCM;
378  int refIndxTop = fRefDCMTop - fMinDCM;
379  int refIndxSide = fRefDCMSide - fMinDCM;
380  TMatrixD A;
381  A.ResizeTo(fNumDCM, fNumDCM);
382  double y[num];
383  for (int x=0; x<num; ++x) y[x]=0.;
384 
385  //see docdb-XXXX, solve a system of linear equations that corresponds to minimizing the chhi squared for the absolute offsets. Fill matrix in this step
386  for(int i=0; i<fNumDCM; ++i){
387  double temp[num];
388  for (int x=0; x<num; ++x) temp[x]=0.;
389  for(int j=0; j<fNumDCM; ++j){
390  if (i==j) continue;
391  if (j<i){
392  temp[i] += weightDCMOffsetArrayRun[j][i]*2.0;
393  temp[j] += -weightDCMOffsetArrayRun[j][i]*2.0;
395  }
396  if (j>i){
397  temp[i] += weightDCMOffsetArrayRun[i][j]*2.0;
398  temp[j] += -weightDCMOffsetArrayRun[i][j]*2.0;
400  }
401  }
402 
403  if (i == refIndxTop){
404  for (int j=0; j<fNumDCM; ++j) temp[j] = 0.0;
405  temp[refIndxTop]=1.0;
406  }
407  if (i == refIndxSide){
408  for (int j=0; j<fNumDCM; ++j) temp[j] = 0.0;
409  temp[refIndxSide]=1.0;
410  }
411 
412  y[i]*=-1.0;
413  TVectorD tempvec;
414  tempvec.Use(fNumDCM,temp);
415  TMatrixDRow(A,i) = tempvec;
416  }
417 
418  y[refIndxTop] = fRefOffsetTop;
419  y[refIndxSide] = fRefOffsetSide;
420  TVectorD yvec;
421  yvec.Use(fNumDCM,y);
422 
423  //calculate uncertainties. Errors are highly correlated, this is an over estimate from the diagonal elements
424  for(int i=0; i<fNumDCM; ++i){
425  for(int j=0; j<fNumDCM; ++j){
426  if (i==j) continue;
428  else if (j>i) absoluteDCMErrorsRun[i] += 2.0*weightDCMOffsetArrayRun[i][j];
429  }
431  }
432  absoluteDCMErrorsRun[refIndxTop] = 0.0;
433  absoluteDCMErrorsRun[refIndxSide] = 0.0;
434 
435  //solve for absolute offsets
436  TDecompLU lu(A);
437  bool inv = true;
438  try{
439  lu.Decompose();
440  }
441  catch(...){
442  std::cout<<"WARNING: Could not perform matrix inversion to calculate absolute DCM offsets for run: "<<run.run()<<". The matrix inversion was performed over the range: "<<fMinDCM<<" to "<<fMaxDCM<<". If this run had missing diblocks in this range the matrix inversion will be invalid. If this is not desired please adjust the DCM range in the FCL file."<<std::endl;
443  inv = false;
444  }
445  if (inv) lu.Invert(A);
446  if (inv) yvec *= A;
447 
448  for (int i=0; i<fNumDCM; ++i){
449  if (!inv) break;
450  caldp::DCMOffset dcmOffset;
451  dcmOffset.SetDCM(i+fMinDCM);
452  dcmOffset.SetOffset(yvec[i]);
453  dcmOffset.SetSigma(absoluteDCMErrorsRun[i]);
454  dcmOffset.SetStats(countDCMUseRun[i]);
455  dcmOffset.SetStartTime(minTimeRun);
456  dcmOffset.SetEndTime(maxTimeRun);
457  runOffsets->push_back(dcmOffset);
458  offsetMap[run.run()].push_back(dcmOffset);
459  }
460 
461  countDCMOffsetArrayRun.clear();
463  sumDCMOffsetArrayRun.clear();
464  meanDCMOffsetArrayRun.clear();
465  weightDCMOffsetArrayRun.clear();
466  absoluteDCMOffsetRun.clear();
467  absoluteDCMErrorsRun.clear();
468  countDCMUseRun.clear();
469 
470  std::unique_ptr<caldp::DCMSummary> dsum(new caldp::DCMSummary(runsum));
471  run.put(std::move(dsum));
472 
473  run.put(std::move(runOffsets));
474 
475  return true;
476  }
477 
478  //---------------------------------------------------------------------------
480  {
481 std::cout<<"I am cleaning up"<<std::endl;
482  //figure out dcms per diblock from geometry
485  const int NDCMDB = cmap->getNumberOfDCM(daqchannelmap::AFULLAFULL_DIBLOCK);
486 
487  //make summary plots of calibration performance
489 
490  const int runs = offsetMap.size();
491  const int num= fNumDCM;
492  TGraphErrors* grAbsDCM[num];
493  TGraphErrors* grAbsRuns[runs];
494  TGraph* grStatsDCM[num];
495  TGraph* grStatsRuns[runs];
496  TGraphErrors* grAbs;
497  TGraphErrors* grAbsDiff;
498  TGraph* grStats;
499  TGraphErrors* grRelDiff;
500  TH2F* fHistoRelMap;
501  TH1F* fHistoRelDiff;
502  TH2F* fHistoRelDiffMap;
503 
504  double runDCM[num][runs];
505  double runDCMOffset[num][runs];
506  double runDCMOffsetE[num][runs];
507  double runDCMStats[num][runs];
508 
509  grRelDiff = tfs->make<TGraphErrors>();
510  gDirectory->Append(grRelDiff);
511  grRelDiff->SetName("grRelDiff");
512  grRelDiff->SetTitle("Relative offset drift vs diblock gap");
513  grRelDiff->Set(55);
514 
515  grAbsDiff = tfs->make<TGraphErrors>();
516  gDirectory->Append(grAbsDiff);
517 
518 
519  grAbs = tfs->make<TGraphErrors>();
520  gDirectory->Append(grAbs);
521  grStats = tfs->make<TGraph>();
522  gDirectory->Append(grStats);
523  fHistoRelMap = tfs->make<TH2F>("fHistoRelMap","Relative offsets (ns), all runs, DCM Y- DCM X;dcm X;dcm Y",
524  fNumDCM,(double)fMinDCM-0.5,(double)fMaxDCM+0.5,fNumDCM,(double)fMinDCM-0.5,(double)fMaxDCM+0.5);
525  fHistoRelDiffMap = tfs->make<TH2F>("fHistoRelDiffMap","Relative offsets differences (ns), all runs, DCM Y- DCM X - expected;dcm X;dcm Y",
526  fNumDCM,(double)fMinDCM-0.5,(double)fMaxDCM+0.5,fNumDCM,(double)fMinDCM-0.5,(double)fMaxDCM+0.5);
527  fHistoRelDiff = tfs->make<TH1F>("fHistoRelDiff","DCM i - DCM j - expected;diff (ns);count",201,-100.5,100.5);
528 
529  for (int i=0; i< runs; ++i){
530  grAbsRuns[i] = tfs->make<TGraphErrors>();
531  gDirectory->Append(grAbsRuns[i]);
532  grStatsRuns[i] = tfs->make<TGraph>();
533  gDirectory->Append(grStatsRuns[i]);
534  }
535  for (int i=0; i< num; ++i){
536  grAbsDCM[i] = tfs->make<TGraphErrors>();
537  gDirectory->Append(grAbsDCM[i]);
538  grStatsDCM[i] = tfs->make<TGraph>();
539  gDirectory->Append(grStatsDCM[i]);
540  }
541 
542  //expected pattern of DCM offsets within each FD diblock based on measured cable delays. This is used to plot accuracy of calibration. If
543  //the expected offsets change this should be updated.
544  double dcmoffsetExpectedFDDiblock[12] = {-203.125,-164.0625,-125.0,-78.125,-39.0625,0.0,0.0,-39.0625,-78.125,-125.0,-164.0625,-203.125};
545 
546  //expected pattern of DCM offsets for all ND DCMs based on measured cable delays. There is not a repeatable diblock pattern so all
547  //DCMs are listed. This is used to plot accuracy of calibration. If the expected offsets change this should be updated.
548  double dcmoffsetExpectedND[14] = {-101.5625,-70.3125,-30.25,0.0,-101.5625,-70.3125,-30.25,0.0,-101.5625,-70.3125,-30.25,0.0,-101.5625,-62.5};
549 
550  //initialize bins of the relative offset map so that empty bins can be distinguished from 0
551  for(int i=0;i<fNumDCM;++i){
552  for(int j=0;j<fNumDCM;++j){
553  fHistoRelMap->SetBinContent(j+1,i+1,-5000.0);
554  }
555  }
556 
557  //calculate mean and standard deviation for all relative offset pairs
558  for(int i=0; i<fNumDCM-1; ++i){
559  for(int j=i+1; j<fNumDCM; ++j){
560  float sum = sumDCMOffsetArray[i][j];
561  int count = countDCMOffsetArray[i][j];
562  float sumSqr = sumSquareDCMOffsetArray[i][j];
563  if (count >= fMinCount){
564  meanDCMOffsetArray[i][j] = sum/count;
565  float sigma = sqrt((count*sumSqr - sum*sum)/(count*(count-1.0)));
566  weightDCMOffsetArray[i][j] = 1.0/(sigma*sigma);
567  fHistoRelMap->SetBinContent(j+1,i+1,sum/count);
568  double expected = 0.0;
569  if (geom->DetId() == novadaq::cnv::kFARDET) expected = dcmoffsetExpectedFDDiblock[i%12]-dcmoffsetExpectedFDDiblock[j%12];
570  else if (geom->DetId() == novadaq::cnv::kNEARDET) expected = dcmoffsetExpectedND[i]-dcmoffsetExpectedND[j];
571  fHistoRelDiffMap->SetBinContent(j+1,i+1,(sum/count)-expected);
572  fHistoRelDiff->Fill((sum/count)-expected);
573  }
574  else {
575  meanDCMOffsetArray[i][j] = 0;
576  weightDCMOffsetArray[i][j] = 0;
577  }
578  }
579  }
580 
581 
582  int refIndxTop = fRefDCMTop - fMinDCM;
583  int refIndxSide = fRefDCMSide - fMinDCM;
584  TMatrixD A;
585  A.ResizeTo(fNumDCM, fNumDCM);
586  double y[num];
587  for (int x=0; x<num; ++x) y[x]=0.;
588 
589  //see docdb-XXXX, solve a system of linear equations that corresponds to minimizing the chhi squared for the absolute offsets. Fill matrix in this step
590  for(int i=0; i<fNumDCM; ++i){
591  double temp[num];
592  for (int x=0; x<num; ++x) temp[x]=0.;
593  for(int j=0; j<fNumDCM; ++j){
594  if (i==j) continue;
595  if (j<i){
596  temp[i] += weightDCMOffsetArray[j][i]*2.0;
597  temp[j] += -weightDCMOffsetArray[j][i]*2.0;
598  y[i] += meanDCMOffsetArray[j][i]*weightDCMOffsetArray[j][i]*2.0;
599  }
600  if (j>i){
601  temp[i] += weightDCMOffsetArray[i][j]*2.0;
602  temp[j] += -weightDCMOffsetArray[i][j]*2.0;
603  y[i] += -meanDCMOffsetArray[i][j]*weightDCMOffsetArray[i][j]*2.0;
604  }
605  }
606 
607  if (i == refIndxTop){
608  for (int j=0; j<fNumDCM; ++j) temp[j] = 0.0;
609  temp[refIndxTop]=1.0;
610  }
611  if (i == refIndxSide){
612  for (int j=0; j<fNumDCM; ++j) temp[j] = 0.0;
613  temp[refIndxSide]=1.0;
614  }
615 
616  y[i]*=-1.0;
617  TVectorD tempvec;
618  tempvec.Use(fNumDCM,temp);
619  TMatrixDRow(A,i) = tempvec;
620  }
621 
622  y[refIndxTop] = fRefOffsetTop;
623  y[refIndxSide] = fRefOffsetSide;
624  TVectorD yvec;
625  yvec.Use(fNumDCM,y);
626 
627  //calculate uncertainties
628  for(int i=0; i<fNumDCM; ++i){
629  for(int j=0; j<fNumDCM; ++j){
630  if (i==j) continue;
631  if (j<i) absoluteDCMErrors[i] += 2.0*weightDCMOffsetArray[j][i];
632  else if (j>i) absoluteDCMErrors[i] += 2.0*weightDCMOffsetArray[i][j];
633  }
635  }
636  absoluteDCMErrors[refIndxTop] = 0.0;
637  absoluteDCMErrors[refIndxSide] = 0.0;
638 
639  //solve for absolute offsets
640  TDecompLU lu(A);
641  try{
642  lu.Decompose();
643  }
644  catch(...){
645  std::cout<<"ERROR: Could not perform matrix inversion to calculate absolute DCM offsets for the processed runs. The matrix inversion was performed over the range: "<<fMinDCM<<" to "<<fMaxDCM<<". If this set of runs had missing diblocks in this range the matrix inversion will be invalid. If this is not desired please adjust the DCM range in the FCL file. The program will exit now."<<std::endl;
646  return;
647  }
648  lu.Invert(A);
649  yvec *= A;
650 
651  //write results to csv file
652  const TString mcStr = fIsMC ? "_mc" : "";
653  FILE* fConstsCSV = fopen("calib_abs_dcm_delay_consts"+mcStr+".csv","w");
654 
655  fprintf(fConstsCSV, "#dcmoffset,sigma,stats,channel,tv\n");
656 
657  //there is a drift in the dcm offsets in diblocks 7-14. The cause is unknown for now offsets of each diblock are adjusted by hand to put them in
658  //the right place. This should not be necessary in the future when the cause in known and can be removed
659  double dbcorrectionFD[14] = {0.0,
660  0.0,
661  0.0,
662  0.0,
663  0.0,
664  0.0,
665  4.06,
666  3.88,
667  3.82,
668  4.24,
669  6.31,
670  7.29,
671  10.23,
672  11.47};
673 
674  //same style correction for the ND
675  double dbcorrectionND[14] = {15.55,
676  8.2,
677  0.0,
678  0.0};
679 
680  //Multipy by -1 since the correction facto is the opposite sign of the offset
681  //store channel as diblock*1000+dcm
682  //make CVS file with constants to put in DB
683  for (int i=0; i<fNumDCM; ++i){
684  int db = ((i+fMinDCM)/NDCMDB)+1;
685  int dcm = (i+fMinDCM)%NDCMDB;
686  if ((i+fMinDCM)%NDCMDB == 0){
687  db--;
688  dcm = NDCMDB;
689  }
690  if (geom->DetId() == novadaq::cnv::kFARDET) absoluteDCMOffset[i] = yvec[i] + dbcorrectionFD[db-1];
691  else if (geom->DetId() == novadaq::cnv::kNEARDET) absoluteDCMOffset[i] = yvec[i] + dbcorrectionND[db-1];
692  int chan = db*1000 + dcm;
693  fprintf(fConstsCSV, "%g,%g,%d,%d,%d\n",
695  }
696 
697  fclose(fConstsCSV);
698 
699  //save summary calibration histograms on run-by-run basis
700  if (fSaveRunCalibration){
701  int runcount= 0;
702  for (std::map<int, std::vector<caldp::DCMOffset> >::iterator itr = offsetMap.begin(); itr != offsetMap.end(); ++itr){
703  int rnum = itr->first;
704  std::vector<caldp::DCMOffset> offsets = itr->second;
705 
706  TString dummy;
707  std::ostringstream convert;
708  convert << rnum;
709  dummy = convert.str();
710  grAbsRuns[runcount]->SetName("grAbsOffset_run_"+dummy);
711  grAbsRuns[runcount]->SetTitle("Absolute DCM timing offsets (ns) for run: "+dummy);
712  grAbsRuns[runcount]->Set(offsets.size());
713  grStatsRuns[runcount]->SetName("grStatsOffset_run_"+dummy);
714  grStatsRuns[runcount]->SetTitle("Number of track segments used by dcm for run: "+dummy);
715  grStatsRuns[runcount]->Set(offsets.size());
716  for (unsigned int i=0; i<offsets.size(); ++i){
717  int tdcm = offsets[i].DCM();
718  double toffset = offsets[i].Offset();
719  double tsigma = offsets[i].Sigma();
720  int tstats = offsets[i].Stats();
721  grStatsRuns[runcount]->SetPoint(i,tdcm,tstats);
722  grAbsRuns[runcount]->SetPoint(i,tdcm,toffset);
723  grAbsRuns[runcount]->SetPointError(i,0.0,tsigma);
724  runDCM[tdcm-fMinDCM][runcount]=rnum;
725  runDCMOffset[tdcm-fMinDCM][runcount]=toffset;
726  runDCMOffsetE[tdcm-fMinDCM][runcount]=tsigma;
727  runDCMStats[tdcm-fMinDCM][runcount]=tstats;
728  }
729  runcount++;
730  }
731 
732  for(int i=0; i<num; ++i){
733  int tdcm = i + fMinDCM;
734  TString dummy;
735  std::ostringstream convert;
736  convert << tdcm;
737  dummy = convert.str();
738  grAbsDCM[i]->SetName("grAbsOffset_dcm_"+dummy);
739  grAbsDCM[i]->SetTitle("Absolute DCM timing offsets (ns) across all run for dcm: "+dummy);
740  grAbsDCM[i]->Set(runs);
741  grStatsDCM[i]->SetName("grStatsOffset_dcm_"+dummy);
742  grStatsDCM[i]->SetTitle("Number of track segments used for all runs by dcm: "+dummy);
743  grStatsDCM[i]->Set(runs);
744  for(int j=0; j<runs; ++j){
745  grAbsDCM[i]->SetPoint(j,runDCM[i][j],runDCMOffset[i][j]);
746  grAbsDCM[i]->SetPointError(j,0.0,runDCMOffsetE[i][j]);
747  grStatsDCM[i]->SetPoint(j,runDCM[i][j],runDCMStats[i][j]);
748  }
749  }
750  }
751 
752  grAbs->SetName("grAbsOffset_total");
753  grAbs->SetTitle("Absolute DCM timing offsets (ns), all runs");
754  grAbs->Set(num);
755  grStats->SetName("grStatsOffset_total");
756  grStats->SetTitle("Number of track segments used by dcm, all runs");
757  grStats->Set(num);
758  grAbsDiff->SetName("grAbsDiff_total");
759  grAbsDiff->SetTitle("Absolute DCM timing offset - expected value (ns), all runs");
760  grAbsDiff->Set(num);
761 
762  for(int i=0; i< fNumDCM; ++i){
763  if (geom->DetId() == novadaq::cnv::kFARDET) grAbsDiff->SetPoint(i,i+fMinDCM,absoluteDCMOffset[i]-dcmoffsetExpectedFDDiblock[i%12]);
764  else if (geom->DetId() == novadaq::cnv::kNEARDET) grAbsDiff->SetPoint(i,i+fMinDCM,absoluteDCMOffset[i]-dcmoffsetExpectedND[i]);
765  grAbsDiff->SetPointError(i,0.0,absoluteDCMErrors[i]);
766  grAbs->SetPoint(i,i+fMinDCM,absoluteDCMOffset[i]);
767  grAbs->SetPointError(i,0.0,absoluteDCMErrors[i]);
768  grStats->SetPoint(i,i+fMinDCM,countDCMUse[i]);
769  }
770 
771 
772  }
773 
774 } //end namespace
775 
776 //-----------------------------------------------------------------------------
777 
DCMTimingOffset(fhicl::ParameterSet const &p)
void SetStartTime(uint32_t aStart)
Definition: DCMSummary.h:40
std::vector< std::vector< float > > meanDCMOffsetArray
virtual unsigned int getNumberOfDCM(DiBlock_TYPE dbt) const =0
How many DCMs does this diblock have?
TH2 * rh
Definition: drawXsec.C:5
std::vector< std::vector< float > > weightDCMOffsetArrayRun
void reconfigure(const fhicl::ParameterSet &pset)
virtual bool endSubRun(art::SubRun &srun)
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::string fDCMLabel
Where to find dcm stats.
int fRefDCMTop
DCM number to use as reference for top chain, default is 6.
void SetDCMUseIndx(int i, int val)
Definition: DCMSummary.h:50
std::vector< float > absoluteDCMOffset
virtual bool endRun(art::Run &r)
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
RunNumber_t run() const
Definition: Run.h:47
int fMinCount
Minimum number of counts for a DCM pair.
void SetDCMOffsetSumIndx(int i, int j, float val)
Definition: DCMSummary.h:48
int fRefDCMSide
DCM number to use as reference for side chain, default is 7.
void SetStartTime(uint32_t aStart)
Set start time offset is valid for.
Definition: DCMOffset.h:43
Definition: Run.h:31
std::vector< float > absoluteDCMOffsetRun
void SetOffset(float aOffset)
Set the timing offset (ns)
Definition: DCMOffset.h:37
void convert(std::string dir="cc_numu/C12")
Definition: convert.C:107
virtual bool beginRun(art::Run &run)
void SetDCMOffsetSumSquareIndx(int i, int j, float val)
Definition: DCMSummary.h:49
fclose(fg1)
Far Detector at Ash River, MN.
CDPStorage service.
int fMinDCM
Number of lowest DCM in calibration set.
void SetEndTime(uint32_t aEnd)
Set end time offset is valid for.
Definition: DCMOffset.h:45
static DAQChannelMap * getInstance(int detID)
float fRefOffsetTop
absolute offset, in ns, to fix refence DCM to for top chain, default is 0.0
std::vector< std::vector< float > > weightDCMOffsetArray
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
std::vector< float > absoluteDCMErrorsRun
Near Detector in the NuMI cavern.
void SetDCM(int aDCM)
Set DCM.
Definition: DCMOffset.h:35
std::vector< float > absoluteDCMErrors
std::vector< std::vector< float > > sumDCMOffsetArray
int fMaxDCM
Number of highest DCM in calibration set.
const double j
Definition: BetheBloch.cxx:29
Identifier for diblocks using a 32/32 configuration.
float fRefOffsetSide
absolute offset, in ns, to fix refence DCM to for side chain, default is 0.0
std::vector< std::vector< int > > countDCMOffsetArrayRun
std::vector< int > countDCMUseRun
double sigma(TH1F *hist, double percentile)
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
std::vector< std::vector< float > > meanDCMOffsetArrayRun
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
float DCMOffsetSumSquareIndx(int i, int j) const
Definition: DCMSummary.h:37
static const double A
Definition: Units.h:82
T * make(ARGS...args) const
static const double sCmPerNsec
bool fOverrideRef
override default of fixing calibration of dcm 6 with an offset of 0
void SetStats(int aStats)
Set number of track segments for dcm in analysis, quality measure.
Definition: DCMOffset.h:41
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
int num
Definition: f2_nu.C:119
std::vector< std::vector< float > > sumDCMOffsetArrayRun
void geom(int which=0)
Definition: geom.C:163
int DCMOffsetCountIndx(int i, int j) const
Definition: DCMSummary.h:35
bool fInitialized
Flag if initialization has happened.
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
void SetDCMOffsetCountIndx(int i, int j, int val)
Definition: DCMSummary.h:47
virtual bool filter(art::Event &evt)
cmap::CMap class source code
Definition: CMap.cxx:17
std::map< int, std::vector< caldp::DCMOffset > > offsetMap
float DCMOffsetSumIndx(int i, int j) const
Definition: DCMSummary.h:36
int fNumDCM
Number of DCMs calibration will be performed on.
std::vector< std::vector< float > > sumSquareDCMOffsetArrayRun
void SetEndTime(uint32_t aEnd)
Set end time offset is valid for.
Definition: DCMSummary.h:42
nova::dbi::DataType GetDataType()
Definition: RunHistory.h:452
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
void SetSigma(float aSigma)
Set uncertaint on offset (ns)
Definition: DCMOffset.h:39
int DCMUseIndx(int i) const
Definition: DCMSummary.h:38
std::vector< std::vector< float > > sumSquareDCMOffsetArray
uint32_t StartTime() const
return start time offset is valid for
Definition: DCMSummary.h:27
list runs
Definition: shutoffs.py:47
Encapsulate the geometry of one entire detector (near, far, ndos)
uint32_t EndTime() const
return end time offset is valid for
Definition: DCMSummary.h:29
fvar< T > inv(const fvar< T > &x)
Definition: inv.hpp:12
std::vector< std::vector< int > > countDCMOffsetArray
enum BeamMode string