DriftCorrection_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
// Class:       DriftCorrection
// Module Type: producer
// File:        DriftCorrection_module.cc
//
// Generated at Mon Nov 19 15:49:27 2012 by Ruth Toner using artmod
// from art v1_02_06.
////////////////////////////////////////////////////////////////////////

//C++ and C:
#include <cmath>
#include <iostream>
#include <fstream>
#include <vector>

//Framework:
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/FileBlock.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Persistency/Provenance/Timestamp.h" 
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "cetlib_except/exception.h"
#include "messagefacility/MessageLogger/MessageLogger.h" 
#include "canvas/Persistency/Common/FindOne.h"

//NovaSoft:
#include "CalibrationDataProducts/DriftResponse.h"
#include "Geometry/Geometry.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "CMap/service/CMapService.h"
#include "CalibrationUtils/CalibUtil.h"
#include "Calibration/CDPStorage.h"

// ROOT includes
#include "TTree.h"
#include "TH1F.h"
#include "TFile.h"
#include "TH2F.h"
#include "TGraphAsymmErrors.h"
#include "TGraph.h"
#include "TVector.h"

#include <iostream>
#include <fstream>


namespace caldp { 
  class DriftResponse;   
}

namespace calib {

  class DriftCorrection : public art::EDProducer {

  public:
    explicit DriftCorrection(fhicl::ParameterSet const & p);
    virtual ~DriftCorrection();

    void produce(art::Event & e) override;
    void reconfigure(fhicl::ParameterSet const & p);
    void beginJob() override;
    void beginRun(art::Run & r) override;
    void respondToCloseOutputFiles(art::FileBlock const& fb) override;
    

  private:

    struct respStruct
    {
      ///Mean response
      double         mean;   
      ///Error on mean response   
      double         meanerr;
      ///Number of hits in mean calculation
      int            nhits;
      ///RMS of hits in mean calculation
      int            rms;     
      ///Run start time
      std::uint32_t  starttime;
      ///Run end time 
      std::uint32_t  endtime;
    };

    //FCL options:
    ///FCL input: Make diagnostic histograms?
    bool        fDoDiagnosticPlot;
    ///FCL input: Label of the DriftResponse collection:
    std::string fDriftRespLabel;
    ///FCL input: Number of days to use in aggregation:
    double      fNDay;
    ///FCL input: Automatically remove low-statistics runs?
    bool        fRemoveLowStat;
    ///FCL input: Run number at which to automatically cut off calibration
    int         fFlushRunNum;
    ///FCL input: Jump size at which to make new validity context
    double      fJumpSize;
    ///FCL input: Write a CSV File?
    bool        fDoCSVFile;
    ///FCL input: file name for output CSV File
    std::string fCSVFile;
    ///FCL input: Size of sliding window determining run quality.
    int         fQualWindow;
    ///FCL input: Is this a normalization job?
    bool        fDoNormalize;

    ///Holds maximum run number for this job
    int fMaxRun; 

    ///Graph holding all runs
    TGraphAsymmErrors *tgraphall; 
    ///Graph holding all good quality runs
    TGraphAsymmErrors *tgraphgood; 
    ///Graph holding all validity context information
    TGraphAsymmErrors *tgraphbreak; 

    ///Map of Channel->(run,response)
    std::map<int, std::map<int, respStruct> > respMap;

    ///Current list of run numbers
    std::vector<int> nRun;

    ///Run quality map: channel->(run,quality)
    std::map<int, std::map<int,bool> > runQual;

    ///Map of validity context start points (run number) to Drift constants.
    std::map<int,caldp::DriftResponse> BreakList;

    ///Function to determine location of validity context breaks and corresponding Drift constants.
    void LookForBreaks(std::map<int, respStruct> &chanlist, int channum);

    struct graphStruct
    {
      ///Mean response
      double mean;      
      ///Error on mean response
      double meanerr;
      ///End time   
      double time;   
      ///End time
      double timeerr;
    };

