ToFPositionRecoAnalysis_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////
/// \file    ToFPositionRecoAnalysis_module.cc
/// \module  analyzer
/// \brief   Experiment with position reconstruction
/// \author  John Rabaey, University of Dallas, <jrabaey@udallas.edu>
/// \date    March 2019
////////////////////////////////////////////////////////////////////////////  

// framework
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art_root_io/TFileService.h"
#include "art_root_io/TFileDirectory.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindMany.h"

// nova
#include "RawData/RawBeamline.h"
#include "BeamlineRecoBase/BeamlineDigit.h"
#include "BeamlineRecoBase/BeamlinePID.h"
#include "BeamlineRecoBase/MuonStackTrack.h"
#include "BeamlineRecoBase/ToF.h"
#include "BeamlineRecoBase/WCTrack.h"
#include "TestBeam/ToFPositionFinder.h"
#include "BeamlineUtils/BeamlineChannelMap.h"

// ROOT
#include "TH1D.h"
#include "TH1I.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TGraph.h"

// stl
#include <iostream>
#include <algorithm>
#include <fstream>

//a whole bunch of histograms?
bool dumpPosns = false;

//main histogram window constants
double MIN_TOF = 40;
double MAX_TOF = 47;
double N_BINS = 100;

double MIN_DT = -5;
double MAX_DT = 5;

double MAX_SPREAD = 4;
double SPREAD_CUTOFF = 4;

bool POSN_ZOOM_IN = true;

//position reconstruction constants
bool VARY_SIDE = false;
double SCINT_SIDE_LENGTH = 2.5; // In ns // use if vary_side is false
double MIN_SCINT_SIDE = 0; //use if vary_side is true
double MAX_SCINT_SIDE = 4.0;
double SCINT_SIDE_STEP = .1;


//position reconstruction dump binning
int POSN_N_BINS = 40;

enum ToFCounter {
  kToFCounterError,
  kToFCounterUS,
  kToFCounterDS,
  kToFCounterDSSiPM
};

//channel numbers
const int nChannels = 4;
std::map<ToFCounter,std::vector<unsigned int>> channels = {
  { kToFCounterDS,    {15,14,13,12} },
  { kToFCounterUS,    {8, 11,10,9 } },
  { kToFCounterDSSiPM,{7, 6, 5, 4 } }
};
std::vector<unsigned int> requiredChannels = {};// channels[kToFCounterDSSiPM];

std::vector<ToFCounter> activeToFCounters = {
  kToFCounterUS,
  kToFCounterDS,
  kToFCounterDSSiPM
};

std::map<ToFCounter,TString> tofCounterStrings = {
  {kToFCounterUS,     "USToF__"},
  {kToFCounterDS,     "DSToF__"},
  {kToFCounterDSSiPM, "SiPMToF"},
  {kToFCounterError,  "ErrorTF"}
};

std::map<ToFCounter,double> tofCounterDelays = {
  { kToFCounterUS,     21.0 },
  { kToFCounterDS,     0    },
  { kToFCounterDSSiPM, 0    }
};

// -----------------------------------------------------------------------
namespace beamlinereco {
  class ToFPositionRecoAnalysis : public art::EDAnalyzer {

  public:

    ToFPositionRecoAnalysis(const fhicl::ParameterSet& pset);

    void reconfigure(const fhicl::ParameterSet& pset);
    void beginJob() override;
    void analyze(const art::Event& evt);
    void endJob() override;
    
  private:

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

    std::unordered_map<unsigned int, std::string> channelMap;

    art::InputTag fRecoToFLabel;
    art::InputTag fToFDigitLabel;

    std::vector<std::pair<unsigned int, std::string>> fChannels;
    double fSamplingInterval;

    TH1D* hToF;
    std::map<int, std::map<int, TH1D*>>     hdtByChannel;
    std::map<ToFCounter,TH1D*>              hFirstHits;
    std::map<ToFCounter,std::vector<TH2D*>> hPosnRecoDump;
    std::map<ToFCounter,std::vector<TH2D*>> hPosn;
    std::map<ToFCounter,std::vector<TH1D*>> hPosnSpread;

