Commit f3effdfc authored by dino's avatar dino
Browse files

The inner representation of the calculated basis is now only float (removed int16 and int32)

Introduced a "Tukey's biweight" type metrics and removed the event-by-event rescaling to MAXNORM of the experimental data before calculating the chi2 loop; now the the calculated trace is rescaled at every access.
Although not yet used, the exponential decay of the transfer function can now be give for each segment (no more only slice-wise)


git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@1045 170316e4-aea8-4b27-aad4-0380ec0519c9
parent dc026a9d
......@@ -36,7 +36,7 @@ PSAFilter::PSAFilter() :
fBasisFile.clear();
fPsaXtalkFile.clear();
fPsaTshiftFile.clear();
//fPsaTshiftFile.clear();
fTraceLengthPSA = defTraceLengthPSA;
fWriteNumTraces = 0;
fPsaMinSegMult = 1;
......@@ -251,10 +251,10 @@ void PSAFilter::GetParameters(UInt_t *error_code)
ok = data.size() > 0;
fPsaXtalkFile = data;
}
else if( stringEq(keyw, "TshiftFile") ) {
ok = data.size() > 0;
fPsaTshiftFile = data;
}
//else if( stringEq(keyw, "TshiftFile") ) {
// ok = data.size() > 0;
// fPsaTshiftFile = data;
//}
else if( stringEq(keyw, "Verbose") ) {
fVerbose = true;
ok = true;
......@@ -328,6 +328,7 @@ void PSAFilter::GetParameters(UInt_t *error_code)
if(fDeadSegment >= CrystalInterface::kNbSegments)
fDeadSegment = -1;
}
// this should be taken from the general calibration file of the detector, PreprocessingFilterPSA.conf
else if( stringEq(keyw, "TauSlice") ) {
ok = 7 == sscanf(data.c_str(), "%f %f %f %f %f %f %f", fTauSlice, fTauSlice+1, fTauSlice+2, fTauSlice+3, fTauSlice+4, fTauSlice+5, fTauSlice+6);
fTauGiven = ok;
......@@ -379,6 +380,7 @@ Int_t PSAFilter::SetInput(int slot)
return 0;
}
// Part of the histogramming could be moved here (e.g. VIGRU)
Int_t PSAFilter::SetOutput(int slot)
{
fFramePSA->Reset();
......
......@@ -215,8 +215,8 @@ protected:
std::string fOdirPrefix;
std::string fBasisFile;
std::string fPsaXtalkFile; // File of cross talk correction coefficients from preprocessing
std::string fPsaTshiftFile; // PreprocessingFilterPSA.conf
std::string fPsaXtalkFile; // File of cross talk correction coefficients from preprocessing (xdir_1325-1340.cal)
//std::string fPsaTshiftFile; // this functionality is now in PreprocessingFilterPSA, which has to provide fully aligned data
std::string fnPsaHits; //
Int_t fTraceLengthPSA;
Int_t fWriteNumTraces;
......
This diff is collapsed.
......@@ -9,8 +9,8 @@
//#include "commonDefs.h"
// Define this symbol to write to file a correction to the timestamp of the events
//#define TSHIFT_CORRECTION
#ifdef TSHIFT_CORRECTION
//#define WRITE_TSHIFT_CORRECTION
#ifdef WRITE_TSHIFT_CORRECTION
#include "dataPairFile.h"
#endif
......@@ -56,7 +56,7 @@ const int WSAMP = 60; // number of samples per channel in the output waves
struct pointFull : public pointExp
{
pointFull() : samp_first(0), samp_last(0), samp_used(0),
isValid(false), isSelected(false), isInitialized(false), isPreSearched(false) {}
isValid(false), isSelected(false), isInitialized(false), isPreSearched(false), baseScale(1.f) {}
int samp_first; // must be a multiple of 4 to use sse intrinsics
int samp_last;
int samp_used;
......@@ -65,11 +65,14 @@ struct pointFull : public pointExp
bool isSelected; // used to select events
bool isInitialized;
bool isPreSearched;
float enerSeg; // the energy of the segment for GridSearch2
float baseScale; // temporarily used to manage the chi2 mapping
PsaOut_t PsaOut[ADF::CrystalInterface::kNbSegments];
void Reset() {
pointExp::Reset();
samp_first = samp_last = samp_used = 0;
isValid = isSelected= isInitialized = isPreSearched = false;
baseScale = 1;
}
};
......@@ -90,7 +93,7 @@ private:
bool gPsaSegCenter;
bool gFullSearch;
Float_t tauSlice[7];
float tauSGCC[37]; // the 36 segments and the CC
#ifdef PSA_LOCALSPECTRA
nDhist<unsigned short> *baseZXY;
......@@ -104,7 +107,7 @@ private:
nDhist<unsigned short> *matrRZE;
#endif //PSA_LOCALSPECTRA
#ifdef TSHIFT_CORRECTION
#ifdef WRITE_TSHIFT_CORRECTION
dataPairFile<ULong64_t, float> *tstmpFileO;
#endif
......@@ -162,14 +165,14 @@ public:
void SetHitSegThreshold(Float_t thres) {fHitSegThreshold = thres;}
Float_t GetHitSegThreshold() {return fHitSegThreshold;}
bool WriteWorkingWave(float *pW, bool append);
bool WritePartialWave(gs_type *pW, float fscale, bool append);
bool WriteWorkingWave(float *pW, bool append);
bool WritePartialWave(float *pW, float fscale, bool append);
int WriteTraces (PsaData *pD);
int WritePsaHits (PsaOut_t *pOut, int index = 0);
void AddToSolution (pointFull &pS, pointPsa &bestPoint, int numsamp, float scaleFact);
void AddToSolution (pointFull &pS, pointPsa &bestPoint1, pointPsa &bestPoint2, float bestFac, int numsamp, float scaleFact);
void SaveTotalTrace(pointFull &pS, PsaData *pD);
Float_t FitT0AfterPSA(pointFull &pS, int tsamp);
Float_t FitT0AfterPSA(pointFull &pS);
private:
void Extend(float *dd, int n1, int n2) {
......@@ -181,37 +184,22 @@ private:
// intergrid final interpolation (still just a uniform random distribution)
// be aware that this depends on the basis and its grid step (assumed 2 mm)
float xInterp() {return (rand()%200 - 100)*0.01f; }
float yInterp() {return (rand()%200 - 100)*0.01f; }
float zInterp() {return (rand()%200 - 125)*0.01f; }
float xInterp() {return 0; }
float yInterp() {return 0; }
float zInterp() {return 0; }
//float xInterp() {return (rand()%200 - 100)*0.01f; }
//float yInterp() {return (rand()%200 - 100)*0.01f; }
//float zInterp() {return (rand()%200 - 125)*0.01f; }
//float xInterp() {return (gCoarseOnly) ? 0 : (rand()%200 - 100)*0.01f; }
//float yInterp() {return (gCoarseOnly) ? 0 : (rand()%200 - 100)*0.01f; }
//float zInterp() {return (gCoarseOnly) ? 0 : (rand()%200 - 125)*0.01f; }
void incrPosMats(int netChSeg, int nhi, float fx, float fy, float fz, float rr);
private:
float Chi2InnerLoop(const float *pReal, const float *pBase);
float Chi2InnerLoop(const int *pReal, const int *pBase);
float Chi2InnerLoop(const short *pReal, const short *pBase);
float Chi2InnerLoop(const float *pReal, const float *pBase1, const float *pBase2);
float Chi2OuterLoop(const float **pAmpliReal, const float **pAmpliBase, const int nss, const float chi2min);
float Chi2OuterLoop(const float **pAmpliReal, const float **pAmpliBase1, const float **pAmpliBase2, const int nss, const float chi2min);
int loopSamples;
void SetLoopSamples(int nn) {loopSamples = nn;}
float baseScale;
void SetBaseScale(float scale) {baseScale = scale;}
float diffScale;
void SetDiffScale(float scale) {diffScale = scale;}
#ifdef USE_SSE_VERSION
__m128 loopFac1, loopFac2;
void SetLoopFactors(float ff1, float ff2) {loopFac1 = _mm_set_ps(ff1, ff1, ff1, ff1); loopFac2 = _mm_set_ps(ff2, ff2, ff2, ff2);}
#else
float loopFac1, loopFac2;
void SetLoopFactors(float ff1, float ff2) {loopFac1 = ff1; loopFac2 = ff2;}
#endif
float Chi2InnerLoop(const float *pReal, const float *pBase, float bScale);
float Chi2InnerLoop(const float *pReal, const float *pBase1, float bScale1, const float *pBase2, float bScale2);
// For the Chi2OuterLoop methods, look to versions before r1039
};
#if defined(__GNUC__)
......@@ -228,12 +216,13 @@ private:
///////////// 1 interaction //////////
//////////////////////////////////////
#define COMPUTE_VALUE0(nn) ( diff = _mm_sub_ps(realTrace[nn], baseTrace[nn]), chis = _mm_add_ps(chis, _mm_mul_ps(diff, diff)) ) // pow( d ,2 )
#define COMPUTE_VALUE1(nn) ( chis = _mm_add_ps(chis, _mm_andnot_ps(masks, _mm_sub_ps(realTrace[nn], baseTrace[nn])) ) ) // pow(|d|,1 )
#define COMPUTE_VALUE2(nn) ( chis = _mm_add_ps(chis, _mm_sqrt_ps(_mm_andnot_ps(masks, _mm_sub_ps(realTrace[nn], baseTrace[nn])) ) ) ) // pow(|d|,1/2)
#define COMPUTE_VALUE3(nn) ( chis = _mm_add_ps(chis, _mm_sqrt_ps(_mm_sqrt_ps(_mm_andnot_ps(masks, _mm_sub_ps(realTrace[nn], baseTrace[nn])) )) ) ) // pow(|d|,1/4)
#define C1_WAVE(nn) ( wave = _mm_mul_ps(loopFac, baseTrace[nn]) )
#define COMPUTE_VALUE0(nn) ( diff = _mm_sub_ps(realTrace[nn], C1_WAVE(nn)), chis = _mm_add_ps(chis, _mm_mul_ps(diff, diff)) ) // pow( d ,2 )
#define COMPUTE_VALUE1(nn) ( chis = _mm_add_ps(chis, _mm_andnot_ps(masks, _mm_sub_ps(realTrace[nn], C1_WAVE(nn))) ) ) // pow(|d|,1 )
#define COMPUTE_VALUE2(nn) ( chis = _mm_add_ps(chis, _mm_sqrt_ps(_mm_andnot_ps(masks, _mm_sub_ps(realTrace[nn], C1_WAVE(nn))) ) ) ) // pow(|d|,1/2)
#define COMPUTE_VALUE3(nn) ( chis = _mm_add_ps(chis, _mm_sqrt_ps(_mm_sqrt_ps(_mm_andnot_ps(masks, _mm_sub_ps(realTrace[nn], C1_WAVE(nn))) )) ) ) // pow(|d|,1/4)
#define COMPUTE_VALUE(nn) COMPUTE_VALUE3(nn)
#define COMPUTE_VALUE(nn) COMPUTE_VALUE2(nn)
// As realTrace is used repeatedly, it will be in the cache after the first point. Therefore it is not worth prefetching it (actually it seems to be a bit slower)
// The basis is used only once ==> we should prefetch it with _MM_HINT_NTA (non temporal with minimal cache pollution). But _MM_HINT_T0 seems to be a bit faster.
......@@ -246,7 +235,7 @@ private:
# define PREFETCH_BASE(nn) ( _mm_prefetch ((const char *)&baseTrace[ nn], _MM_HINT_NTA) )
#endif
// To unroll the innermost loop using the preprocessor, SFIXEDLEN is assumed to be at least 32 or 36 or 40 or 44 or 48 (--> SS_NUM = [8 9 10 11 12] )
// To unroll the innermost loop using the preprocessor, LOOP_SAMP is assumed to be at least 32 or 36 or 40 or 44 or 48 (--> SS_NUM = [8 9 10 11 12] )
#if SSE_NUM > 8
#define PREFETCH_REAL_08_ PREFETCH_REAL( 8)
......@@ -285,16 +274,21 @@ private:
#define COMPUTE_VALUE_11_
#endif
// the number of used samples is fixed to SFIXEDLEN
INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const float *pBase)
// the number of used samples is fixed to LOOP_SAMP
INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const float *pBase, float bScale)
{
const __m128 masks = _mm_set1_ps(-0.0f);
const __m128 zeros = _mm_setzero_ps();
__m128* realTrace = (__m128*)pReal;
__m128* baseTrace = (__m128*)pBase;
__m128 diff = _mm_setzero_ps();
__m128 chis = _mm_setzero_ps();
bScale *= RESCALE;
__m128 loopFac = _mm_set_ps(bScale, bScale, bScale, bScale);
__m128* realTrace = (__m128*)pReal;
__m128* baseTrace = (__m128*)pBase;
__m128 wave = _mm_setzero_ps();
__m128 diff = _mm_setzero_ps();
__m128 chis = _mm_setzero_ps();
//// intermix the prefetch and compute lines to get the desired prefetch lag
//PREFETCH_REAL( 0); PREFETCH_REAL( 1); PREFETCH_REAL( 2); PREFETCH_REAL( 3);
......@@ -319,7 +313,6 @@ INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const
return chi2;
}
///////////////////////////////////////
///////////// 2 interactions //////////
///////////////////////////////////////
......@@ -377,16 +370,23 @@ INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const
#define COMPUTE2_VALUE_11_
#endif
// the number of used samples is fixed to SFIXEDLEN
// the number of used samples is fixed to LOOP_SAMP
// the weight of the two components must be set by a preceeding call to SetLoopFactors(ff1, ff2)
INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const float *pBase1, const float *pBase2)
INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const float *pBase1, float bScale1, const float *pBase2, float bScale2)
{
const __m128 masks = _mm_set1_ps(-0.0f);
const __m128 zeros = _mm_setzero_ps();
bScale1 *= RESCALE;
__m128 loopFac1 = _mm_set_ps(bScale1, bScale1, bScale1, bScale1);
bScale2 *= RESCALE;
__m128 loopFac2 = _mm_set_ps(bScale2, bScale2, bScale2, bScale2);
__m128* realTrace = (__m128*)pReal;
__m128* baseTrace1 = (__m128*)pBase1;
__m128* baseTrace2 = (__m128*)pBase2;
__m128 wave = _mm_setzero_ps();
__m128 diff = _mm_setzero_ps();
__m128 chis = _mm_setzero_ps();
......@@ -412,144 +412,60 @@ INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const
return chi2;
}
INLINE_ALWAYS float PSAFilterGridSearch::Chi2OuterLoop(const float **pAmpliReal, const float **pAmpliBase, const int nss, const float chi2min)
{
const __m128 masks = _mm_set1_ps(-0.0f);
const __m128 zeros = _mm_setzero_ps();
__m128 diff = _mm_setzero_ps();
__m128 chis = _mm_setzero_ps();
float chi2 = 0;
for(int iss = 0; iss < nss; iss++) {
__m128* realTrace = (__m128*)pAmpliReal[iss];
__m128* baseTrace = (__m128*)pAmpliBase[iss];
COMPUTE_VALUE( 0); COMPUTE_VALUE( 1); COMPUTE_VALUE( 2); COMPUTE_VALUE( 3);
COMPUTE_VALUE( 4); COMPUTE_VALUE( 5); COMPUTE_VALUE( 6); COMPUTE_VALUE( 7);
COMPUTE_VALUE_08_; COMPUTE_VALUE_09_; COMPUTE_VALUE_10_; COMPUTE_VALUE_11_;
chi2 = _mm_cvtss_f32(_mm_hadd_ps(_mm_hadd_ps(chis, zeros), zeros)); // sum the 4 values add to chi2
if(chi2 > chi2min)
break;
}
return chi2;
}
#else // #ifdef USE_SSE_VERSION
// the standard, non-vectorized versions
// loopSamples is defined by a preceeding call to SetLoopSamples(uSamples)
// the two multiplications baseScale and diffScale cost about 5% in terms of speed
INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const float *pBase)
// The multiplication * baseScale costs ~ 5% in terms of speed. It could be avoided by pre-scaling the experimental trace
// in such a way that bScale == 1. In this case the mapped-metrics vector is applied with a varying scale and loses its
// statistical meaning (if it has one like with the biweight). Of course we could rescale in accessing it but it wouled
INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const float *pBase, float bScale)
{
const float *metrics = fMetrics + RMETRIC;
float chi2 = 0;
for(int nn = 0; nn < loopSamples; nn++) {
float diff = pReal[nn] - /*baseScale * */pBase[nn];
bScale *= RESCALE;
for(int nn = 0; nn < LOOP_SAMP; nn++) {
float fdiff = pReal[nn] - bScale * pBase[nn];
int idiff = int(fdiff);
#ifdef SPECMETRICS
// This block is to check that the size of "metrics", is big enough to address all
// possible values of "diff". If enabled, it slows very much down the execution.
// As nDhist is not thread safe, be sure that you are NOT USING THREADS.
int idiff = int(diff + RMETRIC);
if(specMetrics) specMetrics->incr(0, idiff);
if(idiff<=0 || idiff>=2*RMETRIC-1) {
int nerrr = 1; // to catch the outliers with the debugger
}
#endif
chi2 += metrics[int(diff/* * diffScale*/)];
//chi2 += pow(fabs(diff), float(METRIC)); // direct use of the pow function --> very slow
//chi2 += diff*diff; // direct calculation of pow( d , 2 )
//chi2 += fabs(diff); // direct calculation of pow(|d|, 1 )
//chi2 += sqrt(fabs(diff)); // direct calculation of pow(|d|, 1/2)
//chi2 += sqrt(sqrt(fabs(diff))); // direct calculation of pow(|d|, 1/4)
}
return chi2;
}
INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const int *pReal, const int *pBase)
{
const float *metrics = fMetrics + RMETRIC;
float chi2 = 0;
for(int nn = 0; nn < loopSamples; nn++) {
int diff = pReal[nn] - pBase[nn];
#ifdef DORESCALE
diff >>= RESCALE;
#endif
#ifdef SPECMETRICS // be sure that specMetrics is defined and that you are NOT USING THREADS
int idiff = int(diff + RMETRIC);
specMetrics->incr(0, idiff);
if(idiff<=0 || idiff>=2*RMETRIC-1) {
int nerrr = 1; // to catch the outliers with the debugger
}
#endif
chi2 += metrics[diff];
}
return chi2;
}
INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const short *pReal, const short *pBase)
{
const float *metrics = fMetrics + RMETRIC;
float chi2 = 0;
for(int nn = 0; nn < loopSamples; nn++) {
short diff = pReal[nn] - pBase[nn];
#ifdef SPECMETRICS // be sure that specMetrics is defined and that you are NOT USING THREADS
int idiff = int(diff + RMETRIC);
specMetrics->incr(0, idiff);
if(idiff<=0 || idiff>=2*RMETRIC-1) {
if(specMetrics) specMetrics->incr(idiff+RMETRIC);
if(idiff<=-RMETRIC+2 || idiff> RMETRIC-2) {
int nerrr = 1; // to catch the outliers with the debugger
}
#endif
chi2 += metrics[diff];
chi2 += metrics[idiff];
//chi2 += pow(fabs(diff), float(PMETRIC)); // direct use of the pow function --> very slow
//chi2 += diff*diff; // direct calculation of pow( d , 2 )
//chi2 += fabs(diff); // direct calculation of pow(|d|, 1 )
//chi2 += sqrt(fabs(diff)); // direct calculation of pow(|d|, 1/2)
//chi2 += sqrt(sqrt(fabs(diff))); // direct calculation of pow(|d|, 1/4)
//const float wwww = -1.f/(2*30.f*30.f);
//chi2 += 1.f - exp(diff*diff*wwww); // direct use of the 1-exp(d^2) function --> very slow
}
return chi2;
}
// loop with 2 hits
// the weight of the two components must be set by a preceeding call to SetLoopFactors(ff1, ff2)
INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const float *pBase1, const float *pBase2)
INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const float *pBase1, float bScale1, const float *pBase2, float bScale2)
{
const float *metrics = fMetrics + RMETRIC;
float chi2 = 0;
for(int nn = 0; nn < loopSamples; nn++) {
float diff = pReal[nn] - (loopFac1*pBase1[nn] + loopFac2*pBase2[nn]);
#ifdef SPECMETRICS // be sure that specMetrics is defined and that you are NOT USING THREADS
int idiff = int(diff + RMETRIC);
specMetrics->incr(0, idiff);
if(idiff<=0 || idiff>=2*RMETRIC-1) {
int nerrr = 1; // to catch the outliers with the debugger
}
#endif
//chi2 += metrics[int(diff)];
//chi2 += sqrt(fabs(diff)); // to test direct calculation of the metrics
chi2 += diff*diff; // direct calculation of pow( d , 2 )
bScale1 *= RESCALE;
bScale2 *= RESCALE;
for(int nn = 0; nn < LOOP_SAMP; nn++) {
float fdiff = pReal[nn] - (bScale1*pBase1[nn] + bScale2*pBase2[nn]);
int idiff = int(fdiff);
chi2 += metrics[idiff];
//chi2 += sqrt(fabs(fdiff)); // to test direct calculation of the metrics
//chi2 += fdiff*fdiff; // direct calculation of pow( d , 2 )
}
return chi2;
}
INLINE_ALWAYS float PSAFilterGridSearch::Chi2OuterLoop(const float **pAmpliReal, const float **pAmpliBase, const int nss, const float chi2min)
{
float chi2 = 0;
for(int iss = 0; iss < nss; iss++) {
const float *pReal = pAmpliReal[iss];
const float *pBase = pAmpliBase[iss];
for(int nn = 0; nn < loopSamples; nn++) {
float diff = pReal[nn] - pBase[nn];
chi2 += diff*diff; // direct calculation of pow( d , 2 )
}
if(chi2 > chi2min)
break;
}
return chi2;
}
#endif // #ifdef USE_SSE_VERSION
#endif // #else // #ifdef USE_SSE_VERSION
#endif // ADF_PSAFILTERGRID_SEARCH_H
......@@ -17,14 +17,14 @@
#endif
#endif
// Represents a single interaction point and is used for the signals of the basis and also for experimental data.
// PointPsa represents a single interaction point and is used for the signals of the basis and also for experimental data.
// Because of this, the experimental traces cannot be longer than NTIME samples
// For each signal we keep TCHAN = NSEGS+3 single traces:
// - the first NSEGS traces are for the segment.
// - the 3nd last trace [at TCHAN-3==INDCC==NSEGS] is the (possibly differentiated) core
// - the 2nd last trace [at TCHAN-2] is a non-differentiated copy of the net charge segment.
// - the last trace [at TCHAN-1] is a non-differentiated copy of the core.
// When used for the PSA basis the traces are usually normalized to MAXNORM and the net charge segment is already differentiated.
// The basis traces are normalized to MAXNORM (normally 1.f) and the net charge segment is already differentiated.
// To use the SSE extensions, the structure is aligned to a 16-byte boundary and NTIME must be a multiple of 4.
class pointPsa
{
......@@ -32,7 +32,7 @@ public:
float x, y, z; // coordinates of the interaction (mm)
int netChargeSegment; // its net-charge segment
MEMALIGN_16 // compiler-specific alignment directive
gs_type amplitude[TCHAN][NTIME]; // 36 segments + CC + non differentiated copies of net-charge segment and the core
float amplitude[TCHAN][NTIME]; // 36 segments + CC + non differentiated copies of net-charge segment and the core
int addr_first; // valid/useful start of the amplitude
int addr_last; // valid/useful last of the amplitude (margin included)
float sumTransients; // sum of all transient signals
......@@ -125,14 +125,14 @@ public:
int idt = (int)dt;
int nt1 = NTIME - idt - 1;
float xdt = dt - idt;
gs_type temporary[NTIME];
gs_type *pD = amplitude[nc] + idt;
float temporary[NTIME];
float *pD = amplitude[nc] + idt;
for(int ii = 0; ii < nt1; ii++, pD++) {
float val = pD[0] + (pD[1]-pD[0])*xdt;
temporary[ii] = gs_type(val);
temporary[ii] = val;
}
memcpy(amplitude[nc], temporary, sizeof(gs_type)*nt1);
gs_type val = temporary[nt1-1];
memcpy(amplitude[nc], temporary, sizeof(float)*nt1);
float val = temporary[nt1-1];
pD = amplitude[nc] + nt1;
for(int ii = nt1; ii < NTIME; ii++)
*pD++ = val;
......@@ -142,14 +142,13 @@ public:
int idt = (int)dt;
float xdt = idt + 1 - dt;
int nt0 = idt + 1;
gs_type temporary[NTIME] = {0}; // initialization to zero not really needed
gs_type *pD = amplitude[nc];
float temporary[NTIME] = {0}; // initialization to zero not really needed
float *pD = amplitude[nc];
for(int ii = nt0; ii < NTIME; ii++, pD++) {
float val = pD[0] + (pD[1]-pD[0])*xdt;
temporary[ii] = gs_type(val);
temporary[ii] = pD[0] + (pD[1]-pD[0])*xdt;
}
memcpy(amplitude[nc]+nt0, temporary, sizeof(gs_type)*(NTIME-nt0));
gs_type val = temporary[nt0];
memcpy(amplitude[nc]+nt0, temporary, sizeof(float)*(NTIME-nt0));
float val = temporary[nt0];
pD = amplitude[nc];
for(int ii = 0 ; ii < nt0; ii++)
*pD++ = val;
......@@ -175,7 +174,7 @@ public:
void diffNetCharge(int lag) {
if(lag > 0) {
gs_type *ptdata = amplitude[netChargeSegment];
float *ptdata = amplitude[netChargeSegment];
// the first lag channels stay as they are as their reference is zero
for(int nn = NTIME-1; nn >= lag; nn--) {
ptdata[nn] -= ptdata[nn - lag];
......@@ -184,13 +183,13 @@ public:
}
void saveNCandCC() {
memcpy(amplitude[INDNC], amplitude[netChargeSegment], sizeof(gs_type)*NTIME);
memcpy(amplitude[INDQC], amplitude[INDCC], sizeof(gs_type)*NTIME);
memcpy(amplitude[INDNC], amplitude[netChargeSegment], sizeof(float)*NTIME);
memcpy(amplitude[INDQC], amplitude[INDCC], sizeof(float)*NTIME);
}
void restoreNCandCC() {
memcpy(amplitude[netChargeSegment], amplitude[INDNC], sizeof(gs_type)*NTIME);
memcpy(amplitude[INDCC], amplitude[INDQC], sizeof(gs_type)*NTIME);
memcpy(amplitude[netChargeSegment], amplitude[INDNC], sizeof(float)*NTIME);
memcpy(amplitude[INDCC], amplitude[INDQC], sizeof(float)*NTIME);
}
// moving average filter
......@@ -198,20 +197,20 @@ public:
if(nn<2 || nn>=NTIME)
return;
float fact = 1.f/nn;
gs_type temp[NTIME];
float temp[NTIME];
for(int ns = 0; ns < NCHAN; ns++) { // including the CC
gs_type *pampl = amplitude[ns];
gs_type ampl0 = pampl[0];
float *pampl = amplitude[ns];
float ampl0 = pampl[0];
float sum = float(nn*ampl0);
int i1 = -nn;
int i2 = 0;
for(; i1 < 0; i1++, i2++) {
sum += pampl[i2] - ampl0;
temp[i2] = gs_type(sum*fact);
temp[i2] = sum*fact;
}
for(; i2 < NTIME; i1++, i2++) {
sum += pampl[i2] - pampl[i1];
temp[i2] = gs_type(sum*fact);
temp[i2] = sum*fact;
}
memcpy(amplitude[ns], temp, sizeof(temp));
}
......@@ -219,7 +218,7 @@ public:
// convolution using a given kernel
void conv(double* kernel, int kernel_len, int chan, int shift) {
gs_type temp[2*NTIME];
float temp[2*NTIME];
for (int i = 0; i < NTIME; i++) {
double accu = 0; // accumulator register
int len = kernel_len > i ? i : kernel_len;
......@@ -227,10 +226,10 @@ public:
for (int j = 0; j < len; j++) {
accu += kernel[j] * amplitude[chan][i - j];
}
temp[i] = (gs_type)accu;
temp[i] = float(accu);
}
for(int nn = NTIME; nn < 2* NTIME; nn++) temp[nn] = temp[NTIME-2];
memcpy(amplitude[chan], temp+shift, sizeof(gs_type)*NTIME);
memcpy(amplitude[chan], temp+shift, sizeof(float)*NTIME);
}
void convsegs(double* kernel, int kernel_len, int shift) {
......@@ -252,7 +251,7 @@ public:
float ss = 0;
for(int nn = 0; nn < len; nn++) {
ss += vv[nn];
amplitude[sg][nn] = (gs_type)ss;
amplitude[sg][nn] = ss;
}
}
......@@ -264,35 +263,36 @@ public:
float vv = 0;
for(int nn = 0; nn < NTIME; nn++) {
vv = gg*amplitude[sg][nn] + aa*vv;
amplitude[sg][nn] = (gs_type)vv;
amplitude[sg][nn] = vv;
}
}
}
// delta->exp, one value for each slice and the CC
// delta->exp, one value for each segment and the CC
void convDeltaToExp(const float *Tau) {
for(int sg = 0; sg < NCHAN; sg++) { // including CC
float tau = (sg==INDCC) ? Tau[6] : Tau[sg%6];
float tau = Tau[sg];
float aa = (float)exp(-TIME_STEP/tau);
float gg = 1.f-aa;
float vv = 0;
gs_type * psg = amplitude[sg];
float * psg = amplitude[sg];