ana::SBNOsc::CRTHitTiming Class Reference

Electron neutrino event selection. More...

Inheritance diagram for ana::SBNOsc::CRTHitTiming:

Public Member Functions
	CRTHitTiming ()
void	Initialize (fhicl::ParameterSet *config=NULL)
void	Finalize ()
int	GetPDG (const std::vector< simb::MCParticle > &mc_particles, int trackID)
bool	ProcessEvent (const gallery::Event &ev, const std::vector< event::Interaction > &truth, std::vector< event::RecoInteraction > &reco)
Public Member Functions inherited from core::SelectionBase
	SelectionBase ()
virtual	~SelectionBase ()
Public Member Functions inherited from core::ProcessorBase
	ProcessorBase ()
virtual	~ProcessorBase ()
virtual void	FillTree ()
virtual void	FillRecoTree ()
virtual void	EventCleanup ()
template<class T >
TBranch *	AddBranch (std::string name, T *obj)
template<class T >
TBranch *	AddRecoBranch (std::string name, T *obj)

Protected Member Functions
void	FillFebMap ()
char	GetAuxDetType (geo::AuxDetGeo const &adgeo)
char	MacToType (int mac)
int	MacToRegion (int mac)
int	MacToAuxDetID (int mac, int chan)
std::array< int, 3 >	GetCRTChannel (int aux_det_id, int aux_det_sens_id)
Protected Member Functions inherited from core::ProcessorBase
virtual void	Initialize (char *config=NULL)
virtual void	Setup (char *config=NULL)
virtual void	Setup (fhicl::ParameterSet *config=NULL)
virtual void	Teardown ()
void	BuildEventTree (gallery::Event &ev)
void	SetupServices (gallery::Event &ev)
void	UpdateSubRuns (gallery::Event &ev)
void	UpdateFileMeta (gallery::Event &ev)

Protected Attributes
std::map< int, std::vector < std::pair< int, int > > >	fFebMap
std::string	fCRTHitTag
bool	fIsICARUS
bool	_verbose
unsigned	event_ind
TH1D *	_true_v_hit_time
TH1D *	_true_v_hit_time_coince
TH1D *	_true_v_hit_time_c
TH1D *	_true_v_hit_time_d
TH1D *	_true_v_hit_time_m
TH1D *	_hitpe_muon
TH1D *	_hitpe_proton
TH1D *	_hitpe_neutron
Protected Attributes inherited from core::ProcessorBase
unsigned long	fEventIndex
	An incrementing index. More...
Experiment	fExperimentID
	Experiment identifier. More...
ProviderManager *	fProviderManager
	Interface for provider access. More...
std::string	fOutputFilename
	The output filename. More...
std::string	fProviderConfig
	A custom provider config fcl file. More...
std::vector< geo::BoxBoundedGeo >	fActiveVolumes
	List of active volumes in configured detector. More...
bool	fWriteTree
	Enable writing of the main tree. More...
TFile *	fOutputFile
	The output ROOT file. More...
TTree *	fTree
	The output ROOT tree. More...
event::Event *	fEvent
	The standard output event data structure. More...
bool	fWriteRecoTree
	Enable writing of the reco tree. More...
TTree *	fRecoTree
	The output reco ROOT tree. More...
event::RecoEvent *	fRecoEvent
	The standard output reco event data structure. More...
TTree *	fSubRunTree
	Subrun output tree. More...
SubRun *	fSubRun
	Standard output subrun structure. More...
TTree *	fFileMetaTree
	File metadata output tree. More...
FileMeta *	fFileMeta
	standard output file metadata structure More...
TParameter< int > *	fExperimentParameter
	Saves value of experiment enum. More...
std::set< std::pair< int, int > >	fSubRunCache
	Cache stored subruns. More...
art::InputTag	fTruthTag
	art tag for MCTruth information More...
art::InputTag	fFluxTag
	art tag for MCFlux information More...
std::vector< art::InputTag >	fWeightTags
	art tag(s) for MCEventWeight information More...
art::InputTag	fMCTrackTag
	art tag for MCTrack More...
art::InputTag	fMCShowerTag
	art tag for MCShower More...
