DriftResponseCalc_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       DriftResponseCalc
// Module Type: producer
// File:        DriftResponseCalc_module.cc
//
// Generated at Mon Nov 19 15:48:48 2012 by Ruth Toner using artmod
// from art v1_02_06.
////////////////////////////////////////////////////////////////////////

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

//Framework:
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Run.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/PCHit.h"
#include "CalibrationDataProducts/DriftResponse.h"
#include "Geometry/Geometry.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "CMap/service/CMapService.h"
#include "Calibrator/Calibrator.h"
#include "CalibrationUtils/CalibUtil.h"

// ROOT includes
#include "TH1F.h"
#include "TH2F.h"

namespace caldp { 
  class PCHit; 
  class DriftResponse; 
}


namespace calib {

  class DriftResponseCalc : public art::EDProducer {
  public:
    explicit DriftResponseCalc(fhicl::ParameterSet const & p);
    virtual ~DriftResponseCalc();
    
    void produce(art::Event & e) override;
    void reconfigure(fhicl::ParameterSet const & p);
    void beginJob() override;
    void beginRun(art::Run & r) override;
    void endRun(art::Run & r) override;
    
  private:
    
    //Function to create ID for individual plane and FEB
    daqchannelmap::dchan GetChanID(daqchannelmap::plane_t plane, daqchannelmap::cell_t cell, novadaq::cnv::DetId det);
    
    //Fill map with hits:
    void FillPeCorrHits(art::Handle<std::vector<caldp::PCHit> > &cells, 
                        std::map< daqchannelmap::dchan, std::vector<double> > &mMeanPe,
                        novadaq::cnv::DetId detid);
    
    //Calculate mean, error on mean, number of bottom 90% hits
    void CalculateMean(std::vector<double> &chanvector,double &fMean,double &fMeanErr,int &fNHits, double &fRMS);
    
    //Draw distribution for bottom 90% hits
    void DrawMean(std::vector<double> &chanvector,TH1F *hpecorr);
    
    //Input string names
    std::string  fPCHitLabel  ;///Label of PCHit data
    std::string  fQualXYName  ;///Instance name, XY qual hits
    std::string  fQualZName   ;///Instance name, Z qual hits
    std::string  fQualAvgName ;///Instance name, Avg qual hits
    std::string  fFillType    ;///Method for calculating means (which quals to use)
    bool         fDrawMeans   ;///Save response histograms to file, true/false?
    int          fDiBlockNum  ;///DiBlock Number to use (set to 0 with fDCMNum=0 to run all)
    int          fDCMNum      ;///DCM Number to use (set to 0 with fDiBlockNum=0 to run all)
        
    //Histograms of PE to get Truncated Mean:
    std::map< daqchannelmap::dchan, std::vector<double> > mMeanPeXY;
    std::map< daqchannelmap::dchan, std::vector<double> > mMeanPeZ;
    std::map< daqchannelmap::dchan, std::vector<double> > mMeanPeAvg;
    
    //Time variables:
    std::vector<double> eventTime;
    std::uint32_t minTime ;
    std::uint32_t maxTime ;
    
    //Data or MC tag
    bool isData;

