Commit 8eeb0dd5 authored by dino's avatar dino
Browse files

The traces of the signal basis are reordered in memory to keep the samples...

The traces of the signal basis are reordered in memory to keep the samples with the natural spacing of 10ns close to each other. 
The 0,10,20... ns values are in the first part of the wave
The 5,15,25... ns values are in the second part.
This increases the efficiency of the memory pre-fetch and caching mechanisms because the grid-search always advances in steps of 10ns.

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@848 170316e4-aea8-4b27-aad4-0380ec0519c9
parent dabd4e9b
......@@ -56,7 +56,8 @@ char hmask[NCHAN][NCHAN+1] = {
"0000220000000000000000000000220000219" // F6
};
const float eScale = 5.f;
const float eScale = 5.f;
const float dScale = 10.f;
const int specLenE = 16*1024;
const int specLenT = 2*1000;
const int matLen = 100;
......@@ -65,7 +66,7 @@ const int matOff = 50;
#ifdef USENETCHARGES
const int diffLag = 0; // use the net-charge signals as they are
#else
const int diffLag = 0; // use them with a delayed-differentiation (in units of signal samples) if !=0
const int diffLag = 5; // use them with a delayed-differentiation (in units of signal samples) if !=0
//# define SHOWDIFFERENTIATED // in the saved traces show the differentiated version of the net-charge segments
#endif
......@@ -139,7 +140,7 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
int rval = 0;
if(fPsaSegCenter) {
rval = fBasis.ReadCenters(fBasisFile);
rval = fBasis.ReadSegmentCenters(fBasisFile);
if(rval != 0)
return rval;
// not much more to do
......@@ -154,6 +155,15 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
return rval;
#endif
fBasis.SmoothBasisSignals(SMOOTH);
if(diffLag > 0) {
// Modify the signal basis because the grid-search uses the differentiated net-charge signals
fBasis.DifferentiateNetChargeSignals(diffLag);
}
fBasis.ReorderSamples();
// prepare the metrics vectors for the chi2 loops
float largest = 0;
for(int i=0; i<NMETRIC; i++) {
......@@ -164,10 +174,7 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
}
}
if(diffLag > 0) {
// Modify the signal basis because the grid-search uses the differentiated net-charge signals
fBasis.DifferentiateNetChargeSignals(diffLag);
}
// repack basis to keep locality of samples
//specEner = new nDhist<unsigned int>(40, specLenE);
//specEner->setFileName(fOdirPrefix+"Psa?Ener.spec");
......@@ -319,9 +326,9 @@ Int_t PSAFilterGridSearch::Process(int slot)
#endif
#ifdef USEADAPTIVE
SearchAdaptive(&pS, netChSeg, localMask, samp_first, samp_last);
int usamp = SearchAdaptive(&pS, netChSeg, localMask, samp_first, samp_last);
#else
SearchFullGrid(&pS, netChSeg, localMask, samp_first, samp_last);
int usamp = SearchFullGrid(&pS, netChSeg, localMask, samp_first, samp_last);
#endif
float fx = fBasis.Pts[netChSeg][pS.bestPt].x;
......@@ -342,7 +349,8 @@ Int_t PSAFilterGridSearch::Process(int slot)
// traces
// The trace of the best point is normalized to the experimental amplitude and accumulated
if(fWriteNumTraces > 0) {
short netcharge[NTIME]; // to reconstruct the netcharge if grid-search is done using the differentiated net-charge signals.
short netcharge[NTIME]; // to reconstruct the netcharge if grid-search is done using the differentiated net-charge signals.
int time_first = TSHIFT; // for the 5-ns part this is NTIME/(SAMP_STEP/TIME_STEP)+TSHIFT
fact = eScale/fact;
for(int segm = 0; segm < NCHAN; segm++) {
float * realTrace = pS.rAmplitude[segm];
......@@ -353,17 +361,15 @@ Int_t PSAFilterGridSearch::Process(int slot)
// If, as here, the original trace is needed we reproduce it here by deconvolution of the differentiated one
if(diffLag > 0 && segm == netChSeg ) {
memcpy(netcharge, baseTrace, sizeof(netcharge));
for(int nn = 2*diffLag; nn < NTIME; nn++)
netcharge[nn] += netcharge[nn-2*diffLag];
for(int nn = diffLag; nn < NTIME/2; nn++)
netcharge[nn] += netcharge[nn-diffLag];
baseTrace = netcharge;
}
#endif
int iTime0 = TSHIFT + pS.bestdt;
int iTime1 = NTIME-1; // plot using the whole signal
int iSamp = samp_first;// + iTime0/(SAMP_STEP/TIME_STEP);
// this trace is produced always in steps of 10 ns
for(int iTime = iTime0; iTime <= iTime1; iTime += SAMP_STEP/TIME_STEP, iSamp++) {
realTrace[iSamp] += fact*baseTrace[iTime];
realTrace += samp_first;
baseTrace += time_first;
for(int nn = 0; nn < usamp; nn++) {
(*realTrace++) += (*baseTrace++)*fact;
}
}
}
......@@ -519,8 +525,32 @@ int PSAFilterGridSearch::WriteTraces(int slot)
}
}
if(chanToCopy < WCHAN) {
// next come the hits
memset(wbuf, 0, sizeof(wbuf));
for(int ns = chanToCopy; ns < WCHAN; ns++) {
PsaData *pD = DD + slot;
int numHits = pD->numHits;
for(int nh = 0; nh < numHits; nh++) {
int pos = nh*10;
if(pos >= (int)toWrite)
break;
ADF::PSAHit *pH = &pD->theHits[nh];
pos++;
wbuf[pos++] = pH->GetID();
wbuf[pos++] = (int)(pH->GetE()*eScale);
wbuf[pos++] = (int)(pH->GetX()*dScale);
wbuf[pos++] = (int)(pH->GetY()*dScale);
wbuf[pos++] = (int)(pH->GetZ()*dScale);
wbuf[pos++] = (int)(dScale*sqrt(pH->GetX()*pH->GetX()+pH->GetY()*pH->GetY()));
}
size_t nwritten = fwrite(wbuf, sizeof(short), toWrite, fpPsaTraces);
if(nwritten != toWrite) {
rstat = 1;
fWriteNumTraces = 1;
goto finish;
}
// now the empty ones
memset(wbuf, 0, sizeof(wbuf));
for(int ns = chanToCopy; ns < WCHAN-1; ns++) {
size_t nwritten = fwrite(wbuf, sizeof(short), toWrite, fpPsaTraces);
if(nwritten != toWrite) {
rstat = 1;
......@@ -541,106 +571,64 @@ finish:
return rstat;
}
void PSAFilterGridSearch::SearchFullGrid(pointExp *pS,
int netChSeg, char *lMask, int samp_first, int samp_last)
int PSAFilterGridSearch::SearchFullGrid(pointExp *pS, int netChSeg, char *lMask, int samp_first, int samp_last)
{
float *metrics = fMetrics + NMETRIC2;
// the limits on the calculated signals
int time_first = TSHIFT;
int time_last = (samp_last-samp_first+1)*(SAMP_STEP/TIME_STEP);
if(time_last >= NTIME-1)
time_last = NTIME-1;
// number of samples to use in the search
int uSamples = samp_last-samp_first+1;
if( uSamples > NTIME/(SAMP_STEP/TIME_STEP)-TSHIFT )
uSamples = NTIME/(SAMP_STEP/TIME_STEP)-TSHIFT;
int iTime0 = 0;
int iTime1 = 0;
int time_first = TSHIFT; // for the 5-ns part this is NTIME/(SAMP_STEP/TIME_STEP)+TSHIFT
int bestPt = 0;
int dt = 0;
int bestdt = 0;
int bestPt = 0;
float chi2min = pS->chi2min;
float chi2 = 0;
//int mSamp = 0;
const int sPts = 1; // step over the base points
const int sSamp = STEP_SAMP; // step over experimental traces (in 10 ns units)
const int sTime = sSamp*(SAMP_STEP/TIME_STEP); // equivalent step on the calculated basis
#ifdef USESUBSAMPLE
const int rt = 5; // range of subsample search is (-rt...rt)*TIME_STEP
// make sure there is enough space
if(time_first < rt) time_first = rt;
if(time_last+rt >= NTIME) time_last = NTIME-1 - rt;
#endif
for(int iPts = 0; iPts < fBasis.numPts[netChSeg]; iPts += sPts) {
#ifdef USESUBSAMPLE
for(dt= -rt; dt <= rt; dt++) { // subsample step search loop
#endif
iTime0 = time_first + dt; // keep it here for the dt loop (if enabled)
iTime1 = time_last + dt;
int nPts = fBasis.numPts[netChSeg];
for(int iPts = 0; iPts < nPts; iPts++) {
chi2 = 0;
for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') {
short * realTrace = pS->sAmplitude[iSegm];
short * baseTrace = fBasis.Pts[netChSeg][iPts].amplitude[iSegm];
int iSamp = samp_first;
int iTime = iTime0;
for( ; iTime <= iTime1; iTime += sTime, iSamp += sSamp) {
chi2 += metrics[realTrace[iSamp] - baseTrace[iTime]];
realTrace += samp_first;
baseTrace += time_first;
for(int nn = 0; nn < uSamples; nn++) {
chi2 += metrics[(*realTrace++) - (*baseTrace++)];
}
//if(iSamp > mSamp) mSamp = iSamp;
if(chi2 > chi2min)
break;
}
} // end loop over the segments
if(chi2 < chi2min) {
bestPt = iPts;
bestdt = dt;
chi2min = chi2;
}
#ifdef USESUBSAMPLE
} // end loop over the subsample step dt
#endif
} // end loop over the base points iPts
pS->bestdt = bestdt;
pS->bestdt = 0;
pS->bestPt = bestPt;
pS->chi2min = chi2min;
return uSamples;
}
void PSAFilterGridSearch::SearchAdaptive(pointExp *pS,
int netChSeg, char *lMask, int samp_first, int samp_last)
int PSAFilterGridSearch::SearchAdaptive(pointExp *pS, int netChSeg, char *lMask, int samp_first, int samp_last)
{
float *metrics = fMetrics + NMETRIC2;
// the limits on the calculated signals
int time_first = TSHIFT;
int time_last = (samp_last-samp_first+1)*(SAMP_STEP/TIME_STEP);
if(time_last >= NTIME-1)
time_last = NTIME-1;
// number of samples to use in the search
int uSamples = samp_last-samp_first+1;
if( uSamples > NTIME/(SAMP_STEP/TIME_STEP)-TSHIFT )
uSamples = NTIME/(SAMP_STEP/TIME_STEP)-TSHIFT;
int iTime0 = 0;
int iTime1 = 0;
int time_first = TSHIFT; // for the 5-ns part this is NTIME/(SAMP_STEP/TIME_STEP)+TSHIFT
int bestPt = 0;
int dt = 0;
int bestdt = 0;
int bestPt = 0;
float chi2min = pS->chi2min;
float chi2 = 0;
//int mSamp = 0;
const int sSamp = STEP_SAMP; // step over experimental traces (in 10 ns units)
const int sTime = sSamp*(SAMP_STEP/TIME_STEP); // equivalent step on the calculated basis
#ifdef USESUBSAMPLE
const int rt = 5; // range of subsample search is (-rt...rt)*TIME_STEP
// make sure there is enough space
if(time_first < rt) time_first = rt;
if(time_last+rt >= NTIME) time_last = NTIME-1 - rt;
#endif
float chi2 = 0;
// the coarse-fine structure of this segment
SignalBasis::cfgrid &sgrid = fBasis.cflist[netChSeg];
......@@ -653,42 +641,32 @@ void PSAFilterGridSearch::SearchAdaptive(pointExp *pS,
int knpts = sgrid.numCoarse;
for(int kloop = 0; kloop < knpts; kloop++) {
int kPts = sgrid.indFine[sgrid.posCoarse[kloop]];
#ifdef USESUBSAMPLE
for(dt= -rt; dt <= rt; dt++) { // subsample step search loop
#endif
iTime0 = time_first + dt; // keep it here for the dt loop (if enabled)
iTime1 = time_last + dt;
chi2 = 0;
for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') {
short * realTrace = pS->sAmplitude[iSegm];
short * baseTrace = fBasis.Pts[netChSeg][kPts].amplitude[iSegm];
int iSamp = samp_first; // here to inspect it after the loop
int iTime = iTime0; // here to inspect it after the loop
for( ; iTime <= iTime1; iTime += sTime, iSamp += sSamp) {
chi2 += metrics[realTrace[iSamp] - baseTrace[iTime]];
realTrace += samp_first;
baseTrace += time_first;
for(int nn = 0; nn < uSamples; nn++) {
chi2 += metrics[(*realTrace++) - (*baseTrace++)];
}
//if(iSamp > mSamp) mSamp = iSamp;
if(chi2 > chi2min)
break;
}
} // end loop over the segments
if(chi2 < chi2min) {
bestPt = kPts;
bestdt = dt;
chi2min = chi2;
kbest = kloop;
}
#ifdef USESUBSAMPLE
} // end loop over the subsample step dt
#endif
} // end loop over the base points iPts
pS->bestdt = bestdt;
pS->bestdt = 0;
pS->bestPt = bestPt;
pS->chi2min = chi2min;
//return;
//return uSamples;
/////////////////////////////
// loop on the fine grid //
......@@ -699,41 +677,32 @@ void PSAFilterGridSearch::SearchAdaptive(pointExp *pS,
int jloop2 = sgrid.posCoarse[kbest+1];
for(int jloop = jloop1; jloop < jloop2; jloop++) {
int jPts = sgrid.indFine[jloop];
#ifdef USESUBSAMPLE
for(dt= -rt; dt <= rt; dt++) { // subsample step search loop
#endif
iTime0 = time_first + dt; // keep it here for the dt loop (if enabled)
iTime1 = time_last + dt;
chi2 = 0;
for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') {
short * realTrace = pS->sAmplitude[iSegm];
short * baseTrace = fBasis.Pts[netChSeg][jPts].amplitude[iSegm];
int iSamp = samp_first; // here to inspect it after the loop
int iTime = iTime0; // here to inspect it after the loop
for( ; iTime <= iTime1; iTime += sTime, iSamp += sSamp) {
chi2 += metrics[realTrace[iSamp] - baseTrace[iTime]];
realTrace += samp_first;
baseTrace += time_first;
for(int nn = 0; nn < uSamples; nn++) {
chi2 += metrics[(*realTrace++) - (*baseTrace++)];
}
//if(iSamp > mSamp) mSamp = iSamp;
if(chi2 > chi2min)
break;
}
} // end loop over the segments
if(chi2 < chi2min) {
bestPt = jPts;
bestdt = dt;
chi2min = chi2;
jbest = jloop;
}
#ifdef USESUBSAMPLE
} // end loop over the subsample step dt
#endif
} // end loop over the base points iPts
pS->bestdt = bestdt;
pS->bestdt = 0;
pS->bestPt = bestPt;
pS->chi2min = chi2min;
return uSamples;
}
Int_t PSAFilterGridSearch::PostProcess(int slot)
......
......@@ -68,8 +68,8 @@ private:
void PrepareEvent (PsaData *pD, pointExp *pS);
void SetToSegCenter(PsaData *pD, pointExp *pS);
void SearchFullGrid(pointExp *pS, int netChSeg, char *lMask, int addr_first, int addr_last);
void SearchAdaptive(pointExp *pS, int netChSeg, char *lMask, int addr_first, int addr_last);
int SearchFullGrid(pointExp *pS, int netChSeg, char *lMask, int addr_first, int addr_last);
int SearchAdaptive(pointExp *pS, int netChSeg, char *lMask, int addr_first, int addr_last);
protected:
// this is written in a thread-safe mode and can be called in parallel, using different data slots
Int_t Process(int slot = 0);
......@@ -89,8 +89,8 @@ public:
private:
// this is needed to have the possibility to write the traces in a thread-safe mode
short slot_fAmplitude[TCOUNT*TMODULO][NCHAN*NTIME]; // locally saved experimental trace
short slot_rAmplitude[TCOUNT*TMODULO][NCHAN*NTIME]; // locally saved "fitted" trace
short slot_fAmplitude[TCOUNT*TMODULO][(NCHAN+1)*NTIME]; // locally saved experimental trace
short slot_rAmplitude[TCOUNT*TMODULO][(NCHAN+1)*NTIME]; // locally saved "fitted" trace
};
......
......@@ -21,11 +21,11 @@ const int TIME_STEP = 5; // (ns) time step of the signal basis ==>
const int NTIME = 100; // number of points of the basis signals (100*5ns = 500ns)
const int TSHIFT = 5; // the "T0" (samples) of the calculated traces
const int SMOOTH = 10; // with of the moving average used to smooth the basis
const int SMOOTH = 10; // with of the moving average used to smooth the basis (== 50 ns)
const float XTCORR = 0.002f; // fraction of CC to add to the segments
const float DECONV = 4600.f; // to deconvolve the main pole
const int NSAMP = defTraceLengthPSA; // length of experimental data (taken as float data)
const int TIME0 = defTriggerSample; // pretrigger of the experimental traces
const int STEP_SAMP = 2; // step over experimental traces (in 10 ns units)
const int STEP_SAMP = 1; // step over experimental traces (in 10 ns units)
#endif // GRIDSEARCHPARAMS_H_INCLUDED
......@@ -59,18 +59,21 @@ int SignalBasis::ReadBasis(std::string fname, bool keep)
if(bVentu && !bBartB) {
printf("Reading the basis of signals ...\n");
return ReadBasisFormatVentu(ofp1, keep);
int readErr = ReadBasisFormatVentu(ofp1, keep);
return readErr;
}
else if (!bVentu && bBartB) {
size_t numSignals = fSize/sizeof(bbhit_t);
printf("Reading Bart\'s basis containing %d signals ...\n", (int)numSignals);
return ReadBasisFormatBartB(ofp1, keep, (int)numSignals);
int numSignals = int(fSize/sizeof(bbhit_t));
printf("Reading Bart\'s basis containing %d signals ...\n", numSignals);
int readErr = ReadBasisFormatBartB(ofp1, keep, numSignals);
return readErr;
}
else {
printf("Cannot decide format of signal basis\n");
fclose(ofp1);
return 130;
}
return 0;
}
int SignalBasis::ReadBasisFormatVentu(FILE * ofp1, bool keep)
......@@ -87,12 +90,10 @@ int SignalBasis::ReadBasisFormatVentu(FILE * ofp1, bool keep)
totlen += npoints;
}
float Xav=0, Yav=0, Zav=0;
int nerr = 0;
for(int ii=0; ii<NSEGS; ii++) {
npoints = numPts[ii];
Pts[ii] = new pointPsa[npoints];
float xav = 0, yav = 0, zav = 0;
for(int jj=0; jj<npoints; jj++) {
size_t rv;
rv = fread(&(Pts[ii][jj].fullChargeCollected), sizeof(bool), 1, ofp1); ///
......@@ -104,47 +105,20 @@ int SignalBasis::ReadBasisFormatVentu(FILE * ofp1, bool keep)
rv = fread(&(Pts[ii][jj].y), sizeof(float), 1, ofp1); ///
rv = fread(&(Pts[ii][jj].z), sizeof(float), 1, ofp1); ///
rv = fread( Pts[ii][jj].amplitude, sizeof(short), NCHAN*NTIME, ofp1); ///
Pts[ii][jj].masksignals[NCHAN] = 0; ///
Pts[ii][jj].smooth(SMOOTH); ///
Pts[ii][jj].masksignals[NCHAN] = 0; // to close the string
int ntseg = Pts[ii][jj].netChargeSegment;
if(ntseg != ii)
nerr++;
xav += Pts[ii][jj].x;
yav += Pts[ii][jj].y;
zav += Pts[ii][jj].z;
averPt[ntseg].x += Pts[ii][jj].x;
averPt[ntseg].y += Pts[ii][jj].y;
averPt[ntseg].z += Pts[ii][jj].z;
}
Xav += xav; Yav += yav; Zav += zav;
xav /= npoints; yav /= npoints; zav /= npoints;
printf("Segment %3d contains %6d points", ii, npoints);
printf(" and is centered at %7.2f %7.2f %7.2f\n", xav, yav, zav);
if(!keep) {
delete [] Pts[ii];
Pts[ii] = NULL;
}
}
Xav /= totlen; Yav /= totlen; Zav /= totlen;
printf("This basis contains %6d points", totlen);
printf(" and is centered at %7.2f %7.2f %7.2f\n", Xav, Yav, Zav);
// now calculate the real center of segments
// could differ from the previous one if Pts[ii][jj].netChargeSegment != ii
for(int ii=0; ii<NSEGS; ii++) {
averPt[NSEGS].x += averPt[ii].x;
averPt[NSEGS].y += averPt[ii].y;
averPt[NSEGS].z += averPt[ii].z;
averPt[ii].x /= numPts[ii];
averPt[ii].y /= numPts[ii];
averPt[ii].z /= numPts[ii];
}
averPt[NSEGS].x /= totlen;
averPt[NSEGS].y /= totlen;
averPt[NSEGS].z /= totlen;
fclose(ofp1);
CalculateSegmentCenters();
bCenters = true;
bFullGrid = keep;
......@@ -212,13 +186,13 @@ char hmask[NCHAN][NCHAN+1] = {
// Read the basis to determine the number of points in the 36 segments
// whichSeg records the net-charge segment of each of the points
while(true) {
if( fread(&bbhit, sizeof(float), 5, ofp1) != 5)
if( fread(&bbhit, sizeof(float), 5, ofp1) != 5)
break;
if( fread(&bbhit.T0, sizeof(int), 1, ofp1) != 1)
if( fread(&bbhit.Trig, sizeof(int), 1, ofp1) != 1)
break;
if( fread(bbhit.Es, sizeof(float), NCHAN, ofp1) != NCHAN)
if( fread(bbhit.Es, sizeof(float), NCHAN, ofp1) != NCHAN)
break;
if( fread(bbhit.Tr, sizeof(float), NCHAN*120, ofp1) != NCHAN*120)
if( fread(bbhit.Tr, sizeof(float), NCHAN*120, ofp1) != NCHAN*120)
break;
// look for perfect signals
......@@ -264,7 +238,7 @@ char hmask[NCHAN][NCHAN+1] = {
}
if(np != 1 || npp) {
printf("warning: minor problem at signal %5d %3d %3d [%7.3f %7.3f %7.3f] %7.5f %7.5f %7.5f %7.5f\n",
printf("warning: problem at signal %5d %3d %3d [%7.3f %7.3f %7.3f] %7.5f %7.5f %7.5f %7.5f\n",
totsigs, nmax, npp, bbhit.Pos[0], bbhit.Pos[1], bbhit.Pos[2], bbhit.Eint, bbhit.Es[0], ep1, ep2);
}
......@@ -291,13 +265,13 @@ char hmask[NCHAN][NCHAN+1] = {
fseek(ofp1, 0, SEEK_SET);
for(int nns = 0; ; nns++) {
if( fread(&bbhit, sizeof(float), 5, ofp1) != 5)
if( fread(&bbhit, sizeof(float), 5, ofp1) != 5)
break;
if( fread(&bbhit.T0, sizeof(int), 1, ofp1) != 1)
if( fread(&bbhit.Trig, sizeof(int), 1, ofp1) != 1)
break;
if( fread(bbhit.Es, sizeof(float), NCHAN, ofp1) != NCHAN)
if( fread(bbhit.Es, sizeof(float), NCHAN, ofp1) != NCHAN)
break;
if( fread(bbhit.Tr, sizeof(float), NCHAN*120, ofp1) != NCHAN*120)
if( fread(bbhit.Tr, sizeof(float), NCHAN*120, ofp1) != NCHAN*120)
break;
int ii = whichSeg[nns];
......@@ -321,34 +295,13 @@ char hmask[NCHAN][NCHAN+1] = {
memcpy(Pts[ii][jj].masksignals, hmask[ii], NCHAN+1);
Pts[ii][jj].addr_first = TSHIFT;
Pts[ii][jj].addr_last = NTIME-1;
Pts[ii][jj].smooth(SMOOTH);
averPt[ii].x += Pts[ii][jj].x;
averPt[ii].y += Pts[ii][jj].y;
averPt[ii].z += Pts[ii][jj].z;
//Pts[ii][jj].smooth(SMOOTH); // this should not be done here
numPts[ii]++;
}
// calculate centers
int totlen = 0;
for(int ns = 0; ns < NSEGS; ns++) {
averPt[NSEGS].x += averPt[ns].x;
averPt[NSEGS].y += averPt[ns].y;
averPt[NSEGS].z += averPt[ns].z;
averPt[ns].x /= numPts[ns];
averPt[ns].y /= numPts[ns];
averPt[ns].z /= numPts[ns];
printf("Segment %3d contains %6d points", ns, numPts[ns]);
printf(" and is centered at %7.2f %7.2f %7.2f\n", averPt[ns].x, averPt[ns].y, averPt[ns].z);
totlen += numPts[ns];
}
averPt[NSEGS].x /= totlen;
averPt[NSEGS].y /= totlen;
averPt[NSEGS].z /= totlen;
printf("This basis contains %6d points", totlen);
printf(" and is centered at %7.2f %7.2f %7.2f\n", averPt[NSEGS].x, averPt[NSEGS].y, averPt[NSEGS].z);
CalculateSegmentCenters();
bCenters = true;
bFullGrid = true;
......@@ -357,7 +310,7 @@ char hmask[NCHAN][NCHAN+1] = {
return 0;
}
int SignalBasis::ReadCenters(std::string fname)
int SignalBasis::ReadSegmentCenters(std::string fname)
{
bCenters = false;
bFullGrid = false;
......@@ -877,3 +830,70 @@ int SignalBasis::DifferentiateNetChargeSignals(int lag)
return count;
}
// passes a moving average window on all signals
void SignalBasis::SmoothBasisSignals(int nsmoo)
{
if(nsmoo >=2) {
int count = 0;
for(int ii=0; ii<NSEGS; ii++) {
int npoints = numPts[ii];
for(int jj=0; jj<npoints; jj++) {
Pts[ii][jj].smooth(nsmoo);