    ///Map of Run->(channel,response)
    std::map<int,std::map<int,float> > mCorrvals;

    std::map<int,std::map<int, graphStruct> > mAllRun;
    std::map<int,std::map<int, graphStruct> > mGoodRun;
    std::map<int,std::map<int, graphStruct> > mBreakRun;
  };
 

  ///Find validity contexts and Drift constants in an individual channel, for all runs in this channel 
  ///yet to be assigned to a validity context.
  void DriftCorrection::LookForBreaks(std::map<int, respStruct> &chanlist, int channum){

    //Clear out the run list:
    nRun.clear();

    //Make a list of runs currently active for this channel
    std::map<int, respStruct >::iterator runIter;
    for (runIter = chanlist.begin(); runIter != chanlist.end(); ++runIter ) nRun.push_back(runIter->first);

    //Make sure that the runs are properly ordered
    std::sort(nRun.begin(),nRun.end());

    //What is the final run on the list?
    int runfinal = nRun[nRun.size()-1];

    int nruni; //Holds a run number.
    double mean; //Hols a mean
    double rms; //Holds an rms
    double nhits; //Holds the number of hits
    double meanerr; //Holds an error on a mean
    std::uint32_t starttime; //Holds a start time
    std::uint32_t endtime = 0; //Holds an end time
    double residual; //Holds a residual for an individual run off a running mean

    //All active run numbers in this channel are "bad" by default
    int nrun;
    for (unsigned int irun=0; irun < nRun.size(); irun++) {
      nrun = nRun[irun];
      if (runQual[channum].find(nrun) == runQual[channum].end())  {
        runQual[channum].insert(std::make_pair(nrun,false));
      } 
    }
 
    //Vector to find median
    std::vector<double> sortvector;

    //Find running median of most recent runs in list.
    int irun=0;
    for (int i=0; i<fQualWindow; i++){
      irun = nRun.size() - fQualWindow + i; //Position in ordered run list
      if (irun < 0) continue;
      nruni = nRun[irun]; //The run number at list position irun
      mean = chanlist[nruni].mean; //Response for run nruni
      sortvector.push_back(mean); //Add the response to the running list
    }
        
    //Get the median entry:
    std::sort(sortvector.begin(),sortvector.end());
    int mentry = (int)(((double)sortvector.size())/2.0);
    double runningMedian = sortvector[mentry];
        
    //Loop back over running list to find good quality entries:
    for (int i=0; i<fQualWindow; i++){
      irun = nRun.size() - fQualWindow + i; //Position in ordered run list
      if (irun < 0) continue;
      nruni = nRun[irun]; //The run number at list position irun
      mean = chanlist[nruni].mean; //Response for run nruni
      meanerr = chanlist[nruni].meanerr; //Error on mean for run nruni
      nhits = chanlist[nruni].nhits; //Number of hits in run nruni
      starttime = chanlist[nruni].starttime; //Start time for run nruni
      endtime = chanlist[nruni].endtime; //End run for nruni

      //The following run types are automatically "bad" quality:
      if (mean<=0) { //Badly defined response
        runQual[channum][nruni] = false;
        continue;
      } else if ((fabs(meanerr/mean)>0.03 || nhits<200 ) && fRemoveLowStat==true ) { //"Low stat" run (if opt to remove automatically)
        runQual[channum][nruni] = false;
        continue;
      } else if (starttime>3.0e9||starttime<1.3e9||endtime<1.3e9){ //Runs with poorly defined times
        runQual[channum][nruni] = false;
        continue;
      } 

      //Calculate the difference between this run's response and the running median
      residual = fabs((mean/runningMedian) - 1.0); //Deviation from median in range.
      if (residual<fJumpSize) {
        runQual[channum][nruni] = true; 
        //Note: run has to be within range of median at ANY point in sliding range to be declared "good"
      }
    }

    //Records for old combined values in validity context
    double old_mean=0;
    double old_meanerr=0;
    double old_rms=0;
    double old_nhits=0;
    
    //Records for recalculated combined values in validity context
    double new_mean=0;
    double new_rms=0;
    double new_nhits=0;
    double new_meanerr=0;

    //Values for the next run
    double next_start=0;
    double next_mean=0;
    double next_err=0;
    double next_run=-1;

    //Counter variables
    double time_tot=0; //Total time in this validity context
    double time_0=0;   //Time at start of validty context
    int    run_0=-1;   //Run at start of validity context
    int    run_lastinvc=-1; //Last run in validity context

    //Boolean: do we flush existing responses at this run automatically?
    bool flushtherun=false;    
    bool atflushrun=false;

    //Look at all quality-assured runs (i.e., not final fQualWindow), unless told to flush at current run.

    int terminalirun; //Run list index (not run number) at which to stop

    if ( fFlushRunNum>=0 && runfinal == fFlushRunNum ) { 
      //Analyze up to final run in current list
      flushtherun=true;
      terminalirun=nRun.size(); 
    } else {
      //Analyze up to final run whose quality has been determined
      //(i.e., not ones currently in initial quality sliding window)
      terminalirun = nRun.size() - fQualWindow;
    }

    //Now loop over all runs up until current terminal window
    for (irun=0; irun < terminalirun; ++irun) {
      
      //This is the Flush run.
      if (irun==(terminalirun-1)&&flushtherun) atflushrun=true;

      nrun = nRun[irun];//Current run number

      //If the run is bad, skip to the next one.
      if (runQual[channum][nrun]) {

        //Current response values:
        mean = chanlist[nrun].mean; //Mean for this run
        meanerr = (mean>0) ? chanlist[nrun].meanerr : 0; //Error for this run
        rms = chanlist[nruni].rms; //RMS at this run
        nhits = chanlist[nrun].nhits; //Total hits for this run
        
        //Start and end times:
        starttime = chanlist[nrun].starttime; //Start time of this run
        endtime = chanlist[nrun].endtime; //End time of this run
        
        if (fDoDiagnosticPlot){
          //-----------DIAGNOSTIC--------------//
          //Add good runs to structure for adding to graph:
          if (mGoodRun[channum].find(nrun) == mGoodRun[channum].end() ){
            graphStruct graphAll;
            graphAll.mean    =  mean;
            graphAll.meanerr =  meanerr;
            graphAll.time    =  ((double)endtime + (double)starttime)/2.0;
            graphAll.timeerr  =  (double)endtime - graphAll.time;
            mGoodRun[channum].insert(std::make_pair(nrun,graphAll));
          }
        }
        
        //----------------------------------//
        
        //For all runs not at end of window:
        if (irun<terminalirun-1) {
          
          //Find the next good run in range:
          int igood=1;
          nruni = nRun[irun+igood];//The very next run

          //Keep counting up until either find the next good quality run, or hit end
          while ((!runQual[channum][nruni])&&(irun+igood<=terminalirun)) {
            igood++;
            nruni = nRun[irun+igood];
          }
          
          //If next good run is in range, find its values:
          if (irun+igood<terminalirun){ 
            next_run = nRun[irun+igood];
            next_start = chanlist[next_run].starttime;
            next_mean = chanlist[next_run].mean;
            next_err = chanlist[next_run].meanerr;
          } else continue; //Nothing new to look at in this range
        } else if (flushtherun) { //This is the last run, but it's the flush run
          next_run = nrun;
          next_start = endtime;
          next_mean = mean;
          next_err = meanerr;
        } else continue;
        
        
        //If this is the first run in the VC, set the "initial" conditions
        if (time_tot==0) {
          time_0 = starttime;
          run_0 = nrun;
        }
        
        //Add current hits to old total to make "new" total
        new_nhits = old_nhits + nhits;
        
        //If there are hits to use, start calculating new values:
        if (new_nhits!=0){
          
          //Combine together current mean response and old one to form "new" mean response
          new_mean = (old_nhits*old_mean + nhits*mean)/(new_nhits);
          
          //Some useful expressions to calculate new RMS
          double sqsum = rms*rms*(nhits-1) + nhits*mean*mean;
          double old_sqsum = old_rms*old_rms*(old_nhits-1) + old_nhits*old_mean*old_mean;
          double new_rms2 = (sqsum + old_sqsum - (new_nhits*new_mean*new_mean))/(new_nhits-1);
          
          //Update the RMS:
          new_rms=sqrt(new_rms2);
          
          //Update the error on mean response:
          new_meanerr = new_rms/sqrt(new_nhits);
          
          //Total elapsed time in this validity context:
          time_tot = endtime - time_0;
          
          //How do the mean and time change between now and next run in list?
          double nextjump = ( new_meanerr > 0 ) ? fabs(next_mean - new_mean)/new_mean : 0;
          double nexterrfrac = ( next_mean > 0 ) ? (fabs(next_err/next_mean)) : 0;
          
          //Search for a new break based on A) Total time in run, B) Time between end of this run and start of next
          //C) Fractional jump between current mean and next run's, D) Whether we're right at the "Flush Run" :
          if (
              time_tot>fNDay*86400.
              ||(next_start-endtime)>fNDay*86400./2.
              ||(nextjump>fJumpSize&&nexterrfrac<0.3)
              ||(atflushrun==true)
              ){
            
            //Add an entry to the Break List for this channel and run!
            BreakList[run_0].AddChannelResponse(channum,new_mean,new_meanerr,new_rms,new_nhits);
            
            //Add Diagnostic entries:
            if (fDoDiagnosticPlot){
              //-----------DIAGNOSTIC--------------//
              //Add entry to graphs for this break validity context:
              if (mBreakRun[channum].find(run_0) == mBreakRun[channum].end() ){
                graphStruct graphAll;
                graphAll.mean    =  new_mean;
                graphAll.meanerr =  new_meanerr;
                graphAll.time    =  ((double)endtime + (double)time_0)/2.0;
                graphAll.timeerr  =  (double)endtime - graphAll.time;
                mBreakRun[channum].insert(std::make_pair(run_0,graphAll));
              }
              //----------------------------------//
            }
            
            //Reset all counters and records:
            time_tot=0;
            new_mean=0;
            new_rms=0;
            new_nhits=0;
            new_meanerr=0;
            old_mean=0;
            old_rms=0;
            old_nhits=0;
            
            //Record the last run in the validity context
            if (nrun>run_lastinvc) run_lastinvc = nrun;
            
          }         
          
          //The updated version will be the "old" version for the next run:
          old_mean=new_mean;
          old_rms=new_rms;
          old_nhits=new_nhits;
          old_meanerr=new_meanerr;
          
        }
      } else if (atflushrun&&old_nhits>0){ //This is a bad run, but we need to flush anyway

        //Use what we have now to make a break:
        BreakList[run_0].AddChannelResponse(channum,old_mean,old_meanerr,old_rms,old_nhits);
        
        //-----------DIAGNOSTIC--------------//
        //Add good runs:
        if (mBreakRun[channum].find(run_0) == mBreakRun[channum].end() ){
          graphStruct graphAll;
            graphAll.mean    =  old_mean;
            graphAll.meanerr =  old_meanerr;
            graphAll.time    =  ((double)endtime + (double)time_0)/2.0;
            graphAll.timeerr  =  (double)endtime - graphAll.time;
            mBreakRun[channum].insert(std::make_pair(run_0,graphAll));
          }
          //----------------------------------//


          time_tot=0;

          new_mean=0;
          new_rms=0;
          new_nhits=0;
          new_meanerr=0;
            
          old_mean=0;
          old_meanerr=0;
          old_rms=0;
          old_nhits=0;

          //Last run to be used in a validity context:
          if (nrun>run_lastinvc) run_lastinvc = nrun;

      }
    }

    std::vector<int> eraselist;
    int run_n;
    //Clear everything in map up to final run in last validity context.
    //First, make a record of thigns that need to get erased:
    for (runIter = chanlist.begin(); runIter != chanlist.end(); ++runIter ) {
      run_n = runIter->first;
      if (run_n <= run_lastinvc ) {
        eraselist.push_back(run_n);
      } else break;
    }
    //Then run back over and do the erasing:
    for (unsigned int i=0; i<eraselist.size(); i++){
      int run_n = eraselist[i];
      chanlist.erase(run_n);
    }

  }


