@@ -91,6 +91,8 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
9191 const std::map<CRTBackTrackerAlg::Category, bool > &spacePointRecoStatusMap = {},
9292 const std::map<int , std::pair<bool , bool >> &trackRecoStatusMap = {});
9393
94+ double SPTimeDelta (std::vector<std::pair<int , double >> &time_set);
95+
9496private:
9597
9698 CRTGeoAlg fCRTGeoAlg ;
@@ -163,8 +165,8 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
163165 // strip hit to select the strip which has ADC above threshold
164166 std::vector<uint32_t > _sh_channel;
165167 std::vector<int16_t > _sh_tagger;
166- std::vector<int64_t > _sh_ts0;
167- std::vector<int64_t > _sh_ts1;
168+ std::vector<double > _sh_ts0;
169+ std::vector<double > _sh_ts1;
168170 std::vector<uint32_t > _sh_unixs;
169171 std::vector<double > _sh_pos;
170172 std::vector<double > _sh_err;
@@ -180,16 +182,16 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
180182 std::vector<double > _sh_truth_time;
181183
182184 // cluster from x-y coincidence for CRTSpacePoint , this is what we normally call a CRT hit
183- std::vector<int64_t > _cl_ts0;
184- std::vector<int64_t > _cl_ts1;
185+ std::vector<double > _cl_ts0;
186+ std::vector<double > _cl_ts1;
185187 std::vector<uint32_t > _cl_unixs;
186188 std::vector<uint16_t > _cl_nhits;
187189 std::vector<int16_t > _cl_tagger;
188190 std::vector<uint8_t > _cl_composition;
189191 std::vector<std::vector<uint32_t >> _cl_channel_set;
190192 std::vector<std::vector<uint16_t >> _cl_adc_set;
191- std::vector<std::vector<int64_t >> _cl_sh_ts0_set; // ! To store t0 from x-y coincidences
192- std::vector<std::vector<int64_t >> _cl_sh_ts1_set; // ! To store t1 from x-y coincidences
193+ std::vector<std::vector<double >> _cl_sh_ts0_set; // ! To store t0 from x-y coincidences
194+ std::vector<std::vector<double >> _cl_sh_ts1_set; // ! To store t1 from x-y coincidences
193195 std::vector<std::vector<uint16_t >> _cl_sh_feb_mac5_set; // ! MAC5 addresses of StripHit FEBs
194196 std::vector<std::vector<double >> _cl_sh_time_walk_set; // ! Time walk correction
195197 std::vector<std::vector<double >> _cl_sh_prop_delay_set; // ! Light propagation correction
@@ -221,6 +223,8 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
221223 std::vector<double > _cl_sp_ets0;
222224 std::vector<double > _cl_sp_ts1;
223225 std::vector<double > _cl_sp_ets1;
226+ std::vector<double > _cl_sp_dts0;
227+ std::vector<double > _cl_sp_dts1;
224228 std::vector<bool > _cl_sp_complete;
225229
226230 // backtrack truth information from reco level
@@ -435,8 +439,8 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
435439
436440 fTree ->Branch (" sh_channel" , " std::vector<uint32_t>" , &_sh_channel);
437441 fTree ->Branch (" sh_tagger" , " std::vector<int16_t>" , &_sh_tagger);
438- fTree ->Branch (" sh_ts0" , " std::vector<int64_t >" , &_sh_ts0);
439- fTree ->Branch (" sh_ts1" , " std::vector<int64_t >" , &_sh_ts1);
442+ fTree ->Branch (" sh_ts0" , " std::vector<double >" , &_sh_ts0);
443+ fTree ->Branch (" sh_ts1" , " std::vector<double >" , &_sh_ts1);
440444 fTree ->Branch (" sh_unixs" , " std::vector<uint32_t>" , &_sh_unixs);
441445 fTree ->Branch (" sh_pos" , " std::vector<double>" , &_sh_pos);
442446 fTree ->Branch (" sh_err" , " std::vector<double>" , &_sh_err);
@@ -454,16 +458,16 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
454458 fTree ->Branch (" sh_truth_time" , " std::vector<double>" , &_sh_truth_time);
455459 }
456460
457- fTree ->Branch (" cl_ts0" , " std::vector<int64_t >" , &_cl_ts0);
458- fTree ->Branch (" cl_ts1" , " std::vector<int64_t >" , &_cl_ts1);
461+ fTree ->Branch (" cl_ts0" , " std::vector<double >" , &_cl_ts0);
462+ fTree ->Branch (" cl_ts1" , " std::vector<double >" , &_cl_ts1);
459463 fTree ->Branch (" cl_unixs" , " std::vector<uint32_t>" , &_cl_unixs);
460464 fTree ->Branch (" cl_nhits" , " std::vector<uint16_t>" , &_cl_nhits);
461465 fTree ->Branch (" cl_tagger" , " std::vector<int16_t>" , &_cl_tagger);
462466 fTree ->Branch (" cl_composition" , " std::vector<uint8_t>" , &_cl_composition);
463467 fTree ->Branch (" cl_channel_set" , " std::vector<std::vector<uint32_t>>" , &_cl_channel_set);
464468 fTree ->Branch (" cl_adc_set" , " std::vector<std::vector<uint16_t>>" , &_cl_adc_set);
465- fTree ->Branch (" cl_sh_ts0_set" , " std::vector<std::vector<int64_t >>" , &_cl_sh_ts0_set);
466- fTree ->Branch (" cl_sh_ts1_set" , " std::vector<std::vector<int64_t >>" , &_cl_sh_ts1_set);
469+ fTree ->Branch (" cl_sh_ts0_set" , " std::vector<std::vector<double >>" , &_cl_sh_ts0_set);
470+ fTree ->Branch (" cl_sh_ts1_set" , " std::vector<std::vector<double >>" , &_cl_sh_ts1_set);
467471 fTree ->Branch (" cl_sh_feb_mac5_set" , " std::vector<std::vector<uint16_t>>" , &_cl_sh_feb_mac5_set);
468472 fTree ->Branch (" cl_sh_time_walk_set" , " std::vector<std::vector<double>>" , &_cl_sh_time_walk_set);
469473 fTree ->Branch (" cl_sh_prop_delay_set" , " std::vector<std::vector<double>>" , &_cl_sh_prop_delay_set);
@@ -498,6 +502,8 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
498502 fTree ->Branch (" cl_sp_ets0" , " std::vector<double>" , &_cl_sp_ets0);
499503 fTree ->Branch (" cl_sp_ts1" , " std::vector<double>" , &_cl_sp_ts1);
500504 fTree ->Branch (" cl_sp_ets1" , " std::vector<double>" , &_cl_sp_ets1);
505+ fTree ->Branch (" cl_sp_dts0" , " std::vector<double>" , &_cl_sp_dts0);
506+ fTree ->Branch (" cl_sp_dts1" , " std::vector<double>" , &_cl_sp_dts1);
501507 fTree ->Branch (" cl_sp_complete" , " std::vector<bool>" , &_cl_sp_complete);
502508
503509 if (!fDataMode && fTruthMatch )
@@ -1233,6 +1239,8 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
12331239 _cl_sp_ets0.resize (nClusters);
12341240 _cl_sp_ts1.resize (nClusters);
12351241 _cl_sp_ets1.resize (nClusters);
1242+ _cl_sp_dts0.resize (nClusters);
1243+ _cl_sp_dts1.resize (nClusters);
12361244 _cl_sp_complete.resize (nClusters);
12371245
12381246 for (unsigned i = 0 ; i < nClusters; ++i)
@@ -1251,12 +1259,14 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
12511259 const auto striphits = clustersToStripHits.at (cluster.key ());
12521260 _cl_channel_set[i].resize (_cl_nhits[i]);
12531261 _cl_adc_set[i].resize (2 * _cl_nhits[i]);
1254- _cl_sh_ts0_set[i].resize (_cl_nhits[i] );
1255- _cl_sh_ts1_set[i].resize (_cl_nhits[i] );
1262+ _cl_sh_ts0_set[i].resize (_cl_nhits[i]);
1263+ _cl_sh_ts1_set[i].resize (_cl_nhits[i]);
12561264 _cl_sh_feb_mac5_set[i].resize (_cl_nhits[i]);
12571265 _cl_sh_time_walk_set[i].resize (_cl_nhits[i]);
12581266 _cl_sh_prop_delay_set[i].resize (_cl_nhits[i]);
12591267
1268+ std::vector<std::pair<int , double >> ts0_set, ts1_set;
1269+
12601270 for (unsigned ii = 0 ; ii < _cl_nhits[i]; ++ii)
12611271 {
12621272 const auto striphit = striphits[ii];
@@ -1281,24 +1291,28 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
12811291 * Because the Ts0(), Ts1() getters invoked in _cl_sp_ts*, _cl_sh_ts*_set are raw T0/1
12821292 * counters, the time walk and propagation delay are saved as explicit branches here.
12831293 */
1284- if (spacepoints.size () == 1 ) { // need unique position of spacepoint
1285- double pe0 = fCRTGeoAlg .GetSiPM ( striphit->Channel () ).gain * striphit->ADC1 ();
1286- double pe1 = fCRTGeoAlg .GetSiPM ( striphit->Channel () + 1 ).gain * striphit->ADC2 ();
1287- double pe = pe0 + pe1;
1288-
1289- double dist = fCRTGeoAlg .DistanceDownStrip ( spacepoints[0 ]->Pos (), striphit->Channel () );
1294+ if (spacepoints.size () == 1 )
1295+ {
1296+ double pe0 = fCRTGeoAlg .GetSiPM ( striphit->Channel () ).gain * striphit->ADC1 ();
1297+ double pe1 = fCRTGeoAlg .GetSiPM ( striphit->Channel () + 1 ).gain * striphit->ADC2 ();
1298+ double pe = pe0 + pe1;
12901299
1291- double corr = std::pow ( dist - fPEAttenuation , 2.0 ) / std::pow ( fPEAttenuation , 2.0 );
1292- double tw_pe = pe * corr;
1300+ double dist = fCRTGeoAlg .DistanceDownStrip ( spacepoints[0 ]->Pos (), striphit->Channel () );
12931301
1294- _cl_sh_time_walk_set[i][ii] = fTimeWalkNorm * std::exp ( - fTimeWalkScale * tw_pe );
1295- _cl_sh_prop_delay_set[i][ii] = fPropDelay * dist ;
1302+ double corr = std::pow ( dist - fPEAttenuation , 2.0 ) / std::pow ( fPEAttenuation , 2.0 );
1303+ double tw_pe = pe * corr ;
12961304
1297- } else { // fill with nonsense
1298- _cl_sh_time_walk_set[i][ii] = -999999 .;
1299- _cl_sh_prop_delay_set[i][ii] = -999999 .;
1300- }
1305+ _cl_sh_time_walk_set[i][ii] = fTimeWalkNorm * std::exp ( -fTimeWalkScale * tw_pe );
1306+ _cl_sh_prop_delay_set[i][ii] = fPropDelay * dist;
13011307
1308+ ts0_set.push_back ({_cl_sh_feb_mac5_set[i][ii], _cl_sh_ts0_set[i][ii] - _cl_sh_time_walk_set[i][ii] - _cl_sh_prop_delay_set[i][ii]});
1309+ ts1_set.push_back ({_cl_sh_feb_mac5_set[i][ii], _cl_sh_ts1_set[i][ii] - _cl_sh_time_walk_set[i][ii] - _cl_sh_prop_delay_set[i][ii]});
1310+ }
1311+ else
1312+ {
1313+ _cl_sh_time_walk_set[i][ii] = -999999 .;
1314+ _cl_sh_prop_delay_set[i][ii] = -999999 .;
1315+ }
13021316 }
13031317
13041318 if (!fDataMode && fTruthMatch )
@@ -1338,6 +1352,8 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
13381352 _cl_sp_ets0[i] = spacepoint->Ts0Err ();
13391353 _cl_sp_ts1[i] = spacepoint->Ts1 ();
13401354 _cl_sp_ets1[i] = spacepoint->Ts1Err ();
1355+ _cl_sp_dts0[i] = SPTimeDelta (ts0_set);
1356+ _cl_sp_dts1[i] = SPTimeDelta (ts1_set);
13411357 _cl_sp_complete[i] = spacepoint->Complete ();
13421358 }
13431359 else
@@ -1354,6 +1370,8 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
13541370 _cl_sp_ets0[i] = -999999 .;
13551371 _cl_sp_ts1[i] = -999999 .;
13561372 _cl_sp_ets1[i] = -999999 .;
1373+ _cl_sp_dts0[i] = -999999 .;
1374+ _cl_sp_dts1[i] = -999999 .;
13571375 _cl_sp_complete[i] = false ;
13581376 }
13591377 }
@@ -1858,4 +1876,22 @@ void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::
18581876 _tpc_tr_score.resize (nActualTracks);
18591877}
18601878
1879+ double sbnd::crt::CRTAnalysis::SPTimeDelta (std::vector<std::pair<int , double >> &time_set)
1880+ {
1881+ std::sort (time_set.begin (), time_set.end (), [](auto &a, auto &b)
1882+ { return a.first < b.first ; });
1883+
1884+ double sum = 0 ;
1885+
1886+ for (unsigned int i = 0 ; i < time_set.size (); ++i)
1887+ {
1888+ for (unsigned int ii = i + 1 ; ii < time_set.size (); ++ii)
1889+ sum += (time_set[i].second - time_set[ii].second );
1890+ }
1891+
1892+ sum /= (time_set.size () - 1 );
1893+
1894+ return sum;
1895+ }
1896+
18611897DEFINE_ART_MODULE (sbnd::crt::CRTAnalysis)
0 commit comments