Commit 56889753 authored by dino's avatar dino
Browse files

Multi-segment events are now analyzed in order of decreasing segment energy....

Multi-segment events are now analyzed in order of decreasing segment energy. The the result of the already-analyzed segments is subtracted from the experimental trace so as to remove/reduce their contribution on the lower energy ones.

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@875 170316e4-aea8-4b27-aad4-0380ec0519c9
parent 52c717e1
......@@ -46,24 +46,13 @@ PSAFilter::PSAFilter() :
fInblock.SetModeIO(ConfAgent::kRead);
fOublock.SetModeIO(ConfAgent::kWrite);
for(int slot = 0; slot < TCOUNT*TMODULO; slot++) {
DD[slot].numHits = 0;
for(UInt_t nh = 0; nh < ADF::CrystalInterface::kNbSegments; nh++) {
DD[slot].theHits[nh].Reset();
}
}
}
PSAFilter::~PSAFilter()
{
for(int slot = 0; slot < TCOUNT*TMODULO; slot++) {
for(int nn=0; nn<CrystalInterface::kNbSegments; nn++){
delete [] DD[slot].fTracesSG[nn];
}
for(int nn = 0; nn < CrystalInterface::kNbCores; nn++) {
delete [] DD[slot].fTracesCC[nn];
}
DD[slot].Destroy();
}
}
......@@ -158,21 +147,7 @@ void PSAFilter::process_initialise (UInt_t *error_code)
fFrameIO.SetStatus(BaseFrameIO::kIdle);
for(int slot = 0; slot < TCOUNT*TMODULO; slot++) {
for(int nn = 0; nn < CrystalInterface::kNbSegments; nn++) {
DD[slot].fTracesSG[nn] = new Float_t[fTraceLengthPSA];
memset(DD[slot].fTracesSG[nn], 0, sizeof(Float_t)*fTraceLengthPSA);
DD[slot].SegE[nn] = 0;
}
for(int nn = 0; nn < CrystalInterface::kNbCores; nn++) {
DD[slot].fTracesCC[nn] = new Float_t[fTraceLengthPSA];
memset(DD[slot].fTracesCC[nn], 0, sizeof(Float_t)*fTraceLengthPSA);
DD[slot].CoreE[nn] = DD[slot].CoreT[nn] = 0;
}
DD[slot].numHits = 0;
for(UInt_t nh = 0; nh < ADF::CrystalInterface::kNbSegments; nh++) {
DD[slot].theHits[nh].Reset();
}
DD[slot].InitInput(fTraceLengthPSA);
}
// version-specific initializations
......@@ -319,8 +294,8 @@ Int_t PSAFilter::SetOutput(int slot)
crystal_status = pD->crystal_status;
CoreE0 = pD->CoreE[0];
CoreE1 = pD->CoreE[1];
CoreT0 = pD->CoreT[0];
CoreT1 = pD->CoreT[1];
CoreT0 = pD->CoreT[0] + pD->shiftT0;
CoreT1 = pD->CoreT[1] + pD->shiftT0;
//myAno.Set(crystal_id, 0);
//myAno.fRealSize = sizeof(crystal_id);
......
......@@ -60,6 +60,106 @@ struct CSlot
};
#endif // TCOUNT > 1
// Data interface to PSAGridSearch.
// Packed in a structure/class to be able to run independent serches in parallel
struct PsaData
{
// event identification
UShort_t crystal_id;
UShort_t crystal_status;
UInt_t evnumber;
ULong64_t timestamp;
Float_t shiftT0;
// input data (length of traces managed by the owner of the object)
Float_t *fTracesSG[ADF::CrystalInterface::kNbSegments];
Float_t *fTracesCC[ADF::CrystalInterface::kNbCores];
Float_t SegE[ADF::CrystalInterface::kNbSegments];
Float_t CoreE[ADF::CrystalInterface::kNbCores];
Float_t CoreT[ADF::CrystalInterface::kNbCores];
// output data (with this algo there cannot be more hits than segments)
UInt_t numHits;
ADF::PSAHit theHits[ADF::CrystalInterface::kNbSegments];
// This is the the resulting "fitted" trace, used to save Psa__xxxx__Traces.samp and to fit T0
// To enable using short or float without including GridSearchParams.h, keep both pointer types
short *sfAmplitude; // locally saved experimental trace (length managed by the owner of the object)
short *srAmplitude; // locally saved "fitted" trace (length managed by the owner of the object)
float *ffAmplitude; // locally saved experimental trace (length managed by the owner of the object)
float *frAmplitude; // locally saved "fitted" trace (length managed by the owner of the object)
Int_t retValue; // used only by the threaded part
bool selectIt; // used to select events
PsaData() : crystal_id(0), crystal_status(0), evnumber(0), timestamp(0), shiftT0(0) {
memset(fTracesSG, 0, sizeof(fTracesSG));
memset(fTracesCC, 0, sizeof(fTracesCC));
memset(SegE, 0, sizeof(SegE));
memset(CoreE, 0, sizeof(CoreE));
memset(CoreT, 0, sizeof(CoreT));
numHits = 0;
for(UInt_t nh = 0; nh < ADF::CrystalInterface::kNbSegments; nh++)
theHits[nh].Reset();
numHits = 0;
sfAmplitude = srAmplitude = NULL;
ffAmplitude = frAmplitude = NULL;
retValue = 0;
selectIt = false;
}
void InitInput(Int_t traceLengthPSA) {
for(int nn = 0; nn < ADF::CrystalInterface::kNbSegments; nn++) {
fTracesSG[nn] = new Float_t[traceLengthPSA];
memset(fTracesSG[nn], 0, sizeof(Float_t)*traceLengthPSA);
SegE[nn] = 0;
}
for(int nn = 0; nn < ADF::CrystalInterface::kNbCores; nn++) {
fTracesCC[nn] = new Float_t[traceLengthPSA];
memset(fTracesCC[nn], 0, sizeof(Float_t)*traceLengthPSA);
CoreE[nn] = CoreT[nn] = 0;
}
}
void InitOutput(Int_t totLengthOut, float) {
ffAmplitude = new float[totLengthOut];
memset(ffAmplitude, 0, sizeof(float)*totLengthOut);
frAmplitude = new float[totLengthOut];
memset(frAmplitude, 0, sizeof(float)*totLengthOut);
}
void InitOutput(Int_t totLengthOut, short) {
sfAmplitude = new short[totLengthOut];
memset(sfAmplitude, 0, sizeof(short)*totLengthOut);
srAmplitude = new short[totLengthOut];
memset(srAmplitude, 0, sizeof(short)*totLengthOut);
}
void Destroy() {
for(int nn=0; nn<ADF::CrystalInterface::kNbSegments; nn++) {
if(fTracesSG)
delete [] fTracesSG[nn];
fTracesSG[nn] = NULL;
}
for(int nn = 0; nn < ADF::CrystalInterface::kNbCores; nn++) {
if(fTracesCC)
delete [] fTracesCC[nn];
fTracesCC[nn] = NULL;
}
if(ffAmplitude)
delete [] ffAmplitude;
ffAmplitude = NULL;
if(frAmplitude)
delete [] frAmplitude;
frAmplitude = NULL;
if(sfAmplitude)
delete [] sfAmplitude;
sfAmplitude = NULL;
if(srAmplitude)
delete [] srAmplitude;
srAmplitude = NULL;
}
float *Get_fAmplitude(float) {return ffAmplitude;}
float *Get_rAmplitude(float) {return frAmplitude;}
short *Get_fAmplitude(short) {return sfAmplitude;}
short *Get_rAmplitude(short) {return srAmplitude;}
};
//! Base class for PSA filters
class PSAFilter : public ADF::NarvalFilter
{
......@@ -89,27 +189,8 @@ protected:
//ADF::Anonymous myAno;
public:
// The data-interface to the GridSearch algorithm
// Packed in a struct so as to be able to run independent serches in parallel
typedef struct {
// event identification
UShort_t crystal_id;
UShort_t crystal_status;
UInt_t evnumber;
ULong64_t timestamp;
// input data
Float_t *fTracesSG[ADF::CrystalInterface::kNbSegments];
Float_t *fTracesCC[ADF::CrystalInterface::kNbCores];
Float_t SegE[ADF::CrystalInterface::kNbSegments];
Float_t CoreE[ADF::CrystalInterface::kNbCores];
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
bool selectIt; // used to select events
} PsaData;
// The data-interface to the GridSearch algorithm
PsaData DD[TCOUNT*TMODULO];
protected:
......
This diff is collapsed.
......@@ -51,8 +51,14 @@ private:
nDhist<unsigned int> *specEner;
nDhist<unsigned int> *specTzero;
nDhist<unsigned int> *specSigma;
nDhist<unsigned int> *matrXYZR;
nDhist<unsigned int> *matrSeg;
nDhist<unsigned short> *matrXYZR;
nDhist<unsigned short> *matrSeg;
nDhist<unsigned short> *matrZE;
nDhist<unsigned short> *matrRE;
nDhist<unsigned short> *matrZRE;
nDhist<unsigned short> *matrRZE;
nDhist<unsigned short> *matrZREcEs;
nDhist<unsigned short> *matrRZEcEs;
#endif //LOCALSPECTRA
std::string fnPsaTraces;
......@@ -68,11 +74,13 @@ private:
TH3F *PSAxyz;
#endif // _FromGRU_
void MakeSegmentMap(int neighbours);
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);
int SearchAdaptive(pointExp *pS, int netChSeg, char *lMask, int addr_first, int addr_last);
void MakeSegmentMap (int neighbours);
void PrepareEvent (PsaData *pD, pointExp *pS);
void SetToSegCenter (PsaData *pD, pointExp *pS);
void MakeLocalMask (char *localMask, pointExp &pS, int netChSeg);
int PreSearchCoarse(pointExp &pS, int sMult);
int SearchFullGrid (pointExp *pS, int netChSeg, char *lMask, int addr_first, int addr_last);
int SearchAdaptive (pointExp *pS, int netChSeg, char *lMask, int addr_first, int addr_last, bool bCoarseOnly = false);
protected:
// this is written in a thread-safe mode and can be called in parallel, using different data slots
Int_t Process(int slot = 0);
......@@ -88,12 +96,11 @@ public:
void SetHitSegThreshold(Float_t thres) {fHitSegThreshold = thres;}
Float_t GetHitSegThreshold() {return fHitSegThreshold;}
int WriteTraces(int slot);
private:
// this is needed to have the possibility to write the traces in a thread-safe mode
gs_type slot_fAmplitude[TCOUNT*TMODULO][NCHAN*WSAMP]; // locally saved experimental trace
gs_type slot_rAmplitude[TCOUNT*TMODULO][NCHAN*WSAMP]; // locally saved "fitted" trace
int WriteTraces(PsaData *pD);
void StorePartialTrace(pointExp &pS, pointPsa *bestPoint, float scaleFact, int samp_first, int usamp);
int WritePartialTrace(pointExp &pS, pointPsa *bestPoint, float scaleFact, int samp_first, int usamp);
void SaveTotalTrace (pointExp &pS, int slot);
Float_t FitT0FromCore (pointExp &pS, int tsamp);
};
......
......@@ -87,9 +87,9 @@ public:
}
}
// convolution using a given kernel
void conv(double* kernel, int kernellen, int chan) {
void conv(double* kernel, int kernellen, int chan, int shift) {
gs_type temp[NTIME];
gs_type temp[2*NTIME];
for (int i = 0; i < NTIME; i++) {
double accu = 0; // accumulator register
int len = kernellen > i ? i : kernellen;
......@@ -99,13 +99,14 @@ public:
}
temp[i] = (gs_type)accu;
}
memcpy(amplitude[chan], temp, sizeof(temp));
for(int nn = NTIME; nn < 2* NTIME; nn++) temp[nn] = temp[NTIME-2];
memcpy(amplitude[chan], temp+shift, sizeof(gs_type)*NTIME);
}
void convsegs(double* kernel, int kernellen) {
for(int ns = 0; ns < NSEGS; ns++) conv(kernel, kernellen, ns);
void convsegs(double* kernel, int kernellen, int shift) {
for(int ns = 0; ns < NSEGS; ns++) conv(kernel, kernellen, ns, shift);
}
void convcore(double* kernel, int kernellen) {
conv(kernel, kernellen, NSEGS);
void convcore(double* kernel, int kernellen, int shift) {
conv(kernel, kernellen, NSEGS, shift);
}
// exp -> step
void deconvExpToStep(float tau, int len, int sg) {
......@@ -173,7 +174,7 @@ public:
int idGe;
float enerGe;
float sumSegm;
int netChargeNumber;
int numNetCharges;
int netChargeSegnum[NCHAN];
float netChargeEnergy[NCHAN];
gs_type oAmplitude[NCHAN][NTIME]; // the original experimental data
......@@ -184,6 +185,11 @@ public:
int bestdt;
float chi2min;
// used to perform search in energy-decreasing order
int segOrder[NSEGS]; // the net-charge segments
float eneOrder[NSEGS]; // their energy
int ptCoarse[NSEGS]; // the best coarse-grid point (if PRECOARSESEARCH)
pointExp() {
Init();
}
......@@ -192,13 +198,16 @@ private:
sumSegm = 0;
idGe = 0;
enerGe = 0;
netChargeNumber = 0;
numNetCharges = 0;
memset(netChargeSegnum, 0, sizeof(netChargeSegnum));
memset(netChargeEnergy, 0, sizeof(netChargeEnergy));
memset(oAmplitude, 0, sizeof(oAmplitude));
memset(fAmplitude, 0, sizeof(fAmplitude));
memset(sAmplitude, 0, sizeof(sAmplitude));
memset(rAmplitude, 0, sizeof(rAmplitude));
memset(segOrder, 0, sizeof(segOrder));
memset(eneOrder, 0, sizeof(eneOrder));
memset(ptCoarse, 0, sizeof(ptCoarse));
}
public:
void Reset() {
......@@ -206,12 +215,12 @@ public:
Init();
}
void MakeSearchWave(float fact, char *mask) {
float yMax = NMETRIC2/2 - 10;
float yMin = -NMETRIC2/2 + 10;
const float yMax = NMETRIC2/2 - 10; // set range limits
const float yMin = -NMETRIC2/2 + 10;
for(int iSegm = 0; iSegm < NSEGS; iSegm++) {
if(mask[iSegm] != '0') {
gs_type *sA = sAmplitude[iSegm];
gs_type *fA = fAmplitude[iSegm];
gs_type *sA = sAmplitude[iSegm]; // the wave to be composed
gs_type *fA = fAmplitude[iSegm]; // taking the data from this
for(int kk=0; kk<NTIME; kk++) {
float yy = fact*fA[kk];
if(yy > yMax)
......@@ -233,9 +242,10 @@ public:
}
void fRemoveFitted(float fact, char *mask, pointPsa *pPsa, int samp_first, int time_first) {
for(int iSegm = 0; iSegm < NSEGS; iSegm++) {
// remove contribution for all segments
//if(mask[iSegm] != '0') {
gs_type *baseTrace = pPsa->amplitude[iSegm];
gs_type *realTrace = fAmplitude[iSegm];
gs_type *baseTrace = pPsa->amplitude[iSegm]; // base point to be subtracted
gs_type *realTrace = fAmplitude[iSegm]; // from the experimental data
if(iSegm == netChargeSegment) {
baseTrace = pPsa->netcharge;
}
......
......@@ -3,8 +3,7 @@
#include "commonDefs.h"
#define VER "PSA code: AGS(Apr2009) by RV"
#define VER "PSA code: AGS(May2010) by RV&DB"
#define GS_FLOAT // signal and signal basis used as float instead of short
......@@ -20,12 +19,16 @@ typedef short gs_type;
#ifdef USENETCHARGES
const int DIFFLAG = 0; // use the net-charge signals as they are
#else
const int DIFFLAG = 5; // use them with a delayed-differentiation (in units of signal samples) if !=0
//const int DIFFLAG = 0; // don't use the net-charge segments
const int DIFFLAG = 5; // use them with a delayed-differentiation (in units of signal samples)
//# define SHOWDIFFERENTIATED // in the saved traces show the differentiated version of the net-charge segments
#endif
#define ORDEREDSEARCH // energy-ordered decomposition with removal of previous signals
//#define SHOWPARTIALS // if using the energy-ordered decomposition write the partial traces
#ifdef ORDEREDSEARCH
# define PRECOARSESEARCH // a coarse-only preliminary search to remove the other net-charges
//# define WRITEPARTIALS // if using the energy-ordered decomposition write the partial traces
#endif
const double METRIC = 0.3; // norm for the figure of merit
const int NMETRIC = 131072; // 2^17
......
......@@ -915,11 +915,28 @@ void SignalBasis::ExperimentalResponse()
double kernelsegs[] = {0.006224, 0.125911, 0.219505, 0.160006, 0.122884, 0.106105, 0.074517, 0.049438, 0.033024, 0.029816, 0.020622, 0.014916, 0.009906, 0.009114, 0.005847, 0.002675, 0.004706, 0.001936, 0.002851};
// from the average of the response function of cores of crystal 2B and 2G
double kernelcore[] = {0.030127, 0.235100, 0.289680, 0.154834, 0.092118, 0.063705, 0.042655, 0.028896, 0.020688, 0.013784, 0.009559, 0.006759, 0.005118, 0.003887, 0.003090};
const int shiftleft = 0;
// gaussian with 35 ns fwhm
//double kernelsegs[] = {0.000004,0.000078,0.000936,0.007179,0.035006,0.108547,0.214045,0.268411,0.214045,0.108547,0.035006,0.007179,0.000936,0.000078,0.000004};
//double kernelcore[] = {0.000004,0.000078,0.000936,0.007179,0.035006,0.108547,0.214045,0.268411,0.214045,0.108547,0.035006,0.007179,0.000936,0.000078,0.000004};
//const int shiftleft = 6;
// gaussian with 40 ns fwhm
//double kernelsegs[] = {0.000004,0.000048,0.000459,0.003086,0.014679,0.049373,0.117430,0.197492,0.234859,0.197492,0.117430,0.049373,0.014679,0.003086,0.000459,0.000048,0.000004};
//double kernelcore[] = {0.000004,0.000048,0.000459,0.003086,0.014679,0.049373,0.117430,0.197492,0.234859,0.197492,0.117430,0.049373,0.014679,0.003086,0.000459,0.000048,0.000004};
//const int shiftleft = 5;
//// gaussian with 45 ns fwhm
//double kernelsegs[] = {0.000003,0.000033,0.000255,0.001510,0.006809,0.023348,0.060882,0.120727,0.182051,0.208764,0.182051,0.120727,0.060882,0.023348,0.006809,0.001510,0.000255,0.000033,0.000003};
//double kernelcore[] = {0.000003,0.000033,0.000255,0.001510,0.006809,0.023348,0.060882,0.120727,0.182051,0.208764,0.182051,0.120727,0.060882,0.023348,0.006809,0.001510,0.000255,0.000033,0.000003};
//const int shiftleft = 6;
//// gaussian with 50 ns fwhm
//double kernelsegs[] = {0.000003,0.000024,0.000155,0.000820,0.003467,0.011743,0.031861,0.069249,0.120570,0.168164,0.187887,0.168164,0.120570,0.069249,0.031861,0.011743,0.003467,0.000820,0.000155,0.000024,0.000003};
//double kernelcore[] = {0.000003,0.000024,0.000155,0.000820,0.003467,0.011743,0.031861,0.069249,0.120570,0.168164,0.187887,0.168164,0.120570,0.069249,0.031861,0.011743,0.003467,0.000820,0.000155,0.000024,0.000003};
//const int shiftleft = 9;
for(int ii=0; ii<NSEGS; ii++) {
int npoints = numPts[ii];
for(int jj=0; jj<npoints; jj++) {
Pts[ii][jj].convsegs(kernelsegs, sizeof(kernelsegs)/sizeof(double));
Pts[ii][jj].convcore(kernelcore, sizeof(kernelcore)/sizeof(double));
Pts[ii][jj].convsegs(kernelsegs, sizeof(kernelsegs)/sizeof(double), shiftleft);
Pts[ii][jj].convcore(kernelcore, sizeof(kernelcore)/sizeof(double), shiftleft);
}
}
}
......
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