art::InputTag	fMCParticleTag
	art tag for MCParticle More...
std::string	fGeneratorProcess
	process_name of process used to run genie. Used to extract subrun/POT information More...

Additional Inherited Members
Public Attributes inherited from core::ProcessorBase
std::vector < event::RecoInteraction > *	fReco
	Reco interaction list. More...

Detailed Description

Electron neutrino event selection.

Definition at line 41 of file CRTHitTiming.cc.

Constructor & Destructor Documentation

ana::SBNOsc::CRTHitTiming::CRTHitTiming ( )

inline

Constructor.

Definition at line 44 of file CRTHitTiming.cc.

{}

Member Function Documentation

void ana::SBNOsc::CRTHitTiming::FillFebMap ( )

inlineprotected

Definition at line 227 of file CRTHitTiming.cc.

  {
    if(!this->fFebMap.empty())
        return;

    std::string fullFileName;
    cet::search_path searchPath("FW_SEARCH_PATH");
    searchPath.find_file("feb_map.txt",fullFileName);
    std::ifstream fin;
    fin.open(fullFileName,std::ios::in);
    if(fin.good()) std::cout << "opened file 'feb_map.txt' for reading..." << std::endl;
    else
        throw cet::exception("FillFebMap") << "Unable to find/open file 'feb_map.txt'" << std::endl;
    std::vector<std::string> row;
    std::string line, word;
    while(getline(fin,line)) {
        row.clear();
        std::stringstream s(line);
        int mod;
        while (std::getline(s, word, ',')) {
            row.push_back(word);
        }
        mod = std::stoi(row[0]);
        (this->fFebMap)[mod].push_back(std::make_pair(std::stoi(row[1]),std::stoi(row[2])));
        if(row.size()>3)
            (this->fFebMap)[mod].push_back(std::make_pair(std::stoi(row[3]),std::stoi(row[4])));
    }
    std::cout << "filled febMap with " << (this->fFebMap).size() << " entries" << std::endl;
    fin.close();
  } //FillFebMap()

void ana::SBNOsc::CRTHitTiming::Finalize ( )

inlinevirtual

Finalize and write objects to the output file.

Implements core::ProcessorBase.

Definition at line 68 of file CRTHitTiming.cc.

                  {
    fOutputFile->cd();
    _true_v_hit_time_coince->Write();
    _true_v_hit_time->Write();
    _true_v_hit_time_c->Write();
    _true_v_hit_time_m->Write();
    _true_v_hit_time_d->Write();
    _hitpe_muon->Write();
    _hitpe_proton->Write();
    _hitpe_neutron->Write();
  }

char ana::SBNOsc::CRTHitTiming::GetAuxDetType ( geo::AuxDetGeo const &  adgeo )

inlineprotected

Definition at line 258 of file CRTHitTiming.cc.

{
  std::string volName(adgeo.TotalVolume()->GetName());
  if (volName.find("MINOS") != std::string::npos) return 'm';
  if (volName.find("CERN")  != std::string::npos) return 'c';
  if (volName.find("DC")    != std::string::npos) return 'd';

  return 'e';
}//GetAuxDetType()

std::array<int, 3> ana::SBNOsc::CRTHitTiming::GetCRTChannel ( int  aux_det_id, int  aux_det_sens_id  )

inlineprotected

