Skip to content

Commit

Permalink
Merge pull request #445 from francescopoppi/NewTimingCRTfix
Browse files Browse the repository at this point in the history 
Major update in the CRT Timing structure + some improvements
  • Loading branch information
@rennney
rennney authored Sep 20, 2022
2 parents 256e7fa + 0750077 commit 3982df1
Show file tree
Hide file tree
Showing 5 changed files with 263 additions and 40 deletions.
18 changes: 11 additions & 7 deletions icaruscode/CRT/CRTDataAnalysis_module.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@
#include "larcorealg/Geometry/GeometryCore.h"
#include "larcorealg/Geometry/AuxDetGeometryCore.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
//#include "icaruscode/Decode/DataProducts/ExtraTriggerInfo.h"
#include "sbnobj/Common/Trigger/ExtraTriggerInfo.h"
//#include "icaruscode/CRT/CRTUtils/CRTHitRecoAlg.h"

Expand Down Expand Up @@ -210,7 +209,7 @@ namespace crt {
float fYErrHit; ///< stat error of CRT hit reco Y (cm)
float fZErrHit; ///< stat error of CRT hit reco Z (cm)
uint64_t fT0Hit; ///< hit time w.r.t. PPS
uint64_t fT1Hit; ///< hit time w.r.t. global event time
Long64_t fT1Hit; ///< hit time w.r.t. global trigger
int fHitReg; ///< region code of CRT hit
int fHitSubSys;
int fNHit; ///< number of CRT hits for this event
Expand Down Expand Up @@ -333,7 +332,7 @@ namespace crt {
fHitNtuple->Branch("yErr", &fYErrHit, "yErr/F");
fHitNtuple->Branch("zErr", &fZErrHit, "zErr/F");
fHitNtuple->Branch("t0", &fT0Hit, "t0/l");
fHitNtuple->Branch("t1", &fT1Hit, "t1/l");
fHitNtuple->Branch("t1", &fT1Hit, "t1/L");
fHitNtuple->Branch("region", &fHitReg, "region/I");
// fHitNtuple->Branch("tagger", &ftagger, "tagger/C");
fHitNtuple->Branch("subSys", &fHitSubSys, "subSys/I");
Expand Down Expand Up @@ -413,7 +412,14 @@ namespace crt {
for(int chan=0; chan<32; chan++) {
std::pair<double,double> const chg_cal = fChannelMap->getSideCRTCalibrationMap((int)crtList[febdat_i]->fMac5,chan);
float pe = (crtList[febdat_i]->fAdc[chan]-chg_cal.second)/chg_cal.first;
if(pe<=fPEThresh) continue;
// In order to have Reset TS1 hits in CRTData from Side CRT, we have to explicitly include them
// The current threshold cut (6.5 PE) was applied to filter out noise, but this also filters out
// Reset events which are random trigger around the pedestal. These Reset hits are removed in
// CRT Hit reconstruction. Top CRT has in internal triggering logic and threshold that screens
// from the noise (hence presel = true for all the hits).
// Please revise this in the future if also T0 Reset hits need to be kept in CRTData.
// To do so, include !0crtList[febdat_i]->IsReference_TS0()
if(pe<=fPEThresh && !crtList[febdat_i]->IsReference_TS1()) continue;
presel = true;
}
}else if ( type == 'c' ) {
Expand Down Expand Up @@ -503,7 +509,6 @@ namespace crt {
fHitReg = fCrtutils->AuxDetRegionNameToNum(fCrtutils->MacToRegion(mactmp));
fHitSubSys = fCrtutils->MacToTypeCode(mactmp);


m_gate_crt_diff = m_gate_start_timestamp - hit.ts0_ns;

auto ittmp = hit.pesmap.find(mactmp);
Expand All @@ -513,8 +518,7 @@ namespace crt {
mf::LogError("CRTDataAnalysis") << "mactmp = " << mactmp << std::endl;
mf::LogError("CRTDataAnalysis") << "could not find mac in pesmap!" << std::endl;
continue;
}

}
int chantmp = (*ittmp).second[0].first;

fHitMod = fCrtutils->MacToAuxDetID(mactmp, chantmp);
Expand Down
6 changes: 3 additions & 3 deletions icaruscode/CRT/CRTDecoder/decodercrtjob.fcl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,12 +41,12 @@ physics:
{
daqCRT: @local::crtdaq_icarus
daqTrigger: @local::decodeTriggerV2

triggerconfig: @local::extractTriggerConfig
}

#define the producer and filter modules for this path, order matters,
#filters reject all following items. see lines starting physics.producers below
daq: [ daqTrigger, daqCRT]
daq: [ triggerconfig, daqTrigger, daqCRT]

#define the output stream, there could be more than one if using filters
stream1: [ out1 ]
Expand All @@ -73,5 +73,5 @@ outputs:
compressionLevel: 1
}
}

physics.producers.daqTrigger.DecoderTool.TrigConfigLabel: triggerconfig
#services.IICARUSChannelMap.ChannelMappingTool: @local::ChannelMappingPostGres
10 changes: 7 additions & 3 deletions icaruscode/CRT/CRTSimHitProducer_module.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,19 +142,23 @@ namespace crt {
art::Handle<sbn::ExtraTriggerInfo> trigger_handle;
event.getByLabel( fTriggerLabel, trigger_handle );
if( trigger_handle.isValid() )
m_trigger_timestamp = trigger_handle->triggerTimestamp;
m_trigger_timestamp = trigger_handle->triggerTimestamp;
else
mf::LogError("CRTSimHitProducer") << "No raw::Trigger associated to label: " << fTriggerLabel.label() << "\n" ;
mf::LogError("CRTSimHitProducer") << "No raw::Trigger associated to label: " << fTriggerLabel.label() << "\n" ;
} else{
std::cout << "Trigger Data product " << fTriggerLabel.label() << " not found!\n" ;
}

mf::LogInfo("CRTSimHitProducer")
<<"Number of SiPM hits = "<<crtList.size();
//m_trigger_timestamp = event.getProduct<sbn::ExtraTriggerInfo>(fTriggerLabel).triggerTimestamp;

mf::LogInfo("CRTSimHitProducer")
<<"Number of SiPM hits = "<<crtList.size();

vector<art::Ptr<CRTData>> crtData = hitAlg.PreselectCRTData(crtList, m_trigger_timestamp);

vector<std::pair<CRTHit, vector<int>>> crtHitPairs = hitAlg.CreateCRTHits(crtData);
vector<std::pair<CRTHit, vector<int>>> crtHitPairs = hitAlg.CreateCRTHits(crtData, m_trigger_timestamp);
//vector<std::pair<CRTHit, vector<int>>> crtHitPairs = hitAlg.CreateCRTHits(crtList);

mf::LogInfo("CRTSimHitProducer")
Expand Down
99 changes: 81 additions & 18 deletions icaruscode/CRT/CRTUtils/CRTHitRecoAlg.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ void CRTHitRecoAlg::reconfigure(const Config& config){
}

//---------------------------------------------------------------------------------------
vector<art::Ptr<CRTData>> CRTHitRecoAlg::PreselectCRTData(vector<art::Ptr<CRTData>> crtList, uint64_t trigger_timestamp){
vector<art::Ptr<CRTData>> CRTHitRecoAlg::PreselectCRTData(const vector<art::Ptr<CRTData>> &crtList, uint64_t trigger_timestamp){
if (fVerbose) mf::LogInfo("CRTHitRecoAlg: ") << "In total " << crtList.size() << " CRTData found in an event" << '\n';
vector<art::Ptr<CRTData>> crtdata;
bool presel = false;
Expand Down Expand Up @@ -85,11 +85,12 @@ vector<art::Ptr<CRTData>> CRTHitRecoAlg::PreselectCRTData(vector<art::Ptr<CRTDat

}
mf::LogInfo("CRTHitRecoAlg:") << "Found " << crtdata.size() << " after preselection "<< '\n';
//std::cout<<trigger_timestamp<<std::endl;
return crtdata;
}

//---------------------------------------------------------------------------------------
vector<pair<sbn::crt::CRTHit, vector<int>>> CRTHitRecoAlg::CreateCRTHits(vector<art::Ptr<CRTData>> crtList) {
vector<pair<sbn::crt::CRTHit, vector<int>>> CRTHitRecoAlg::CreateCRTHits(vector<art::Ptr<CRTData>> crtList, uint64_t trigger_timestamp) {

vector<pair<CRTHit, vector<int>>> returnHits;
vector<int> dataIds;
Expand All @@ -103,22 +104,46 @@ vector<pair<sbn::crt::CRTHit, vector<int>>> CRTHitRecoAlg::CreateCRTHits(vector<

// sort by the time
std::sort(crtList.begin(), crtList.end(), compareBytime);


//Load Delays map for Top CRT
CRT_delay_map const FEB_delay_map = LoadFEBMap();
std::vector<std::pair<int,ULong64_t>> CRTReset;
ULong64_t TriggerArray[305]={0};
for (size_t crtdat_i=0; crtdat_i<crtList.size(); crtdat_i++) {
uint8_t mac = crtList[crtdat_i]->fMac5;
int adid = fCrtutils->MacToAuxDetID(mac,0);
char type = fCrtutils->GetAuxDetType(adid);
//For the time being, Only Top CRT delays are loaded, nothing to do for Side CRT yet
if (type == 'c' && crtList[crtdat_i]->IsReference_TS1()) {
ULong64_t Ts0T1ResetEvent = crtList[crtdat_i]->fTs0 + FEB_delay_map.at((int)mac+73).T0_delay - FEB_delay_map.at((int)mac+73).T1_delay;
TriggerArray[(int) mac]=Ts0T1ResetEvent;
CRTReset.emplace_back((int) mac,Ts0T1ResetEvent);
}
}
const int trigger_offset= 60; //Average distance between Global Trigger and Trigger_timestamp (ns)
ULong64_t GlobalTrigger= trigger_timestamp;
if (!CRTReset.empty()) GlobalTrigger = GetMode(CRTReset);
//Add average difference between trigger_timestamp and Global trigger
else GlobalTrigger=GlobalTrigger-trigger_offset;// In this event, the T1 Reset was probably "vetoed" by the T0 Reset
for (int i=0; i<304; i++){
if (TriggerArray[i]==0) TriggerArray[i]=GlobalTrigger;
}
//std::cout<<"Global Trigger "<<GlobalTrigger<<std::endl;
//loop over time-ordered CRTData
for (size_t febdat_i=0; febdat_i<crtList.size(); febdat_i++) {

uint8_t mac = crtList[febdat_i]->fMac5;
int adid = fCrtutils->MacToAuxDetID(mac,0); //module ID

string region = fCrtutils->GetAuxDetRegion(adid);
char type = fCrtutils->GetAuxDetType(adid);
CRTHit hit;

dataIds.clear();

//CERN modules (intramodule coincidence)
if ( type == 'c' ) {
hit = MakeTopHit(crtList[febdat_i]);
hit = MakeTopHit(crtList[febdat_i], TriggerArray);
if(IsEmptyHit(hit))
nMissC++;
else {
Expand Down Expand Up @@ -211,7 +236,7 @@ vector<pair<sbn::crt::CRTHit, vector<int>>> CRTHitRecoAlg::CreateCRTHits(vector<
mf::LogInfo("CRTHitRecoAlg: ") << "attempting to produce MINOS hit from " << coinData.size()
<< " data products..." << '\n';

CRTHit hit = MakeSideHit(coinData);
CRTHit hit = MakeSideHit(coinData, TriggerArray);

if(IsEmptyHit(hit)){
unusedDataIndex.push_back(indices[index_i]);
Expand Down Expand Up @@ -257,7 +282,7 @@ vector<pair<sbn::crt::CRTHit, vector<int>>> CRTHitRecoAlg::CreateCRTHits(vector<
//--------------------------------------------------------------------------------------------
// Function to make filling a CRTHit a bit faster
sbn::crt::CRTHit CRTHitRecoAlg::FillCRTHit(vector<uint8_t> tfeb_id, map<uint8_t,vector<pair<int,float>>> tpesmap,
float peshit, uint64_t time0, uint64_t time1, int plane,
float peshit, uint64_t time0, Long64_t time1, int plane,
double x, double ex, double y, double ey, double z, double ez, string tagger){
CRTHit crtHit;
crtHit.feb_id = tfeb_id;
Expand All @@ -266,7 +291,7 @@ sbn::crt::CRTHit CRTHitRecoAlg::FillCRTHit(vector<uint8_t> tfeb_id, map<uint8_t,
crtHit.ts0_s_corr = time0 / 1'000'000'000;
crtHit.ts0_ns = time0 % 1'000'000'000;
crtHit.ts0_ns_corr = time0;
crtHit.ts1_ns = time1 % 1'000'000'000;
crtHit.ts1_ns = time1 /*% 1'000'000'000*/; //TODO: Update the CRTHit data product /sbnobj/common/CRT . Discussion with SBND people needed
crtHit.ts0_s = time0 / 1'000'000'000;
crtHit.plane = plane;
crtHit.x_pos = x;
Expand All @@ -282,7 +307,7 @@ sbn::crt::CRTHit CRTHitRecoAlg::FillCRTHit(vector<uint8_t> tfeb_id, map<uint8_t,
} // CRTHitRecoAlg::FillCRTHit()

//------------------------------------------------------------------------------------------
sbn::crt::CRTHit CRTHitRecoAlg::MakeTopHit(art::Ptr<CRTData> data){
sbn::crt::CRTHit CRTHitRecoAlg::MakeTopHit(art::Ptr<CRTData> data, ULong64_t GlobalTrigger[305]){

uint8_t mac = data->fMac5;
if(fCrtutils->MacToType(mac)!='c')
Expand All @@ -301,12 +326,11 @@ sbn::crt::CRTHit CRTHitRecoAlg::MakeTopHit(art::Ptr<CRTData> data){
TVector3 postrig;
bool findx = false, findz = false;
int maxx=0, maxz=0;

double sum=0;
for(int chan=0; chan<32; chan++) {
// std::cout<<"Top Gain: "<<ftopGain<<" Top Pedestal: "<<ftopPed<< '\n';
sum=sum+data->fAdc[chan];
float pe = (data->fAdc[chan]-ftopPed)/ftopGain;
// float pe = (data->fAdc[chan]-fQPed)/fQSlope;
// if(pe<=fPEThresh) continue;
// if(pe<=fPEThresh) continue;
nabove++;
int adsid = fCrtutils->ChannelToAuxDetSensitiveID(mac,chan);
petot += pe;
Expand Down Expand Up @@ -368,7 +392,7 @@ sbn::crt::CRTHit CRTHitRecoAlg::MakeTopHit(art::Ptr<CRTData> data){

auto const& adsGeo = adGeo.SensitiveVolume(adsid_max); //trigger strip
uint64_t thit = data->fTs0;
uint64_t thit1 = data->fTs1;
Long64_t thit1 = data->fTs1;

if(adsid_max<8){
thit -= (uint64_t)round(abs((92+hitpos.X())*fPropDelay));
Expand All @@ -388,6 +412,10 @@ sbn::crt::CRTHit CRTHitRecoAlg::MakeTopHit(art::Ptr<CRTData> data){
hitpointerr[0] = adsGeo.HalfWidth1()*2/sqrt(12);
hitpointerr[1] = adGeo.HalfHeight();
hitpointerr[2] = adsGeo.HalfWidth1()*2/sqrt(12);
thit1 = (Long64_t)(thit-GlobalTrigger[(int)mac+73]);

//Remove T1 Reset event not correctly flagged, remove T1 reset events, remove T0 reset events
if((sum<10000 && thit1<2'001'000 && thit1>2'000'000)||data->IsReference_TS1() || data->IsReference_TS0()) return FillCRTHit({},{},0,0,0,0,0,0,0,0,0,0,"");

CRTHit hit = FillCRTHit({mac},pesmap,petot,thit,thit1,plane,hitpoint[0],hitpointerr[0],
hitpoint[1],hitpointerr[1],hitpoint[2],hitpointerr[2],region);
Expand Down Expand Up @@ -463,7 +491,7 @@ sbn::crt::CRTHit CRTHitRecoAlg::MakeBottomHit(art::Ptr<CRTData> data){
} // CRTHitRecoAlg::MakeBottomHit

//-----------------------------------------------------------------------------------
sbn::crt::CRTHit CRTHitRecoAlg::MakeSideHit(vector<art::Ptr<CRTData>> coinData) {
sbn::crt::CRTHit CRTHitRecoAlg::MakeSideHit(vector<art::Ptr<CRTData>> coinData, ULong64_t GlobalTrigger[305]) {

vector<uint8_t> macs;
map< uint8_t, vector< pair<int,float> > > pesmap;
Expand Down Expand Up @@ -1031,10 +1059,11 @@ sbn::crt::CRTHit CRTHitRecoAlg::MakeSideHit(vector<art::Ptr<CRTData>> coinData)
hitpointerr[1] = adsGeo.HalfWidth1()*2/sqrt(12);
hitpointerr[2] = (zmax-zmin)/sqrt(12);
}

Long64_t thit1=(Long64_t)(thit-GlobalTrigger[(int)macs.at(0)]);


//generate hit
CRTHit hit = FillCRTHit(macs,pesmap,petot,thit,t1hit,plane,hitpoint[0],hitpointerr[0],
CRTHit hit = FillCRTHit(macs,pesmap,petot,thit,thit1,plane,hitpoint[0],hitpointerr[0],
hitpoint[1],hitpointerr[1],hitpoint[2],hitpointerr[2],region);

return hit;
Expand All @@ -1051,3 +1080,37 @@ bool CRTHitRecoAlg::IsEmptyHit(CRTHit hit) {

return false;
}

//-----------------------------------------------------------------------------
ULong64_t icarus::crt::GetMode(std::vector<std::pair<int, ULong64_t>> vector) {

sort(vector.begin(), vector.end(), icarus::crt::sortbytime);

int modecounter = 0;
int isnewmodecounter = 0;
ULong64_t Mode = 0;
ULong64_t isnewMode = 0;
bool isFirst = true;
for (auto i : vector) {
if (!isFirst) {
if (i.second == Mode) modecounter++;
else if (i.second !=isnewMode) {
isnewMode = i.second;
isnewmodecounter = 1;
}
else if (i.second == isnewMode) {
isnewmodecounter++;
if (isnewmodecounter > modecounter) {
Mode = isnewMode;
modecounter = isnewmodecounter;
}
}
}
else {
isFirst = false;
Mode = i.second;
modecounter++;
}
}
return Mode;
}
Loading

0 comments on commit 3982df1

Please sign in to comment.