15 #include "Utilities/func/MathUtil.h" 19 #include "NovaDAQConventions/DAQConventions.h" 20 #include "DAQChannelMap/DAQChannelMap.h" 32 #include "Utilities/AssociationUtil.h" 75 std::map<int, std::vector<art::Ptr<caldp::PCHit> > >
fDCMHitMap;
149 produces< std::vector<caldp::DCMStat> >();
150 produces< art::Assns<caldp::DCMStat, rb::Track> >();
184 fOutTree = tfs->
make<TTree>(
"TimingCalib",
"Timing Calibaration");
232 switch(geom->
DetId()){
243 assert(0 &&
"Unknown detector");
257 std::unique_ptr< std::vector<caldp::DCMStat> > dcmstats(
new std::vector<caldp::DCMStat>);
272 std::vector< art::Ptr<rb::Track> > tracks;
283 for (
unsigned int i=0;
i<tracks.size(); ++
i){
291 std::vector<art::Ptr<caldp::TCTrack> > tct = fmtc.at(
i);
293 if (!(fmtc.isValid()))
continue;
295 if (tct.size() != 1)
continue;
324 for(
unsigned int j=0;
j<tct[0]->NPCHit(); ++
j){
333 std::vector<art::Ptr<caldp::PCHit> > dcmcells = itr->second;
338 std::vector<float> cellnumtemp;
339 std::vector<float> planenumtemp;
340 std::vector<float> dcmtimetemp;
341 std::vector<float> dcmstimetemp;
342 std::vector<float> dcmrdtemp;
343 std::vector<float> dcmpetemp;
344 std::vector<float> dcmpltemp;
346 std::vector<double>
sigma;
347 for(
unsigned int k=0; k<dcmcells.size(); ++k){
348 double tns = dcmcells[k]->TNS();
349 double rd = dcmcells[k]->ReadDist();
350 double pl = dcmcells[k]->FlightLen();
351 if (
fCorrectTimewalk) tns += tct[0]->ReadoutDistCorrection(dcmcells[k]->ReadDist(), detId);
355 if (
fAdjustSpeed) tmpFiberSpeed = tct[0]->CalcFiberVelocity(rd, detId);
356 dcmstimetemp.push_back(tns - pl/sCmPerNsec - rd/tmpFiberSpeed);
357 dcmtimetemp.push_back(tns);
358 dcmrdtemp.push_back(rd);
359 dcmpetemp.push_back(dcmcells[k]->PE());
360 dcmpltemp.push_back(pl);
361 cellnumtemp.push_back(dcmcells[k]->Cell());
362 planenumtemp.push_back(dcmcells[k]->Plane());
363 double stdev = tct[0]->TNSUncertainty(dcmcells[k]->PE(), dcmcells[k]->GoodTime(), detId);
364 sigma.push_back(stdev);
365 double w = 1.0/(stdev*
stdev);
366 weights.push_back(w);
372 for(
unsigned int k=0; k<dcmstimetemp.size(); ++k){
373 st += dcmstimetemp[k]*weights[k];
374 weight += weights[k];
378 unsigned int maxpos = 0;
379 for(
unsigned int k=0; k<dcmstimetemp.size(); ++k){
380 float recresid = st - dcmstimetemp[k];
381 if (
fabs(recresid/sigma[k]) >
fabs(maxresid)){
382 maxresid =
fabs(recresid/sigma[k]);
386 if (dcmcells.size() <= 2)
break;
389 dcmcells.erase(dcmcells.begin()+maxpos);
390 cellnumtemp.erase(cellnumtemp.begin()+maxpos);
391 planenumtemp.erase(planenumtemp.begin()+maxpos);
392 dcmtimetemp.erase(dcmtimetemp.begin()+maxpos);
393 dcmstimetemp.erase(dcmstimetemp.begin()+maxpos);
394 dcmrdtemp.erase(dcmrdtemp.begin()+maxpos);
395 dcmpetemp.erase(dcmpetemp.begin()+maxpos);
396 dcmpltemp.erase(dcmpltemp.begin()+maxpos);
397 weights.erase(weights.begin()+maxpos);
398 sigma.erase(sigma.begin()+maxpos);
403 dcmnum.push_back(itr->first);
404 dcmhits.push_back(dcmstimetemp.size());
407 dcmTime.
dcmNum = itr->first;
408 dcmTime.
dcmHits = dcmcells.size();
410 for(
unsigned int k=0; k<dcmstimetemp.size(); ++k){
411 double recresid = st - dcmstimetemp[k];
413 sum += weights[k]*(dcmstimetemp[k] - st)*(dcmstimetemp[k]-st);
414 cellnum.push_back(cellnumtemp[k]);
415 planenum.push_back(planenumtemp[k]);
416 dcmtime.push_back(dcmtimetemp[k]);
417 dcmstime.push_back(dcmstimetemp[k]);
418 dcmrd.push_back(dcmrdtemp[k]);
419 dcmpe.push_back(dcmpetemp[k]);
420 dcmpl.push_back(dcmpltemp[k]);
423 dcms.
dcmList.push_back(dcmTime);
431 dcmstats->push_back(dcms);
447 evt.
put(std::move(dcmstats));
448 evt.
put(std::move(dcmAssns));
SubRunNumber_t subRun() const
back track the reconstruction to the simulation
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
int fMinDCMHits
minimum number of hits on a dcm to keep results
std::vector< float > dcmrd
distance to readout
virtual unsigned int getNumberOfDCM(DiBlock_TYPE dbt) const =0
How many DCMs does this diblock have?
std::vector< float > planenum
plane number
float sigma
uncertainty in average hit time
fvar< T > fabs(const fvar< T > &x)
int fMinDCM
minimum number of good dcms on a track
bool fDebug
Add exta debugging info to ntuple.
TimingCalibration(fhicl::ParameterSet const &p)
std::vector< float > dcmtime
uncorrected times of hits
constexpr std::uint32_t timeHigh() const
std::vector< float > cellnum
cell number
Vertical planes which measure X.
bool fCorrectTimewalk
adjust hit times based on timewalk effect for distance to readout
void reconfigure(const fhicl::ParameterSet &pset)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
void beginRun(art::Run &run) override
Horizontal planes which measure Y.
virtual ~TimingCalibration()
std::map< int, std::vector< art::Ptr< caldp::PCHit > > > fDCMHitMap
fhicl::ParameterSet fPSetTB
ProductID put(std::unique_ptr< PROD > &&product)
std::string fTCTrackLabel
Where to find calibration tracks.
std::string fTrackLabel
Where to find reco Tracks.
int dcmHits
number of hits in dcm
fhicl::ParameterSet fPSetFD
Far Detector at Ash River, MN.
std::vector< float > dcmsimtime
simultaneous hit time in dcm
static DAQChannelMap * getInstance(int detID)
int fRemoveThreshold
number of sigma away from mean to remove hits
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Near Detector in the NuMI cavern.
float dcmTime
weighted average time of hits in dcm, correcting for time-of-flight and path lenght ...
std::vector< float > dcmstime
times with distance to readout, time-of-flight, timewalk corrections
std::vector< float > dcmnum
dcm index
std::vector< float > dcmpl
path length along track
Identifier for diblocks using a 32/32 configuration.
int Diblock() const
Return diblock value.
EventNumber_t event() const
T * make(ARGS...args) const
static const double sCmPerNsec
bool fRemoveOutlier
remove outlier hits in DCMs based on timing resolution
double GetTimingOffset(unsigned int const &plane, unsigned int const &cell, bool const &isData)
Get the timing offset for a given plane, cell. Useful downstream to check calibration.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
fhicl::ParameterSet fPSetND
std::vector< float > dcmsimresid
difference between simultaneous hit time and average
assert(nhit_max >=nhit_nbins)
std::vector< float > dcmhits
number of hits in dcm
bool fAdjustSpeed
adjust fiberspeed based on distance to readout instead of assuming flat speed
cmap::CMap class source code
std::vector< DCMTime > dcmList
int DCM() const
Return dcm value.
float fFiberSpeed
fixed fiberspeed
TVector3 Stop() const
Position of the final trajectory point.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
void produce(art::Event &e) override
std::vector< float > dcmpe
number of photoelectrons
Encapsulate the geometry of one entire detector (near, far, ndos)
int dcmNum
Number of the dcm.