  void DriftCorrection::reconfigure(fhicl::ParameterSet const & p)
  {
    fDriftRespLabel       =   p.get< std::string >("DriftRespLabel"); //Label of DriftResponse module
    fNDay                 =   p.get< double      >("NDay"); //Number of days to employ
    fRemoveLowStat        =   p.get< bool        >("RemoveLowStat"); //Number of days to employ
    fFlushRunNum          =   p.get< int         >("FlushRunNum"); //Run number at which to automatically flush runs
    fDoDiagnosticPlot     =   p.get< bool        >("DoDiagnosticPlot"); //Make file with Drift Constants and data vs time
    fJumpSize             =   p.get< double      >("JumpSize"); //Fractional jump needed to declare new correction
    fDoCSVFile            =   p.get< bool        >("DoCSVFile"); //Write the results to a CSV file
    fCSVFile              =   p.get< std::string >("CSVFile"); //CSV File Name
    fQualWindow           =   p.get< int         >("QualWindow"); //Sliding quality window
    fDoNormalize          =   p.get< bool        >("DoNormalize"); //Make file with Drift Constants and data vs time
  }
  
  DriftCorrection::DriftCorrection(fhicl::ParameterSet const & p)
  {
    reconfigure(p);

    //DriftCorrection output produced every Run
    produces< std::vector<caldp::DriftResponse >, art::InRun >();


  }

  DriftCorrection::~DriftCorrection()
  {
  }

  void DriftCorrection::produce(art::Event & e)
  {
  }

  void DriftCorrection::beginJob()
  {
    fMaxRun=-1;
  }

  void DriftCorrection::beginRun(art::Run & r)
  {

    //Geometry service:
    art::ServiceHandle<geo::Geometry> fGeom;
    //Channel map service:
    art::ServiceHandle<cmap::CMap> cmap;

    //Collection of mean responses:
    art::Handle<std::vector<caldp::DriftResponse> > meancol;
    r.getByLabel(fDriftRespLabel, meancol);

    int           runnum;
    std::uint32_t timestart;
    std::uint32_t timeend;
    std::uint32_t channum;

    //Run over collection of mean responses:
    for (size_t i = 0; i < meancol->size(); ++i){

      //Grab drift response:
      const caldp::DriftResponse *driftmean = &meancol->at(i);
      
      //Grab map of channel to response for this run:
      std::map< int, caldp::MeanStruct > chanMap = driftmean->ReturnMap();

      //Grab run/timing information for this entry:
      runnum = driftmean->Run();
      timestart = driftmean->StartTime();
      timeend = driftmean->EndTime();

      //Poorly defined start or end time to run.  Skip this run.
      if (timestart>3.0e9||timestart<1.3e9||timeend<1.3e9) {
          LOG_VERBATIM("DriftCorrection") << "Run has poorly defined start or end time.  Skipping." ;
        continue;
      }
      //Check that runs are in order
      if (runnum>fMaxRun) {
        fMaxRun=runnum;
      } else throw cet::exception("Rethrow") << "ERROR: Runs have been improperly ordered.  Aborting.\n";

      //Add a new entry to the output map
      if (BreakList.find(runnum) == BreakList.end())
        BreakList.insert(std::pair<int, caldp::DriftResponse >(runnum,caldp::DriftResponse(runnum,timestart,timeend)));

      //Iterate over channels make map of channels->(run,response):
      std::map< int, caldp::MeanStruct >::iterator chanIter;
      for (chanIter = chanMap.begin(); chanIter != chanMap.end(); ++chanIter) {
 
        //Channel number:
        channum = chanIter->first;

        //Structure with response info
        respStruct response;
        response.mean    = chanIter->second.Mean();
        response.rms    = chanIter->second.RMS();
        response.meanerr = chanIter->second.MeanErr();
        response.nhits   = chanIter->second.NHits();
        response.starttime = timestart;
        response.endtime   = timeend;

        //If haven't added empty map for this channel yet, do so:
        if (respMap.find(channum) == respMap.end())  {
          respMap.insert(std::pair<int, std::map<int, respStruct> >(channum,std::map<int, respStruct>()));
          //Run quality for an individual channel and run:
          runQual.insert(std::pair<int, std::map<int, bool> >(channum,std::map<int, bool>()));
        }

        //Add unique entry for this run:
        if (respMap[channum].find(runnum) == respMap[channum].end() ){
          respMap[channum].insert(std::make_pair(runnum,response));
        } else throw cet::exception("Rethrow") << "ERROR: run was already added to the map; runs not ordered correctly.  Aborting.\n";

        //---------------DIAGNOSTIC----------------//
        //Fill structure for later diagnostic plots:
        graphStruct graphAll;
        graphAll.mean    =  chanIter->second.Mean();
        graphAll.meanerr =  chanIter->second.MeanErr();
        graphAll.time    =  ((double)timeend + (double)timestart)/2.0;
        graphAll.timeerr  =  (double)timeend - graphAll.time;

        //Add entries to all maps:
        if (mAllRun.find(channum) == mAllRun.end()){
          mAllRun.insert(std::pair<int, std::map<int, graphStruct> >(channum,std::map<int, graphStruct>()));
          mGoodRun.insert(std::pair<int, std::map<int, graphStruct> >(channum,std::map<int, graphStruct>()));
          mBreakRun.insert(std::pair<int, std::map<int, graphStruct> >(channum,std::map<int, graphStruct>()));
        }
        //Add "total" run:
        if (mAllRun[channum].find(runnum) == mAllRun[channum].end() ){
          mAllRun[channum].insert(std::make_pair(runnum,graphAll));
        } else throw cet::exception("Rethrow") << "ERROR: run was never added to the map; something strange has happened.  Aborting.\n";
        //------------------------------------------//

        //Start looking for breaks after 9 runs have been accumulated!
        //Do this individually for each run, looking over runs in order
        if ( (int)respMap[channum].size() > (fQualWindow-1) || runnum == fFlushRunNum )
          LookForBreaks(respMap[channum],channum);

      }//End of channel loop

    }//End of meancol loop

  }