    //Output histograms for use later for debugging
    std::map<daqchannelmap::dchan, TH1F*> srHistprint;
    TH2F* hMeanDet;
    TH1F *hRMS;
    TH1F *hMeanFrac;
    TH1F *hNHits;
    
  };

DriftResponseCalc::DriftResponseCalc(fhicl::ParameterSet const & p)
{
  reconfigure(p);

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

}

DriftResponseCalc::~DriftResponseCalc()
{
}


void DriftResponseCalc::reconfigure(fhicl::ParameterSet const & pset)
{
      fPCHitLabel       =   pset.get< std::string >("PCHitLabel"); //Label of PCHitList module
      fQualXYName       =   pset.get< std::string >("QualXYName"); //Instance label, "XY" quality hits
      fQualZName        =   pset.get< std::string >("QualZName");  //Instance label, "Z" quality hits
      fQualAvgName      =   pset.get< std::string >("QualAvgName");//Instance label, "Avg" quality hits
      fDrawMeans        =   pset.get< bool        >("DrawMeans");  //Do or do not save mean histos
      fFillType         =   pset.get< std::string >("FillType");   //How to calc means
      fDiBlockNum       =   pset.get< int         >("DiBlockNum"); //Number of DiBlock to run (if = 0, run all)
      fDCMNum           =   pset.get< int         >("DCMNum");     //Number of DiBlock to run (if = 0, run all)

}


daqchannelmap::dchan DriftResponseCalc::GetChanID(daqchannelmap::plane_t plane, daqchannelmap::cell_t cell, novadaq::cnv::DetId det){
  ///Returns unique ID, channel number for pixel 0 in corresponding FEB

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

  //---Find and assign a unique ID number to the current FEB---//
  //Get the DAQ channel ID for this cell/plane/det:
  daqchannelmap::lchan logchannel = cmap->Map()->encodeLChan(det, plane, cell);
  daqchannelmap::dchan daqchannel = cmap->Map()->encodeDChan(logchannel);
  //Get the DiBlock number for this DAQ channel:
  daqchannelmap::diblock_t diblock =  cmap->Map()->getDiBlock(daqchannel);
  //Get the DMC for this DAQ channel:
  daqchannelmap::dcm_id_t dcm =  cmap->Map()->getDCM(daqchannel);
  //Get the FEB position for this DAQ channel:
  daqchannelmap::feb_t dcmlink =  cmap->Map()->getFEB(daqchannel);

  //Pixel 0 is the unique channel ID:
  daqchannelmap::dchan chanid = cmap->Map()->encodeDChan(det, diblock, dcm, dcmlink, 0);

 //Return unique ID for FEB:
 return chanid;

}

void DriftResponseCalc::FillPeCorrHits(art::Handle<std::vector<caldp::PCHit> > &cells, std::map< daqchannelmap::dchan, std::vector<double> > &mMeanPe, novadaq::cnv::DetId detid){
  ///Adds PECorr to vector for correct channel.

  daqchannelmap::dchan chanid; //Channel ID
  double fPeCorr=-1; //PECorr value to insert
  double attencorr = 1.0;
  size_t i=0;

  //Calibrator necessary to get Atten information
  art::ServiceHandle<calib::Calibrator> calibrator;

  //Loop over PCHits:
  for (i = 0; i < cells->size(); ++i){

    const caldp::PCHit *cell = &cells->at(i);

    //Go to the next hit if something strange happened with path calculation:
    if (cell->Path()<=0) continue;

    //Get the unique channel ID:
    chanid = GetChanID(cell->Plane(),cell->Cell(),detid);

    //Check if channel was added to map:
    if (mMeanPe.find(chanid) != mMeanPe.end() ) {
     
      //Make mock cellhit for use in Attenuation Calibrator
      rb::CellHit cellhit;
      cellhit.SetPlane(cell->Plane());
      cellhit.SetCell(cell->Cell());

      //Set hit as data or MC
      if (!isData) cellhit.SetMC(); 

      //Attenuation constant for this hit:
      attencorr = calibrator->GetAttenScale(cellhit, cell->W());

      //Skip this hit if cell is not calibrated.
      if (attencorr<=0) continue;

      //Correct the PE hit for Attenuation and Path length and add to vector:
      fPeCorr = cell->PE()*attencorr/cell->Path();
      mMeanPe[chanid].push_back(fPeCorr);

    }//End "if channel found in map"
 
  }//End cell loop

}

void DriftResponseCalc::CalculateMean(std::vector<double> &chanvector,double &fMean,double &fMeanErr,int &fNHits, double &fRMS){
  ///Function returns truncated mean (fMean), error on mean (fMeanErr), number of hits in truncated mean (fNHits), and RMS of
  ///response (fRMS) for vector of PECorrs for a single channel.

  //Initialize all to 0:
  fNHits   = 0;
  fMean    = 0;
  fMeanErr = 0;
  
  //Sort the list:
  std::sort(chanvector.begin(),chanvector.end());
  
  //Find boundary of the smallest 90% of hits
  fNHits = (int)(chanvector.size()*9/10);

  //Return if there are too few hits in the final list
  if (fNHits<=1) return;

  //Loop over the vector and calculate sum of hits
  double dsum=0;
  int i=0;
  for (i=0; i<fNHits; ++i) dsum += chanvector[i];

  //Also number of hits, for use as double in math
  double dnhits = (double)fNHits;

  //Calculate the mean:
  fMean=dsum/dnhits;

  //Sum up difference from mean squared:
  double dsumvar2=0;
  for (i=0; i<fNHits; ++i) dsumvar2 += (chanvector[i]-fMean)*(chanvector[i]-fMean);
  
  //Unbiased estimate of true variance:
  double dstdvar = sqrt(dsumvar2/(dnhits - 1));

  //Estimate of uncertainty on mean:
  fMeanErr = dstdvar/sqrt(dnhits);
  fRMS = dstdvar;
}

void DriftResponseCalc::DrawMean(std::vector<double> &chanvector,  TH1F *hpecorr){

  //Sort the list and grab bottom 90% of events:
  std::sort(chanvector.begin(),chanvector.end());
  int nhits = (int)(chanvector.size()*9/10);

  //Fill the histogram
  for (int i=0; i<nhits; ++i) hpecorr->Fill(chanvector[i]);

}


void DriftResponseCalc::beginJob()
{
}

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