Definition at line 326 of file CRTHitTiming.cc.

                                                                      {
      // const char auxDetType, ) {
      const int adid = aux_det_id; //adsc.AuxDetID(); //CRT module ID number (from gdml)
      const int adsid = aux_det_sens_id; //adsc.AuxDetSensitiveID(); //CRT strip ID number (from gdml)
      const geo::AuxDetGeo& adGeo = fProviderManager->GetAuxDetGeometryProvider()->AuxDet(adid);
      const char auxDetType = GetAuxDetType(adGeo); //CRT module type (c, d, or m)

      int channel0ID=-1, channel1ID=-1;
      int mac5;
      int mac5dual;
      switch (auxDetType){
          case 'c' :
              mac5 =fFebMap[adid][0].first;
              channel0ID = 2 * adsid + 0;
              channel1ID = 2 * adsid + 1;
              if(mac5<107||mac5>230)
                  std::cout << "WARNING: mac5 out of bounds for c-type!" << std::endl;
              if(channel0ID<0 || channel0ID > 31 || channel1ID<0 || channel1ID>31)
                  std::cout << "WARNING: channel out of bounds for c-type!" << std::endl;
              break;
          case 'd' : 
              mac5 = fFebMap[adid][0].first;
              channel0ID = adsid;
              if(mac5<93||mac5>106)
                  std::cout << "WARNING: mac5 out of bounds for d-type!" << std::endl;
              if(channel0ID<0 || channel0ID > 63)
                  std::cout << "WARNING: channel out of bounds for d-type!" << std::endl;
              break;
          case 'm' : 
              mac5 = fFebMap[adid][0].first;
              channel0ID = adsid/2 + 10*(fFebMap[adid][0].second-1);
              if(mac5<1||mac5>92)
                  std::cout << "WARNING: mac5 out of bounds for m-type!" << std::endl;
              if(channel0ID<0 || channel0ID > 31)
                  std::cout << "WARNING: channel out of bounds for m-type!" << std::endl;
              if (fFebMap[adid].size()==2)  {
                  mac5dual = fFebMap[adid][1].first;
                  if(mac5dual<1||mac5dual>92)
                      std::cout << "WARNING: mac5dual out of bounds for m-type!" << std::endl;
              }
              break;
      }
      return {mac5, channel0ID, channel1ID};
  }

int ana::SBNOsc::CRTHitTiming::GetPDG ( const std::vector< simb::MCParticle > &  mc_particles, int  trackID  )

inline

Definition at line 80 of file CRTHitTiming.cc.

                                                                         {
    for (const simb::MCParticle &part: mc_particles) {
      if (part.TrackId() == trackID) return part.PdgCode();
    }
    return -1;
  }

void ana::SBNOsc::CRTHitTiming::Initialize ( fhicl::ParameterSet *  config = NULL )

inlinevirtual

Initialization.

Parameters

config

A configuration, as a FHiCL ParameterSet object

Implements core::ProcessorBase.

Definition at line 51 of file CRTHitTiming.cc.

                                                  {
    fCRTHitTag = config ? config->get<std::string>("CRTHitTag", "crtsimhit") : "crtsimhit";
    fIsICARUS = config ? config->get<bool>("IsICARUS", true) : true;
    event_ind = 0;
    _verbose = false;
    _true_v_hit_time = new TH1D("true_v_hit_time", "true_v_hit_time", 200, -100., 100.);
    _true_v_hit_time_c = new TH1D("true_v_hit_time_c", "true_v_hit_time_c", 200, -100., 100.);
    _true_v_hit_time_d = new TH1D("true_v_hit_time_d", "true_v_hit_time_d", 200, -100., 100.);
    _true_v_hit_time_m = new TH1D("true_v_hit_time_m", "true_v_hit_time_m", 200, -100., 100.);
    _true_v_hit_time_coince = new TH1D("true_v_hit_time_coince", "true_v_hit_time_coince", 200, -100., 100.);
    _hitpe_muon = new TH1D("hitpe_muon", "hitpe_muon", 1000, 0., 1000.);
    _hitpe_proton = new TH1D("hitpe_proton", "hitpe_proton", 1000, 0., 1000.);
    _hitpe_neutron = new TH1D("hitpe_neutron", "hitpe_neutron", 1000, 0., 1000.);
    FillFebMap();
  }

int ana::SBNOsc::CRTHitTiming::MacToAuxDetID ( int  mac, int  chan  )

inlineprotected

Definition at line 300 of file CRTHitTiming.cc.

  {
      char type = MacToType(mac);
      if (type == 'e') return INT_MAX;

      int pos=1;
      if(type=='m')
          pos = chan/10 + 1;

      if((this->fFebMap).empty())
          std::cout << "ERROR in MacToAuxDetID: FEBMap is empty!" << std::endl;

      for(const auto&  p : this->fFebMap) {
          if(p.second[0].first == mac && p.second[0].second==pos)
              return (uint32_t)p.first;
          if(p.second.size()==2)
              if(p.second[1].first==mac && p.second[1].second==pos)
                  return (uint32_t)p.first;
      }


    std::cout << "ERROR in MacToAuxDetID: auxDetID not set!" << std::endl;
    return INT_MAX;
  }

