23 #include "art_root_io/TFileService.h" 25 #include "cetlib_except/exception.h" 32 #include "NovaDAQConventions/DAQConventions.h" 43 #include "TGraphAsymmErrors.h" 121 std::map<int, std::map<int, respStruct> >
respMap;
133 void LookForBreaks(std::map<int, respStruct> &chanlist,
int channum);
150 std::map<int,std::map<int, graphStruct> >
mAllRun;
151 std::map<int,std::map<int, graphStruct> >
mGoodRun;
158 void DriftCorrection::LookForBreaks(std::map<int, respStruct> &chanlist,
int channum){
164 std::map<int, respStruct >::iterator runIter;
165 for (runIter = chanlist.begin(); runIter != chanlist.end(); ++runIter ) nRun.push_back(runIter->first);
168 std::sort(nRun.begin(),nRun.end());
171 int runfinal = nRun[nRun.size()-1];
184 for (
unsigned int irun=0; irun < nRun.size(); irun++) {
186 if (runQual[channum].find(nrun) == runQual[channum].end()) {
192 std::vector<double> sortvector;
196 for (
int i=0;
i<fQualWindow;
i++){
197 irun = nRun.size() - fQualWindow +
i;
198 if (irun < 0)
continue;
200 mean = chanlist[nruni].mean;
201 sortvector.push_back(mean);
205 std::sort(sortvector.begin(),sortvector.end());
206 int mentry = (
int)(((
double)sortvector.size())/2.0);
207 double runningMedian = sortvector[mentry];
210 for (
int i=0;
i<fQualWindow;
i++){
211 irun = nRun.size() - fQualWindow +
i;
212 if (irun < 0)
continue;
214 mean = chanlist[nruni].mean;
215 meanerr = chanlist[nruni].meanerr;
216 nhits = chanlist[nruni].nhits;
217 starttime = chanlist[nruni].starttime;
218 endtime = chanlist[nruni].endtime;
222 runQual[channum][nruni] =
false;
224 }
else if ((fabs(meanerr/mean)>0.03 || nhits<200 ) && fRemoveLowStat==
true ) {
225 runQual[channum][nruni] =
false;
227 }
else if (starttime>3.0e9||starttime<1.3e9||endtime<1.3e9){
228 runQual[channum][nruni] =
false;
233 residual = fabs((mean/runningMedian) - 1.0);
234 if (residual<fJumpSize) {
235 runQual[channum][nruni] =
true;
242 double old_meanerr=0;
250 double new_meanerr=0;
265 bool flushtherun=
false;
266 bool atflushrun=
false;
272 if ( fFlushRunNum>=0 && runfinal == fFlushRunNum ) {
275 terminalirun=nRun.size();
279 terminalirun = nRun.size() - fQualWindow;
283 for (irun=0; irun < terminalirun; ++irun) {
286 if (irun==(terminalirun-1)&&flushtherun) atflushrun=
true;
291 if (runQual[channum][nrun]) {
294 mean = chanlist[
nrun].mean;
295 meanerr = (mean>0) ? chanlist[nrun].meanerr : 0;
296 rms = chanlist[nruni].rms;
297 nhits = chanlist[
nrun].nhits;
300 starttime = chanlist[
nrun].starttime;
301 endtime = chanlist[
nrun].endtime;
303 if (fDoDiagnosticPlot){
306 if (mGoodRun[channum].find(nrun) == mGoodRun[channum].end() ){
311 graphAll.
timeerr = (double)endtime - graphAll.
time;
319 if (irun<terminalirun-1) {
323 nruni = nRun[irun+igood];
326 while ((!runQual[channum][nruni])&&(irun+igood<=terminalirun)) {
328 nruni = nRun[irun+igood];
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;
338 }
else if (flushtherun) {
353 new_nhits = old_nhits + nhits;
359 new_mean = (old_nhits*old_mean + nhits*
mean)/(new_nhits);
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);
367 new_rms=
sqrt(new_rms2);
370 new_meanerr = new_rms/
sqrt(new_nhits);
373 time_tot = endtime - time_0;
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;
382 time_tot>fNDay*86400.
383 ||(next_start-endtime)>fNDay*86400./2.
384 ||(nextjump>fJumpSize&&nexterrfrac<0.3)
389 BreakList[run_0].AddChannelResponse(channum,new_mean,new_meanerr,new_rms,new_nhits);
392 if (fDoDiagnosticPlot){
395 if (mBreakRun[channum].find(run_0) == mBreakRun[channum].
end() ){
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;
417 if (nrun>run_lastinvc) run_lastinvc =
nrun;
425 old_meanerr=new_meanerr;
428 }
else if (atflushrun&&old_nhits>0){
431 BreakList[run_0].AddChannelResponse(channum,old_mean,old_meanerr,old_rms,old_nhits);
435 if (mBreakRun[channum].find(run_0) == mBreakRun[channum].
end() ){
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;
459 if (nrun>run_lastinvc) run_lastinvc =
nrun;
464 std::vector<int> eraselist;
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);
475 for (
unsigned int i=0;
i<eraselist.size();
i++){
476 int run_n = eraselist[
i];
477 chanlist.erase(run_n);
486 fNDay = p.
get<
double >(
"NDay");
487 fRemoveLowStat = p.
get<
bool >(
"RemoveLowStat");
488 fFlushRunNum = p.
get<
int >(
"FlushRunNum");
489 fDoDiagnosticPlot = p.
get<
bool >(
"DoDiagnosticPlot");
490 fJumpSize = p.
get<
double >(
"JumpSize");
491 fDoCSVFile = p.
get<
bool >(
"DoCSVFile");
493 fQualWindow = p.
get<
int >(
"QualWindow");
494 fDoNormalize = p.
get<
bool >(
"DoNormalize");
503 produces< std::vector<caldp::DriftResponse >,
art::InRun >();
534 std::uint32_t timestart;
535 std::uint32_t timeend;
536 std::uint32_t channum;
539 for (
size_t i = 0;
i < meancol->size(); ++
i){
545 std::map< int, caldp::MeanStruct > chanMap = driftmean->
ReturnMap();
548 runnum = driftmean->
Run();
550 timeend = driftmean->
EndTime();
553 if (timestart>3.0e9||timestart<1.3e9||timeend<1.3e9) {
554 MF_LOG_VERBATIM(
"DriftCorrection") <<
"Run has poorly defined start or end time. Skipping." ;
560 }
else throw cet::exception(
"Rethrow") <<
"ERROR: Runs have been improperly ordered. Aborting.\n";
567 std::map< int, caldp::MeanStruct >::iterator chanIter;
568 for (chanIter = chanMap.begin(); chanIter != chanMap.end(); ++chanIter) {
571 channum = chanIter->first;
575 response.
mean = chanIter->second.Mean();
576 response.
rms = chanIter->second.RMS();
577 response.
meanerr = chanIter->second.MeanErr();
578 response.
nhits = chanIter->second.NHits();
584 respMap.insert(std::pair<
int, std::map<int, respStruct> >(channum,std::map<int, respStruct>()));
586 runQual.insert(std::pair<
int, std::map<int, bool> >(channum,std::map<int, bool>()));
592 }
else throw cet::exception(
"Rethrow") <<
"ERROR: run was already added to the map; runs not ordered correctly. Aborting.\n";
597 graphAll.
mean = chanIter->second.Mean();
598 graphAll.
meanerr = chanIter->second.MeanErr();
599 graphAll.
time = ((double)timeend + (
double)timestart)/2.0;
600 graphAll.
timeerr = (double)timeend - graphAll.
time;
604 mAllRun.insert(std::pair<
int, std::map<int, graphStruct> >(channum,std::map<int, graphStruct>()));
605 mGoodRun.insert(std::pair<
int, std::map<int, graphStruct> >(channum,std::map<int, graphStruct>()));
606 mBreakRun.insert(std::pair<
int, std::map<int, graphStruct> >(channum,std::map<int, graphStruct>()));
611 }
else throw cet::exception(
"Rethrow") <<
"ERROR: run was never added to the map; something strange has happened. Aborting.\n";
632 std::ofstream filecsv;
637 std::vector<caldp::DriftResponse> driftvect;
640 std::map< int, std::vector< std::pair<double,int> > > meanrecord;
643 std::map<int,caldp::DriftResponse>::iterator brkIter;
652 std::map< int, caldp::MeanStruct > chanMap = driftresp.
ReturnMap();
656 std::map< int, caldp::MeanStruct >::iterator chanIter;
657 for (chanIter = chanMap.begin(); chanIter != chanMap.end(); ++chanIter) {
658 MF_LOG_DEBUG(
"DriftCorrection") <<
" channel " << chanIter->first
659 <<
" : mean=" << chanIter->second.Mean()
660 <<
" : meanerr=" << chanIter->second.MeanErr()
661 <<
" : rms=" << chanIter->second.RMS()
662 <<
" : nhits=" << chanIter->second.NHits();
665 MF_LOG_VERBATIM(
"DriftCorrection") <<
"Push back chan=" << chanIter->first <<
" with " << chanIter->second.Mean() <<
"/" << chanIter->second.NHits() ;
669 filecsv << chanIter->second.Mean() <<
"," << chanIter->second.MeanErr() <<
"," << chanIter->first <<
"," << driftresp.
StartTime() <<
"\n";
673 meanrecord[chanIter->first].push_back(
std::make_pair(chanIter->second.Mean(),chanIter->second.NHits()));
690 double meansum,nsum,newmean;
693 std::map< int, std::vector< std::pair<double,int> > >::iterator meanIter;
694 for (meanIter = meanrecord.begin(); meanIter != meanrecord.end(); ++meanIter) {
696 chanid = meanIter->first;
697 std::vector< std::pair<double,int> > meanlist = meanIter->second;
704 for (
unsigned int ival=0;
ival<meanlist.size(); ++
ival){
707 nsum += meanlist[
ival].second;
711 newmean = (nsum>0) ? meansum/nsum : -1;
718 filecsv << newmean <<
"," << -1 <<
"," << chanid <<
"," << 0 <<
"\n";
723 driftvect.push_back(normresp);
740 std::map<int,std::map<int, graphStruct> >::iterator graphIter;
743 for (graphIter =
mAllRun.begin(); graphIter !=
mAllRun.end(); ++graphIter) {
746 channum = graphIter->first;
761 std::map<int, graphStruct>::iterator runIter;
762 for (runIter =
mAllRun[channum].begin(); runIter !=
mAllRun[channum].end(); ++runIter) {
763 gstruct=runIter->second;
764 gall_y(iv) = gstruct.
mean;
765 gall_yerr(iv) = gstruct.
meanerr;
766 gall_x(iv) = gstruct.
time;
767 gall_xerr(iv) = gstruct.
timeerr;
771 tgraphall = tfs->make<TGraphAsymmErrors>(gall_x,gall_y,gall_xerr,gall_xerr,gall_yerr,gall_yerr);
772 tgraphall->SetName(Form(
"graph_all_%d_%d_%d_%d",det,diblock,dcm,feb));
784 for (runIter =
mGoodRun[channum].begin(); runIter !=
mGoodRun[channum].end(); ++runIter) {
785 gstruct=runIter->second;
786 ggood_y(iv) = gstruct.
mean;
787 ggood_yerr(iv) = gstruct.
meanerr;
788 ggood_x(iv) = gstruct.
time;
789 ggood_xerr(iv) = gstruct.
timeerr;
793 tgraphgood = tfs->make<TGraphAsymmErrors>(ggood_x,ggood_y,ggood_xerr,ggood_xerr,ggood_yerr,ggood_yerr);
794 tgraphgood->SetName(Form(
"graph_good_%d_%d_%d_%d",det,diblock,dcm,feb));
807 for (runIter =
mBreakRun[channum].begin(); runIter !=
mBreakRun[channum].end(); ++runIter) {
808 gstruct=runIter->second;
809 gbreak_y(iv) = gstruct.
mean;
810 gbreak_yerr(iv) = gstruct.
meanerr;
811 gbreak_x(iv) = gstruct.
time;
812 gbreak_xerr(iv) = gstruct.
timeerr;
816 tgraphbreak = tfs->make<TGraphAsymmErrors>(gbreak_x,gbreak_y,gbreak_xerr,gbreak_xerr,gbreak_yerr,gbreak_yerr);
817 tgraphbreak->SetName(Form(
"graph_break_%d_%d_%d_%d",det,diblock,dcm,feb));
827 MF_LOG_VERBATIM(
"DriftCorrection") <<
"ADDING DRIFTVECT SIZE " << driftvect.size() ;
double meanerr
Error on mean response.
void produce(art::Event &e) override
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.
diblock
print "ROW IS " print row
void AddChannelResponse(int fOffChan, double fMean, double fMeanErr, double fRMS, int fNHits)
Return subrun end time.
bool fDoCSVFile
FCL input: Write a CSV File?
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
::xsd::cxx::tree::exception< char > exception
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)
std::map< int, std::map< int, float > > mCorrvals
Map of Run->(channel,response)
DEFINE_ART_MODULE(TestTMapFile)
std::uint32_t endtime
Run end time.
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.
void respondToCloseOutputFiles(art::FileBlock const &fb) override
void LookForBreaks(std::map< int, respStruct > &chanlist, int channum)
Function to determine location of validity context breaks and corresponding Drift constants...
std::map< int, std::map< int, graphStruct > > mBreakRun
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
std::string fDriftRespLabel
FCL input: Label of the DriftResponse collection:
T get(std::string const &key) const
std::map< int, std::map< int, respStruct > > respMap
Map of Channel->(run,response)
TGraphAsymmErrors * tgraphall
Graph holding all runs.
void reconfigure(fhicl::ParameterSet const &p)
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.
uint32_t StartTime() const
Return run number.
std::map< int, MeanStruct > const & ReturnMap() const
std::string fCSVFile
FCL input: file name for output CSV File.
std::uint32_t starttime
Run start time.
void DriftResponseToStore(std::vector< caldp::DriftResponse > &drvec)
virtual ~DriftCorrection()
#define MF_LOG_VERBATIM(category)
int nhits
Number of hits in mean calculation.
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.
void beginRun(art::Run &r) override
int rms
RMS of hits in mean calculation.
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.
double mean
Mean response.
diblock_t getDiBlock(dchan daqchan) const
Decode the diblock ID from a dchan.
double mean
Mean response.
Encapsulate the geometry of one entire detector (near, far, ndos)
bool fDoDiagnosticPlot
FCL input: Make diagnostic histograms?
feb_t getFEB(dchan daqchan) const
Decode the feb id from a dchan.