  void DriftCorrection::respondToCloseOutputFiles(art::FileBlock const&)
  {

    //Channel map:
    art::ServiceHandle<cmap::CMap> cmap;

    std::ofstream filecsv;
    if (fDoCSVFile) 
      filecsv.open(fCSVFile.c_str());
    
    //Vector containing drift responses for later storage
    std::vector<caldp::DriftResponse> driftvect;

    //Vector for use in normalization code:
    std::map< int, std::vector< std::pair<double,int> > > meanrecord;

    //Iterator for break list
    std::map<int,caldp::DriftResponse>::iterator brkIter;

    //Loop back over break list
    for (brkIter = BreakList.begin(); brkIter != BreakList.end(); ++brkIter) {

      //Drift Response object:
      caldp::DriftResponse driftresp = brkIter->second;
  
      //Grab map of channel to response for this run, just to print:
      std::map< int, caldp::MeanStruct > chanMap = driftresp.ReturnMap();

      //Loop for debugging and records:
      LOG_VERBATIM("DriftCorrection") << "Start=" << driftresp.StartTime() << " End=" << driftresp.EndTime() ;
      std::map< int, caldp::MeanStruct >::iterator chanIter;
      for (chanIter = chanMap.begin(); chanIter != chanMap.end(); ++chanIter) {
        LOG_DEBUG("DriftCorrection") << "   channel " << chanIter->first 
                                     << " : mean=" << chanIter->second.Mean()
                                     << " : meanerr=" << chanIter->second.MeanErr()
                                     << " : rms=" << chanIter->second.RMS()
                                     << " : nhits=" << chanIter->second.NHits();
        
        //Print out the breaks which were made:
        LOG_VERBATIM("DriftCorrection") << "Push back chan=" << chanIter->first << " with " << chanIter->second.Mean() << "/" << chanIter->second.NHits() ;
        
        //If requested, and not normalizing, turn output into a csv file as well:
        if (fDoCSVFile&&!fDoNormalize) 
          filecsv << chanIter->second.Mean() << "," << chanIter->second.MeanErr() << "," << chanIter->first << "," <<  driftresp.StartTime() << "\n";
        
        //If we are normalizing, save the mean and nhits for this run for the channel
        if (fDoNormalize)
          meanrecord[chanIter->first].push_back(std::make_pair(chanIter->second.Mean(),chanIter->second.NHits()));

      }//end channel loop

        //Store drift response for later, if not normalizing:
      if (!fDoNormalize) driftvect.push_back(driftresp);

    }//end break run loop


    //If we are doing normalization, calculate the mean per channel:
    if (fDoNormalize){

      //Response at time "zero" for normalization
      caldp::DriftResponse normresp(0,0,0); 

      int chanid=-1; //channel
      double meansum,nsum,newmean; //iterators and means

      //Run back over channels and their stored response and nhit info:
      std::map< int, std::vector< std::pair<double,int> > >::iterator meanIter;
      for (meanIter = meanrecord.begin(); meanIter != meanrecord.end(); ++meanIter) {

        chanid = meanIter->first; //channel number
        std::vector< std::pair<double,int> > meanlist = meanIter->second; //list of responses
        
        //Reset counters:
        meansum=0;
        nsum=0;
        
        //Loop over response values to calculate the mean for this channel:
        for (unsigned int ival=0; ival<meanlist.size(); ++ival){
          //Sum of mean and number of hits:
          meansum += (meanlist[ival].first)*(meanlist[ival].second);
          nsum += meanlist[ival].second;
        }
        
        //Calculate the mean:
        newmean = (nsum>0) ? meansum/nsum : -1;
        
        //Add this as a DriftResponse:
        normresp.AddChannelResponse(chanid,newmean,-1,-1,nsum);
        
        //Print to CSV file if requested:
        if (fDoCSVFile)
          filecsv << newmean << "," << -1 << "," << chanid << "," << 0 << "\n";
        
      }//end channel loop
      
      //Save the "normalized" version to the drift vector
      driftvect.push_back(normresp);
      
    }//end normalize if statement
    
    //Close the CSV file:
    if (fDoCSVFile) filecsv.close();

    //-------------DIAGNOSTIC--------------//
    //Make plots of the runs and breaks vs time, if requested
    if (fDoDiagnosticPlot){
      int det, diblock, dcm, feb;
      daqchannelmap::dchan channum;

      //File Service:
      art::ServiceHandle<art::TFileService> tfs;

      //Iterator for a map containing graph points:
      std::map<int,std::map<int, graphStruct> >::iterator graphIter;

      //Loop over all channels to make plots of Data and Fit vs Time.
      for (graphIter = mAllRun.begin(); graphIter != mAllRun.end(); ++graphIter) {
        
        //Channel number:
        channum = graphIter->first;

        //Get "human" values using channel service:
        det = cmap->Map()->getDetector(channum);
        diblock = cmap->Map()->getDiBlock(channum);
        dcm = cmap->Map()->getDCM(channum);
        feb = cmap->Map()->getFEB(channum);

        //----Save ALL Runs to graph----------------------------//      
        TVector gall_y(mAllRun[channum].size());
        TVector gall_yerr(mAllRun[channum].size());
        TVector gall_x(mAllRun[channum].size());
        TVector gall_xerr(mAllRun[channum].size());
        int iv=0;
        graphStruct gstruct;
        std::map<int, graphStruct>::iterator runIter;
        for (runIter = mAllRun[channum].begin(); runIter != mAllRun[channum].end(); ++runIter) {
          gstruct=runIter->second;
          gall_y(iv)    = gstruct.mean;
          gall_yerr(iv) = gstruct.meanerr;
          gall_x(iv)    = gstruct.time;
          gall_xerr(iv) = gstruct.timeerr;
          iv++;
        }       
        //Make the TGraphAsymm here and write it out
        tgraphall  = tfs->make<TGraphAsymmErrors>(gall_x,gall_y,gall_xerr,gall_xerr,gall_yerr,gall_yerr);
        tgraphall->SetName(Form("graph_all_%d_%d_%d_%d",det,diblock,dcm,feb));
        tgraphall->Write();
        //-------------------------------------------------------//
        
        //----Do the same for Good Quality Runs--------------------------//
        //CHECK FOR CHANNEL'S EXISTENCE:
        iv=0;
        if (mGoodRun.find(channum) != mGoodRun.end())  {
          TVector ggood_y(mGoodRun[channum].size());
          TVector ggood_yerr(mGoodRun[channum].size());
          TVector ggood_x(mGoodRun[channum].size());
          TVector ggood_xerr(mGoodRun[channum].size());
          for (runIter = mGoodRun[channum].begin(); runIter != mGoodRun[channum].end(); ++runIter) {
            gstruct=runIter->second;
            ggood_y(iv)    = gstruct.mean;
            ggood_yerr(iv) = gstruct.meanerr;
            ggood_x(iv)    = gstruct.time;
            ggood_xerr(iv) = gstruct.timeerr;
            iv++;           
          }       
          //Make the TGraphAsymm here and write it out
          tgraphgood  = tfs->make<TGraphAsymmErrors>(ggood_x,ggood_y,ggood_xerr,ggood_xerr,ggood_yerr,ggood_yerr);
          tgraphgood->SetName(Form("graph_good_%d_%d_%d_%d",det,diblock,dcm,feb));
          tgraphgood->Write();
        }
        //-------------------------------------------------------//
        
        //----Do the same for Break Runs-------------------------//
        //CHECK FOR CHANNEL'S EXISTENCE
        iv=0;
        if (mBreakRun.find(channum) != mBreakRun.end())  {
          TVector gbreak_y(mBreakRun[channum].size());
          TVector gbreak_yerr(mBreakRun[channum].size());
          TVector gbreak_x(mBreakRun[channum].size());
          TVector gbreak_xerr(mBreakRun[channum].size());
          for (runIter = mBreakRun[channum].begin(); runIter != mBreakRun[channum].end(); ++runIter) {
            gstruct=runIter->second;
            gbreak_y(iv)    = gstruct.mean;
            gbreak_yerr(iv) = gstruct.meanerr;
            gbreak_x(iv)    = gstruct.time;
            gbreak_xerr(iv) = gstruct.timeerr;
            iv++;
          }
          //Make the TGraphAsymm here and write it out
          tgraphbreak  = tfs->make<TGraphAsymmErrors>(gbreak_x,gbreak_y,gbreak_xerr,gbreak_xerr,gbreak_yerr,gbreak_yerr);
          tgraphbreak->SetName(Form("graph_break_%d_%d_%d_%d",det,diblock,dcm,feb));
          tgraphbreak->Write();
        }
        //-------------------------------------------------------//
      }//End graph iteration loop
    }//--------END DIAGNOSTIC-------------// 


    //Write to file 
    art::ServiceHandle<calib::CDPStorage> cdpstorage;
    LOG_VERBATIM("DriftCorrection") << "ADDING DRIFTVECT SIZE " << driftvect.size() ;
    cdpstorage->DriftResponseToStore(driftvect);
   
 
  }
  
} 
DEFINE_ART_MODULE(calib::DriftCorrection)