int ana::SBNOsc::CRTHitTiming::MacToRegion ( int  mac )

inlineprotected

Definition at line 279 of file CRTHitTiming.cc.

                          {

      if(mac>=107 && mac<=190) return 30; //top
      if(mac>=191 && mac<=204) return 31; //rim west
      if(mac>=205 && mac<=218) return 32; //rim east
      if(mac>=219 && mac<=224) return 33; //rim south
      if(mac>=225 && mac<=230) return 34; //rim north
      if(            mac<=12 ) return 40; //west side, south stack
      if(mac>=13  && mac<=24 ) return 41; //west side, center stack
      if(mac>=25  && mac<=36 ) return 42; //west side, north stack
      if(mac>=37  && mac<=48 ) return 43; //east side, south stack
      if(mac>=49  && mac<=60 ) return 44; //east side, center stack
      if(mac>=61  && mac<=72 ) return 45; //east side, north stack
      if(mac>=73  && mac<=84 ) return 46; //south
      if(mac>=85  && mac<=92 ) return 47; //north
      if(mac>=93 && mac<=106) return 50; //bottom

      std::cout << "ERROR in MacToRegion: unknown mac address " << mac << std::endl;
      return 0;
  }

char ana::SBNOsc::CRTHitTiming::MacToType ( int  mac )

inlineprotected

Definition at line 268 of file CRTHitTiming.cc.

  {

      int reg = MacToRegion(mac);
      if(reg>=30&&reg<40) return 'c';
      if(reg>=40&&reg<50) return 'm';
      if(reg==50) return 'd';
      std::cout << "ERROR in MacToType: type not set!" << std::endl;
      return 'e';
  }

bool ana::SBNOsc::CRTHitTiming::ProcessEvent ( const gallery::Event &  ev, const std::vector< event::Interaction > &  truth, std::vector< event::RecoInteraction > &  reco  )

inlinevirtual

Process one event.

Parameters

ev	A single event, as a gallery::Event
Reconstructed	interactions

Returns

True to keep event

Implements core::ProcessorBase.

