Commit d3519f8c authored by dino's avatar dino
Browse files

Better treatment of "broken" segments in the PSA. These segments are...

Better treatment of "broken" segments in the PSA. These segments are recognized by the fact that the diagonal element of the cross-talk matrix is zero. There should probably be also a dedicated keyword in PSAFilter

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@946 170316e4-aea8-4b27-aad4-0380ec0519c9
parent 8009a8fd
......@@ -108,13 +108,15 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
cout << endl;
MakeSegmentMap(2);
SegmentMapMake(2);
if(fPsaXtalkFile.length() > 0) {
string fname = fConfPath + fPsaXtalkFile;
bool readOK = fBasis.ReadXtalkCoeffs(fname, fVerbose);
if(!readOK)
return 111;
SegmentMapTest(); // check if broken segments
}
if(fPsaTshiftFile.length() > 0) {
......@@ -267,7 +269,7 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
}
// Prepare the map of neighbours to be used in the PSA
void PSAFilterGridSearch::MakeSegmentMap(int neighbours)
void PSAFilterGridSearch::SegmentMapMake(int neighbours)
{
neighbours = abs(neighbours);
if(neighbours > 5) neighbours = 5;
......@@ -318,6 +320,18 @@ void PSAFilterGridSearch::MakeSegmentMap(int neighbours)
}
}
// exclude elements corresponding to broken segments
void PSAFilterGridSearch::SegmentMapTest() {
for(int ns = 0; ns < 36; ns++) {
if(fBasis.xTalkProp[ns][ns] == 0) {
for(int nn = 0; nn < 36; nn++) {
hmask[ns][nn] = '0';
hmask[nn][ns] = '0';
}
}
}
}
// modified to process nslots events for each call
// this complication is to call the GPU implementation with locEvts (full) events
......@@ -571,6 +585,8 @@ void PSAFilterGridSearch::PrepareEvent(PsaData *pD, pointFull *pS)
//find hit segments
pS->sumSegm = 0;
Int_t numsegs = 0;
for(Int_t isg=0; isg<NCHAN; isg++) pS->masksignals[isg] = '0';
pS->masksignals[NCHAN] = 0;
for(Int_t isg=0; isg<NSEGS; isg++) {
if(pD->SegE[isg] > fHitSegThreshold) {
pS->netChargeSegnum[numsegs] = isg;
......@@ -600,12 +616,18 @@ void PSAFilterGridSearch::PrepareEvent(PsaData *pD, pointFull *pS)
for(int ii = 0; ii < ntocopy; ii++)
pS->fAmplitude[atseg][ii] = (gs_type)pD->fTracesSG[atseg][ii];
}
for(int ii = 0; ii < ntocopy; ii++)
pS->tAmplitude[atseg][ii] = (gs_type)pD->fTracesSG[atseg][ii];
}
// copy also the CC
for(int ii = 0; ii < ntocopy; ii++)
for(int ii = 0; ii < ntocopy; ii++) {
pS->fAmplitude[NSEGS][ii] = (gs_type)pD->fTracesCC[0][ii];
}
for(int ii = 0; ii < ntocopy; ii++) {
pS->tAmplitude[NSEGS][ii] = (gs_type)pD->fTracesCC[0][ii];
}
// and save the whole trace into oAmplitude for any possible use
// save fAmplitude into oAmplitude for any possible use
memcpy(pS->oAmplitude, pS->fAmplitude, sizeof(pS->fAmplitude));
}
......@@ -698,6 +720,8 @@ void PSAFilterGridSearch::StorePartialTrace(pointFull &pS, pointPsa *bestPoint,
int time_first = TSHIFT-PRETRIG;
float fact1 = 1.f/scaleFact; // to get the right vertical gain
for(int segm = 0; segm < NCHAN; segm++) {
if(fBasis.segBroken[segm])
continue;
realTrace = pS.rAmplitude[segm];
baseTrace = bestPoint->amplitude[segm];
#ifndef SHOWDIFFERENTIATED
......@@ -763,7 +787,11 @@ void PSAFilterGridSearch::SaveTotalTrace(pointFull &pS, PsaData *pD)
gs_type *sdata;
sdata = pD->Get_fAmplitude(gs_type(0));
for(int segm = 0; segm < NCHAN; segm++) {
#ifdef SHOWFULLINPUTTRACE
gs_type *fdata = pS.tAmplitude[segm];
#else
gs_type *fdata = pS.oAmplitude[segm];
#endif
for(int samp = 0; samp < wnsamp; samp++)
*sdata++ = gs_type(fdata[samp]*sScale);
}
......
......@@ -90,7 +90,8 @@ private:
TH3F *PSAxyz;
#endif // PSA_FromGRU_
void MakeSegmentMap (int neighbours);
void SegmentMapMake (int neighbours); // based on Manhattan distance
void SegmentMapTest (); // exclude elements related to broken segments
void PrepareEvent (PsaData *pD, pointFull *pS);
void SetToSegCenter (PsaData *pD, pointFull *pS);
int ProcessEvent (pointFull *n_pS, int num);
......
......@@ -251,8 +251,9 @@ public:
int numNetCharges;
int netChargeSegnum[NCHAN];
float netChargeEnergy[NCHAN];
gs_type oAmplitude[NCHAN][NTIME]; // the original experimental data
gs_type fAmplitude[NCHAN][NTIME]; // the experimental data used by the GS
gs_type tAmplitude[NCHAN][NTIME]; // the total original experimental data
gs_type fAmplitude[NCHAN][NTIME]; // the part of experimental data used by the GS
gs_type oAmplitude[NCHAN][NTIME]; // a copy of fAmplitude to simplify analysis
gs_type sAmplitude[NCHAN][NTIME]; // used by the search loop
gs_type rAmplitude[NCHAN][NTIME]; // "fitted" trace
int bestPt;
......@@ -275,8 +276,9 @@ private:
numNetCharges = 0;
memset(netChargeSegnum, 0, sizeof(netChargeSegnum));
memset(netChargeEnergy, 0, sizeof(netChargeEnergy));
memset(oAmplitude, 0, sizeof(oAmplitude));
memset(tAmplitude, 0, sizeof(tAmplitude));
memset(fAmplitude, 0, sizeof(fAmplitude));
memset(oAmplitude, 0, sizeof(oAmplitude));
memset(sAmplitude, 0, sizeof(sAmplitude));
memset(rAmplitude, 0, sizeof(rAmplitude));
memset(segOrder, 0, sizeof(segOrder));
......
......@@ -26,6 +26,8 @@ typedef short gs_type;
//# define SHOWDIFFERENTIATED // in the saved traces show the differentiated version of the net-charge segments
#endif
#define SHOWFULLINPUTTRACE // comment this to see only the part used by the grid search
#define ORDEREDSEARCH // energy-ordered decomposition with removal of previous signals
#ifdef ORDEREDSEARCH
//# define PRECOARSESEARCH // a coarse-only preliminary search to remove the other net-charges
......
......@@ -1137,6 +1137,10 @@ bool SignalBasis::ReadXtalkCoeffs(string cname, bool verbose)
}
cout << "Found all " << 36*36 << " expected values" << endl;
// should also be possible to "break" a segment from the command line
for(int nn = 0; nn < 36; nn++)
segBroken[nn] = xTalkProp[nn][nn] ? false : true;
// average of non-neighbours and not in the same sector segments
for(int n1 = 0; n1 < 36; n1++) {
float vv = 0;
......@@ -1150,6 +1154,8 @@ bool SignalBasis::ReadXtalkCoeffs(string cname, bool verbose)
if(sector2 == sector1-1 || sector2 == sector1+1 ) continue;
if(sector2+6 == sector1-1 || sector2+6 == sector1+1 ) continue;
if(sector2 == sector1+6-1 || sector2 == sector1+6+1 ) continue;
if(segBroken[n2])
continue; // exclude broken segments
vv += xTalkProp[n1][n2];
nn++;
}
......@@ -1162,12 +1168,13 @@ bool SignalBasis::ReadXtalkCoeffs(string cname, bool verbose)
// gain factor for the correction has to be adjusted
for(int n1 = 0; n1 < 36; n1++) {
for(int n2 = 0; n2 < 36; n2++) {
if(n2==n1) {
if(n2==n1 || segBroken[n2]==0) {
xTalkDiff[n1][n2] = 0;
}
else {
float xdiff = xTalkAver[n1]-xTalkProp[n1][n2];
if(xdiff < 0) xdiff = 0;
if(xdiff < 0)
xdiff = 0;
xTalkDiff[n1][n2] = xdiff*xFactor;
}
}
......
......@@ -50,6 +50,7 @@ public:
float xTalkProp[NSEGS][NSEGS];
float xTalkDiff[NSEGS][NSEGS];
float xTalkAver[NSEGS];
bool segBroken[NSEGS];
bool tShiftsCorr;
float tShifts[NCHAN+1];
......@@ -89,6 +90,7 @@ public:
memset(numPts, 0, sizeof(numPts));
memset(Pts, 0, sizeof(Pts));
memset(cflist, 0, sizeof(cflist));
memset(segBroken, false, sizeof(segBroken));
}
virtual ~SignalBasis() {
for(int ns=0; ns<NSEGS; ns++) {
......
......@@ -842,6 +842,9 @@ bool PreprocessingFilterPSA::readXtalkCoeffs(string cname)
}
cout << "Found all " << 36*36 << " expected values" << endl;
// SINCE THE DIFFERENTIAL CROSS TALK HAS BEEN MOVED TO THE PSA
// xTalkAver AND xTalkDiff ARE UNUSED AND SHOULD BE REMOVED
// average of non-neighbours and not in the same sector segments
for(int n1 = 0; n1 < 36; n1++) {
float vv = 0;
......@@ -855,6 +858,8 @@ bool PreprocessingFilterPSA::readXtalkCoeffs(string cname)
if(sector2==sector1-1 || sector2==sector1+1) continue;
if(sector2+6==sector1-1 || sector2+6==sector1+1) continue;
if(sector2==sector1+6-1 || sector2==sector1+6+1) continue;
if(!CC.xTalkProp[n2][n2])
continue; // exclude broken segments
vv += CC.xTalkProp[n1][n2];
nn++;
}
......@@ -867,7 +872,7 @@ bool PreprocessingFilterPSA::readXtalkCoeffs(string cname)
// gain factor for the correction has to be adjusted
for(int n1 = 0; n1 < 36; n1++) {
for(int n2 = 0; n2 < 36; n2++) {
if(n2==n1)
if(n2==n1 || CC.xTalkProp[n2][n2]==0)
CC.xTalkDiff[n1][n2] = 0;
else
CC.xTalkDiff[n1][n2] = (CC.xTalkAver[n1]-CC.xTalkProp[n1][n2])*xFactor;
......
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