    std::map<ToFCounter,std::vector<int>>   nPhysicalPosns;
    std::map<ToFCounter,std::vector<int>>   nNonphysicalPosns;
    std::map<ToFCounter,int>                nFourHitEvents;
    std::map<ToFCounter,int>                nEvents;
    std::map<ToFCounter,int>                nHighSpreadPosns;

    ToFPositionFinder* posnfinder;

    // services
    art::ServiceHandle<beamlineutil::BeamlineChannelMap> fChannelMap;
  };

  DEFINE_ART_MODULE(ToFPositionRecoAnalysis)

}

// -----------------------------------------------------------------------
beamlinereco::ToFPositionRecoAnalysis::ToFPositionRecoAnalysis(const fhicl::ParameterSet& pset) : EDAnalyzer(pset) {
  this->reconfigure(pset);
}

// -----------------------------------------------------------------------
void beamlinereco::ToFPositionRecoAnalysis::reconfigure(const fhicl::ParameterSet& pset) {
  fRecoToFLabel             = pset.get<art::InputTag>("RecoToFLabel");
  fToFDigitLabel            = pset.get<art::InputTag>("ToFDigitLabel");


  fChannels                 = pset.get<std::vector<std::pair<unsigned int, std::string>>>("Channels");
  fSamplingInterval         = pset.get<double>("TimeSamplingInterval");

  channelMap.clear();
  for (auto ch : fChannels) { channelMap[ch.first] = ch.second; }
}

// -----------------------------------------------------------------------
void beamlinereco::ToFPositionRecoAnalysis::beginJob() {
  
  //set up directories, main ToF distributions
  this->hToF = tfs->make<TH1D>("hToF", "Time of Flight", N_BINS, MIN_TOF, MAX_TOF);

  //Do everything three times, once for each ToF
  for(auto tofc : activeToFCounters) {
    art::TFileDirectory tofDir = tfs->mkdir(tofCounterStrings[tofc].Data());

    //set up dt distributions
    art::TFileDirectory dtDir = tofDir.mkdir("dt_distributions");
  
    for(int i = 0; i < nChannels; i++) {
      for(int j = i+1; j < nChannels; j++) {
        hdtByChannel[channels[tofc][i]][channels[tofc][j]] = 
          dtDir.make<TH1D>(Form("h_dt_ch%d_%d", channels[tofc][i], channels[tofc][j]),
                           Form("dt for channels %d and %d", channels[tofc][i], channels[tofc][j]),
                           N_BINS, MIN_DT, MAX_DT);                              
      }//for j
    }//for i

    //set up plot for the first channel in each hit group
    hFirstHits[tofc] = tofDir.make<TH1D>("first_hits", "Number of times each channel is first in a cluster",4,0,4);
  
    //set up position distribution
    if(VARY_SIDE) {
      //better put things in a folder
      art::TFileDirectory posnstuff = tofDir.mkdir("hPosn");
      //vector needs to contain one double per side length
      for(double d = MIN_SCINT_SIDE; d < MAX_SCINT_SIDE; d += SCINT_SIDE_STEP) {
        this->hPosn[tofc].push_back(posnstuff.make<TH2D>(Form("hPosnS%.3f",d), 
                                                         Form("Hit position with side length constant %.3f", d), 
                                                         N_BINS, 
                                                         (POSN_ZOOM_IN ? 0 :  -d),
                                                         (POSN_ZOOM_IN ? d : 2*d),
                                                         N_BINS,
                                                         (POSN_ZOOM_IN ? 0 :  -d),
                                                         (POSN_ZOOM_IN ? d : 2*d)));
        this->hPosnSpread[tofc].push_back(posnstuff.make<TH1D>(Form("hPosnSpreadS%.3f", d),
                                                               Form("Position spread with side length constant %.3f", d),
                                                               N_BINS, 0, MAX_SPREAD));
        nPhysicalPosns[tofc].push_back(0);
        nNonphysicalPosns[tofc].push_back(0);
      } //for d
    } //if
    
    else {//(!VARY_SIDE)
      //vector only needs one thing
      this->hPosn[tofc].push_back(tofDir.make<TH2D>("hPosn", "Position", 
                                                  N_BINS, 
                                                  (POSN_ZOOM_IN ? 0                 :  -SCINT_SIDE_LENGTH),
                                                  (POSN_ZOOM_IN ? SCINT_SIDE_LENGTH : 2*SCINT_SIDE_LENGTH),
                                                  N_BINS,
                                                  (POSN_ZOOM_IN ? 0                 :  -SCINT_SIDE_LENGTH),
                                                  (POSN_ZOOM_IN ? SCINT_SIDE_LENGTH : 2*SCINT_SIDE_LENGTH)));
      this->hPosnSpread[tofc].push_back(tofDir.make<TH1D>("hPosnSpread", "Position spread",
                                                          N_BINS, 0, MAX_SPREAD)); 
      nPhysicalPosns[tofc].push_back(0);
      nNonphysicalPosns[tofc].push_back(0);
    } //else
    
    //set draw options for 2d histograms
    for(TH2D* posnHist : hPosn[tofc]) {
      posnHist->SetOption("colz");
    }
  
    //set the useful counter
    nEvents[tofc] = 0;
    nHighSpreadPosns[tofc] = 0;
    nFourHitEvents[tofc] = 0;
    
  }//for tofc

  posnfinder = new ToFPositionFinder();
}