Definition at line 94 of file CRTHitTiming.cc.

                                                                                                                             {
      const std::vector<simb::MCParticle> &mc_particles = *ev.getValidHandle<std::vector<simb::MCParticle>>("largeant");

    if (fIsICARUS) {
      auto const &crt_hits_handle = ev.getValidHandle<std::vector<sbn::crt::CRTHit>>(fCRTHitTag);;
      const std::vector<sbn::crt::CRTHit> &crt_hits = *crt_hits_handle;
      const std::vector<sim::AuxDetSimChannel> &crt_sim_channels = *ev.getValidHandle<std::vector<sim::AuxDetSimChannel>>("largeant");
      art::FindManyP<icarus::crt::CRTData, void> hits_to_data(crt_hits_handle, ev, fCRTHitTag);

      // TODO: fix when ICARUS has the necessary associations
      // collect the mapping of track ID to IDE's
      std::map<int, std::vector<std::pair<std::pair<int, int>, sim::AuxDetIDE>>> tracks_to_ides;
      for (const sim::AuxDetSimChannel &simchannel: crt_sim_channels) {
        for (const sim::AuxDetIDE &ide: simchannel.AuxDetIDEs()) {
          tracks_to_ides[ide.trackID].push_back({{simchannel.AuxDetID(), simchannel.AuxDetSensitiveID()}, ide});
        }
      }

      // Now for each CRT hit, try to figure out the IDE(s) that made it
      for (unsigned i = 0; i < crt_hits.size(); i++) {
        const sbn::crt::CRTHit &hit = crt_hits[i];
        std::vector<art::Ptr<icarus::crt::CRTData>> datas = hits_to_data.at(i);
        assert(datas.size() == 1 || datas.size() == 2);
        bool has_coincedence = datas.size() == 2;
        _verbose = false;

        // map to auxdet
        int auxdetid_0 = MacToAuxDetID(datas[0]->Mac5(), datas[0]->ChanTrig());
        // and to type
        const geo::AuxDetGeo& adGeo = fProviderManager->GetAuxDetGeometryProvider()->AuxDet(auxdetid_0);
        const char auxDetType = GetAuxDetType(adGeo); //CRT module type (c, d, or m)

        if (has_coincedence) {
          int auxdetid_1 = MacToAuxDetID(datas[1]->Mac5(), datas[1]->ChanTrig());
          const geo::AuxDetGeo& adGeo = fProviderManager->GetAuxDetGeometryProvider()->AuxDet(auxdetid_1);
          assert(auxDetType == GetAuxDetType(adGeo));
        }

        double hit_time = (int)hit.ts0_ns;

        if (_verbose) std::cout << "CRT Hit at: " << hit.ts1_ns << " " << (int)hit.ts0_ns << " " << (int)hit.ts0_ns_corr << std::endl;
        if (_verbose) std::cout << "Has Coincedence: " << has_coincedence << std::endl;

        double true_time_sum = 0.;
        unsigned n_true_time = 0;

        bool has_muon = false;
        bool has_proton = false;
        bool has_neutron = false;
        for (unsigned j = 0; j < datas.size(); j++) {
          const icarus::crt::CRTData &data = *datas[j];
          if (_verbose) std::cout << "CRT Data at: " << data.TTrig() << " Mac5: " << data.Mac5() << " chan: " << data.ChanTrig() << std::endl;
          const std::vector<icarus::crt::CRTChannelData> &channel_datas = data.ChanData();

        
          for (const icarus::crt::CRTChannelData &chan: channel_datas) {
            const std::vector<int> track_ids = chan.TrackID();
            if (_verbose) std::cout << "Channel Data at: " << chan.T0() << " chan: " << chan.Channel() << std::endl; 
            for (int track_id: track_ids) {
              int pdg_match = abs(GetPDG(mc_particles, track_id));
              if (abs(pdg_match) == 13) has_muon = true;
              else if (abs(pdg_match) == 2212) has_proton = true;
              else if (abs(pdg_match) == 2112) has_neutron = true;
              if (_verbose) std::cout << "True Particle: " << track_id << std::endl;
              for (const auto &ide_pair: tracks_to_ides.at(track_id)) {
                const sim::AuxDetIDE &ide = ide_pair.second;
                std::pair<int, int> channel = ide_pair.first;
                double true_time = ide.entryT - 1.1e6;
                if (_verbose) std::cout << "True Hit at: " << true_time << " X: "  << ide.entryX << " Y: " << ide.entryY <<" Z: " << ide.entryZ << " Chan: " << channel.first << " " << channel.second << std::endl;
                std::array<int, 3> crt_channel = GetCRTChannel(channel.first, channel.second);
                if (_verbose) std::cout << "Mapped FEB: " << crt_channel[0] << "Mapped channel: " << crt_channel[1] << " " << crt_channel[2] << std::endl;

                if (data.Mac5() == crt_channel[0] && (data.ChanTrig() == crt_channel[1] || data.ChanTrig() == crt_channel[2])) {
                  true_time_sum += true_time;
                  n_true_time += 1;
                }
              }
            }
          }
        }
        if (has_coincedence) _true_v_hit_time_coince->Fill((true_time_sum/n_true_time) - hit_time);
        else {
          std::cout << "Hit time: " << hit_time << " avg true time: " << (true_time_sum/n_true_time) << std::endl;
          _true_v_hit_time->Fill((true_time_sum/n_true_time) - hit_time);
        }
        if (auxDetType == 'm') _true_v_hit_time_m->Fill((true_time_sum/n_true_time) - hit_time);
        else if (auxDetType == 'd') _true_v_hit_time_d->Fill((true_time_sum/n_true_time) - hit_time);
        else if (auxDetType == 'c') _true_v_hit_time_c->Fill((true_time_sum/n_true_time) - hit_time);
        else assert(false);

        double pes = (hit.peshit - 2*63.6) / (70.*2);
        if (has_muon) _hitpe_muon->Fill(pes);
        else if (has_proton) _hitpe_proton->Fill(pes);
        else if (has_neutron) _hitpe_neutron->Fill(pes);
        
      }
      // fOutputFile->cd();
    }
    else {
      auto const &crt_hits_handle = ev.getValidHandle<std::vector<sbn::crt::CRTHit>>(fCRTHitTag);;
      const std::vector<sbn::crt::CRTHit> &crt_hits = *crt_hits_handle;
      const std::vector<sim::AuxDetSimChannel> &crt_sim_channels = *ev.getValidHandle<std::vector<sim::AuxDetSimChannel>>("largeant");
      art::FindManyP<sbnd::crt::CRTData, void> hits_to_data(crt_hits_handle, ev, fCRTHitTag);
      auto const &crt_data_handle = ev.getValidHandle<std::vector<sbnd::crt::CRTData>>("crt");
      for (unsigned i = 0; i < crt_hits.size(); i++) {
        int track_id = -1;
        const sbn::crt::CRTHit &hit = crt_hits[i];
        std::vector<art::Ptr<sbnd::crt::CRTData>> datas = hits_to_data.at(i);
        art::FindManyP<sim::AuxDetIDE, void> sims(datas, ev, "crt");
        bool has_muon = false;
        bool has_proton = false;
        bool has_neutron = false;
        for (unsigned j = 0; j < datas.size(); j++) {
          for (unsigned k = 0; k < sims.at(j).size(); k++) {
            track_id = sims.at(j).at(k)->trackID;
            int pdg_match = GetPDG(mc_particles, track_id);
            if (abs(pdg_match) == 13) has_muon = true;
            else if (abs(pdg_match) == 2212) has_proton = true;
            else if (abs(pdg_match) == 2112) has_neutron = true;
          }
        }
        double pes = hit.peshit;
        if (has_muon) _hitpe_muon->Fill(pes);
        else if (has_proton) _hitpe_proton->Fill(pes);
        else if (has_neutron) _hitpe_neutron->Fill(pes);
      }
    }
    event_ind += 1;

    return false; 
  }

