Commit 395dca73 authored by dino's avatar dino
Browse files

Removed obsolete parts from PreprocessingFilterPSA and PSAFilterGridSearch

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@1076 170316e4-aea8-4b27-aad4-0380ec0519c9
parent f96d695a
...@@ -1027,10 +1027,6 @@ ...@@ -1027,10 +1027,6 @@
RelativePath="..\common\cycleServer.h" RelativePath="..\common\cycleServer.h"
> >
</File> </File>
<File
RelativePath="..\common\dataPairFile.h"
>
</File>
<File <File
RelativePath="..\common\misc.cpp" RelativePath="..\common\misc.cpp"
> >
......
...@@ -456,11 +456,19 @@ Int_t PSAFilter::SetOutput(int slot) ...@@ -456,11 +456,19 @@ Int_t PSAFilter::SetOutput(int slot)
CoreE0 = pSlot->CoreE[0]; CoreE0 = pSlot->CoreE[0];
CoreE1 = pSlot->CoreE[1]; CoreE1 = pSlot->CoreE[1];
CoreT0 = pSlot->CoreT[0];
CoreT1 = pSlot->CoreT[1];
//#define ADDPSATIME
#ifdef ADDPSATIME
// Add the T0 correction as determined by the PSA. --> The overall Ge-Ge timing becomes ~10% worse, with larger tails. // Add the T0 correction as determined by the PSA. --> The overall Ge-Ge timing becomes ~10% worse, with larger tails.
// If the T0 shift is subtracted the overall ge-ge timing becomes ~2 worse (but very symmetric, meaning added noise). // If the T0 shift is subtracted the overall ge-ge timing becomes ~2 worse (but very symmetric, meaning added noise).
// Therefore, it is better not to correct. However, one should also question the T0 realignment done in the PSA! // The above seems not to be true if the number of samples used to determine Tzero in PsaFilterGridSearch::FitT0AfterPSA
CoreT0 = pSlot->CoreT[0] /*+ pSlot->shiftT0*/; // is smal: e.g. at nsamps=10 we now get an improvement by 10%. To be tested with a wider dinamycal range.
CoreT1 = pSlot->CoreT[1] /*+ pSlot->shiftT0*/; // Therefore, it is better not to correct. MOreover, one should also question the T0 realignment done in the PSA!
CoreT0 += pSlot->shiftT0;
CoreT1 += pSlot->shiftT0;
#endif
//myAno.Set(crystal_id, 0); //myAno.Set(crystal_id, 0);
//myAno.fRealSize = sizeof(crystal_id); //myAno.fRealSize = sizeof(crystal_id);
......
...@@ -48,24 +48,27 @@ const float lowSegEnergy = 50.f; // segment-energies smaller than this are deco ...@@ -48,24 +48,27 @@ const float lowSegEnergy = 50.f; // segment-energies smaller than this are deco
// The output file debug__0-37-56-F__waves.fsamp is written in pwd and may or may not be rewound at the beginning of every event. // The output file debug__0-37-56-F__waves.fsamp is written in pwd and may or may not be rewound at the beginning of every event.
// They produce a large amount of output and are normally used in Debug mode as a graphical inspection tool. // They produce a large amount of output and are normally used in Debug mode as a graphical inspection tool.
// They are not thread-safe and, due to the fixed-name ouput file, incompatible with multiple crystal analyses. // They are not thread-safe and, due to the fixed-name ouput file, incompatible with multiple crystal analyses.
#if 1 && defined(_DEBUG) // set to 1 to enable writing the intermediate waves while debugging #if 0 && defined(_DEBUG) // set to 1 to enable writing the intermediate waves while debugging
# define WRITEWORKINGWAVE(p1,p2) WriteWorkingWave(p1, p2) # define WRITEWORKINGWAVE(p1,p2) WriteWorkingWave(p1, p2)
# define WRITEPARTIALWAVE(p1,p2,p3) WritePartialWave(p1, p2, p3) # define WRITEPARTIALWAVE(p1,p2,p3) WritePartialWave(p1, p2, p3)
#else #else
# define WRITEWORKINGWAVE(p1,p2) # define WRITEWORKINGWAVE(p1,p2)
# define WRITEPARTIALWAVE(p1,p2,p3) # define WRITEPARTIALWAVE(p1,p2,p3)
//# define WRITEFITTEDWAVES
# ifdef WRITEFITTEDWAVES
# define WRITEWORKINGWAVE1(p1,p2) WriteWorkingWave(p1, p2)
# endif
#endif #endif
PSAFilterGridSearch::PSAFilterGridSearch() PSAFilterGridSearch::PSAFilterGridSearch()
{ {
fpPsaTraces = NULL; fpPsaTraces = NULL;
fpPsaHits = NULL; fpPsaHits = NULL;
fPsaIndex = 0; fPsaIndex = 0;
#ifdef WRITE_TSHIFT_CORRECTION
tstmpFileO = NULL;
#endif
gUseAdaptive = true; gUseAdaptive = true;
gTryTwoHits = false; gTryTwoHits = false;
gCoarseOnly = false; gCoarseOnly = false;
...@@ -254,6 +257,56 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise() ...@@ -254,6 +257,56 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
} }
#endif #endif
#if 0 && defined(PSA_MULTIHIST)
{
// Write the complete basis of a segment a <float> spectrum. Normalization is 1000.
// write only NSEGS+CC as they are going to be used (i.e. with differentiated netChargeSeg and CC)
// To get the non-differentiated version uncomment //#define USE_NETCHARGES in FridSearchParameters.h
// To get the whole internal basis replace NCHAN with TCHAN in the definition of numTraces
const int fineStep = fstep;
const int detRadius = 40;
const int detLength = 90;
const int numTraces = NCHAN;
const int theSegment = 3;
int numPts = fBasis.numPts[theSegment];
char chsnum[20]; sprintf(chsnum, "%2.2d", theSegment);
MultiHist<float> *baseXYZ = new MultiHist<float>(numPts, numTraces, NTIME);
baseXYZ->setFileName(fOdirPrefix+"Psa?BaseXYZ.spec_"+std::string(chsnum));
baseXYZ->setComment("the PSA base of segment "+std::string(chsnum));
int errcount = 0;
float bb[TCHAN][NTIME];
float fact = 1000.f/MAXNORM;
int indSpe = 0;
for(int ns = theSegment; ns <=theSegment; ns++) {
for(int np = 0; np < fBasis.numPts[ns]; np++) {
pointPsa *bpt = &fBasis.Pts[ns][np];
memcpy(bb, bpt->amplitude, sizeof(bb));
// copy exactly as is
//for(int ii = 0; ii < NTIME; ii++) {
// float tmp = bb[ns][ii];
// bb[ns][ii] = bb[INDNC][ii];
// bb[INDNC][ii] = tmp;
// tmp = bb[INDCC][ii];
// bb[INDCC][ii] = bb[INDQC][ii];
// bb[INDQC][ii] = tmp;
//}
float *pd = baseXYZ->getData(indSpe++, 0);
if(pd) {
for (int is = 0; is < numTraces; is++) {
for(int it = 0; it < NTIME; it++) {
*pd++ = bb[is][it]*fact;
}
}
}
else
errcount++;
}
}
baseXYZ->write();
delete baseXYZ;
}
#endif
#if 0 && defined(PSA_MULTIHIST) #if 0 && defined(PSA_MULTIHIST)
{ {
// 3D distribution sum of segment signals (2 mm packed) // 3D distribution sum of segment signals (2 mm packed)
...@@ -366,12 +419,6 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise() ...@@ -366,12 +419,6 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
#endif #endif
} // else { // if(gPsaSegCenter) { } // else { // if(gPsaSegCenter) {
#ifdef WRITE_TSHIFT_CORRECTION
if(doFit_T0) {
tstmpFileO = new dataPairFile<ULong64_t, float>(fOdirPrefix+"Tadjust.fdat", false, true, 1000);
}
#endif
#ifdef PSA_MULTIHIST #ifdef PSA_MULTIHIST
if(fUseMultiHist) { if(fUseMultiHist) {
PsaSpecEner = new MultiHist<unsigned int>(40, specLenE); PsaSpecEner = new MultiHist<unsigned int>(40, specLenE);
...@@ -655,10 +702,10 @@ Int_t PSAFilterGridSearch::Process(int slot, int count) ...@@ -655,10 +702,10 @@ Int_t PSAFilterGridSearch::Process(int slot, int count)
return 0; // nothing more to do return 0; // nothing more to do
} }
// T0 of the event is calculated after the GridSearch by minimizing the Chi2 of the // T0 of the event is calculated after the GridSearch by minimizing the Chi2 of
// Core+NetChargeSegments signal as a function of the time-zero position of the experimental trace. // the 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... // 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, // The problem with this approach is that every iteration perform a full search,
// ignoring the result of the previous one. // ignoring the result of the previous one.
// More logical (and probably much faster) would be that after the first turn, // More logical (and probably much faster) would be that after the first turn,
// a local search (fit?) is performed, starting from the previous positions. // a local search (fit?) is performed, starting from the previous positions.
...@@ -684,6 +731,12 @@ Int_t PSAFilterGridSearch::Process(int slot, int count) ...@@ -684,6 +731,12 @@ Int_t PSAFilterGridSearch::Process(int slot, int count)
PrepareEvent(pSlot, PF, (rpt==0), true); // reset PF, identify the hits and prepare the traces, considering also shifT0 PrepareEvent(pSlot, PF, (rpt==0), true); // reset PF, identify the hits and prepare the traces, considering also shifT0
WRITEWORKINGWAVE((float *)PF.tAmplitude, (rpt!=0)); // false --> also rewind the file WRITEWORKINGWAVE((float *)PF.tAmplitude, (rpt!=0)); // false --> also rewind the file
#ifdef WRITEFITTEDWAVES
if(PF.numNetCharges != 1 || PF.netChargeSegnum[0] != 0) // only if segment A1
break;
WRITEWORKINGWAVE1((float *)PF.tAmplitude, true); // false --> also rewind the file
#endif
if(!PF.isValid) { if(!PF.isValid) {
pSlot->crystal_status = 1; // not to be used (anyway there are no hits) pSlot->crystal_status = 1; // not to be used (anyway there are no hits)
pSlot->numNetSegs = 0; pSlot->numNetSegs = 0;
...@@ -711,7 +764,9 @@ Int_t PSAFilterGridSearch::Process(int slot, int count) ...@@ -711,7 +764,9 @@ Int_t PSAFilterGridSearch::Process(int slot, int count)
chiDone += abs(rv); chiDone += abs(rv);
nCalled++; nCalled++;
WRITEWORKINGWAVE((float *)PF.rAmplitude, true); WRITEWORKINGWAVE((float *)PF.rAmplitude, true);
#ifdef WRITEFITTEDWAVES
WRITEWORKINGWAVE1((float *)PF.rAmplitude, true);
#endif
// 3333333333333333333333333333 move the search-results from PF to pSlot // 3333333333333333333333333333 move the search-results from PF to pSlot
pSlot->numNetSegs = PF.numNetCharges; pSlot->numNetSegs = PF.numNetCharges;
pSlot->isSelected = PF.isSelected; pSlot->isSelected = PF.isSelected;
...@@ -770,13 +825,8 @@ int PSAFilterGridSearch::ProcessTheEvent(pointFull &PF, bool doFull, bool twoHit ...@@ -770,13 +825,8 @@ int PSAFilterGridSearch::ProcessTheEvent(pointFull &PF, bool doFull, bool twoHit
PF.AddBaseTrace(bestPoint, scaleFact, PF.samp_first, TIME0); PF.AddBaseTrace(bestPoint, scaleFact, PF.samp_first, TIME0);
WRITEWORKINGWAVE((float *)PF.wAmplitude, true); WRITEWORKINGWAVE((float *)PF.wAmplitude, true);
// select the segments to use for the chi2 loop // prepare sAmplitude from the subset of active segments
MakeLocalMask(PF); MakeSearchWave(PF, 1.f, netChEner);
// Prepare sAmplitude from the subset of PF.wAmplitude given by localMask
// MakeSearchWave applies also the delayed differentiation of DIFFLAG samples (if defined)
// With a precalculated metric it must also limit the y-values to avoid overflowing the vector.
PF.MakeSearchWave(1.f, PF.localMask, netChEner);
WRITEPARTIALWAVE((float *)PF.sAmplitude, 1.f/RESCALE, true); // scaled down to get it in natural units WRITEPARTIALWAVE((float *)PF.sAmplitude, 1.f/RESCALE, true); // scaled down to get it in natural units
PF.chi2min = float(1.e30); PF.chi2min = float(1.e30);
...@@ -963,10 +1013,8 @@ int PSAFilterGridSearch::ProcessTwoHits(pointFull &PF) ...@@ -963,10 +1013,8 @@ int PSAFilterGridSearch::ProcessTwoHits(pointFull &PF)
PF.AddBaseTrace(bestPoint, scaleFact, PF.samp_first, TIME0); PF.AddBaseTrace(bestPoint, scaleFact, PF.samp_first, TIME0);
WRITEWORKINGWAVE((float *)PF.wAmplitude, true); WRITEWORKINGWAVE((float *)PF.wAmplitude, true);
// Prepare sAmplitude from the subset of PF.wAmplitude given by localMask // prepare sAmplitude from the subset of active segments
// Apply also the delayed differentiation of DIFFLAG samples (if defined) MakeSearchWave(PF, 1.f, netChEner);
MakeLocalMask(PF);
PF.MakeSearchWave(1.f, PF.localMask, netChEner);
WRITEPARTIALWAVE((float *)PF.sAmplitude, 1.f/RESCALE, true); WRITEPARTIALWAVE((float *)PF.sAmplitude, 1.f/RESCALE, true);
bool onlyCoarse = gCoarseOnly; bool onlyCoarse = gCoarseOnly;
...@@ -1667,7 +1715,7 @@ int PSAFilterGridSearch::SearchFullGrid(pointFull &PF, int netChSeg) ...@@ -1667,7 +1715,7 @@ int PSAFilterGridSearch::SearchFullGrid(pointFull &PF, int netChSeg)
chi2 = 0; chi2 = 0;
for(int iSegm=0; iSegm<NCHAN; iSegm++) { for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') { if(lMask[iSegm] != '0') {
float * realTrace = PF.sAmplitude[iSegm] + samp_first; float * realTrace = PF.sAmplitude[iSegm] + samp_first;
float * baseTrace = fBasis.Pts[netChSeg][iPts].amplitude[iSegm] + time_first; float * baseTrace = fBasis.Pts[netChSeg][iPts].amplitude[iSegm] + time_first;
chi2 += Chi2InnerLoop(realTrace, baseTrace, PF.baseScale*RESCALE); chi2 += Chi2InnerLoop(realTrace, baseTrace, PF.baseScale*RESCALE);
#ifndef MATRCHI2 // don't interrupt the loop over the segments while producing MatChi2 #ifndef MATRCHI2 // don't interrupt the loop over the segments while producing MatChi2
...@@ -1703,7 +1751,7 @@ int PSAFilterGridSearch::SearchFullGrid(pointFull &PF, int netChSeg) ...@@ -1703,7 +1751,7 @@ int PSAFilterGridSearch::SearchFullGrid(pointFull &PF, int netChSeg)
chi2 = 0; chi2 = 0;
for(int iSegm=0; iSegm<NCHAN; iSegm++) { for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') { if(lMask[iSegm] != '0') {
float * realTrace = PF.sAmplitude[iSegm] + samp_first; float * realTrace = PF.sAmplitude[iSegm] + samp_first;
float * baseTrace1 = fBasis.Pts[netChSeg][iPts1].amplitude[iSegm] + time_first; float * baseTrace1 = fBasis.Pts[netChSeg][iPts1].amplitude[iSegm] + time_first;
float * baseTrace2 = fBasis.Pts[netChSeg][iPts2].amplitude[iSegm] + time_first; float * baseTrace2 = fBasis.Pts[netChSeg][iPts2].amplitude[iSegm] + time_first;
chi2 += Chi2InnerLoop(realTrace, baseTrace1, bScale1*RESCALE, baseTrace2, bScale2*RESCALE); chi2 += Chi2InnerLoop(realTrace, baseTrace1, bScale1*RESCALE, baseTrace2, bScale2*RESCALE);
...@@ -1770,7 +1818,7 @@ int PSAFilterGridSearch::SearchAdaptive1(pointFull &PF, int netChSeg, bool bCoar ...@@ -1770,7 +1818,7 @@ int PSAFilterGridSearch::SearchAdaptive1(pointFull &PF, int netChSeg, bool bCoar
for(int iSegm=0; iSegm<NCHAN; iSegm++) { for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') { if(lMask[iSegm] != '0') {
nInner++; nInner++;
float * realTrace = PF.sAmplitude[iSegm] + samp_first; float * realTrace = PF.sAmplitude[iSegm] + samp_first;
float * baseTrace = pPtSeg[kPts].amplitude[iSegm] + time_first; float * baseTrace = pPtSeg[kPts].amplitude[iSegm] + time_first;
chi2 += Chi2InnerLoop(realTrace, baseTrace, baseScale*RESCALE); chi2 += Chi2InnerLoop(realTrace, baseTrace, baseScale*RESCALE);
if(chi2 > chi2min) if(chi2 > chi2min)
...@@ -1816,7 +1864,7 @@ int PSAFilterGridSearch::SearchAdaptive1(pointFull &PF, int netChSeg, bool bCoar ...@@ -1816,7 +1864,7 @@ int PSAFilterGridSearch::SearchAdaptive1(pointFull &PF, int netChSeg, bool bCoar
for(int iSegm=0; iSegm<NCHAN; iSegm++) { for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') { if(lMask[iSegm] != '0') {
nInner++; nInner++;
float * realTrace = PF.sAmplitude[iSegm] + samp_first; float * realTrace = PF.sAmplitude[iSegm] + samp_first;
float * baseTrace = pPtSeg[jPts].amplitude[iSegm] + time_first; float * baseTrace = pPtSeg[jPts].amplitude[iSegm] + time_first;
chi2 += Chi2InnerLoop(realTrace, baseTrace, baseScale*RESCALE); chi2 += Chi2InnerLoop(realTrace, baseTrace, baseScale*RESCALE);
if(chi2 > chi2min) if(chi2 > chi2min)
...@@ -1999,7 +2047,7 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar ...@@ -1999,7 +2047,7 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
for(int iSegm=0; iSegm<NCHAN; iSegm++) { for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') { if(lMask[iSegm] != '0') {
nInner++; nInner++;
float * realTrace = PF.sAmplitude[iSegm] + samp_first; float * realTrace = PF.sAmplitude[iSegm] + samp_first;
float * baseTrace1 = pPtSeg[kPts1].amplitude[iSegm] + time_first; float * baseTrace1 = pPtSeg[kPts1].amplitude[iSegm] + time_first;
float * baseTrace2 = pPtSeg[kPts2].amplitude[iSegm] + time_first; float * baseTrace2 = pPtSeg[kPts2].amplitude[iSegm] + time_first;
chi2 += Chi2InnerLoop(realTrace, baseTrace1, bScale1*RESCALE, baseTrace2, bScale2*RESCALE); chi2 += Chi2InnerLoop(realTrace, baseTrace1, bScale1*RESCALE, baseTrace2, bScale2*RESCALE);
...@@ -2047,7 +2095,7 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar ...@@ -2047,7 +2095,7 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
for(int iSegm=0; iSegm<NCHAN; iSegm++) { for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') { if(lMask[iSegm] != '0') {
nInner++; nInner++;
float * realTrace = PF.sAmplitude[iSegm] + samp_first; float * realTrace = PF.sAmplitude[iSegm] + samp_first;
float * baseTrace1 = pPtSeg[jPts].amplitude[iSegm] + time_first; float * baseTrace1 = pPtSeg[jPts].amplitude[iSegm] + time_first;
float * baseTrace2 = pPtSeg[fxPt].amplitude[iSegm] + time_first; float * baseTrace2 = pPtSeg[fxPt].amplitude[iSegm] + time_first;
chi2 += Chi2InnerLoop(realTrace, baseTrace1, bScale1*RESCALE, baseTrace2, bScale2*RESCALE); chi2 += Chi2InnerLoop(realTrace, baseTrace1, bScale1*RESCALE, baseTrace2, bScale2*RESCALE);
...@@ -2088,7 +2136,7 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar ...@@ -2088,7 +2136,7 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
for(int iSegm=0; iSegm<NCHAN; iSegm++) { for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') { if(lMask[iSegm] != '0') {
nInner++; nInner++;
float * realTrace = PF.sAmplitude[iSegm] + samp_first; float * realTrace = PF.sAmplitude[iSegm] + samp_first;
float * baseTrace1 = pPtSeg[fxPt].amplitude[iSegm] + time_first; float * baseTrace1 = pPtSeg[fxPt].amplitude[iSegm] + time_first;
float * baseTrace2 = pPtSeg[jPts].amplitude[iSegm] + time_first; float * baseTrace2 = pPtSeg[jPts].amplitude[iSegm] + time_first;
chi2 += Chi2InnerLoop(realTrace, baseTrace1, bScale1*RESCALE, baseTrace2, bScale2*RESCALE); chi2 += Chi2InnerLoop(realTrace, baseTrace1, bScale1*RESCALE, baseTrace2, bScale2*RESCALE);
...@@ -2142,7 +2190,7 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar ...@@ -2142,7 +2190,7 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
for(int iSegm=0; iSegm<NCHAN; iSegm++) { for(int iSegm=0; iSegm<NCHAN; iSegm++) {
if(lMask[iSegm] != '0') { if(lMask[iSegm] != '0') {
nInner++; nInner++;
float * realTrace = PF.sAmplitude[iSegm] + samp_first; float * realTrace = PF.sAmplitude[iSegm] + samp_first;
float * baseTrace1 = pPtSeg[kPts1].amplitude[iSegm] + time_first; float * baseTrace1 = pPtSeg[kPts1].amplitude[iSegm] + time_first;
float * baseTrace2 = pPtSeg[kPts2].amplitude[iSegm] + time_first; float * baseTrace2 = pPtSeg[kPts2].amplitude[iSegm] + time_first;
chi2 += Chi2InnerLoop(realTrace, baseTrace1, bScale1*RESCALE, baseTrace2, bScale2*RESCALE); chi2 += Chi2InnerLoop(realTrace, baseTrace1, bScale1*RESCALE, baseTrace2, bScale2*RESCALE);
...@@ -2196,6 +2244,33 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar ...@@ -2196,6 +2244,33 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
return chiDone; return chiDone;
} }
void PSAFilterGridSearch::MakeSearchWave(pointFull &PF, float fact, float maxVal)
{
// select the active segments
MakeLocalMask(PF);
// calculate sAmplitude
PF.MakeSearchWave(fact, PF.localMask);
#ifndef USE_SSE_VERSION
// if using the mapped-metric vector, limit result to avoid overruns
float yMax = RMETRIC - maxVal * RESCALE - 10.f; // some extra margin
float yMin = -RMETRIC + maxVal * RESCALE + 10.f; // some extra margin
for(int iSegm = 0; iSegm < NCHAN; iSegm++) {
if(PF.localMask[iSegm] != '0') {
float *sA = PF.sAmplitude[iSegm];
for(int kk=0; kk<NTIME; kk++) {
float yy = sA[kk];
if(yy > yMax)
sA[kk] = yMax;
else if(yy < yMin)
sA[kk] = yMin;
}
}
}
#endif
}
void PSAFilterGridSearch::MakeLocalMask(pointFull &PF) void PSAFilterGridSearch::MakeLocalMask(pointFull &PF)
{ {
// set mask of segments for the search loop // set mask of segments for the search loop
...@@ -2266,9 +2341,8 @@ void PSAFilterGridSearch::PreSearchCoarse(pointFull &PF) ...@@ -2266,9 +2341,8 @@ void PSAFilterGridSearch::PreSearchCoarse(pointFull &PF)
PF.AddBaseTrace(bestPoint, scaleFact, PF.samp_first, TIME0); PF.AddBaseTrace(bestPoint, scaleFact, PF.samp_first, TIME0);
WRITEWORKINGWAVE((float *)PF.wAmplitude, true); WRITEWORKINGWAVE((float *)PF.wAmplitude, true);
// prepare sAmplitude from the subset of wAmplitude given by localMask // prepare sAmplitude from the subset of active segments
MakeLocalMask(PF); MakeSearchWave(PF, 1.f, netChEner);
PF.MakeSearchWave(1.f, PF.localMask, netChEner);
WRITEPARTIALWAVE((float *)PF.sAmplitude, scaleFact, true); WRITEPARTIALWAVE((float *)PF.sAmplitude, scaleFact, true);
PF.chi2min = float(1.e30); PF.chi2min = float(1.e30);
...@@ -2320,11 +2394,6 @@ Int_t PSAFilterGridSearch::AlgoSpecificPostProcess(int slot) ...@@ -2320,11 +2394,6 @@ Int_t PSAFilterGridSearch::AlgoSpecificPostProcess(int slot)
fPsaIndex = (fPsaIndex+1)%16; fPsaIndex = (fPsaIndex+1)%16;
} }
#ifdef WRITE_TSHIFT_CORRECTION
if(doFit_T0 && tstmpFileO)
tstmpFileO->put(pSlot->timestamp, pSlot->shiftT0);
#endif
// no histogramming requested // no histogramming requested
if(!pSlot->isSelected) if(!pSlot->isSelected)
return 0; return 0;
...@@ -2662,19 +2731,28 @@ void PSAFilterGridSearch::MatChi2(float chi, pointPsa* pPsa) ...@@ -2662,19 +2731,28 @@ void PSAFilterGridSearch::MatChi2(float chi, pointPsa* pPsa)
} }
#endif // #ifdef MATRCHI2 #endif // #ifdef MATRCHI2
// use the complete trace // Use all segments and core to compare the "fitted" trace to the experimental one
// The comparison runs from SAMP0/2 (=-4) to nsams-SAMP0/2 of the fitted trace, while the
// experimental one is shifted for up to +-nsteps
// returns the time shift with respect to SAMP0 in units of samples // returns the time shift with respect to SAMP0 in units of samples
Float_t PSAFilterGridSearch::FitT0AfterPSA(pointFull &PF) Float_t PSAFilterGridSearch::FitT0AfterPSA(pointFull &PF)
{ {
// calculate chi2 using SAMP0/2 of pretrigger and 20-SAMP0/2 of risetime // calculate chi2 using SAMP0/2 of pretrigger and 20-SAMP0/2 of risetime
const int nsteps = 3; // <= SAMP0/2 const int nsteps = 4; // <= SAMP0/2
const int nsamps = 20; const int nsamps = 10;
float x2[2*nsteps+1]; memset(x2, 0, sizeof(x2));
const int x2size = 2*nsteps+1;
for(int ll = 0; ll <= 2*nsteps; ll++) { float x2[x2size];
for(int ll = 0; ll < x2size; ll++)
x2[ll] = -1;
int ssmin = -1;
float x2min = -1.f;
int ll = nsteps; // dl = ll-nsteps
while(ll>=0 && ll<=2*nsteps) {
float chi2 = 0; float chi2 = 0;
for(int iSegm=0; iSegm<NCHAN; iSegm++) { for(int iSegm=0; iSegm<NCHAN; iSegm++) {
float * tData = PF.tAmplitude[iSegm] + SAMP0/2 - nsteps + ll; float * tData = PF.tAmplitude[iSegm] + SAMP0/2 + ll-nsteps;
float * rData = PF.rAmplitude[iSegm] + SAMP0/2; float * rData = PF.rAmplitude[iSegm] + SAMP0/2;
for(int ii = 0; ii < nsamps; ii++) { // ii = nsteps for(int ii = 0; ii < nsamps; ii++) { // ii = nsteps
float diff = tData[ii] - rData[ii]; float diff = tData[ii] - rData[ii];
...@@ -2682,22 +2760,20 @@ Float_t PSAFilterGridSearch::FitT0AfterPSA(pointFull &PF) ...@@ -2682,22 +2760,20 @@ Float_t PSAFilterGridSearch::FitT0AfterPSA(pointFull &PF)
} }
} // end loop over the segments } // end loop over the segments
x2[ll] = chi2; x2[ll] = chi2;
} if(ll==nsteps || chi2 < x2min) {
// find the minimum
int ssmin = nsteps;
float x2min = x2[nsteps];
for(int ll = 0; ll <= 2*nsteps; ll++) {
if(x2[ll] < x2min) {
ssmin = ll; ssmin = ll;
x2min = x2[ll]; x2min = chi2;
} }
if(ssmin==0 || ssmin==2*nsteps)
return float(ssmin-nsteps); // if at the edge
if(x2[ssmin+1] < 0)
ll = ssmin+1; // calculate ssmin+1
else if (x2[ssmin-1] < 0)
ll = ssmin-1; // calculate ssmin-1
else
break; // both sides of ssmin already available
} }
// if at the edge
if(ssmin == 0 || ssmin == 2*nsteps)
return float(ssmin-nsteps);
// parabolic interpolation around the minimum // parabolic interpolation around the minimum
float ym = x2[ssmin-1]; float ym = x2[ssmin-1];
float ys = x2[ssmin ]; float ys = x2[ssmin ];
......
#ifndef PSAFILTERGRIDSEARCH_H_INCLUDED #ifndef PSAFILTERGRIDSEARCH_H_INCLUDED
#define PSAFILTERGRIDSEARCH_H_INCLUDED #define PSAFILTERGRIDSEARCH_H_INCLUDED
// Define this symbol to write to file a correction to the timestamp of the events
//#define WRITE_TSHIFT_CORRECTION
#ifdef WRITE_TSHIFT_CORRECTION
#include "dataPairFile.h"
#endif
#include "PSAFilter.h" #include "PSAFilter.h"
#include "GridSearchParams.h" #include "GridSearchParams.h"
#include "SignalBasis.h" #include "SignalBasis.h"
...@@ -103,10 +97,6 @@ private: ...@@ -103,10 +97,6 @@ private:
MultiHist<unsigned short> *PsaMatrSeg; MultiHist<unsigned short> *PsaMatrSeg;
#endif //PSA_MULTIHIST #endif //PSA_MULTIHIST
#ifdef WRITE_TSHIFT_CORRECTION
dataPairFile<ULong64_t, float> *tstmpFileO;
#endif
Int_t gWriteNumTraces; Int_t gWriteNumTraces;
Bool_t gWritePartialTraces; Bool_t gWritePartialTraces;
...@@ -129,9 +119,10 @@ private: ...@@ -129,9 +119,10 @@ private:
int SearchFullGrid (pointFull &PF, int netChSeg); int SearchFullGrid (pointFull &PF, int netChSeg);
int SearchAdaptive1(pointFull &PF, int netChSeg, bool bCoarseOnly = false, bool bFirstOnly = false); // 1 hit int SearchAdaptive1(pointFull &PF, int netChSeg, bool bCoarseOnly = false, bool bFirstOnly = false); // 1 hit
int SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoarseOnly = false, bool bFirstOnly = false); // 2 hits int SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoarseOnly = false, bool bFirstOnly = false); // 2 hits
void MakeSearchWave (pointFull &PF, float fact, float maxVal);
void MakeLocalMask (pointFull &PF); void MakeLocalMask (pointFull &PF);
float CalcChi2GridPts(int segNum, int pdNum, pointFull &PF);
float CalcChi2Residue(pointFull &PF); float CalcChi2Residue(pointFull &PF);
float CalcChi2GridPts(int segNum, int pdNum, pointFull &PF);
bool WriteBaseAverPt(float scale, std::string ext); bool WriteBaseAverPt(float scale, std::string ext);
bool WriteBasePoint (float scale, pointPsa &pt, FILE *FP); bool WriteBasePoint (float scale, pointPsa &pt, FILE *FP);
......
...@@ -333,7 +333,6 @@ public: ...@@ -333,7 +333,6 @@ public:
} }
} }
// possible overflow problems with GS_TYPE_INT16
void add(pointPsa &pt) { void add(pointPsa &pt) {
x += pt.x; x += pt.x;
y += pt.y; y += pt.y;
...@@ -388,9 +387,9 @@ public: ...@@ -388,9 +387,9 @@ public:
} }
}; };
// augmented by members used by the Grid-search algorithm // Augmented by members used by the Grid-search algorithm
// describ