  //Automatic assumption = we are running MC
  isData = false;

  //We want to be able to handle any time possible.
  minTime =  std::numeric_limits<uint32_t>::max()-1;
  maxTime =  std::numeric_limits<uint32_t>::min();

  //Check if map of XY hits is filled; if not, define vectors for all maps.
  //NOTE: Geometry service features dictate this must go in beginRun instead of beginJob.
  if (mMeanPeXY.size()==0){

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

    //Detector identity: 
    novadaq::cnv::DetId detid = geom->DetId();
    
    //Plane, cell, channel ID output:
    unsigned int planeid;
    unsigned int cellid;
    daqchannelmap::dchan chanid;
    
    //---Must make vector (for PCHits) for each FEB in detector.
    
    //Loop over X and Y views:
    for (int view = geo::kX; view <= geo::kY; ++view){

      //Grab set of planes in that view from geometry
      const std::set<unsigned int> planeSet = 
        (view == geo::kX) ? geom->GetPlanesByView(geo::kX) : geom->GetPlanesByView(geo::kY);
      
      //Iterate over std:set of planes:
      std::set<unsigned int >::iterator planeIter;
      for (planeIter = planeSet.begin(); planeIter != planeSet.end(); ++planeIter){

        //Grab plane number for this entry
        planeid = *planeIter;
        
        //Loop over cells in that plane:
        for (unsigned int icx = 0; icx < geom->Plane(planeid)->Ncells(); ++icx){
         
          //Cell number is loop entry
          cellid = icx;
 
          //Find FEB and create ID describing FEB and plane:
          chanid = GetChanID(planeid,cellid,detid);

          //Check to see if running this diblock and dcm:
          daqchannelmap::diblock_t diblock = cmap->Map()->getDiBlock(chanid);
          daqchannelmap::dcm_id_t dcmnum =  cmap->Map()->getDCM(chanid);

          //Are we supposed to run this DCM and DiBlock, according to fcl settings?
          if ( ((fDiBlockNum>=0&&diblock==(unsigned int)fDiBlockNum)||fDiBlockNum<0) &&
               ((fDCMNum>=0&&dcmnum==(unsigned int)fDCMNum)||fDCMNum<0) ){

            //Define new vector of hits for each map if FEB entry doesn't already exist.
            if (mMeanPeXY.find(chanid) == mMeanPeXY.end()) {
              
              mMeanPeXY.insert(std::pair<daqchannelmap::dchan,std::vector<double> >(chanid,std::vector<double>()));
              mMeanPeZ.insert(std::pair<daqchannelmap::dchan,std::vector<double> >(chanid,std::vector<double>()));
              mMeanPeAvg.insert(std::pair<daqchannelmap::dchan,std::vector<double> >(chanid,std::vector<double>()));
            }
          }
          
        }//end cell loop
      }//end plane loop
    }//end view loop
    
  } else {

    //Clear out vectors for individual channels to start with:
    
    daqchannelmap::dchan channum; 
    std::map< daqchannelmap::dchan, std::vector<double> >::iterator chanIter;
    for (chanIter = mMeanPeXY.begin(); chanIter != mMeanPeXY.end(); ++chanIter) {
      channum = chanIter->first;
      mMeanPeXY[channum].clear();
      mMeanPeZ[channum].clear();
      mMeanPeAvg[channum].clear();
    }
    
  }
}