Member Data Documentation

TH1D* ana::SBNOsc::CRTHitTiming::_hitpe_muon

protected

Definition at line 381 of file CRTHitTiming.cc.

TH1D* ana::SBNOsc::CRTHitTiming::_hitpe_neutron

protected

Definition at line 383 of file CRTHitTiming.cc.

TH1D* ana::SBNOsc::CRTHitTiming::_hitpe_proton

protected

Definition at line 382 of file CRTHitTiming.cc.

TH1D* ana::SBNOsc::CRTHitTiming::_true_v_hit_time

protected

Definition at line 376 of file CRTHitTiming.cc.

TH1D* ana::SBNOsc::CRTHitTiming::_true_v_hit_time_c

protected

Definition at line 378 of file CRTHitTiming.cc.

TH1D* ana::SBNOsc::CRTHitTiming::_true_v_hit_time_coince

protected

Definition at line 377 of file CRTHitTiming.cc.

TH1D* ana::SBNOsc::CRTHitTiming::_true_v_hit_time_d

protected

Definition at line 379 of file CRTHitTiming.cc.

TH1D* ana::SBNOsc::CRTHitTiming::_true_v_hit_time_m

protected

Definition at line 380 of file CRTHitTiming.cc.

bool ana::SBNOsc::CRTHitTiming::_verbose

protected

Definition at line 374 of file CRTHitTiming.cc.

unsigned ana::SBNOsc::CRTHitTiming::event_ind

protected

Definition at line 375 of file CRTHitTiming.cc.

std::string ana::SBNOsc::CRTHitTiming::fCRTHitTag

protected

Definition at line 372 of file CRTHitTiming.cc.

std::map<int,std::vector<std::pair<int,int> > > ana::SBNOsc::CRTHitTiming::fFebMap

protected

Definition at line 371 of file CRTHitTiming.cc.

bool ana::SBNOsc::CRTHitTiming::fIsICARUS

protected

Definition at line 373 of file CRTHitTiming.cc.

---

The documentation for this class was generated from the following file:

• CRTHitTiming.cc