DriftCorrection_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 // Class: DriftCorrection
3 // Module Type: producer
4 // File: DriftCorrection_module.cc
5 //
6 // Generated at Mon Nov 19 15:49:27 2012 by Ruth Toner using artmod
7 // from art v1_02_06.
8 ////////////////////////////////////////////////////////////////////////
9 
10 //C++ and C:
11 #include <cmath>
12 #include <iostream>
13 #include <fstream>
14 #include <vector>
15 
16 //Framework:
25 #include "cetlib_except/exception.h"
28 
29 //NovaSoft:
31 #include "Geometry/Geometry.h"
32 #include "NovaDAQConventions/DAQConventions.h"
33 #include "RecoBase/CellHit.h"
36 #include "Calibration/CDPStorage.h"
37 
38 // ROOT includes
39 #include "TTree.h"
40 #include "TH1F.h"
41 #include "TFile.h"
42 #include "TH2F.h"
43 #include "TGraphAsymmErrors.h"
44 #include "TGraph.h"
45 #include "TVector.h"
46 
47 #include <iostream>
48 #include <fstream>
49 
50 
51 namespace caldp {
52  class DriftResponse;
53 }
54 
55 namespace calib {
56 
58 
59  public:
60  explicit DriftCorrection(fhicl::ParameterSet const & p);
61  virtual ~DriftCorrection();
62 
63  void produce(art::Event & e) override;
64  void reconfigure(fhicl::ParameterSet const & p);
65  void beginJob() override;
66  void beginRun(art::Run & r) override;
67  void respondToCloseOutputFiles(art::FileBlock const& fb) override;
68 
69 
70  private:
71 
72  struct respStruct
73  {
74  ///Mean response
75  double mean;
76  ///Error on mean response
77  double meanerr;
78  ///Number of hits in mean calculation
79  int nhits;
80  ///RMS of hits in mean calculation
81  int rms;
82  ///Run start time
83  std::uint32_t starttime;
84  ///Run end time
85  std::uint32_t endtime;
86  };
87 
88  //FCL options:
89  ///FCL input: Make diagnostic histograms?
91  ///FCL input: Label of the DriftResponse collection:
93  ///FCL input: Number of days to use in aggregation:
94  double fNDay;
95  ///FCL input: Automatically remove low-statistics runs?
97  ///FCL input: Run number at which to automatically cut off calibration
99  ///FCL input: Jump size at which to make new validity context
100  double fJumpSize;
101  ///FCL input: Write a CSV File?
103  ///FCL input: file name for output CSV File
105  ///FCL input: Size of sliding window determining run quality.
107  ///FCL input: Is this a normalization job?
109 
110  ///Holds maximum run number for this job
111  int fMaxRun;
112 
113  ///Graph holding all runs
114  TGraphAsymmErrors *tgraphall;
115  ///Graph holding all good quality runs
116  TGraphAsymmErrors *tgraphgood;
117  ///Graph holding all validity context information
118  TGraphAsymmErrors *tgraphbreak;
119 
120  ///Map of Channel->(run,response)
121  std::map<int, std::map<int, respStruct> > respMap;
122 
123  ///Current list of run numbers
124  std::vector<int> nRun;
125 
126  ///Run quality map: channel->(run,quality)
127  std::map<int, std::map<int,bool> > runQual;
128 
129  ///Map of validity context start points (run number) to Drift constants.
130  std::map<int,caldp::DriftResponse> BreakList;
131 
132  ///Function to determine location of validity context breaks and corresponding Drift constants.
133  void LookForBreaks(std::map<int, respStruct> &chanlist, int channum);
134 
135  struct graphStruct
136  {
137  ///Mean response
138  double mean;
139  ///Error on mean response
140  double meanerr;
141  ///End time
142  double time;
143  ///End time
144  double timeerr;
145  };
146 
147  ///Map of Run->(channel,response)
148  std::map<int,std::map<int,float> > mCorrvals;
149 
150  std::map<int,std::map<int, graphStruct> > mAllRun;
151  std::map<int,std::map<int, graphStruct> > mGoodRun;
152  std::map<int,std::map<int, graphStruct> > mBreakRun;
153  };
154 
155 
156  ///Find validity contexts and Drift constants in an individual channel, for all runs in this channel
157  ///yet to be assigned to a validity context.
158  void DriftCorrection::LookForBreaks(std::map<int, respStruct> &chanlist, int channum){
159 
160  //Clear out the run list:
161  nRun.clear();
162 
163  //Make a list of runs currently active for this channel
164  std::map<int, respStruct >::iterator runIter;
165  for (runIter = chanlist.begin(); runIter != chanlist.end(); ++runIter ) nRun.push_back(runIter->first);
166 
167  //Make sure that the runs are properly ordered
168  std::sort(nRun.begin(),nRun.end());
169 
170  //What is the final run on the list?
171  int runfinal = nRun[nRun.size()-1];
172 
173  int nruni; //Holds a run number.
174  double mean; //Hols a mean
175  double rms; //Holds an rms
176  double nhits; //Holds the number of hits
177  double meanerr; //Holds an error on a mean
178  std::uint32_t starttime; //Holds a start time
179  std::uint32_t endtime = 0; //Holds an end time
180  double residual; //Holds a residual for an individual run off a running mean
181 
182  //All active run numbers in this channel are "bad" by default
183  int nrun;
184  for (unsigned int irun=0; irun < nRun.size(); irun++) {
185  nrun = nRun[irun];
186  if (runQual[channum].find(nrun) == runQual[channum].end()) {
187  runQual[channum].insert(std::make_pair(nrun,false));
188  }
189  }
190 
191  //Vector to find median
192  std::vector<double> sortvector;
193 
194  //Find running median of most recent runs in list.
195  int irun=0;
196  for (int i=0; i<fQualWindow; i++){
197  irun = nRun.size() - fQualWindow + i; //Position in ordered run list
198  if (irun < 0) continue;
199  nruni = nRun[irun]; //The run number at list position irun
200  mean = chanlist[nruni].mean; //Response for run nruni
201  sortvector.push_back(mean); //Add the response to the running list
202  }
203 
204  //Get the median entry:
205  std::sort(sortvector.begin(),sortvector.end());
206  int mentry = (int)(((double)sortvector.size())/2.0);
207  double runningMedian = sortvector[mentry];
208 
209  //Loop back over running list to find good quality entries:
210  for (int i=0; i<fQualWindow; i++){
211  irun = nRun.size() - fQualWindow + i; //Position in ordered run list
212  if (irun < 0) continue;
213  nruni = nRun[irun]; //The run number at list position irun
214  mean = chanlist[nruni].mean; //Response for run nruni
215  meanerr = chanlist[nruni].meanerr; //Error on mean for run nruni
216  nhits = chanlist[nruni].nhits; //Number of hits in run nruni
217  starttime = chanlist[nruni].starttime; //Start time for run nruni
218  endtime = chanlist[nruni].endtime; //End run for nruni
219 
220  //The following run types are automatically "bad" quality:
221  if (mean<=0) { //Badly defined response
222  runQual[channum][nruni] = false;
223  continue;
224  } else if ((fabs(meanerr/mean)>0.03 || nhits<200 ) && fRemoveLowStat==true ) { //"Low stat" run (if opt to remove automatically)
225  runQual[channum][nruni] = false;
226  continue;
227  } else if (starttime>3.0e9||starttime<1.3e9||endtime<1.3e9){ //Runs with poorly defined times
228  runQual[channum][nruni] = false;
229  continue;
230  }
231 
232  //Calculate the difference between this run's response and the running median
233  residual = fabs((mean/runningMedian) - 1.0); //Deviation from median in range.
234  if (residual<fJumpSize) {
235  runQual[channum][nruni] = true;
236  //Note: run has to be within range of median at ANY point in sliding range to be declared "good"
237  }
238  }
239 
240  //Records for old combined values in validity context
241  double old_mean=0;
242  double old_meanerr=0;
243  double old_rms=0;
244  double old_nhits=0;
245 
246  //Records for recalculated combined values in validity context
247  double new_mean=0;
248  double new_rms=0;
249  double new_nhits=0;
250  double new_meanerr=0;
251 
252  //Values for the next run
253  double next_start=0;
254  double next_mean=0;
255  double next_err=0;
256  double next_run=-1;
257 
258  //Counter variables
259  double time_tot=0; //Total time in this validity context
260  double time_0=0; //Time at start of validty context
261  int run_0=-1; //Run at start of validity context
262  int run_lastinvc=-1; //Last run in validity context
263 
264  //Boolean: do we flush existing responses at this run automatically?
265  bool flushtherun=false;
266  bool atflushrun=false;
267 
268  //Look at all quality-assured runs (i.e., not final fQualWindow), unless told to flush at current run.
269 
270  int terminalirun; //Run list index (not run number) at which to stop
271 
272  if ( fFlushRunNum>=0 && runfinal == fFlushRunNum ) {
273  //Analyze up to final run in current list
274  flushtherun=true;
275  terminalirun=nRun.size();
276  } else {
277  //Analyze up to final run whose quality has been determined
278  //(i.e., not ones currently in initial quality sliding window)
279  terminalirun = nRun.size() - fQualWindow;
280  }
281 
282  //Now loop over all runs up until current terminal window
283  for (irun=0; irun < terminalirun; ++irun) {
284 
285  //This is the Flush run.
286  if (irun==(terminalirun-1)&&flushtherun) atflushrun=true;
287 
288  nrun = nRun[irun];//Current run number
289 
290  //If the run is bad, skip to the next one.
291  if (runQual[channum][nrun]) {
292 
293  //Current response values:
294  mean = chanlist[nrun].mean; //Mean for this run
295  meanerr = (mean>0) ? chanlist[nrun].meanerr : 0; //Error for this run
296  rms = chanlist[nruni].rms; //RMS at this run
297  nhits = chanlist[nrun].nhits; //Total hits for this run
298 
299  //Start and end times:
300  starttime = chanlist[nrun].starttime; //Start time of this run
301  endtime = chanlist[nrun].endtime; //End time of this run
302 
303  if (fDoDiagnosticPlot){
304  //-----------DIAGNOSTIC--------------//
305  //Add good runs to structure for adding to graph:
306  if (mGoodRun[channum].find(nrun) == mGoodRun[channum].end() ){
307  graphStruct graphAll;
308  graphAll.mean = mean;
309  graphAll.meanerr = meanerr;
310  graphAll.time = ((double)endtime + (double)starttime)/2.0;
311  graphAll.timeerr = (double)endtime - graphAll.time;
312  mGoodRun[channum].insert(std::make_pair(nrun,graphAll));
313  }
314  }
315 
316  //----------------------------------//
317 
318  //For all runs not at end of window:
319  if (irun<terminalirun-1) {
320 
321  //Find the next good run in range:
322  int igood=1;
323  nruni = nRun[irun+igood];//The very next run
324 
325  //Keep counting up until either find the next good quality run, or hit end
326  while ((!runQual[channum][nruni])&&(irun+igood<=terminalirun)) {
327  igood++;
328  nruni = nRun[irun+igood];
329  }
330 
331  //If next good run is in range, find its values:
332  if (irun+igood<terminalirun){
333  next_run = nRun[irun+igood];
334  next_start = chanlist[next_run].starttime;
335  next_mean = chanlist[next_run].mean;
336  next_err = chanlist[next_run].meanerr;
337  } else continue; //Nothing new to look at in this range
338  } else if (flushtherun) { //This is the last run, but it's the flush run
339  next_run = nrun;
340  next_start = endtime;
341  next_mean = mean;
342  next_err = meanerr;
343  } else continue;
344 
345 
346  //If this is the first run in the VC, set the "initial" conditions
347  if (time_tot==0) {
348  time_0 = starttime;
349  run_0 = nrun;
350  }
351 
352  //Add current hits to old total to make "new" total
353  new_nhits = old_nhits + nhits;
354 
355  //If there are hits to use, start calculating new values:
356  if (new_nhits!=0){
357 
358  //Combine together current mean response and old one to form "new" mean response
359  new_mean = (old_nhits*old_mean + nhits*mean)/(new_nhits);
360 
361  //Some useful expressions to calculate new RMS
362  double sqsum = rms*rms*(nhits-1) + nhits*mean*mean;
363  double old_sqsum = old_rms*old_rms*(old_nhits-1) + old_nhits*old_mean*old_mean;
364  double new_rms2 = (sqsum + old_sqsum - (new_nhits*new_mean*new_mean))/(new_nhits-1);
365 
366  //Update the RMS:
367  new_rms=sqrt(new_rms2);
368 
369  //Update the error on mean response:
370  new_meanerr = new_rms/sqrt(new_nhits);
371 
372  //Total elapsed time in this validity context:
373  time_tot = endtime - time_0;
374 
375  //How do the mean and time change between now and next run in list?
376  double nextjump = ( new_meanerr > 0 ) ? fabs(next_mean - new_mean)/new_mean : 0;
377  double nexterrfrac = ( next_mean > 0 ) ? (fabs(next_err/next_mean)) : 0;
378 
379  //Search for a new break based on A) Total time in run, B) Time between end of this run and start of next
380  //C) Fractional jump between current mean and next run's, D) Whether we're right at the "Flush Run" :
381  if (
382  time_tot>fNDay*86400.
383  ||(next_start-endtime)>fNDay*86400./2.
384  ||(nextjump>fJumpSize&&nexterrfrac<0.3)
385  ||(atflushrun==true)
386  ){
387 
388  //Add an entry to the Break List for this channel and run!
389  BreakList[run_0].AddChannelResponse(channum,new_mean,new_meanerr,new_rms,new_nhits);
390 
391  //Add Diagnostic entries:
392  if (fDoDiagnosticPlot){
393  //-----------DIAGNOSTIC--------------//
394  //Add entry to graphs for this break validity context:
395  if (mBreakRun[channum].find(run_0) == mBreakRun[channum].end() ){
396  graphStruct graphAll;
397  graphAll.mean = new_mean;
398  graphAll.meanerr = new_meanerr;
399  graphAll.time = ((double)endtime + (double)time_0)/2.0;
400  graphAll.timeerr = (double)endtime - graphAll.time;
401  mBreakRun[channum].insert(std::make_pair(run_0,graphAll));
402  }
403  //----------------------------------//
404  }
405 
406  //Reset all counters and records:
407  time_tot=0;
408  new_mean=0;
409  new_rms=0;
410  new_nhits=0;
411  new_meanerr=0;
412  old_mean=0;
413  old_rms=0;
414  old_nhits=0;
415 
416  //Record the last run in the validity context
417  if (nrun>run_lastinvc) run_lastinvc = nrun;
418 
419  }
420 
421  //The updated version will be the "old" version for the next run:
422  old_mean=new_mean;
423  old_rms=new_rms;
424  old_nhits=new_nhits;
425  old_meanerr=new_meanerr;
426 
427  }
428  } else if (atflushrun&&old_nhits>0){ //This is a bad run, but we need to flush anyway
429 
430  //Use what we have now to make a break:
431  BreakList[run_0].AddChannelResponse(channum,old_mean,old_meanerr,old_rms,old_nhits);
432 
433  //-----------DIAGNOSTIC--------------//
434  //Add good runs:
435  if (mBreakRun[channum].find(run_0) == mBreakRun[channum].end() ){
436  graphStruct graphAll;
437  graphAll.mean = old_mean;
438  graphAll.meanerr = old_meanerr;
439  graphAll.time = ((double)endtime + (double)time_0)/2.0;
440  graphAll.timeerr = (double)endtime - graphAll.time;
441  mBreakRun[channum].insert(std::make_pair(run_0,graphAll));
442  }
443  //----------------------------------//
444 
445 
446  time_tot=0;
447 
448  new_mean=0;
449  new_rms=0;
450  new_nhits=0;
451  new_meanerr=0;
452 
453  old_mean=0;
454  old_meanerr=0;
455  old_rms=0;
456  old_nhits=0;
457 
458  //Last run to be used in a validity context:
459  if (nrun>run_lastinvc) run_lastinvc = nrun;
460 
461  }
462  }
463 
464  std::vector<int> eraselist;
465  int run_n;
466  //Clear everything in map up to final run in last validity context.
467  //First, make a record of thigns that need to get erased:
468  for (runIter = chanlist.begin(); runIter != chanlist.end(); ++runIter ) {
469  run_n = runIter->first;
470  if (run_n <= run_lastinvc ) {
471  eraselist.push_back(run_n);
472  } else break;
473  }
474  //Then run back over and do the erasing:
475  for (unsigned int i=0; i<eraselist.size(); i++){
476  int run_n = eraselist[i];
477  chanlist.erase(run_n);
478  }
479 
480  }
481 
482 
483  void DriftCorrection::reconfigure(fhicl::ParameterSet const & p)
484  {
485  fDriftRespLabel = p.get< std::string >("DriftRespLabel"); //Label of DriftResponse module
486  fNDay = p.get< double >("NDay"); //Number of days to employ
487  fRemoveLowStat = p.get< bool >("RemoveLowStat"); //Number of days to employ
488  fFlushRunNum = p.get< int >("FlushRunNum"); //Run number at which to automatically flush runs
489  fDoDiagnosticPlot = p.get< bool >("DoDiagnosticPlot"); //Make file with Drift Constants and data vs time
490  fJumpSize = p.get< double >("JumpSize"); //Fractional jump needed to declare new correction
491  fDoCSVFile = p.get< bool >("DoCSVFile"); //Write the results to a CSV file
492  fCSVFile = p.get< std::string >("CSVFile"); //CSV File Name
493  fQualWindow = p.get< int >("QualWindow"); //Sliding quality window
494  fDoNormalize = p.get< bool >("DoNormalize"); //Make file with Drift Constants and data vs time
495  }
496 
497  DriftCorrection::DriftCorrection(fhicl::ParameterSet const & p)
498  {
499  reconfigure(p);
500 
501  //DriftCorrection output produced every Run
502  produces< std::vector<caldp::DriftResponse >, art::InRun >();
503 
504 
505  }
506 
507  DriftCorrection::~DriftCorrection()
508  {
509  }
510 
511  void DriftCorrection::produce(art::Event & e)
512  {
513  }
514 
516  {
517  fMaxRun=-1;
518  }
519 
520  void DriftCorrection::beginRun(art::Run & r)
521  {
522 
523  //Geometry service:
525  //Channel map service:
527 
528  //Collection of mean responses:
530  r.getByLabel(fDriftRespLabel, meancol);
531 
532  int runnum;
533  std::uint32_t timestart;
534  std::uint32_t timeend;
535  std::uint32_t channum;
536 
537  //Run over collection of mean responses:
538  for (size_t i = 0; i < meancol->size(); ++i){
539 
540  //Grab drift response:
541  const caldp::DriftResponse *driftmean = &meancol->at(i);
542 
543  //Grab map of channel to response for this run:
544  std::map< int, caldp::MeanStruct > chanMap = driftmean->ReturnMap();
545 
546  //Grab run/timing information for this entry:
547  runnum = driftmean->Run();
548  timestart = driftmean->StartTime();
549  timeend = driftmean->EndTime();
550 
551  //Poorly defined start or end time to run. Skip this run.
552  if (timestart>3.0e9||timestart<1.3e9||timeend<1.3e9) {
553  LOG_VERBATIM("DriftCorrection") << "Run has poorly defined start or end time. Skipping." ;
554  continue;
555  }
556  //Check that runs are in order
557  if (runnum>fMaxRun) {
558  fMaxRun=runnum;
559  } else throw cet::exception("Rethrow") << "ERROR: Runs have been improperly ordered. Aborting.\n";
560 
561  //Add a new entry to the output map
562  if (BreakList.find(runnum) == BreakList.end())
563  BreakList.insert(std::pair<int, caldp::DriftResponse >(runnum,caldp::DriftResponse(runnum,timestart,timeend)));
564 
565  //Iterate over channels make map of channels->(run,response):
566  std::map< int, caldp::MeanStruct >::iterator chanIter;
567  for (chanIter = chanMap.begin(); chanIter != chanMap.end(); ++chanIter) {
568 
569  //Channel number:
570  channum = chanIter->first;
571 
572  //Structure with response info
573  respStruct response;
574  response.mean = chanIter->second.Mean();
575  response.rms = chanIter->second.RMS();
576  response.meanerr = chanIter->second.MeanErr();
577  response.nhits = chanIter->second.NHits();
578  response.starttime = timestart;
579  response.endtime = timeend;
580 
581  //If haven't added empty map for this channel yet, do so:
582  if (respMap.find(channum) == respMap.end()) {
583  respMap.insert(std::pair<int, std::map<int, respStruct> >(channum,std::map<int, respStruct>()));
584  //Run quality for an individual channel and run:
585  runQual.insert(std::pair<int, std::map<int, bool> >(channum,std::map<int, bool>()));
586  }
587 
588  //Add unique entry for this run:
589  if (respMap[channum].find(runnum) == respMap[channum].end() ){
590  respMap[channum].insert(std::make_pair(runnum,response));
591  } else throw cet::exception("Rethrow") << "ERROR: run was already added to the map; runs not ordered correctly. Aborting.\n";
592 
593  //---------------DIAGNOSTIC----------------//
594  //Fill structure for later diagnostic plots:
595  graphStruct graphAll;
596  graphAll.mean = chanIter->second.Mean();
597  graphAll.meanerr = chanIter->second.MeanErr();
598  graphAll.time = ((double)timeend + (double)timestart)/2.0;
599  graphAll.timeerr = (double)timeend - graphAll.time;
600 
601  //Add entries to all maps:
602  if (mAllRun.find(channum) == mAllRun.end()){
603  mAllRun.insert(std::pair<int, std::map<int, graphStruct> >(channum,std::map<int, graphStruct>()));
604  mGoodRun.insert(std::pair<int, std::map<int, graphStruct> >(channum,std::map<int, graphStruct>()));
605  mBreakRun.insert(std::pair<int, std::map<int, graphStruct> >(channum,std::map<int, graphStruct>()));
606  }
607  //Add "total" run:
608  if (mAllRun[channum].find(runnum) == mAllRun[channum].end() ){
609  mAllRun[channum].insert(std::make_pair(runnum,graphAll));
610  } else throw cet::exception("Rethrow") << "ERROR: run was never added to the map; something strange has happened. Aborting.\n";
611  //------------------------------------------//
612 
613  //Start looking for breaks after 9 runs have been accumulated!
614  //Do this individually for each run, looking over runs in order
615  if ( (int)respMap[channum].size() > (fQualWindow-1) || runnum == fFlushRunNum )
616  LookForBreaks(respMap[channum],channum);
617 
618  }//End of channel loop
619 
620  }//End of meancol loop
621 
622  }
623 
624 
625  void DriftCorrection::respondToCloseOutputFiles(art::FileBlock const&)
626  {
627 
628  //Channel map:
630 
631  std::ofstream filecsv;
632  if (fDoCSVFile)
633  filecsv.open(fCSVFile.c_str());
634 
635  //Vector containing drift responses for later storage
636  std::vector<caldp::DriftResponse> driftvect;
637 
638  //Vector for use in normalization code:
639  std::map< int, std::vector< std::pair<double,int> > > meanrecord;
640 
641  //Iterator for break list
642  std::map<int,caldp::DriftResponse>::iterator brkIter;
643 
644  //Loop back over break list
645  for (brkIter = BreakList.begin(); brkIter != BreakList.end(); ++brkIter) {
646 
647  //Drift Response object:
648  caldp::DriftResponse driftresp = brkIter->second;
649 
650  //Grab map of channel to response for this run, just to print:
651  std::map< int, caldp::MeanStruct > chanMap = driftresp.ReturnMap();
652 
653  //Loop for debugging and records:
654  LOG_VERBATIM("DriftCorrection") << "Start=" << driftresp.StartTime() << " End=" << driftresp.EndTime() ;
655  std::map< int, caldp::MeanStruct >::iterator chanIter;
656  for (chanIter = chanMap.begin(); chanIter != chanMap.end(); ++chanIter) {
657  LOG_DEBUG("DriftCorrection") << " channel " << chanIter->first
658  << " : mean=" << chanIter->second.Mean()
659  << " : meanerr=" << chanIter->second.MeanErr()
660  << " : rms=" << chanIter->second.RMS()
661  << " : nhits=" << chanIter->second.NHits();
662 
663  //Print out the breaks which were made:
664  LOG_VERBATIM("DriftCorrection") << "Push back chan=" << chanIter->first << " with " << chanIter->second.Mean() << "/" << chanIter->second.NHits() ;
665 
666  //If requested, and not normalizing, turn output into a csv file as well:
667  if (fDoCSVFile&&!fDoNormalize)
668  filecsv << chanIter->second.Mean() << "," << chanIter->second.MeanErr() << "," << chanIter->first << "," << driftresp.StartTime() << "\n";
669 
670  //If we are normalizing, save the mean and nhits for this run for the channel
671  if (fDoNormalize)
672  meanrecord[chanIter->first].push_back(std::make_pair(chanIter->second.Mean(),chanIter->second.NHits()));
673 
674  }//end channel loop
675 
676  //Store drift response for later, if not normalizing:
677  if (!fDoNormalize) driftvect.push_back(driftresp);
678 
679  }//end break run loop
680 
681 
682  //If we are doing normalization, calculate the mean per channel:
683  if (fDoNormalize){
684 
685  //Response at time "zero" for normalization
686  caldp::DriftResponse normresp(0,0,0);
687 
688  int chanid=-1; //channel
689  double meansum,nsum,newmean; //iterators and means
690 
691  //Run back over channels and their stored response and nhit info:
692  std::map< int, std::vector< std::pair<double,int> > >::iterator meanIter;
693  for (meanIter = meanrecord.begin(); meanIter != meanrecord.end(); ++meanIter) {
694 
695  chanid = meanIter->first; //channel number
696  std::vector< std::pair<double,int> > meanlist = meanIter->second; //list of responses
697 
698  //Reset counters:
699  meansum=0;
700  nsum=0;
701 
702  //Loop over response values to calculate the mean for this channel:
703  for (unsigned int ival=0; ival<meanlist.size(); ++ival){
704  //Sum of mean and number of hits:
705  meansum += (meanlist[ival].first)*(meanlist[ival].second);
706  nsum += meanlist[ival].second;
707  }
708 
709  //Calculate the mean:
710  newmean = (nsum>0) ? meansum/nsum : -1;
711 
712  //Add this as a DriftResponse:
713  normresp.AddChannelResponse(chanid,newmean,-1,-1,nsum);
714 
715  //Print to CSV file if requested:
716  if (fDoCSVFile)
717  filecsv << newmean << "," << -1 << "," << chanid << "," << 0 << "\n";
718 
719  }//end channel loop
720 
721  //Save the "normalized" version to the drift vector
722  driftvect.push_back(normresp);
723 
724  }//end normalize if statement
725 
726  //Close the CSV file:
727  if (fDoCSVFile) filecsv.close();
728 
729  //-------------DIAGNOSTIC--------------//
730  //Make plots of the runs and breaks vs time, if requested
731  if (fDoDiagnosticPlot){
732  int det, diblock, dcm, feb;
733  daqchannelmap::dchan channum;
734 
735  //File Service:
737 
738  //Iterator for a map containing graph points:
739  std::map<int,std::map<int, graphStruct> >::iterator graphIter;
740 
741  //Loop over all channels to make plots of Data and Fit vs Time.
742  for (graphIter = mAllRun.begin(); graphIter != mAllRun.end(); ++graphIter) {
743 
744  //Channel number:
745  channum = graphIter->first;
746 
747  //Get "human" values using channel service:
748  det = cmap->Map()->getDetector(channum);
749  diblock = cmap->Map()->getDiBlock(channum);
750  dcm = cmap->Map()->getDCM(channum);
751  feb = cmap->Map()->getFEB(channum);
752 
753  //----Save ALL Runs to graph----------------------------//
754  TVector gall_y(mAllRun[channum].size());
755  TVector gall_yerr(mAllRun[channum].size());
756  TVector gall_x(mAllRun[channum].size());
757  TVector gall_xerr(mAllRun[channum].size());
758  int iv=0;
759  graphStruct gstruct;
760  std::map<int, graphStruct>::iterator runIter;
761  for (runIter = mAllRun[channum].begin(); runIter != mAllRun[channum].end(); ++runIter) {
762  gstruct=runIter->second;
763  gall_y(iv) = gstruct.mean;
764  gall_yerr(iv) = gstruct.meanerr;
765  gall_x(iv) = gstruct.time;
766  gall_xerr(iv) = gstruct.timeerr;
767  iv++;
768  }
769  //Make the TGraphAsymm here and write it out
770  tgraphall = tfs->make<TGraphAsymmErrors>(gall_x,gall_y,gall_xerr,gall_xerr,gall_yerr,gall_yerr);
771  tgraphall->SetName(Form("graph_all_%d_%d_%d_%d",det,diblock,dcm,feb));
772  tgraphall->Write();
773  //-------------------------------------------------------//
774 
775  //----Do the same for Good Quality Runs--------------------------//
776  //CHECK FOR CHANNEL'S EXISTENCE:
777  iv=0;
778  if (mGoodRun.find(channum) != mGoodRun.end()) {
779  TVector ggood_y(mGoodRun[channum].size());
780  TVector ggood_yerr(mGoodRun[channum].size());
781  TVector ggood_x(mGoodRun[channum].size());
782  TVector ggood_xerr(mGoodRun[channum].size());
783  for (runIter = mGoodRun[channum].begin(); runIter != mGoodRun[channum].end(); ++runIter) {
784  gstruct=runIter->second;
785  ggood_y(iv) = gstruct.mean;
786  ggood_yerr(iv) = gstruct.meanerr;
787  ggood_x(iv) = gstruct.time;
788  ggood_xerr(iv) = gstruct.timeerr;
789  iv++;
790  }
791  //Make the TGraphAsymm here and write it out
792  tgraphgood = tfs->make<TGraphAsymmErrors>(ggood_x,ggood_y,ggood_xerr,ggood_xerr,ggood_yerr,ggood_yerr);
793  tgraphgood->SetName(Form("graph_good_%d_%d_%d_%d",det,diblock,dcm,feb));
794  tgraphgood->Write();
795  }
796  //-------------------------------------------------------//
797 
798  //----Do the same for Break Runs-------------------------//
799  //CHECK FOR CHANNEL'S EXISTENCE
800  iv=0;
801  if (mBreakRun.find(channum) != mBreakRun.end()) {
802  TVector gbreak_y(mBreakRun[channum].size());
803  TVector gbreak_yerr(mBreakRun[channum].size());
804  TVector gbreak_x(mBreakRun[channum].size());
805  TVector gbreak_xerr(mBreakRun[channum].size());
806  for (runIter = mBreakRun[channum].begin(); runIter != mBreakRun[channum].end(); ++runIter) {
807  gstruct=runIter->second;
808  gbreak_y(iv) = gstruct.mean;
809  gbreak_yerr(iv) = gstruct.meanerr;
810  gbreak_x(iv) = gstruct.time;
811  gbreak_xerr(iv) = gstruct.timeerr;
812  iv++;
813  }
814  //Make the TGraphAsymm here and write it out
815  tgraphbreak = tfs->make<TGraphAsymmErrors>(gbreak_x,gbreak_y,gbreak_xerr,gbreak_xerr,gbreak_yerr,gbreak_yerr);
816  tgraphbreak->SetName(Form("graph_break_%d_%d_%d_%d",det,diblock,dcm,feb));
817  tgraphbreak->Write();
818  }
819  //-------------------------------------------------------//
820  }//End graph iteration loop
821  }//--------END DIAGNOSTIC-------------//
822 
823 
824  //Write to file
826  LOG_VERBATIM("DriftCorrection") << "ADDING DRIFTVECT SIZE " << driftvect.size() ;
827  cdpstorage->DriftResponseToStore(driftvect);
828 
829 
830  }
831 
832 }
834 
double meanerr
Error on mean response.
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
int getDetector(uint32_t anychan) const
TGraphAsymmErrors * tgraphbreak
Graph holding all validity context information.
Structure to hold drift response in a single channel at a single time.
Definition: DriftResponse.h:36
diblock
print "ROW IS " print row
Definition: geo2elec.py:31
void AddChannelResponse(int fOffChan, double fMean, double fMeanErr, double fRMS, int fNHits)
Return subrun end time.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
bool fDoCSVFile
FCL input: Write a CSV File?
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
bool fDoNormalize
FCL input: Is this a normalization job?
std::map< int, std::map< int, graphStruct > > mGoodRun
TGraphAsymmErrors * tgraphgood
Graph holding all good quality runs.
const daqchannelmap::DAQChannelMap * Map() const
Definition: CMap.h:57
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
double meanerr
Error on mean response.
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
std::map< int, std::map< int, float > > mCorrvals
Map of Run->(channel,response)
DEFINE_ART_MODULE(TestTMapFile)
Definition: Run.h:31
bool fRemoveLowStat
FCL input: Automatically remove low-statistics runs?
std::map< int, std::map< int, bool > > runQual
Run quality map: channel->(run,quality)
int fQualWindow
FCL input: Size of sliding window determining run quality.
std::map< int, std::map< int, graphStruct > > mBreakRun
CDPStorage service.
std::string fDriftRespLabel
FCL input: Label of the DriftResponse collection:
if(dump)
void beginJob()
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::map< int, std::map< int, respStruct > > respMap
Map of Channel->(run,response)
TGraphAsymmErrors * tgraphall
Graph holding all runs.
std::vector< int > nRun
Current list of run numbers.
int fMaxRun
Holds maximum run number for this job.
Histograms used by attenuation calibration.
uint32_t EndTime() const
Return subrun start time.
Definition: DriftResponse.h:48
uint32_t StartTime() const
Return run number.
Definition: DriftResponse.h:47
std::map< int, MeanStruct > const & ReturnMap() const
Definition: DriftResponse.h:61
std::string fCSVFile
FCL input: file name for output CSV File.
std::uint32_t starttime
Run start time.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
void DriftResponseToStore(std::vector< caldp::DriftResponse > &drvec)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
int nhits
Number of hits in mean calculation.
TRandom3 r(0)
std::map< int, caldp::DriftResponse > BreakList
Map of validity context start points (run number) to Drift constants.
std::map< int, std::map< int, graphStruct > > mAllRun
dcm_id_t getDCM(dchan daqchan) const
Decode the dcm ID from a dchan.
#define LOG_VERBATIM(category)
int rms
RMS of hits in mean calculation.
Float_t e
Definition: plot.C:35
uint32_t dchan
< DAQ Channel Map Package
double fJumpSize
FCL input: Jump size at which to make new validity context.
double fNDay
FCL input: Number of days to use in aggregation:
int fFlushRunNum
FCL input: Run number at which to automatically cut off calibration.
diblock_t getDiBlock(dchan daqchan) const
Decode the diblock ID from a dchan.
Encapsulate the geometry of one entire detector (near, far, ndos)
bool fDoDiagnosticPlot
FCL input: Make diagnostic histograms?
ival
Definition: test.py:10
feb_t getFEB(dchan daqchan) const
Decode the feb id from a dchan.