Commit fe0f953d authored by dino's avatar dino
Browse files

Preliminary version of differential cross-talk correction and various modifications of SignalBasis

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@852 170316e4-aea8-4b27-aad4-0380ec0519c9
parent 59518583
......@@ -103,8 +103,9 @@ public:
Float_t CoreT[ADF::CrystalInterface::kNbCores];
// output data (with this algo there cannot be more hits than segments)
ADF::PSAHit theHits[ADF::CrystalInterface::kNbSegments];
UInt_t numHits;
Int_t retValue; // used only by the threaded part
UInt_t numHits;
Int_t retValue; // used only by the threaded part
bool selectIt; // used to select events
} PsaData;
PsaData DD[TCOUNT*TMODULO];
......
......@@ -14,47 +14,7 @@ using namespace std;
// of the individual segments but should be excluded in real PSA.
//
char hmask[NCHAN][NCHAN+1] = {
// 0 1 2 3
// 0123456789012345678901234567890123456
// A B C D E F I
"1200002200002000002000002000002200009", // A1
"2120002220000000000000000000002220009", // A2
"0212000222000000000000000000000222009", // A3
"0021200022200000000000000000000022209", // A4
"0002120002220000000000000000000002229", // A5
"0000210000220000000000000000000000229", // A6
"2200001200002200002000002000002000009", // B1
"2220002120002220000000000000000000009", // B2
"0222000212000222000000000000000000009", // B3
"0022200021200022200000000000000000009", // B4
"0002220002120002220000000000000000009", // B5
"0000220000210000220000000000000000009", // B6
"2000002200001200002200002000002000009", // C1
"0000002220002120002220000000000000009", // C2
"0000000222000212000222000000000000009", // C3
"0000000022200021200022200000000000009", // C4
"0000000002220002120002220000000000009", // C5
"0000000000220000210000220000000000009", // C6
"2000002000002200001200002200002000009", // D1
"0000000000002220002120002220000000009", // D2
"0000000000000222000212000222000000009", // D3
"0000000000000022200021200022200000009", // D4
"0000000000000002220002120002220000009", // D5
"0000000000000000220000210000220000009", // D6
"2000002000002000002200001200002200009", // E1
"0000000000000000002220002120002220009", // E2
"0000000000000000000222000212000222009", // E3
"0000000000000000000022200021200022209", // E4
"0000000000000000000002220002120002229", // E5
"0000000000000000000000220000210000229", // E6
"2200002000002000002000002200001200009", // F1
"2220000000000000000000002220002120009", // F2
"0222000000000000000000000222000212009", // F3
"0022200000000000000000000022200021209", // F4
"0002220000000000000000000002220002129", // F5
"0000220000000000000000000000220000219" // F6
};
char hmask[NCHAN][NCHAN+1] = {0};
const float eScale = 5.f;
const float dScale = 10.f;
......@@ -63,6 +23,8 @@ const int specLenT = 2*1000;
const int matLen = 100;
const int matOff = 50;
const bool doSelect = false; // the selection is programmed manually in Process()
#ifdef USENETCHARGES
const int diffLag = 0; // use the net-charge signals as they are
#else
......@@ -138,15 +100,17 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
cout << endl;
MakeSegmentMap();
int rval = 0;
if(fPsaSegCenter) {
rval = fBasis.ReadSegmentCenters(fBasisFile);
rval = fBasis.ReadSegmentCenters(fBasisFile, hmask);
if(rval != 0)
return rval;
// not much more to do
}
else {
rval = fBasis.ReadBasis(fBasisFile);
rval = fBasis.ReadBasis(fBasisFile, hmask);
if(rval != 0)
return rval;
#ifdef USEADAPTIVE
......@@ -155,8 +119,8 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
return rval;
#endif
// fBasis.SmoothBasisSignals(SMOOTH);
fBasis.ConvBasisSignals();
//fBasis.SmoothBasisSignals(SMOOTH); // moving average, superseeded by the
fBasis.ExponentialResponse(float(SMOOTH)); // exponential convolution or the experimental one
if(diffLag > 0) {
// Modify the signal basis because the grid-search uses the differentiated net-charge signals
......@@ -165,6 +129,10 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
fBasis.ReorderSamples();
// the convolution derived from the pulser must be executed after having the samples at 10 ns.
// It will not do good for the second part, but, as this is not yet used, id doesn't matter.
//fBasis.ConvBasisSignals();
// prepare the metrics vectors for the chi2 loops
float largest = 0;
for(int i=0; i<NMETRIC; i++) {
......@@ -237,12 +205,54 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
return 0;
}
void PSAFilterGridSearch::MakeSegmentMap()
{
const int neighbours = 2;
char hh[36][37];
for(int i1 = 0; i1 < 36; i1++) {
memset(hh[i1], '0', 36);
hh[i1][36] = 0;
}
for(int i1 = 0; i1 < 36; i1++) {
int sector1 = i1/6;
int slice1 = i1%6;
for(int nsli = -neighbours; nsli <= neighbours; nsli++) {
int slice2 = slice1 + nsli;
if(slice2 < 0 || slice2 > 5) continue;
int mxsec = neighbours-abs(nsli);
for(int nsec = -mxsec; nsec <= mxsec; nsec++) {
int sector2 = (sector1 + nsec + 6)%6;
int i2 = sector2*6 + slice2;
hh[i1][i2] = (i2==i1) ? '1' : '2';
}
// add all front segments as neighbours
}
if(slice1==0) {
for(int nsec = 1; nsec < 6; nsec++) {
int sector2 = (sector1 + nsec + 6)%6;
int i2 = sector2*6;
hh[i1][i2] = '2';
}
}
hh[i1][36] = 0;
}
for(int i1 = 0; i1 < 36; i1++) {
for(int i2 = 0; i2 < 36; i2++)
hmask[i1][i2] = hh[i1][i2];
hmask[i1][36] = '9';
hmask[i1][37] = 0;
}
}
Int_t PSAFilterGridSearch::Process(int slot)
{
PsaData *pD = DD + slot;
pointExp pS;
pD->numHits = 0;
pD->numHits = 0;
pD->selectIt = doSelect ? false : true;
// not a good event
if(pD->crystal_status)
......@@ -304,7 +314,7 @@ Int_t PSAFilterGridSearch::Process(int slot)
pS.Normalize(fact, localMask);
if(diffLag > 0) {
// apply a delayed differentiation of diffLag samples
// apply a delayed differentiation of diffLag samples to the net charge segment
short *ptdata = pS.sAmplitude[netChSeg];
for(int nn = NTIME-1; nn >= diffLag; nn--) {
ptdata[nn] -= ptdata[nn - diffLag];
......@@ -347,6 +357,12 @@ Int_t PSAFilterGridSearch::Process(int slot)
pD->theHits[pD->numHits].SetStatus((UShort_t)pS.bestPt);
pD->numHits++;
// manual selection of events for waves and histograms
if(doSelect) {
if((netChSeg%6)==0 && fx>0 && fx<6 && fy>-5 && fy<5)
pD->selectIt = true;
}
// traces
// The trace of the best point is normalized to the experimental amplitude and accumulated
if(fWriteNumTraces > 0) {
......@@ -382,13 +398,13 @@ Int_t PSAFilterGridSearch::Process(int slot)
sdata = slot_fAmplitude[slot];
for(int segm = 0; segm < NCHAN; segm++) {
float *fdata = pS.fAmplitude[segm];
for(int samp = 0; samp < NTIME; samp++)
for(int samp = 0; samp < WSAMP; samp++)
*sdata++ = (short)(eScale*fdata[samp]);
}
sdata = slot_rAmplitude[slot];
for(int segm = 0; segm < NCHAN; segm++) {
float *fdata = pS.rAmplitude[segm];
for(int samp = 0; samp < NTIME; samp++)
for(int samp = 0; samp < WSAMP; samp++)
*sdata++ = (short)( fdata[samp]); // escale already considered in the production
}
if(diffLag > 0) {
......@@ -396,8 +412,8 @@ Int_t PSAFilterGridSearch::Process(int slot)
// reapply the delayed differentiation to the experimental net-charge segments
for(int snum = 0; snum < sMult; snum++) {
int netChSeg = pS.netChargeSegnum[snum];
short *ptd = slot_fAmplitude[slot] + NTIME*netChSeg;
for(int nn = NTIME-1; nn >= diffLag; nn--) {
short *ptd = slot_fAmplitude[slot] + WSAMP*netChSeg;
for(int nn = WSAMP-1; nn >= diffLag; nn--) {
ptd[nn] -= ptd[nn - diffLag];
}
}
......@@ -444,20 +460,15 @@ void PSAFilterGridSearch::PrepareEvent(PsaData *pD, pointExp *pS)
for(Int_t atseg=0; atseg<NSEGS; atseg++) {
if(pS->masksignals[atseg] != '0') {
memcpy(pS->fAmplitude[atseg], pD->fTracesSG[atseg], ntocopy*sizeof(float));
// not good if done here --> should be done properly in the preprocessing
// in conjunction with the xtalk correction applied later, it increases too much the
// amplitude of the signals (positive and negative)
//if(DECONV)
//// not good if done here --> should be done properly in the preprocessing
//// in conjunction with the xtalk correction applied later, it increases too much the
//// amplitude of the signals (positive and negative)
// pS->deconv(DECONV, ntocopy, atseg); // and should deconvolve fAmplitude not amplitude
}
}
// copy also the CC
memcpy(pS->fAmplitude[NSEGS], pD->fTracesCC[0], ntocopy*sizeof(float));
if(XTCORR) {
// also this should be done properly in the preprocessing
pS->CoreXtalk(ntocopy);
}
}
void PSAFilterGridSearch::SetToSegCenter(PsaData *pD, pointExp *pS)
......@@ -485,6 +496,10 @@ int PSAFilterGridSearch::WriteTraces(int slot)
if(!fpPsaTraces)
return 0;
PsaData *pD = DD + slot;
if(!pD->selectIt)
return 0;
short wbuf[WSAMP];
int rstat = 0;
......@@ -528,7 +543,6 @@ int PSAFilterGridSearch::WriteTraces(int slot)
if(chanToCopy < WCHAN) {
// next come the hits
memset(wbuf, 0, sizeof(wbuf));
PsaData *pD = DD + slot;
int numHits = pD->numHits;
for(int nh = 0; nh < numHits; nh++) {
int pos = nh*10;
......@@ -709,6 +723,8 @@ int PSAFilterGridSearch::SearchAdaptive(pointExp *pS, int netChSeg, char *lMask,
Int_t PSAFilterGridSearch::PostProcess(int slot)
{
PsaData *pD = DD + slot;
if(!pD->selectIt)
goto skipSpectra;
int numHits = pD->numHits;
if(numHits < 1)
......@@ -779,12 +795,14 @@ Int_t PSAFilterGridSearch::PostProcess(int slot)
//
//if(specSigma && maxItime>0) {
// //overall chi2 (needs room in the spectrum)
// float chi2 = pS.ChiSquare(maxItime);
// float chi2 = pS.ChiSquare(maxItime, hmask[]); // ???
// int isigma = (int)sqrt(chi2);
// specSigma->incr(0, isigma);
// specSigma->incr(sMult, isigma);
//}
skipSpectra:
cServer.Exec(pD->timestamp);
#ifdef WRITE_PSA_HITS
......
......@@ -34,7 +34,7 @@
*/
const int WCHAN = 42; // number of channels to write in the output waves
const int WSAMP = 100; // number of samples per channel
const int WSAMP = 60; // number of samples per channel
class PSAFilterGridSearch : public PSAFilter
{
......@@ -66,6 +66,7 @@ private:
TH3F *PSAxyz;
#endif // _FromGRU_
void MakeSegmentMap();
void PrepareEvent (PsaData *pD, pointExp *pS);
void SetToSegCenter(PsaData *pD, pointExp *pS);
int SearchFullGrid(pointExp *pS, int netChSeg, char *lMask, int addr_first, int addr_last);
......@@ -89,8 +90,8 @@ public:
private:
// this is needed to have the possibility to write the traces in a thread-safe mode
short slot_fAmplitude[TCOUNT*TMODULO][(NCHAN+1)*NTIME]; // locally saved experimental trace
short slot_rAmplitude[TCOUNT*TMODULO][(NCHAN+1)*NTIME]; // locally saved "fitted" trace
short slot_fAmplitude[TCOUNT*TMODULO][(NCHAN+1)*WSAMP]; // locally saved experimental trace
short slot_rAmplitude[TCOUNT*TMODULO][(NCHAN+1)*WSAMP]; // locally saved "fitted" trace
};
......
......@@ -34,7 +34,7 @@ public:
addr_first = TSHIFT;
addr_last = NTIME-1;
memset(amplitude, 0, sizeof(amplitude));
memset(masksignals, '0', NCHAN+1); // including the closing zero
memset(masksignals, '0', NCHAN); masksignals[NCHAN] = 0;
fullChargeCollected = false;
}
// moving average filter
......@@ -79,7 +79,7 @@ public:
}
// exp -> step
void deconv(float tau, int len, int sg) {
void deconvExpToStep(float tau, int len, int sg) {
float aa = (float)exp(-1./tau);
float vv[NTIME];
for(int nn = 1; nn < len; nn++) {
......@@ -92,6 +92,18 @@ public:
amplitude[sg][nn] = (short)ss;
}
}
// delta->exp
void convDeltaToExp(float tau) {
float aa = (float)exp(-1./tau);
float gg = 1.f-aa;
for(int sg = 0; sg < NCHAN; sg++) {
float vv = 0; // accumulator used as float to keep precision
for(int nn = 0; nn < NTIME; nn++) {
vv = gg*amplitude[sg][nn] + aa*vv;
amplitude[sg][nn] = short(vv);
}
}
}
float distance(pointPsa *pp) {
float dx = x - pp->x;
float dy = y - pp->y;
......@@ -149,18 +161,6 @@ public:
}
}
}
void CoreXtalk(int ntocopy) {
for(int ii = 0; ii < ntocopy; ii++)
xCoreFract[ii] = XTCORR*fAmplitude[NSEGS][ii];
for(int ns = 0; ns < NSEGS; ns++) {
if(masksignals[ns] != '0') {
float *fA = fAmplitude[ns];
for(int ii = 0; ii < ntocopy; ii++) {
fA[ii] += xCoreFract[ii];
}
}
}
}
#if 0
// to be verified
void SearchLimits(int netChSeg, float fract, float thresh, int &first, int &last) {
......
......@@ -21,11 +21,11 @@ const int TIME_STEP = 5; // (ns) time step of the signal basis ==>
const int NTIME = 100; // number of points of the basis signals (100*5ns = 500ns)
const int TSHIFT = 5; // the "T0" (samples) of the calculated traces
const int SMOOTH = 10; // with of the moving average used to smooth the basis (== 50 ns)
const float XTCORR = 0.002f; // fraction of CC to add to the segments
const int SMOOTH = 8; // with of the moving average and exponential-response used to smooth the basis (== 40 ns)
const float XTCORR = 0.002f; // fraction of CC to add to the segments (but PSA should not act on the experimental data!!)
const float DECONV = 4600.f; // to deconvolve the main pole
const int NSAMP = defTraceLengthPSA; // length of experimental data (taken as float data)
const int TIME0 = defTriggerSample; // pretrigger of the experimental traces
const int TIME0 = defTriggerSample; // pretrigger of the experimental traces
const int STEP_SAMP = 1; // step over experimental traces (in 10 ns units)
#endif // GRIDSEARCHPARAMS_H_INCLUDED
......@@ -7,7 +7,7 @@
using namespace std;
int SignalBasis::ReadBasis(std::string fname, bool keep)
int SignalBasis::ReadBasis(std::string fname, char hmask[NCHAN][NCHAN+1], bool keep)
{
if(fname.size() < 1) {
printf("Missing signal basis\n");
......@@ -36,7 +36,7 @@ int SignalBasis::ReadBasis(std::string fname, bool keep)
if(numpts[nn] > 0 && numpts[nn]< 100000)
ngood++;
}
bool bVentu = ngood == 36;
bool bVentu = (ngood == 36);
// find size of file (and rewind it)
fseek(ofp1, 0, SEEK_END);
......@@ -59,13 +59,13 @@ int SignalBasis::ReadBasis(std::string fname, bool keep)
if(bVentu && !bBartB) {
printf("Reading the basis of signals ...\n");
int readErr = ReadBasisFormatVentu(ofp1, keep);
int readErr = ReadBasisFormatVentu(ofp1, hmask, keep);
return readErr;
}
else if (!bVentu && bBartB) {
int numSignals = int(fSize/sizeof(bbhit_t));
printf("Reading Bart\'s basis containing %d signals ...\n", numSignals);
int readErr = ReadBasisFormatBartB(ofp1, keep, numSignals);
int readErr = ReadBasisFormatBartB(ofp1, hmask, keep, numSignals);
return readErr;
}
else {
......@@ -76,7 +76,7 @@ int SignalBasis::ReadBasis(std::string fname, bool keep)
return 0;
}
int SignalBasis::ReadBasisFormatVentu(FILE * ofp1, bool keep)
int SignalBasis::ReadBasisFormatVentu(FILE * ofp1, char hmask[NCHAN][NCHAN+1], bool keep)
{
bCenters = false;
bFullGrid = false;
......@@ -106,6 +106,8 @@ int SignalBasis::ReadBasisFormatVentu(FILE * ofp1, bool keep)
rv = fread(&(Pts[ii][jj].z), sizeof(float), 1, ofp1); ///
rv = fread( Pts[ii][jj].amplitude, sizeof(short), NCHAN*NTIME, ofp1); ///
Pts[ii][jj].masksignals[NCHAN] = 0; // to close the string
// replace the mask from the data file with the one provided by the caller
memcpy(Pts[ii][jj].masksignals, hmask[ii], NCHAN+1);
int ntseg = Pts[ii][jj].netChargeSegment;
if(ntseg != ii)
nerr++;
......@@ -125,50 +127,8 @@ int SignalBasis::ReadBasisFormatVentu(FILE * ofp1, bool keep)
return 0;
}
int SignalBasis::ReadBasisFormatBartB(FILE * ofp1, bool keep, int numSignals)
int SignalBasis::ReadBasisFormatBartB(FILE * ofp1, char hmask[NCHAN][NCHAN+1], bool keep, int numSignals)
{
char hmask[NCHAN][NCHAN+1] = {
// 0 1 2 3
// 0123456789012345678901234567890123456
// A B C D E F I
"1200002200002000002000002000002200009", // A1
"2120002220000000000000000000002220009", // A2
"0212000222000000000000000000000222009", // A3
"0021200022200000000000000000000022209", // A4
"0002120002220000000000000000000002229", // A5
"0000210000220000000000000000000000229", // A6
"2200001200002200002000002000002000009", // B1
"2220002120002220000000000000000000009", // B2
"0222000212000222000000000000000000009", // B3
"0022200021200022200000000000000000009", // B4
"0002220002120002220000000000000000009", // B5
"0000220000210000220000000000000000009", // B6
"2000002200001200002200002000002000009", // C1
"0000002220002120002220000000000000009", // C2
"0000000222000212000222000000000000009", // C3
"0000000022200021200022200000000000009", // C4
"0000000002220002120002220000000000009", // C5
"0000000000220000210000220000000000009", // C6
"2000002000002200001200002200002000009", // D1
"0000000000002220002120002220000000009", // D2
"0000000000000222000212000222000000009", // D3
"0000000000000022200021200022200000009", // D4
"0000000000000002220002120002220000009", // D5
"0000000000000000220000210000220000009", // D6
"2000002000002000002200001200002200009", // E1
"0000000000000000002220002120002220009", // E2
"0000000000000000000222000212000222009", // E3
"0000000000000000000022200021200022209", // E4
"0000000000000000000002220002120002229", // E5
"0000000000000000000000220000210000229", // E6
"2200002000002000002000002200001200009", // F1
"2220000000000000000000002220002120009", // F2
"0222000000000000000000000222000212009", // F3
"0022200000000000000000000022200021209", // F4
"0002220000000000000000000002220002129", // F5
"0000220000000000000000000000220000219" // F6
};
bCenters = false;
bFullGrid = false;
bCoarseFine = false;
......@@ -283,7 +243,7 @@ char hmask[NCHAN][NCHAN+1] = {
Pts[ii][jj].netChargeSegment = ii;
Pts[ii][jj].x = bbhit.Pos[0];
Pts[ii][jj].y = bbhit.Pos[1];
Pts[ii][jj].z = bbhit.Pos[2];
Pts[ii][jj].z = bbhit.Pos[2]-1.f; // offset due to SimIon
float * ptf = bbhit.Tr + 10 - TSHIFT; // it seems that bart's basis starts at 10
for(nn = 0; nn < NTIME; nn++)
Pts[ii][jj].amplitude[NSEGS][nn] = short(ptf[nn]*MAXNORM);
......@@ -295,7 +255,6 @@ char hmask[NCHAN][NCHAN+1] = {
memcpy(Pts[ii][jj].masksignals, hmask[ii], NCHAN+1);
Pts[ii][jj].addr_first = TSHIFT;
Pts[ii][jj].addr_last = NTIME-1;
//Pts[ii][jj].smooth(SMOOTH); // this should not be done here
numPts[ii]++;
......@@ -310,7 +269,7 @@ char hmask[NCHAN][NCHAN+1] = {
return 0;
}
int SignalBasis::ReadSegmentCenters(std::string fname)
int SignalBasis::ReadSegmentCenters(std::string fname, char hmask[NCHAN][NCHAN+1])
{
bCenters = false;
bFullGrid = false;
......@@ -328,7 +287,7 @@ int SignalBasis::ReadSegmentCenters(std::string fname)
if( (ofp1 = fopen(fSegCenters.c_str(), "r")) == NULL) {
printf("File %s not found\n", fSegCenters.c_str());
printf("Try reading the full basis ...\n");
return ReadBasis(theFileName, false);
return ReadBasis(theFileName, hmask, false);
}
printf("Reading the center of segments ...\n");
......@@ -846,6 +805,19 @@ void SignalBasis::SmoothBasisSignals(int nsmoo)
}
}
// filters all basis signals with an exponential response of tau (samples of 5 ns)
void SignalBasis::ExponentialResponse(float tau)
{
int count = 0;
for(int ii=0; ii<NSEGS; ii++) {
int npoints = numPts[ii];
for(int jj=0; jj<npoints; jj++) {
Pts[ii][jj].convDeltaToExp(tau);
count++;
}
}
}
// mw convolution on all signals
void SignalBasis::ConvBasisSignals()
{
......
......@@ -75,10 +75,10 @@ typedef struct {
}
}
}
int ReadBasis(std::string fname, bool keep = true);
int ReadBasisFormatVentu(FILE * fp, bool keep);
int ReadBasisFormatBartB(FILE * fp, bool keep, int numSignals);
int ReadSegmentCenters(std::string fname);
int ReadBasis(std::string fname, char hmask[NCHAN][NCHAN+1], bool keep = true);
int ReadBasisFormatVentu(FILE * fp, char hmask[NCHAN][NCHAN+1], bool keep);
int ReadBasisFormatBartB(FILE * fp, char hmask[NCHAN][NCHAN+1], bool keep, int numSignals);
int ReadSegmentCenters(std::string fname, char hmask[NCHAN][NCHAN+1]);
int MakeCoarseFineList(bool verbose);
bool MakeOneSegment(int ns, bool verbose);
bool AssignIt(float px, float py, float pz, float cx, float cy, float cz);
......@@ -87,6 +87,7 @@ typedef struct {
int NearestCoarse(float px, float py, float pz, cgrid *cgr, int numC, float& dist);
void SmoothBasisSignals(int nsmoo);
void ConvBasisSignals();
void ExponentialResponse(float tau);
void CalculateSegmentCenters();
int DifferentiateNetChargeSignals(int lag);
void ReorderSamples();
......
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