Commit 51fc4ea3 authored by dino's avatar dino
Browse files

Default metric in PSAFilterGridSearch is back to 0.3

The signal basis can now be convoluted with the main decay pole of the preamplifier (keyword tauDecay)
Various changes in the flat-binary  spectra and matrices

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@1138 170316e4-aea8-4b27-aad4-0380ec0519c9
parent 29fb21a8
......@@ -179,7 +179,7 @@ const int defTriggerSample = 10; // 60-10=50 samples passed to the PSA
////// TrackingFilter //////
////////////////////////////
//#define TRF_ROOTTREE // enable the root tree in this actor (ENABLED IN THE MAKEFILE)
#define TRF_ROOTTREE // enable the root tree in this actor (ENABLED IN THE MAKEFILE)
//#define TRF_WRITE_TRACKED // enable writing the tracked gammas into the Oft_TrackedEnergies ascii file
//#define TRF_WRITE_GSORT // enable writing input and tracked data in gasp format
#define TFR_WRITE_HITS_MGT // enable writing the input hits in mgt format
......
......@@ -50,7 +50,10 @@ PSAFilter::PSAFilter() :
fDeadSegment = -1; // no one
memset(fTauSlice, 0, sizeof(fTauSlice));
fTauGiven = false;
fTauSliceGiven = false;
memset(fTauSegment, 0, sizeof(fTauSegment));
fTauSegmentGiven = false;
fTauDecay = 0;
theSlots = NULL;
#ifdef PSA_THREADED
......@@ -375,10 +378,24 @@ void PSAFilter::GetParameters(UInt_t *error_code)
// this should be taken from the general calibration file of the detector, PreprocessingFilterPSA.conf
KPAR("TauSlice", "%f %f %f %f %f %f %f", "timing response of the 6 slices and the CC");
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;
fTauSliceGiven = ok;
KEND;
KPAR("EnergyGain", "%f", "scaling factor for binary energy spectra");
KPAR("TauSegment", "%d %f", "timing response of the given segment (can be given multiple times, 36==CC");
int ns = -1;
float ts = 0;
ok = 2 == sscanf(data.c_str(), "%d %f", &ns, &ts);
if(ok && ns >=0 && ns <= CrystalInterface::kNbSegments) {
fTauSegment[ns] = ts;
fTauSegmentGiven = ok;
}
KEND;
KPAR("TauDecay", "%f", "fall time in ns of the preamplifiers (default is not to correct basis for this) ");
ok = 1 == sscanf(data.c_str(), "%f", &fTauDecay);
KEND;
KPAR("EnergyGain", "%f", "scaling factor for energy spectra");
ok = 1 == sscanf(data.c_str(), "%f", &fEnergyGain);
KEND;
......
......@@ -233,7 +233,11 @@ public:
Float_t fDistMetric;
Float_t fTauSlice[7];
Bool_t fTauGiven;
Bool_t fTauSliceGiven;
Float_t fTauSegment[ADF::CrystalInterface::kNbSegments+1]; // last one is CC
Bool_t fTauSegmentGiven;
Float_t fTauDecay;
protected:
ADF::GeSegment *seg;
......
......@@ -187,7 +187,7 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
//WriteBaseAverPt(1000.f/MAXNORM, "aver");
// values from the configuration file
if(fTauGiven) {
if(fTauSliceGiven) {
for(int sector = 0; sector < 6; sector++) {
for(int slice = 0; slice < 6; slice++) {
tauSGCC[6*sector+slice] = fTauSlice[slice];
......@@ -195,9 +195,16 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
}
tauSGCC[36] = fTauSlice[6];
}
// values from the configuration file
if(fTauSegmentGiven) {
for(int ns = 0; ns <= NSEGS; ns++) {
if(fTauSegment[ns] > 0)
tauSGCC[ns] = fTauSegment[ns];
}
}
// convolution with an exponential function
fBasis.ExponentialResponse(tauSGCC);
fBasis.ExponentialResponse(tauSGCC, fTauDecay); // ns
// convolution with the experimental response derived from the pulser
//fBasis.ExperimentalResponse();
......@@ -227,7 +234,7 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
const int numTraces = NCHAN;
MultiHist<short> *baseXYZ = new MultiHist<short>((2*detRadius)/fineStep, (2*detRadius)/fineStep, detLength/fineStep, numTraces, NTIME);
baseXYZ->setFileName(fOdirPrefix+"Psa?BaseXYZ.spec");
baseXYZ->setComment("the whole PSA base");
baseXYZ->setComment("the whole PSA signal basis");
int errcount = 0;
float bb[TCHAN][NTIME];
float fact = 30000.f/MAXNORM;
......@@ -404,6 +411,7 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
#endif // #else // #if 0
}
#if defined(USE_SSE_VERSION)
// to avoid useless work, this test should be done before reading the basis
if(gDistMetric != SSE_PMETRIC) {
cout << "\nERROR message from ";
DEBUG_LINE << endl;
......@@ -861,6 +869,7 @@ int PSAFilterGridSearch::ProcessTheEvent(pointFull &PF, bool doFull, bool twoHit
//if(fx>-12 && fx<-10 && fy>22 && fy<24 && netChSeg==13)
//if((netChSeg%6)==2 && fz>34 && fz<37)
//if(fz>16.f && fz<19.f && fx>6.f && fx<12.f)
//if(fz < 2.f && fy < -28.f/* && fx < -16.f*/)
if(fz < 2.f)
PF.isSelected = true; // the whole (possibly multihit) event is selected
//else
......@@ -1114,6 +1123,18 @@ void PSAFilterGridSearch::PrepareEvent(PsaSlot *pSlot, pointFull &PF, bool doRes
PF.isValid = false;
return;
}
#if 0
// hand selection of events
bool bOK = false;
for(int nOK = 0; nOK < numsegs; nOK++) {
if((segOrder[nOK]%6) == 0)
bOK = true; // first slice
}
if(!bOK) {
PF.isValid = false;
return;
}
#endif
PF.numNetCharges = numsegs;
PF.netSeg = -1;
//PF.netChargeSegment = PF.netChargeSegnum[0];
......@@ -1416,6 +1437,20 @@ int PSAFilterGridSearch::WriteTraces(PsaSlot *pSlot)
ptd[nn] -= aver;
}
}
// same thing for the core
{
float *ptd = pSlot->fAmplitude + NSEGS*WSAMP;
float aver = 0;
for(int nn = 0; nn < DIFFLAG; nn++)
aver += ptd[nn];
aver /= DIFFLAG;
for(int nn = wnsamp-1; nn >= DIFFLAG; nn--) {
ptd[nn] -= ptd[nn - DIFFLAG];
}
for(int nn = 0; nn < DIFFLAG; nn++) {
ptd[nn] -= aver;
}
}
}
#endif // SHOWDIFFERENTIATED
......@@ -2634,14 +2669,16 @@ void PSAFilterGridSearch::incrPosMats(int netChSeg, int nhi, float fx, float fy,
if(netChSeg < 0) {
if(PsaMatrXYZ) {
int slice90 = 90 - netChSeg; // negative value passed is -((netChSeg%6)+1
int slice90 = 90 - netChSeg; // negative value passed is -((netChSeg%6)+1)
PsaMatrXYZ->Incr(0, zAddrU(fz), xAddrU(fx), yAddrU(fy) );
PsaMatrXYZ->Incr(0, 0, xAddrU(fx), yAddrU(fy) ); // total projection
PsaMatrXYZ->Incr(0, slice90, xAddrU(fx), yAddrU(fy) );
PsaMatrXYZ->Incr(0, zAddrU(fz), xAddrU(fx), yAddrU(fy) ); // the slices in 91...96
PsaMatrXYZ->Incr(1, 0, xAddrU(fx), zAddrU(fz) ); // total projection
PsaMatrXYZ->Incr(0, slice90, xAddrU(fx), yAddrU(fy) ); // the slices in 91...96
PsaMatrXYZ->Incr(1, yAddrU(fy), xAddrU(fx), zAddrU(fz) );
PsaMatrXYZ->Incr(2, 0, yAddrU(fy), zAddrU(fz) ); // total projection
PsaMatrXYZ->Incr(1, 0, xAddrU(fx), zAddrU(fz) ); // total projection
PsaMatrXYZ->Incr(1, slice90, xAddrU(fx), zAddrU(fz) ); // the slices in 91...96
PsaMatrXYZ->Incr(2, xAddrU(fx), yAddrU(fy), zAddrU(fz) );
PsaMatrXYZ->Incr(2, 0, yAddrU(fy), zAddrU(fz) ); // total projection
PsaMatrXYZ->Incr(2, slice90, yAddrU(fy), zAddrU(fz) ); // the slices in 91...96
}
return;
......
......@@ -382,7 +382,7 @@ INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const
#if STANDARD_METRIC == STANDARD_METRIC_SQUARE
chi2 += fdiff*fdiff;
#else
int xdiff = int(fdiff);
int idiff = int(fdiff);
chi2 += metric[idiff];
#endif
}
......
......@@ -305,18 +305,60 @@ public:
}
}
// // 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 = Tau[sg];
// float aa = (float)exp(-TIME_STEP/tau);
// float gg = 1.f-aa;
// float vv = 0;
// float * psg = amplitude[sg];
// for(int nn = 0; nn < NTIME; nn++) {
// vv = gg*amplitude[sg][nn] + aa*vv;
// amplitude[sg][nn] = vv;
// }
//#if 1
// // gaussian smoothing 1-8-1 --> att. ~60% at 50 MHz
// // assuming continuation of the edge values
// float v2 = amplitude[sg][0];
// float v1 = 0;
// float v0 = 0;
// psg = amplitude[sg];
// for(int ii = 1; ii < NTIME; ii++) {
// v0 = v1; v1 = v2; v2 = psg[ii];
// psg[ii-1] = 0.1f * (v0 + 8.f*v1 + v2);
// }
// v0 = v1; v1 = v2;
// psg[NTIME-1] = 0.1f * (v0 + 8.f*v1 + v2);
//#endif
// }
// }
// delta->exp, one value for each segment and the CC
void convDeltaToExp(const float *Tau) {
void convDeltaToExp(const float *Tau, float tDecayLong = 0) {
for(int sg = 0; sg < NCHAN; sg++) { // including CC
float tau = Tau[sg];
float aa = (float)exp(-TIME_STEP/tau);
float gg = 1.f-aa;
float vv = 0;
float * psg = amplitude[sg];
//for(int nn = 0; nn < 10; nn++) psg[nn]=0.f;
//for(int nn = 10; nn < NTIME; nn++) psg[nn]=1000.f;
for(int nn = 0; nn < NTIME; nn++) {
vv = gg*amplitude[sg][nn] + aa*vv;
amplitude[sg][nn] = vv;
}
if(tDecayLong > 0) {
float tau2 = tDecayLong;
float aa2 = (float)exp(-TIME_STEP/tau2);
float vv1 = 0;
float vv2 = 0;
for(int nn = 0; nn < NTIME; nn++) {
vv2 = (amplitude[sg][nn] - vv1) + aa2*vv2;
vv1 = amplitude[sg][nn];
amplitude[sg][nn] = vv2;
}
}
#if 1
// gaussian smoothing 1-8-1 --> att. ~60% at 50 MHz
// assuming continuation of the edge values
......
......@@ -21,7 +21,7 @@
#define STANDARD_METRIC_2SQRT 4
#define STANDARD_METRIC_TUKEYS 5
#define STANDARD_METRIC STANDARD_METRIC_SQUARE
#define STANDARD_METRIC STANDARD_METRIC_NONE
#ifdef STANDARD_METRIC
# if STANDARD_METRIC == STANDARD_METRIC_ABSVAL // abs(v)
......@@ -35,12 +35,12 @@
# elif STANDARD_METRIC == STANDARD_METRIC_TUKEYS // Tuckey's biweight-like (1-exp(d^2))
const float PMETRIC = 0;
# elif STANDARD_METRIC == STANDARD_METRIC_NONE // generic, may be changed at program start, from the file PSAFilter.conf
const float PMETRIC = 0.3f; // value optimized by Roberto Venturelli
const float PMETRIC = 0.3f; // value optimized by Roberto Venturelli
# else // other values of STANDARD_METRIC
const float PMETRIC = 0.3f; // value optimized by Roberto Venturelli
const float PMETRIC = 0.3f; // value optimized by Roberto Venturelli
# endif //
#else // #ifdef STANDARD_METRIC // generic, may be changed at program start, from the file PSAFilter.conf
const float PMETRIC = 0.3f; // value optimized by Roberto Venturelli
const float PMETRIC = 0.3f; // value optimized by Roberto Venturelli
#endif // #ifdef STANDARD_METRIC
const int RMETRIC = 256*1024; // max tabulated deviation
......@@ -64,7 +64,7 @@ const float MAXNORM = 1000.f; // normalization of the signal basis
#endif
#endif
// the vectorized version only for the standard values .._MODULO, .._SQUARE, .._1SQRT, .._2SQRT
// the vectorized version only for the standard values .._ABSVAL, .._SQUARE, .._1SQRT, .._2SQRT
#if defined(USE_SSE_VERSION)
# if(STANDARD_METRIC != STANDARD_METRIC_ABSVAL) &&\
(STANDARD_METRIC != STANDARD_METRIC_SQUARE) &&\
......@@ -74,14 +74,14 @@ const float MAXNORM = 1000.f; // normalization of the signal basis
# endif
#endif
#define USE_NETCHARGES // to use the net charge signals in the grid search
#define USE_NETCHARGES // comment this line to use delay-differentiated net charge segments
#define USE_CORETRACE // to use also the core signal in the grid search
//#define SIMPLIFY_SMALL_ENERGIES // simplified search for low energy hits if segment multiplicity > 2
#ifdef USE_NETCHARGES
const int DIFFLAG = 0; // use the net-charge signals as they are
#else
// const int DIFFLAG = 0; // use the net-charge signals as they are
//#else
const int DIFFLAG = 8; // use the net-charge segments with a delayed-differentiation (in units of signal samples)
#endif
......
......@@ -1078,19 +1078,19 @@ void SignalBasis::SmoothBasisSignals(int nsmoo)
}
// Filter all basis signals with an exponential response of tau.
// TauSGCC contains the decay constants for the 36 slices and the core
void SignalBasis::ExponentialResponse(const float *TauSGCC)
// TauSGCC contains the decay constants for the 36 segments and the core
void SignalBasis::ExponentialResponse(const float *TauSGCC, const float TauDecay)
{
int count = 0;
for(int ii=0; ii<NSEGS; ii++) {
int npoints = numPts[ii];
for(int jj=0; jj<npoints; jj++) {
segPts[ii][jj].convDeltaToExp(TauSGCC);
segPts[ii][jj].convDeltaToExp(TauSGCC, TauDecay);
count++;
}
averPt[ii].convDeltaToExp(TauSGCC);
averPt[ii].convDeltaToExp(TauSGCC, TauDecay);
}
averPt[INDCC].convDeltaToExp(TauSGCC);
averPt[INDCC].convDeltaToExp(TauSGCC, TauDecay);
}
// Convolution on all signals with an experimental response function
......
......@@ -139,7 +139,7 @@ public:
int NearestCoarse(float px, float py, float pz, cgrid *cgr, int numC, float& dist);
void SmoothBasisSignals(int nsmoo);
void ExperimentalResponse();
void ExponentialResponse(const float *TauSGCC);
void ExponentialResponse(const float *TauSGCC, const float TauDecay = 0);
void CalculateSegmentCenters();
void FindNeighbours();
int StoreNetCharges();
......
......@@ -686,13 +686,15 @@ Int_t PostPSAFilter::Process()
float fy = pLoc->fy;
float fz = pLoc->fz;
int slice90 = 91 + (pLoc->Sg%6);
PostMatrXYZ->Incr(0, 0, xAddrU(fx), yAddrU(fy) ); // total projection
PostMatrXYZ->Incr(0, zAddrU(fz), xAddrU(fx), yAddrU(fy) );
PostMatrXYZ->Incr(0, 0, xAddrU(fx), yAddrU(fy) ); // total projection
PostMatrXYZ->Incr(0, slice90, xAddrU(fx), yAddrU(fy) ); // the 6 slices in 91...96
PostMatrXYZ->Incr(1, 0, xAddrU(fx), zAddrU(fz) ); // total projection
PostMatrXYZ->Incr(1, yAddrU(fy), xAddrU(fx), zAddrU(fz) );
PostMatrXYZ->Incr(2, 0, yAddrU(fy), zAddrU(fz) ); // total projection
PostMatrXYZ->Incr(1, 0, xAddrU(fx), zAddrU(fz) ); // total projection
PostMatrXYZ->Incr(1, slice90, xAddrU(fx), zAddrU(fz) ); // the 6 slices in 91...96
PostMatrXYZ->Incr(2, xAddrU(fx), yAddrU(fy), zAddrU(fz) );
PostMatrXYZ->Incr(2, 0, yAddrU(fy), zAddrU(fz) ); // total projection
PostMatrXYZ->Incr(2, slice90, yAddrU(fy), zAddrU(fz) ); // the 6 slices in 91...96
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment