25 #include "cetlib_except/exception.h" 33 #include "NovaDAQConventions/DAQConventions.h" 68 void FillPeCorrHits(
art::Handle<std::vector<caldp::PCHit> > &cells,
73 void CalculateMean(std::vector<double> &chanvector,
double &fMean,
double &fMeanErr,
int &fNHits,
double &fRMS);
76 void DrawMean(std::vector<double> &chanvector,TH1F *hpecorr);
89 std::map< daqchannelmap::dchan, std::vector<double> >
mMeanPeXY;
90 std::map< daqchannelmap::dchan, std::vector<double> >
mMeanPeZ;
91 std::map< daqchannelmap::dchan, std::vector<double> >
mMeanPeAvg;
115 produces< std::vector<caldp::DriftResponse >,
art::InRun >();
119 DriftResponseCalc::~DriftResponseCalc()
130 fDrawMeans = pset.
get<
bool >(
"DrawMeans");
132 fDiBlockNum = pset.
get<
int >(
"DiBlockNum");
133 fDCMNum = pset.
get<
int >(
"DCMNum");
168 double attencorr = 1.0;
175 for (i = 0; i < cells->size(); ++
i){
180 if (cell->
Path()<=0)
continue;
183 chanid = GetChanID(cell->
Plane(),cell->
Cell(),detid);
186 if (mMeanPe.find(chanid) != mMeanPe.end() ) {
194 if (!isData) cellhit.
SetMC();
200 if (attencorr<=0)
continue;
203 fPeCorr = cell->
PE()*attencorr/cell->
Path();
204 mMeanPe[chanid].push_back(fPeCorr);
212 void DriftResponseCalc::CalculateMean(std::vector<double> &chanvector,
double &fMean,
double &fMeanErr,
int &fNHits,
double &fRMS){
222 std::sort(chanvector.begin(),chanvector.end());
225 fNHits = (
int)(chanvector.size()*9/10);
228 if (fNHits<=1)
return;
233 for (i=0; i<fNHits; ++
i) dsum += chanvector[i];
236 double dnhits = (double)fNHits;
243 for (i=0; i<fNHits; ++
i) dsumvar2 += (chanvector[i]-fMean)*(chanvector[
i]-fMean);
246 double dstdvar =
sqrt(dsumvar2/(dnhits - 1));
249 fMeanErr = dstdvar/
sqrt(dnhits);
253 void DriftResponseCalc::DrawMean(std::vector<double> &chanvector, TH1F *hpecorr){
256 std::sort(chanvector.begin(),chanvector.end());
257 int nhits = (
int)(chanvector.size()*9/10);
260 for (
int i=0;
i<nhits; ++
i) hpecorr->Fill(chanvector[
i]);
281 if (mMeanPeXY.size()==0){
293 unsigned int planeid;
303 const std::set<unsigned int> planeSet =
307 std::set<unsigned int >::iterator planeIter;
308 for (planeIter = planeSet.begin(); planeIter != planeSet.end(); ++planeIter){
311 planeid = *planeIter;
314 for (
unsigned int icx = 0; icx < geom->
Plane(planeid)->
Ncells(); ++icx){
320 chanid = GetChanID(planeid,cellid,detid);
327 if ( ((fDiBlockNum>=0&&diblock==(
unsigned int)fDiBlockNum)||fDiBlockNum<0) &&
328 ((fDCMNum>=0&&dcmnum==(
unsigned int)fDCMNum)||fDCMNum<0) ){
331 if (mMeanPeXY.find(chanid) == mMeanPeXY.end()) {
333 mMeanPeXY.insert(std::pair<
daqchannelmap::dchan,std::vector<double> >(chanid,std::vector<double>()));
334 mMeanPeZ.insert(std::pair<
daqchannelmap::dchan,std::vector<double> >(chanid,std::vector<double>()));
335 mMeanPeAvg.insert(std::pair<
daqchannelmap::dchan,std::vector<double> >(chanid,std::vector<double>()));
348 std::map< daqchannelmap::dchan, std::vector<double> >::iterator chanIter;
349 for (chanIter = mMeanPeXY.begin(); chanIter != mMeanPeXY.end(); ++chanIter) {
350 channum = chanIter->first;
351 mMeanPeXY[channum].clear();
352 mMeanPeZ[channum].clear();
353 mMeanPeAvg[channum].clear();
375 e.
getByLabel(fPCHitLabel, fQualXYName, cellsxy);
378 e.
getByLabel(fPCHitLabel, fQualZName, cellsz);
381 e.
getByLabel(fPCHitLabel, fQualAvgName, cellsavg);
388 eventTime.push_back(evtTimeS);
389 if (evtTimeS>maxTime) maxTime = evtTimeS;
390 if (evtTimeS<minTime) minTime = evtTimeS;
393 FillPeCorrHits(cellsxy,mMeanPeXY,detid);
394 FillPeCorrHits(cellsz,mMeanPeZ,detid);
395 FillPeCorrHits(cellsavg,mMeanPeAvg,detid);
413 int nplanes = planeSetX.size() + planeSetY.size();
419 hMeanDet = tfs->
make<TH2F>(
"hMeanDet",
420 "Response Plane vs. FEB; Plane Number; FEB",
421 nplanes+1,0,nplanes+1,3,0,3.);
423 hRMS = tfs->
make<TH1F>(
"hRMS",
424 "RMS of Truncated Response; RMS of Truncated Response (PE); Modules",
427 hMeanFrac = tfs->
make<TH1F>(
"hMeanFrac",
428 "Fractional Error on Mean; Frac Error on Mean; Modules",
431 hNHits = tfs->
make<TH1F>(
"hNHits",
432 "Number of Hits in FEB; Number of Hits in FEB; Modules",
438 std::unique_ptr<std::vector<caldp::DriftResponse > >fDriftResponse(
new std::vector<caldp::DriftResponse >);
458 std::ofstream myfile;
460 myfile.open (
"chaninfo.txt");
466 std::map< daqchannelmap::dchan, std::vector<double> >::iterator chanIter;
467 for (chanIter = mMeanPeXY.begin(); chanIter != mMeanPeXY.end(); ++chanIter) {
470 fChanId=chanIter->first;
471 std::vector<double> chanvector;
474 if (fFillType==
"Successive"||fFillType==
"Exclusive"){
477 chanvector=chanIter->second;
481 CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);
484 if (fMean==0||(fMean>0&&fMeanErr/fMean>0.03)){
485 if (fFillType==
"Exclusive") chanvector.clear();
486 chanvector.insert(chanvector.end(),
487 (mMeanPeZ.at(fChanId)).begin(),
488 (mMeanPeZ.at(fChanId)).
end());
492 CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);
495 if (fMean==0||(fMean>0&&fMeanErr/fMean>0.03)){
496 if (fFillType==
"Exclusive") chanvector.clear();
497 chanvector.insert(chanvector.end(),
498 (mMeanPeAvg.at(fChanId)).begin(),
499 (mMeanPeAvg.at(fChanId)).
end());
503 CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);
510 }
else if (fFillType==
"XY"){
513 chanvector=chanIter->second;
517 CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);
519 }
else if (fFillType==
"Z"){
522 chanvector.insert(chanvector.end(),
523 (mMeanPeZ.at(fChanId)).begin(),
524 (mMeanPeZ.at(fChanId)).
end());
528 CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);
530 }
else if (fFillType==
"Avg"){
533 chanvector.insert(chanvector.end(),
534 (mMeanPeAvg.at(fChanId)).begin(),
535 (mMeanPeAvg.at(fChanId)).
end());
539 CalculateMean(chanvector,fMean,fMeanErr,fNHits,fRMS);
552 myfile <<
"//--------Statistics for ChanId=" << fChanId <<
" (initial size=" << chanvector.size() <<
")\n";
553 myfile <<
"// Mean=" << fMean <<
"+/-" << fMeanErr <<
" RMS=" << fRMS <<
" n=" << fNHits <<
" qual=" << fPathQual <<
"\n";
555 std::string histname = Form(
"h_r%d_id%d",fRun,fChanId);
557 srHistprint[fChanId] = tfs->
make<TH1F>(Form(
"%s",histname.c_str()),
558 "Truncated PECorr;PECorr;Number of Hits",
560 DrawMean(chanvector,srHistprint[fChanId]);
565 int planefill = (
int)(fChanId/1e6);
566 if (planefill==prevplane){
568 }
else { FEBcount=0; }
571 hMeanDet->SetBinContent(planefill+1,FEBcount+1,fMean);
572 hMeanFrac->Fill(fMeanErr/fMean);
574 hNHits->Fill(fNHits);
583 fDriftResponse->push_back(driftresp);
586 r.
put(std::move(fDriftResponse));
587 if (fDrawMeans) myfile.close();
std::map< daqchannelmap::dchan, std::vector< double > > mMeanPeZ
Structure to hold drift response in a single channel at a single time.
diblock
print "ROW IS " print row
void AddChannelResponse(int fOffChan, double fMean, double fMeanErr, double fRMS, int fNHits)
Return subrun end time.
std::string fFillType
Instance name, Avg qual hits.
int fDiBlockNum
Save response histograms to file, true/false?
constexpr std::uint32_t timeHigh() const
Vertical planes which measure X.
unsigned int Ncells() const
Number of cells in this plane.
int Cell() const
Return cell value.
std::string fQualXYName
Label of PCHit data.
const daqchannelmap::DAQChannelMap * Map() const
int fDCMNum
DiBlock Number to use (set to 0 with fDCMNum=0 to run all)
const PlaneGeo * Plane(unsigned int i) const
art::ProductID put(std::unique_ptr< PROD > &&)
DEFINE_ART_MODULE(TestTMapFile)
lchan encodeLChan(int detId, plane_t plane, cell_t cell) const
Horizontal planes which measure Y.
"Pre-calibration hit". Common input to calibration procedures
float Path() const
Return Path value.
void SetPlane(unsigned short plane)
float W() const
Return W value.
void SetCell(unsigned short cell)
std::map< daqchannelmap::dchan, std::vector< double > > mMeanPeAvg
float PE() const
Return PE value.
double GetAttenScale(rb::CellHit const &cellhit, double w)
for PE->PECorr conversion
std::vector< double > eventTime
T get(std::string const &key) const
std::map< daqchannelmap::dchan, std::vector< double > > mMeanPeXY
DCM Number to use (set to 0 with fDiBlockNum=0 to run all)
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
uint32_t feb_t
Type for DCM link port. Counts from 0.
EPathQuality
Methods used to estimate path length of a track through a cell.
std::string fQualAvgName
Instance name, Z qual hits.
int Plane() const
Return plane value.
bool fDrawMeans
Method for calculating means (which quals to use)
std::map< daqchannelmap::dchan, TH1F * > srHistprint
uint32_t dcm_id_t
Type for DCM number, counts from 1.
void SetMC(bool isMC=true)
static float min(const float a, const float b, const float c)
Histograms used by attenuation calibration.
bool CheckForChannel(int fOffChan) const
Return number of hits for given offline channel.
A rawdata::RawDigit with channel information decoded.
T * make(ARGS...args) const
uint32_t plane_t
Type for plane.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
std::string fQualZName
Instance name, XY qual hits.
const std::set< unsigned int > & GetPlanesByView(View_t v=kXorY) const
uint32_t diblock_t
Type for diblocks and blocks.
dcm_id_t getDCM(dchan daqchan) const
Decode the dcm ID from a dchan.
dchan encodeDChan(int detID, diblock_t diblock, dcm_id_t dcm, feb_t feb, pixel_t pixel) const
uint32_t dchan
< DAQ Channel Map Package
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
diblock_t getDiBlock(dchan daqchan) const
Decode the diblock ID from a dchan.
Encapsulate the geometry of one entire detector (near, far, ndos)
feb_t getFEB(dchan daqchan) const
Decode the feb id from a dchan.