Commit b1de6255 authored by dino's avatar dino
Browse files

T0Fit of T0 in PSAFilterGridSearch done on the trace of all segments+core

Number of fitting cycles set to 1 as a default

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@1068 170316e4-aea8-4b27-aad4-0380ec0519c9
parent d2fe761a
......@@ -28,7 +28,9 @@ const int numXYZ = 10; // sub PsaMatrXYZR for up to numXYZ-1 hits
const size_t nPsaHitsValues = 16; // number of parameters written in Psa__0-nn-F__Hits.fdat
const int numFit_T0 = 4; // max number of GridSearch--FitT0AfterPSA loops (0 == no fit)
const int numFit_T0 = 1; // max number of GridSearch--FitT0AfterPSA loops (0 == no T0 fit); 1 should be sufficient
const float shiftDTmin = 0.20f; // samples; stop looping if dT0 smaller than this
const float shiftT0max = 4.00f; // samples; stop looping if TotalShift lager than this
const float tauDefault = 35.f; // default time constant of the preamp response function for segments and core (in nanoseconds)
......@@ -623,7 +625,7 @@ void PSAFilterGridSearch::SegmentMapTest() {
}
}
// Process() handles "count" (~100) events for each call to reduce the overhead of managing parallel theads of execution.
// Handles "count" (~100) events for each call to reduce the overhead of managing parallel theads of execution.
// The logical structure of Process is a loop of Preparation-Execution-Analysis on the "count" events.
Int_t PSAFilterGridSearch::Process(int slot, int count)
......@@ -652,28 +654,27 @@ Int_t PSAFilterGridSearch::Process(int slot, int count)
}
// T0 of the event is calculated after the GridSearch by minimizing the Chi2 of the
// core signal as a function of the time-zero position of the experimental trace.
// Core+NetChargeSegments signal as a function of the time-zero position of the experimental trace.
// If the T0 change is significant (but not out of some limits), GridSearch is repeated...
// The problem with this approach is that all iteration perform a full search,
// ignoring the result of the previous one.
// More logical (and probably much faster) would be that after the first turn,
// a local search (fit?) is performed, starting from the previous positions.
const int rptmax = max(1, numFit_T0); // controls looping search -- fitT0
const float dT0min = 0.2f;
const float shiftT0max = 4.00f; // samples
for(int evn = 0; evn < count; evn++) { // loop on the events
int nCalled = 0;
int chiDone = 0;
bool extraLoop = false;
bool doTwoHits = false;
if(rptmax<2 && gTryTwoHits && gUseAdaptive) {
bool extraLoop = false; // the final loop in which also the search for two hits can be done
bool doTwoHits = false; // looking for two hits in a segment is done in the extraloop
if(numFit_T0<1) { // no T0 fit --> execute the final loop rightaway
extraLoop = true;
doTwoHits = true;
if(gTryTwoHits && gUseAdaptive)
doTwoHits = true; // looking also for two hits, if requested
}
for(int rpt = 0; (rpt<rptmax) || extraLoop; rpt++) { // loop on T0, but also the extra loop for the 2-hits search
for(int rpt = 0; rpt<numFit_T0 || extraLoop; rpt++) { // loop on T0, but also the extra loop for the 2-hits search
// 1111111111111111111111111111 prepare one event
int iSlot = slot + evn;
......@@ -704,7 +705,7 @@ Int_t PSAFilterGridSearch::Process(int slot, int count)
WRITEWORKINGWAVE((float *)PF.wAmplitude, true);
// 2222222222222222222222222222 call the grid-search algo for one event
int rv = ProcessTheEvent(PF, doTwoHits);
int rv = ProcessTheEvent(PF, (rpt==0), doTwoHits);
chiDone += abs(rv);
nCalled++;
WRITEWORKINGWAVE((float *)PF.rAmplitude, true);
......@@ -725,36 +726,27 @@ Int_t PSAFilterGridSearch::Process(int slot, int count)
if(extraLoop)
break;
// 4444444444444444444444444444 adjust T0, if requested
if(numFit_T0 == 0) {
if(gTryTwoHits && gUseAdaptive) {
extraLoop = true;
doTwoHits = true;
continue;
}
else
break;
}
// 4444444444444444444444444444 adjust T0
float dT0 = FitT0AfterPSA(PF);
pSlot->shiftT0 += dT0; // Units of samples. Accumulated if different runs
// stop looping if shift too small or total-shift too large
if((rpt==rptmax-1) || abs(dT0)<dT0min || abs(pSlot->shiftT0)>shiftT0max) {
float DT = FitT0AfterPSA(PF);
pSlot->shiftT0 += DT; // Units of samples. Accumulated if different runs
// stop looping if all requested loops done or shift too small or total-shift too large
if( rpt>=numFit_T0-1 || abs(DT)<shiftDTmin || abs(pSlot->shiftT0)>shiftT0max ) {
extraLoop = true;
if(gTryTwoHits && gUseAdaptive) {
doTwoHits = true;
}
continue;
if(gTryTwoHits && gUseAdaptive)
doTwoHits = true; // looking also for two hits, if requested
}
} // for(int rpt = 0; (rpt<rptmax) || extraLoop; rpt++)
} // for(int rpt = 0; (rpt<numFit_T0) || extraLoop; rpt++)
} // for(int evn = 0; evn < count; evn++)
return 0;
}
int PSAFilterGridSearch::ProcessTheEvent(pointFull &PF, bool twoHits)
// doFull==false to make a local search around the previous solution instead of a real coarse-fine search
// (to be done)
int PSAFilterGridSearch::ProcessTheEvent(pointFull &PF, bool doFull, bool twoHits)
{
int nCalled = 0;
int chiDone = 0;
......@@ -800,6 +792,7 @@ int PSAFilterGridSearch::ProcessTheEvent(pointFull &PF, bool twoHits)
// this must become a setup parameter from the configuration file
const bool onlyFirst = false; // just to see the effect of placing the results at the center of the segment
int rv = SearchAdaptive1(PF, netChSeg, onlyCoarse, onlyFirst); // coarse-fine or coarse-only with one hit
chiDone += abs(rv);
nCalled++;
//if(simplified || gCoarseOnly) {
......@@ -1048,7 +1041,7 @@ int PSAFilterGridSearch::ProcessTwoHits(pointFull &PF)
return chiDone;
}
void PSAFilterGridSearch::PrepareEvent(PsaSlot *pSlot, pointFull &PF, bool reset, bool traces)
void PSAFilterGridSearch::PrepareEvent(PsaSlot *pSlot, pointFull &PF, bool doReset, bool useTraces)
{
if(pSlot->crystal_status) {
// not a good event (is this still possible ?)
......@@ -1070,58 +1063,65 @@ void PSAFilterGridSearch::PrepareEvent(PsaSlot *pSlot, pointFull &PF, bool reset
}
#endif
if(reset)
PF.Reset(); // clear everything
PF.enerGe = pSlot->CoreE[0]; // copy Core energy
if(doReset) {
PF.Reset(); // clear everything
PF.enerGe = pSlot->CoreE[0]; // copy Core energy
#ifdef SELECTEVENTS
PF.isSelected = false; // set to false if selection foreseen
PF.isSelected = false; // set to false if selection foreseen
#else
PF.isSelected = true; // otherwise assume selected
PF.isSelected = true; // otherwise assume selected
#endif
//find hit segments
PF.sumSegm = 0;
Int_t numsegs = 0;
int *segOrder = PF.netChargeSegnum;
float *eneOrder = PF.netChargeEnergy;
for(Int_t isg=0; isg<NSEGS; isg++) {
if(pSlot->SegE[isg] > fHitSegThreshold) {
segOrder[numsegs] = isg;
eneOrder[numsegs] = pSlot->SegE[isg];
numsegs++;
PF.sumSegm += pSlot->SegE[isg];
//find hit segments
PF.sumSegm = 0;
Int_t numsegs = 0;
int *segOrder = PF.netChargeSegnum;
float *eneOrder = PF.netChargeEnergy;
for(Int_t isg=0; isg<NSEGS; isg++) {
if(pSlot->SegE[isg] > fHitSegThreshold) {
segOrder[numsegs] = isg;
eneOrder[numsegs] = pSlot->SegE[isg];
numsegs++;
PF.sumSegm += pSlot->SegE[isg];
}
}
if(numsegs < 1 || numsegs < fPsaMinSegMult || numsegs > fPsaMaxSegMult) {
PF.isValid = false;
return;
}
PF.numNetCharges = numsegs;
PF.netSeg = -1;
//PF.netChargeSegment = PF.netChargeSegnum[0];
//Core = Total net-seg ??? --> check what this test means for low energies
if( (PF.sumSegm > 0.9f*PF.enerGe) && (PF.sumSegm < 1.1f*PF.enerGe) )
PF.isFullCharge = true;
else
PF.isFullCharge = false;
PF.isValid = true;
if(useTraces && numsegs>1) {
// order the net-charge segments according to the released energy (not done for GridSearch==SegCenter)
for(int n1 = numsegs-1; n1 >=0; n1--) {
for(int n2 = 0; n2 < n1; n2++) {
if(eneOrder[n2] < eneOrder[n2+1]) {
swap(segOrder[n2], segOrder[n2+1]);
swap(eneOrder[n2], eneOrder[n2+1]);
}
}
}
}
}
if(numsegs < 1 || numsegs < fPsaMinSegMult || numsegs > fPsaMaxSegMult) {
PF.isValid = false;
return;
else { // most of setup already done in the previous turn
#ifdef SELECTEVENTS
PF.isSelected = false; // set to false if selection foreseen
#else
PF.isSelected = true; // otherwise assume selected
#endif
PF.netSeg = -1;
PF.isValid = true;
}
PF.numNetCharges = numsegs;
PF.netSeg = -1;
//PF.netChargeSegment = PF.netChargeSegnum[0];
//Core = Total net-seg ??? --> check what this test means for low energies
if( (PF.sumSegm > 0.9f*PF.enerGe) && (PF.sumSegm < 1.1f*PF.enerGe) )
PF.isFullCharge = true;
else
PF.isFullCharge = false;
PF.isValid = true;
if(!traces)
if(!useTraces)
return;
// order the net-charge segments according to the released energy (not done for GridSearch==SegCenter)
for(int n1 = numsegs-1; n1 >=0; n1--) {
for(int n2 = 0; n2 < n1; n2++) {
if(eneOrder[n2] < eneOrder[n2+1]) {
swap(segOrder[n2], segOrder[n2+1]);
swap(eneOrder[n2], eneOrder[n2+1]);
}
}
}
// Now fill tAmplitude with up to NTIME samples (including the pretrigger)
// and considering the possible time shift given by pSlot->shiftT0
......@@ -1736,7 +1736,7 @@ int PSAFilterGridSearch::SearchFullGrid(pointFull &PF, int netChSeg)
int PSAFilterGridSearch::SearchAdaptive1(pointFull &PF, int netChSeg, bool bCoarseOnly, bool bFirstOnly)
{
int time_first = TIME0; // start of basis
int time_first = TIME0; // start of basis
int samp_first = PF.samp_first; // start of signal
float baseScale = PF.baseScale; // scaling signals to data
char *lMask = PF.localMask; // selection of segments
......@@ -2658,44 +2658,27 @@ void PSAFilterGridSearch::MatChi2(float chi, pointPsa* pPsa)
}
#endif // #ifdef MATRCHI2
// not (yet?) working well
// use the complete trace
// returns the time shift with respect to SAMP0 in units of samples
Float_t PSAFilterGridSearch::FitT0AfterPSA(pointFull &PF)
{
float tData[NTIME];
float tRfit[NTIME];
// core
float * data = PF.tAmplitude[NSEGS];
float * rfit = PF.rAmplitude[NSEGS];
for(int ii = 0; ii < NTIME; ii++) {
tData[ii] = data[ii];
tRfit[ii] = rfit[ii];
}
// plus net-charge segments
// !!! does not work if SHOWDIFFERENTIATED (but we could extend rAmplitude to contain a non-differentiated copy of the net-charge)
for(int ns = 0; ns < PF.numNetCharges; ns++) {
int nss = PF.netChargeSegnum[ns];
data = PF.tAmplitude[nss];
rfit = PF.rAmplitude[nss];
for(int ii = 0; ii < NTIME; ii++) {
tData[ii] += data[ii];
tRfit[ii] += rfit[ii];
}
}
rfit = tRfit + SAMP0/2;
const int nsteps = 3;
const int nsamps = SAMP0/2+20; // calculate chi2 using SAMP0/2 of pretrigger and 20 of risetime
float x2[2*nsteps+1]; //memset(x2, 0, sizeof(x2));
{
// calculate chi2 using SAMP0/2 of pretrigger and 20-SAMP0/2 of risetime
const int nsteps = 3; // <= SAMP0/2
const int nsamps = 20;
float x2[2*nsteps+1]; memset(x2, 0, sizeof(x2));
for(int ll = 0; ll <= 2*nsteps; ll++) {
float *data = tData + SAMP0/2 - nsteps + ll;
float x2t = 0;
for(int ii = nsteps; ii < nsamps; ii++) {
float diff = data[ii] - rfit[ii];
x2t += diff*diff;
}
x2[ll] = x2t;
float chi2 = 0;
for(int iSegm=0; iSegm<NCHAN; iSegm++) {
float * tData = PF.tAmplitude[iSegm] + SAMP0/2 - nsteps + ll;
float * rData = PF.rAmplitude[iSegm] + SAMP0/2;
for(int ii = 0; ii < nsamps; ii++) { // ii = nsteps
float diff = tData[ii] - rData[ii];
chi2 += diff*diff;
}
} // end loop over the segments
x2[ll] = chi2;
}
// find the minimum
......@@ -2718,7 +2701,7 @@ Float_t PSAFilterGridSearch::FitT0AfterPSA(pointFull &PF)
float yp = x2[ssmin+1];
float yz = yp+ym-2.f*ys;
const float flim = 1.e-6f;
const float flim = 1.f; // determined by inspection
if(abs(yz) < flim) // to avoid dividing by 0
return float(ssmin-nsteps);
......
......@@ -121,7 +121,7 @@ private:
void PrepareEvent (PsaSlot *pSlot, pointFull &PF, bool reset, bool traces);
void SetToSegCenter (PsaSlot *pSlot, pointFull &PF);
void ShiftMoveTrace (float *data, int nsamp, float delta, float *accu);
int ProcessTheEvent(pointFull &PF, bool twoHits = false);
int ProcessTheEvent(pointFull &PF, bool doFull, bool twoHits);
int ProcessTwoHits (pointFull &PF);
void InitialSolution(pointFull &PF, bool keep = false);
void PreSearchCoarse(pointFull &PF);
......
......@@ -24,7 +24,7 @@
// - 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.
// The basis traces are normalized to MAXNORM (normally 1.f) and the net charge segment is already differentiated.
// The basis traces are normalized to MAXNORM (normally 1000.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
......
......@@ -37,10 +37,10 @@
# 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
# 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
......
......@@ -885,7 +885,9 @@ bool SignalBasis::CoarseFineForOneSegment(int ns, bool both, bool verbose)
cflist[ns].indFine[indF++] = jcc;
nAss[jcc]++;
for(int jj = 0; jj < npts; jj++) {
if(jj == jcc) continue; // pivot already inserted
if(jj == jcc) {
continue; // pivot already inserted
}
if( AssignIt(ppts[jj].x, ppts[jj].y, ppts[jj].z, cpx, cpy, cpz) ) {
if( !IsCoarse(jj, cgr, numC) ) {
cflist[ns].indFine[indF++] = jj;
......@@ -912,7 +914,9 @@ bool SignalBasis::CoarseFineForOneSegment(int ns, bool both, bool verbose)
cflist[ns].indFine[indF++] = jcc;
nAss[jcc] = -1;
for(int jj = 0; jj < npts; jj++) {
if(jj == jcc) continue; // pivot already inserted
if(jj == jcc) {
continue; // pivot already inserted
}
if( nAss[jj] == indC ) {
if( !IsCoarse(jj, cgr, numC) ) {
cflist[ns].indFine[indF++] = jj;
......
......@@ -162,7 +162,7 @@ inline bool SignalBasis::AssignIt(float px, float py, float pz, float cx, float
float d2 = dx*dx + dy*dy + dz*dz;
bool inside = ( d2 < (cstep+surface)*(cstep+surface) );
#else
// all points inside a centered on this coarse grid point
// all points inside a "cube" centered on this coarse grid point
float dx = px - cx; if(dx < 0) dx = -dx;
float dy = py - cy; if(dy < 0) dy = -dy;
float dz = pz - cz; if(dz < 0) dz = -dz;
......
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