// -----------------------------------------------------------------------
void beamlinereco::ToFPositionRecoAnalysis::analyze(const art::Event& e) {
  
  // Get ToFs for each reco version
  art::Handle<std::vector<brb::ToF>> tofHandle;
  std::vector<art::Ptr<brb::ToF>> tofs;
  if (e.getByLabel(fRecoToFLabel, tofHandle)) {
    art::fill_ptr_vector(tofs, tofHandle);
  }
  unsigned int nToFs = tofs.size();
  
  // Get time of flight for each ToF
  std::vector<double> flightTimes; 
  for(unsigned int i = 0; i < nToFs; i++) {
    flightTimes.push_back(tofs[i]->Time());
  }

  //fill tof histogram
  for(unsigned int i = 0; i < nToFs; i++) {
    hToF->Fill(flightTimes[i]);
  }
  
  
  // Get associated digits for each reco version
  art::Handle<std::vector<brb::BeamlineDigit>> tofDigitHandle;
  std::vector<art::Ptr<brb::BeamlineDigit>> tofDigits;
  if (e.getByLabel(fToFDigitLabel, tofDigitHandle)) {
    art::fill_ptr_vector(tofDigits, tofDigitHandle);
  }

  // Retrieve start times by channel
  std::unordered_map<int,double> startTimeByChannel;
  //first, initialize all required channels to -1. This flag will allow us to skip the event 
  //  if the channel doesn't get a good start time value
  for(auto tofc : activeToFCounters) {
    for(int ch : channels[tofc]) {
      startTimeByChannel[ch] = -1;
    }//for ch
  }//for tofc

  for(auto d : tofDigits){
    //get start time by channel
    startTimeByChannel[fChannelMap->OnlineDigitChannel(d->ChannelID())] = d->StartTimeInNanoSec();
  }//for
    
  //skip the event if it doesn't have all required channels
  for(int ch : requiredChannels) {
    if(startTimeByChannel[ch] <= 0) {
      return;
    }//if
  }//for
  
  //debugging
  //std::cout << "Got start times for run " << e.run() << " subrun " << e.subRun() << " event " << e.id().event() << "\n";

  //everything else needs to be wrapped in a big for loop and executed once for each tof counter
  for(auto tofc : activeToFCounters) {
    //*****************THE DT DISTRIBUTIONS*******************
    for (int i=0; i<4; i++) {
      for (int j=i+1; j<4; j++) {       
        if(startTimeByChannel[channels[tofc][i]] > 0 && startTimeByChannel[channels[tofc][j]] > 0) {
          double dt = startTimeByChannel[channels[tofc][i]] - startTimeByChannel[channels[tofc][j]];
          hdtByChannel[channels[tofc][i]][channels[tofc][j]]->Fill(dt);
        }//if
      }//for i
    }//for j

    //debugging
    //std::cout << "Made dt distributions for run " << e.run() << " subrun " << e.subRun() << " event " << e.id().event() << "\n";

    //************THE FIRST CHANNEL HIT DISTRIBUTIONS*************
    bool fourHits = true;
    for(int i=0; i<4; i++) {
      fourHits = fourHits && (startTimeByChannel[channels[tofc][i]] > 0);
    }
    int firstBin = 0;
    for(int i=0; i<4; i++) {
      firstBin = startTimeByChannel[channels[tofc][i]] < startTimeByChannel[channels[tofc][firstBin]] ? i : firstBin;
    } //for i
    if(fourHits) {
      hFirstHits[tofc]->Fill(firstBin);
    }

    //debugging
    //std::cout << "Made first channel distributions for run " << e.run() << " subrun " << e.subRun() << " event " << e.id().event() << "\n";

    //******************POSITION RECONSTRUCTION DUMP**************
    //only if hits recorded on all four corners
    if(fourHits) {
      if(dumpPosns) {
        //make histograms
        TH2D* eventPosnDump
          = new TH2D(posnfinder->GetSummaryHist(Form("run%d_subrun%d_event%d_posndump", e.run(), e.subRun(), e.id().event()),
                                                Form("run%d_subrun%d_event%d_posndump", e.run(), e.subRun(), e.id().event()),
                                                startTimeByChannel[channels[tofc][0]],
                                                startTimeByChannel[channels[tofc][1]],
                                                startTimeByChannel[channels[tofc][2]],
                                                startTimeByChannel[channels[tofc][3]],
                                                SCINT_SIDE_LENGTH));
        hPosnRecoDump[tofc].push_back(eventPosnDump);
      }
      
      if(VARY_SIDE) {
        int i = 0; //got to keep track in two places...
        for(double s = MIN_SCINT_SIDE; s < MAX_SCINT_SIDE; s += SCINT_SIDE_STEP) {
          double spread = posnfinder->Get4CornerSpread(startTimeByChannel[channels[tofc][0]],
                                                       startTimeByChannel[channels[tofc][1]],
                                                       startTimeByChannel[channels[tofc][2]],
                                                       startTimeByChannel[channels[tofc][3]],
                                                       s);
          hPosnSpread[tofc][i]->Fill(spread);
      
          posn_t loc = posnfinder->Get4CornerPosition(startTimeByChannel[channels[tofc][0]],
                                                      startTimeByChannel[channels[tofc][1]],
                                                      startTimeByChannel[channels[tofc][2]],
                                                      startTimeByChannel[channels[tofc][3]],
                                                      s);
          if(spread < SPREAD_CUTOFF) {
            hPosn[tofc][i]->Fill(loc.x,loc.y);
            bool phys = posnfinder->IsPhysical(loc,s);
            nPhysicalPosns[tofc][i] += phys ? 1 : 0;
            nNonphysicalPosns[tofc][i] += phys ? 0 : 1;
          }//if
          i++;
      }//for
      }//if
      
      else { //(!VARY_SIDE)
        double t1 = startTimeByChannel[channels[tofc][0]];
        double t2 = startTimeByChannel[channels[tofc][1]];
        double t3 = startTimeByChannel[channels[tofc][2]];
        double t4 = startTimeByChannel[channels[tofc][3]];
        double spread = posnfinder->Get4CornerSpread(t1,t2,t3,t4,SCINT_SIDE_LENGTH);
        hPosnSpread[tofc][0]->Fill(spread);
        
        posn_t loc = posnfinder->Get4CornerPosition(t1,t2,t3,t4,SCINT_SIDE_LENGTH);
        if(spread < SPREAD_CUTOFF) {
          hPosn[tofc][0]->Fill(loc.x,loc.y);
          
          bool phys = posnfinder->IsPhysical(loc,SCINT_SIDE_LENGTH);
          nPhysicalPosns[tofc][0] += phys ? 1 : 0;
          nNonphysicalPosns[tofc][0] += phys ? 0 : 1;
          
        }//if
        else {
          std::cout << Form("High spread for run %d subrun %d event %d with hits at %.2f,%.2f,%.2f,%.2f and reco posn %.2f,%.2f: %.2f\n", e.run(), e.subRun(), e.id().event(), t1, t2, t3, t4, loc.x, loc.y, spread);
          nHighSpreadPosns[tofc]++;
        }//else
      }//else
      
      nFourHitEvents[tofc]++;
    }//if fourhits
    
    
    //keep count  
    nEvents[tofc]++;
    
    
    //debugging
    //std::cout << "Made position histograms for run " << e.run() << " subrun " << e.subRun() << " event " << e.id().event() << "\n";
  }//for tofc
}//analyze(event)

void beamlinereco::ToFPositionRecoAnalysis::endJob() {
  if(dumpPosns) {
    art::TFileDirectory posnDumpDir = tfs->mkdir("Position_reco_dump");
    for(auto tofc : activeToFCounters) {
      //send the posn reco dump to a folder
      TString name = "Posns_" + tofCounterStrings[tofc];
      art::TFileDirectory posnRecoDumpDir = posnDumpDir.mkdir(name.Data());
      for(auto eventHist : hPosnRecoDump[tofc]) {
        posnRecoDumpDir.make<TH2D>(*eventHist);
      }//for eventHist
    }//for tofc
  }//if dumpPosns
  
  //prep file
  std::ofstream dtInfo; 
  dtInfo.open("dtInfo.txt");
  for(auto tofc : activeToFCounters) {
    // Write all ToF dt information to a text file for easy reading and viewing
    for (int i=0; i<4; i++) {
      for (int j=i+1; j<4; j++) {       
        //calculate values
        double mean = hdtByChannel[channels[tofc][i]][channels[tofc][j]]->GetMean(1); //mean of the histogram along the x=1 axis
        double rms =  hdtByChannel[channels[tofc][i]][channels[tofc][j]]->GetRMS(1); //standard deviation
      
        //write to file
        dtInfo << std::string(40,'-') << "\n";
        dtInfo << Form("| %2d | %2d | %7.4f | %6.4f |\n", channels[tofc][i], channels[tofc][j], mean, rms);
      } //for i
    } //for j
  } //for tofc
  dtInfo.close();

  for(auto tofc : activeToFCounters) {
    //print the number of physical and all posns
    if(VARY_SIDE) {
      int i = 0;
      for(double s = MIN_SCINT_SIDE; s < MAX_SCINT_SIDE; s += SCINT_SIDE_STEP) {
        std::cout << Form("In " + tofCounterStrings[tofc] + " There were %d physical posns ", nPhysicalPosns[tofc][i])
                  << Form("for side = %.3f\n", s);
        std::cout << Form("In " + tofCounterStrings[tofc] + " There were %d nonphysical posns ", nNonphysicalPosns[tofc][i])
                  << Form("for side = %.3f\n", s);
        i++;
      }//for
    }//if
    else { //!VARY_SIDE
      std::cout << Form("In " + tofCounterStrings[tofc] + " There were %d physical posns.\n", nPhysicalPosns[tofc][0]);
      std::cout << Form("In " + tofCounterStrings[tofc] + " There were %d nonphysical posns.\n", nNonphysicalPosns[tofc][0]);
    }//else
    
    std::cout << Form("In total, for " + tofCounterStrings[tofc] + " there were %d events with spread above %.3f.\n", nHighSpreadPosns[tofc], SPREAD_CUTOFF);
    std::cout << Form("In total, for " + tofCounterStrings[tofc] + " there were %d events with four hits.\n", nFourHitEvents[tofc]);
    std::cout << Form("In total, for " + tofCounterStrings[tofc] + " there were %d events.\n", nEvents[tofc]);
    
    //redo the labels for the first hits histograms
    //and set a nice range
    for(int i = 0; i < 4; i++) {
      hFirstHits[tofc]->GetXaxis()->SetBinLabel(i+1, Form("ch%d",channels[tofc][i]));
      hFirstHits[tofc]->GetYaxis()->SetRangeUser(0, 1.2*hFirstHits[tofc]->GetBinContent(hFirstHits[tofc]->GetMaximumBin()));
      hFirstHits[tofc]->SetOption("bar1");
    }//for i
    
  }//for tofc
}//endJob