void DriftResponseCalc::produce(art::Event & e)
{

  //Check if this is Data:
  //if (e.isRealData()) isData = true;

  //Geometry service
  art::ServiceHandle<geo::Geometry> geom;

  //Grab detector ID from geometry
  novadaq::cnv::DetId detid = geom->DetId();

  //---Grab collections of PCHits for each path quality--//
  //Highest XY quality events:
  art::Handle<std::vector<caldp::PCHit> > cellsxy;
  e.getByLabel(fPCHitLabel, fQualXYName, cellsxy);
  //Medium Z quality events:
  art::Handle<std::vector<caldp::PCHit> > cellsz;
  e.getByLabel(fPCHitLabel, fQualZName, cellsz);
  //Lowest Avg quality events:
  art::Handle<std::vector<caldp::PCHit> > cellsavg;
  e.getByLabel(fPCHitLabel, fQualAvgName, cellsavg);

  //Uncomment to access track/hit association:
  //art::FindOne<rb::Track> fa(cellsxy, e, art::InputTag("caliblist", fQualXYName));

  //Update min and max times.  This should later be edited with Jon Paley's info.
  std::uint32_t evtTimeS  = e.time().timeHigh();
  eventTime.push_back(evtTimeS);
  if (evtTimeS>maxTime) maxTime = evtTimeS;
  if (evtTimeS<minTime) minTime = evtTimeS;

  //Fill all cell types!
  FillPeCorrHits(cellsxy,mMeanPeXY,detid);
  FillPeCorrHits(cellsz,mMeanPeZ,detid);
  FillPeCorrHits(cellsavg,mMeanPeAvg,detid);
 
  //Uncomment to access association information, add own code.
  //    if(fa.at(i).isValid()){
  //      std::cout << "Track cells=" << fa.at(i).ref().NCell() << std::endl;
  //  }
  
}//End produce function

void DriftResponseCalc::endRun(art::Run & r)
{
  
  //Geometry service:
  art::ServiceHandle<geo::Geometry> geom;

  //Some information for histograms:
  const std::set<unsigned int> planeSetX =  geom->GetPlanesByView(geo::kX);
  const std::set<unsigned int> planeSetY =  geom->GetPlanesByView(geo::kY);
  int nplanes = planeSetX.size() + planeSetY.size();

  //---Summary histograms for this run--------------//

  art::ServiceHandle<art::TFileService> tfs;

  hMeanDet = tfs->make<TH2F>("hMeanDet",
                             "Response Plane vs. FEB; Plane Number; FEB",
                             nplanes+1,0,nplanes+1,3,0,3.);

  hRMS = tfs->make<TH1F>("hRMS",
                         "RMS of Truncated Response; RMS of Truncated Response (PE); Modules",
                         100,0,100.);

  hMeanFrac = tfs->make<TH1F>("hMeanFrac",
                             "Fractional Error on Mean; Frac Error on Mean; Modules",
                              100,0,0.10);

  hNHits = tfs->make<TH1F>("hNHits",
                             "Number of Hits in FEB; Number of Hits in FEB; Modules",
                              1000,0,50000.);

  //-----------------------------------------------//

  //Ptr to DriftResp information:
  std::unique_ptr<std::vector<caldp::DriftResponse > >fDriftResponse(new std::vector<caldp::DriftResponse >);

  //Run-specific values for Drift response
  int fRun = r.run();

  //To be changed with Jon Paley's updated info
  //art::Timestamp fStart = r.beginTime();
  //art::Timestamp fEnd   = r.endTime();

  caldp::DriftResponse driftresp(fRun,minTime,maxTime);
    
  //Remaining DriftResponse values to be calculated for each FEB:
  double fMean   =0;
  double fRMS    =0;
  double fMeanErr=0;
  int    fNHits  =0;
  int    fChanId =0;

  int FEBcount=0;
  int prevplane=-1;
  std::ofstream myfile;
  if (fDrawMeans) {
    myfile.open ("chaninfo.txt");
  }

  calib::EPathQuality fPathQual = calib::kXYAdjacents;

  //Iterate over map of channels -> XY hits to start with.
  std::map< daqchannelmap::dchan, std::vector<double> >::iterator chanIter;
  for (chanIter = mMeanPeXY.begin(); chanIter != mMeanPeXY.end(); ++chanIter) {

    //Grab ID and vector for XY:
    fChanId=chanIter->first;
    std::vector<double> chanvector;


    if (fFillType=="Successive"||fFillType=="Exclusive"){

      //vector of XY hits:
      chanvector=chanIter->second;

      //Calculate mean, error on mean, rms, and number of hits for XY hits
      fPathQual = calib::kXYAdjacents;
      CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);

      //Not enough statistics.  Add the Z hits.
      if (fMean==0||(fMean>0&&fMeanErr/fMean>0.03)){    
        if (fFillType=="Exclusive") chanvector.clear();
        chanvector.insert(chanvector.end(),
                          (mMeanPeZ.at(fChanId)).begin(),
                          (mMeanPeZ.at(fChanId)).end());
        
        //Re-calculate mean, etc. with Z hits added
        fPathQual = calib::kZAdjacents;
        CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);
        
        //Not enough statistics.  Add the Avg hits.
        if (fMean==0||(fMean>0&&fMeanErr/fMean>0.03)){
          if (fFillType=="Exclusive") chanvector.clear();
          chanvector.insert(chanvector.end(),
                            (mMeanPeAvg.at(fChanId)).begin(),
                            (mMeanPeAvg.at(fChanId)).end());
          
          //Re-calculate mean, etc. with Avg hits added
          fPathQual = calib::kAverage;
          CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);
         
          
        }//(End calc with added Avg hits)
        
      }//(End calc with added Z hits)

    } else if (fFillType=="XY"){

      //Vector from iterator:
     chanvector=chanIter->second;

     //Calculate mean, error on mean, and number of hits for XY hits:
     fPathQual = calib::kXYAdjacents;
     CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);

    } else if (fFillType=="Z"){

      //Z vector from other map:
      chanvector.insert(chanvector.end(),
                        (mMeanPeZ.at(fChanId)).begin(),
                        (mMeanPeZ.at(fChanId)).end());

     //Calculate mean, error on mean, and number of hits for Z hits:
      fPathQual = calib::kZAdjacents;
      CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);      

    } else if (fFillType=="Avg"){

      //Avg vector from other map:
      chanvector.insert(chanvector.end(),
                            (mMeanPeAvg.at(fChanId)).begin(),
                            (mMeanPeAvg.at(fChanId)).end());

     //Calculate mean, error on mean, and number of hits for Avg hits:
      fPathQual = calib::kAverage;
      CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);
      
    }

    //Final calibrated value is not good enough:
    if (fNHits==0){
      fMean=-1;
      fMeanErr=-1;
      fRMS=-1;
    }//(End failed calibration)

    if (fDrawMeans) {

      myfile << "//--------Statistics for ChanId=" << fChanId << " (initial size=" << chanvector.size() << ")\n";
      myfile << "//        Mean=" << fMean << "+/-" << fMeanErr << " RMS=" << fRMS << " n=" << fNHits << " qual=" << fPathQual <<"\n";

      std::string histname = Form("h_r%d_id%d",fRun,fChanId);
      
      srHistprint[fChanId] = tfs->make<TH1F>(Form("%s",histname.c_str()),
                             "Truncated PECorr;PECorr;Number of Hits",
                             200,0.,200.);
      DrawMean(chanvector,srHistprint[fChanId]);

    }

    //Fill 2D response plot for (plane vs FEB)
    int planefill = (int)(fChanId/1e6);
    if (planefill==prevplane){
      FEBcount++;
    } else { FEBcount=0; }
    prevplane=planefill;
    if (fMean>0) {
      hMeanDet->SetBinContent(planefill+1,FEBcount+1,fMean);
      hMeanFrac->Fill(fMeanErr/fMean);
      hRMS->Fill(fRMS);
      hNHits->Fill(fNHits);
    }

    //---Insert channel response-----//
    if (!driftresp.CheckForChannel(fChanId))
      driftresp.AddChannelResponse(fChanId,fMean,fMeanErr,fRMS,fNHits);

  }//(End calc over entries in XY hit map)

  fDriftResponse->push_back(driftresp);

  //Add the DriftResponse structure.
  r.put(std::move(fDriftResponse));
    if (fDrawMeans)  myfile.close();

}

}//end namespace
DEFINE_ART_MODULE(calib::DriftResponseCalc)