Commit 7dbf5f95 authored by dino's avatar dino
Browse files

Added M=1 spectra in CrystalProducerATCA (useful for calibrations)

Added fit of T0 after PSA

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@1031 170316e4-aea8-4b27-aad4-0380ec0519c9
parent 691a53fe
......@@ -21,9 +21,9 @@ MESSAGE("Root: ${ROOT_INCLUDE_DIR} ${ROOT_LIBRARY_DIR}")
FIND_PACKAGE(Boost)
MESSAGE("Boost: ${Boost_INCLUDE_DIR} ${Boost_LIBRARY_DIRS}")
# Enable the vectorized version of the PSA FOM-loop
ADD_DEFINITIONS(-DUSE_SSE_VERSION)
ADD_DEFINITIONS(-msse4)
## Enable the vectorized version of the PSA FOM-loop
#ADD_DEFINITIONS(-DUSE_SSE_VERSION)
#ADD_DEFINITIONS(-msse4)
# The sources to build femul
SET(ADFDIR "../myADF0.2")
......
......@@ -32,13 +32,15 @@ const int numXYZ = 10; // sub matrXYZ for up to numXYZ-1 hits
const int nPsaHitsValues = 16; // number of parameters written in Psa__-nn-F__Hits.fdat
const bool doFit_T0 = false;
const bool doFit_T0 = true;
const float SMOOTH[7] = {35,35,35,35,35,35,35}; // default decay constant of the preamp response function for the 6 slices and the core (in nanoseconds)
const float minSegEnergy = 10.f; // energy threshold for a segment to be part of the decomposition
const float lowSegEnergy = 50.f; // segment-energies smaller than this are decomposed in a simplified way (and histogrammed separately)
//#define SHOWDIFFERENTIATED // in the saved traces show the differentiated version of the net-charge segments
//#define SELECTEVENTS // "manual" selection of events in ProcessEvents()
//#define SPECMETRICS // to produce specMetrics in SearchFullGrid and SearchAdaptive ==> slow and not thread-safe
//#define TRACESXYZ // WriteTraces() sorts on position files xx_yy_zz_xyz.dat
//#define MATCHI2 // 3D matrix with distribution of chi2; meaningful only for segmult=1; done only in SearchFullGrid
......@@ -83,7 +85,7 @@ PSAFilterGridSearch::PSAFilterGridSearch()
baseZXY = NULL; specEner = NULL; specStat = NULL; specTzero = NULL; specSigma = NULL; specMetrics = NULL;
matrXYZR = NULL; matrSeg = NULL; matrRZE = NULL;
#endif //PSA_LOCALSPECTRA
fHitSegThreshold = 10; // keV
fHitSegThreshold = minSegEnergy; // keV
#ifdef PSA_FromGRU_
//PSACCE_raw = new TH1I ("PSACCE_raw","PSACCE_raw",32768,0,32767);
......@@ -334,7 +336,7 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
// prepare the metrics vectors for the chi2 loops
float largest = 0;
for(int i=0; i<NMETRIC; i++) {
int iv = i >= RMETRIC ? i-RMETRIC : RMETRIC-i;
int iv = (i >= RMETRIC) ? (i-RMETRIC) : (RMETRIC-i); // abs()
float fv = (float)pow(double(iv), METRIC);
fMetrics[i] = fv;
if(fv > largest) largest = fv;
......@@ -359,12 +361,14 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
hGroup.add(specStat);
if(!gPsaSegCenter) {
//specTzero = new nDhist<unsigned int>(NSEGS+4, specLenT);
//specTzero->setFileName(fOdirPrefix+"Psa?Tzero.spec");
//specTzero->setComment ("T90 and best dt (if relevant)");
//hGroup.add(specTzero);
if(doFit_T0) {
specTzero = new nDhist<unsigned int>(NSEGS+4, specLenT);
specTzero->setFileName(fOdirPrefix+"Psa?Tzero.spec");
specTzero->setComment ("T90 and best dt (if relevant)");
hGroup.add(specTzero);
}
//specSigma = new nDhist<unsigned int>(NSEGS+4, specLenT);
//specSigma = new nDhist<unsigned int>(NSEGS+4, 10000);
//specSigma->setFileName(fOdirPrefix+"Psa?Sigma.spec");
//specSigma->setComment ("chisquare");
//hGroup.add(specSigma);
......@@ -678,7 +682,7 @@ Int_t PSAFilterGridSearch::Process(int slot, int nslots)
memcpy(pD->PsaOut + nn, pS.PsaOut + nn, sizeof(PsaOut_t));
}
if(doFit_T0) {
float dt0 = FitT0FromCore(pS, pS.samp_first + pS.samp_used);
float dt0 = FitT0AfterPSA(pS, pS.samp_first + pS.samp_used);
pD->shiftT0 = dt0; // samples
}
// While writing the traces, save the results in the proper slot
......@@ -733,14 +737,21 @@ int PSAFilterGridSearch::ProcessTheEvent(pointFull &pS)
pS.bestPt1 = -1;
bool onlyCoarse = gCoarseOnly;
bool simplified = false;
#ifdef SIMPLIFY_SMALL_ENERGIES
if(snn>=2 && netChEner<100.f) // if mult>2 and the energy of this point is small, make only the coarse search
onlyCoarse = true; // the condition should be refined but this saves a lot of useless calculations
if(snum>1 && netChEner<lowSegEnergy) { // if mult>1 and the energy of this point is small, make only the coarse search
onlyCoarse = true; // the condition should be refined but this saves a lot of useless calculations
simplified = true;
}
#endif
if(gUseAdaptive || onlyCoarse) {
// 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
pS.samp_used = SearchAdaptive1(&pS, netChSeg, onlyCoarse, onlyFirst); // coarse-fine or coarse-only with one hit
if(simplified || gCoarseOnly) {
int tmp = pS.bestPt1; pS.bestPt1 = pS.randPt1; pS.randPt1 = tmp;
// REMEMBER TO DO IT ALSO IN ProcessTwoHits
}
pS.bestFac = 1.f; // second point not present
pS.bestPt2 = -1;
}
......@@ -760,7 +771,7 @@ int PSAFilterGridSearch::ProcessTheEvent(pointFull &pS)
pOut->fy1 = bestPoint1.y + yInterp();
pOut->fz1 = bestPoint1.z + zInterp();
pOut->dTns = 0;
pOut->chiSq = pS.chi2min; // not clear what kind of information is given by this number
pOut->chiSq = pS.chi2min;
pOut->netChargeSeg = netChSeg;
pOut->bestPoint1 = pS.bestPt1;
pOut->bestPoint2 = -1;
......@@ -907,7 +918,7 @@ int PSAFilterGridSearch::ProcessTwoHits(pointFull &pS)
pOut->fy2 = bestPoint2.y + yInterp();
pOut->fz2 = bestPoint2.z + zInterp();
pOut->dTns = 0;
pOut->chiSq = pS.chi2min; // not clear what kind of information is given by this number
pOut->chiSq = pS.chi2min;
pOut->netChargeSeg = netChSeg;
pOut->bestPoint1 = pS.bestPt1;
pOut->bestPoint2 = pS.bestPt2;
......@@ -1574,8 +1585,14 @@ int PSAFilterGridSearch::SearchAdaptive1(pointFull *pS, int netChSeg, bool bCoar
pS->bestPt1 = bestPt;
pS->chi2min = chi2min;
if(bCoarseOnly)
if(bCoarseOnly) {
// in case the caller wants to smooth the data, select at random a point of the associated fine grid
int jloop1 = sgrid.posCoarse[kbest];
int jloop2 = sgrid.posCoarse[kbest+1];
int jloop = jloop1 + rand()%(jloop2-jloop1);
pS->randPt1 = sgrid.indFine[jloop];
return uSamples;
}
/////////////////////////////
// loop on the fine grid //
......@@ -1935,8 +1952,8 @@ void PSAFilterGridSearch::PreSearchCoarse(pointFull &pS)
bool onlyCenter = false;
#ifdef SIMPLIFY_SMALL_ENERGIES
if(snn>=2 && netChEner<100.f) // if mult>2 and the energy of this point is small, make only the coarse search
onlyCenter = true; // the condition should be refined but this saves a lot of useless calculations
if(snum>1 && netChEner<lowSegEnergy) // if mult>1 and the energy of this point is small, make only the coarse search
onlyCenter = true; // the condition should be refined but this saves a lot of useless calculations
#endif
int samp_used = SearchAdaptive1(&pS, netChSeg, true, onlyCenter);
......@@ -2026,7 +2043,7 @@ Int_t PSAFilterGridSearch::PostProcess(int slot)
if(pD->isSelected) {
PsaOut_t * pOut = pD->PsaOut;
for(int nh = 0; nh < numNetSegs; nh++, pOut++) {
int netChSeg = pOut->netChargeSeg;
int netChSeg = pOut->netChargeSeg;
float fx = pOut->fx1;
float fy = pOut->fy1;
float fz = pOut->fz1;
......@@ -2039,38 +2056,111 @@ Int_t PSAFilterGridSearch::PostProcess(int slot)
#endif
// positions inside the crystal
if(matrXYZR) {
int sz = netChSeg%6;
int nz = (int)(fz/10);
matrXYZR->Incr(0, 0, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(0, 1, int(fx+matOff), int(fz ));
matrXYZR->Incr(0, 2, int(fy+matOff), int(fz ));
matrXYZR->Incr(0, 3, int(rr+matOff), int(fz ));
matrXYZR->Incr(0, 4+sz, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(0,10+nz, int(fx+matOff), int(fy+matOff));
nz = (int)(fz/2);
matrXYZR->Incr(0,20+nz, int(fx+matOff), int(fy+matOff));
if(nh < numXYZ-1) {
int nz = (int)(fz/10);
matrXYZR->Incr(nh+1, 0, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nh+1, 1, int(fx+matOff), int(fz ));
matrXYZR->Incr(nh+1, 2, int(fy+matOff), int(fz ));
matrXYZR->Incr(nh+1, 3, int(rr+matOff), int(fz ));
matrXYZR->Incr(nh+1, 4+sz, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nh+1,10+nz, int(fx+matOff), int(fy+matOff));
nz = (int)(fz/2);
matrXYZR->Incr(nh+1,20+nz, int(fx+matOff), int(fy+matOff));
int sz = netChSeg%6;
int nz = (int)(fz/10);
int nw = (int)(fz/2 );
if(pOut->enerSG >= lowSegEnergy) {
// not too-low energy
{ // 0 --> all hits
matrXYZR->Incr(0, 0, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(0, 1, int(fx+matOff), int(fz ));
matrXYZR->Incr(0, 2, int(fy+matOff), int(fz ));
matrXYZR->Incr(0, 3, int(rr+matOff), int(fz ));
matrXYZR->Incr(0, 4+sz, int(fx+matOff), int(fy+matOff)); // the slices
matrXYZR->Incr(0,10+nz, int(fx+matOff), int(fy+matOff)); // every cm
matrXYZR->Incr(0,20+nw, int(fx+matOff), int(fy+matOff)); // every 2 mm
}
if(nh < numXYZ-3) {
// 1--6 --> 1st...6th hit
int nhi = nh + 1;
matrXYZR->Incr(nhi, 0, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi, 1, int(fx+matOff), int(fz ));
matrXYZR->Incr(nhi, 2, int(fy+matOff), int(fz ));
matrXYZR->Incr(nhi, 3, int(rr+matOff), int(fz ));
matrXYZR->Incr(nhi, 4+sz, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi,10+nz, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi,20+nw, int(fx+matOff), int(fy+matOff));
}
}
else {
// low energy
{ // 7 --> all hits
int nhi = numXYZ-3;
matrXYZR->Incr(nhi, 0, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi, 1, int(fx+matOff), int(fz ));
matrXYZR->Incr(nhi, 2, int(fy+matOff), int(fz ));
matrXYZR->Incr(nhi, 3, int(rr+matOff), int(fz ));
matrXYZR->Incr(nhi, 4+sz, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi,10+nz, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi,20+nw, int(fx+matOff), int(fy+matOff));
}
if(nh==0) {
// 8 --> first hit
int nhi = numXYZ-2;
matrXYZR->Incr(nhi, 0, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi, 1, int(fx+matOff), int(fz ));
matrXYZR->Incr(nhi, 2, int(fy+matOff), int(fz ));
matrXYZR->Incr(nhi, 3, int(rr+matOff), int(fz ));
matrXYZR->Incr(nhi, 4+sz, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi,10+nz, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi,20+nw, int(fx+matOff), int(fy+matOff));
}
if(nh==numNetSegs-1) {
// 9 --> last hit for low energy
int nhi = numXYZ-1;
matrXYZR->Incr(nhi, 0, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi, 1, int(fx+matOff), int(fz ));
matrXYZR->Incr(nhi, 2, int(fy+matOff), int(fz ));
matrXYZR->Incr(nhi, 3, int(rr+matOff), int(fz ));
matrXYZR->Incr(nhi, 4+sz, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi,10+nz, int(fx+matOff), int(fy+matOff));
matrXYZR->Incr(nhi,20+nw, int(fx+matOff), int(fy+matOff));
}
}
}
if(matrSeg) {
matrSeg->Incr(0,netChSeg, 0, int(fx+matOff), int(fy+matOff));
matrSeg->Incr(0,netChSeg, 1, int(fx+matOff), int(fz ));
matrSeg->Incr(0,netChSeg, 2, int(fy+matOff), int(fz ));
matrSeg->Incr(0,netChSeg, 3, int(rr+matOff), int(fz ));
if(nh < numXYZ-1) {
matrSeg->Incr(nh+1,netChSeg, 0, int(fx+matOff), int(fy+matOff));
matrSeg->Incr(nh+1,netChSeg, 1, int(fx+matOff), int(fz ));
matrSeg->Incr(nh+1,netChSeg, 2, int(fy+matOff), int(fz ));
matrSeg->Incr(nh+1,netChSeg, 3, int(rr+matOff), int(fz ));
if(pOut->enerSG >= lowSegEnergy) {
// not too-low energy
{ // 0 --> all hits
matrSeg->Incr(0,netChSeg, 0, int(fx+matOff), int(fy+matOff));
matrSeg->Incr(0,netChSeg, 1, int(fx+matOff), int(fz ));
matrSeg->Incr(0,netChSeg, 2, int(fy+matOff), int(fz ));
matrSeg->Incr(0,netChSeg, 3, int(rr+matOff), int(fz ));
}
if(nh < numXYZ-3) {
// 1--6 --> 1st...6th hit
int nhi = nh + 1;
matrSeg->Incr(nhi,netChSeg, 0, int(fx+matOff), int(fy+matOff));
matrSeg->Incr(nhi,netChSeg, 1, int(fx+matOff), int(fz ));
matrSeg->Incr(nhi,netChSeg, 2, int(fy+matOff), int(fz ));
matrSeg->Incr(nhi,netChSeg, 3, int(rr+matOff), int(fz ));
}
}
else {
// low energy
{ // 7 --> all hits
int nhi = numXYZ-3;
matrSeg->Incr(nhi,netChSeg, 0, int(fx+matOff), int(fy+matOff));
matrSeg->Incr(nhi,netChSeg, 1, int(fx+matOff), int(fz ));
matrSeg->Incr(nhi,netChSeg, 2, int(fy+matOff), int(fz ));
matrSeg->Incr(nhi,netChSeg, 3, int(rr+matOff), int(fz ));
}
if(nh==0) {
// 8 --> first hit
int nhi = numXYZ-2;
matrSeg->Incr(nhi,netChSeg, 0, int(fx+matOff), int(fy+matOff));
matrSeg->Incr(nhi,netChSeg, 1, int(fx+matOff), int(fz ));
matrSeg->Incr(nhi,netChSeg, 2, int(fy+matOff), int(fz ));
matrSeg->Incr(nhi,netChSeg, 3, int(rr+matOff), int(fz ));
}
if(nh==numNetSegs-1) {
// 9 --> last hit for low energy
int nhi = numXYZ-1;
matrSeg->Incr(nhi,netChSeg, 0, int(fx+matOff), int(fy+matOff));
matrSeg->Incr(nhi,netChSeg, 1, int(fx+matOff), int(fz ));
matrSeg->Incr(nhi,netChSeg, 2, int(fy+matOff), int(fz ));
matrSeg->Incr(nhi,netChSeg, 3, int(rr+matOff), int(fz ));
}
}
}
if(matrRZE && numNetSegs==1) {
......@@ -2116,14 +2206,14 @@ Int_t PSAFilterGridSearch::PostProcess(int slot)
specStat->Incr(numNetSegs, netChSeg); // distribution of segments for n-hits
specStat->Incr(36, netChSeg); // total distribution of segments
}
if(specTzero) {
// t0 after grid-search (possibly not up to date)
int bestdt = (int)(pOut->dTns+specLenT/2); // this is already in ns
specTzero->incr(netChSeg, bestdt);
}
//if(specTzero) {
// // t0 after grid-search (possibly not up to date) ==> indeed dTns is always 0
// int bestdt = (int)(pOut->dTns+specLenT/2); // this is already in ns
// specTzero->incr(netChSeg, bestdt);
//}
if(specSigma) {
// chisquare of search
int chi2min = (int)(pOut->chiSq);
int chi2min = (int)(pOut->chiSq)>>8; // ad-hoc adjustment to the spectrum range
specSigma->Incr(netChSeg, chi2min);
specSigma->Incr(NSEGS, chi2min); // OR of all
}
......@@ -2144,11 +2234,11 @@ Int_t PSAFilterGridSearch::PostProcess(int slot)
if(specTzero) {
// These times are not very meaningful as they are referred to the beginning of the trace
int bestT0 = int( pD->CoreT[0] *SAMP_STEP);
specTzero->incr(NSEGS, bestT0); // original distribution
specTzero->incr(NSEGS, bestT0); // original distribution
int bestT1 = int((pD->CoreT[0]+pD->shiftT0)*SAMP_STEP);
specTzero->incr(NSEGS+1, bestT1); // corrected distribution
specTzero->incr(NSEGS+1, bestT1); // corrected distribution
int bestT2 = int((pD->CoreT[0]-pD->shiftT0)*SAMP_STEP);
specTzero->incr(NSEGS+2, bestT2); // corrected distribution
specTzero->incr(NSEGS+2, bestT2); // corrected distribution (with reversed sign)
int bestDT = int(pD->shiftT0*SAMP_STEP + specLenT/2);
specTzero->incr(NSEGS+3, bestDT); // the correction
}
......@@ -2402,7 +2492,7 @@ void PSAFilterGridSearch::MatChi2(float chi, pointPsa* pPsa)
// not (yet?) working well
// returns the time shift with respect to SAMP0 in units of samples
Float_t PSAFilterGridSearch::FitT0FromCore(pointFull &pS, int tsamp)
Float_t PSAFilterGridSearch::FitT0AfterPSA(pointFull &pS, int tsamp)
{
float tData[NTIME];
float tRfit[NTIME];
......@@ -2433,7 +2523,7 @@ Float_t PSAFilterGridSearch::FitT0FromCore(pointFull &pS, int tsamp)
for(int ll = 0; ll <= 2*nsteps; ll++) {
float *tfit = rfit + nsteps - ll;
float x2t = 0;
// calculate chi2 using 5=SAMP0-3 samples of pretrigger and 20 of risetime
// calculate chi2 using SAMP0-nsteps (8-3=5) samples of pretrigger and 20 of risetime
for(int ii = nsteps; ii < 25; ii++) {
float diff = float(data[ii] - tfit[ii]);
x2t += diff*diff;
......@@ -2441,6 +2531,7 @@ Float_t PSAFilterGridSearch::FitT0FromCore(pointFull &pS, int tsamp)
x2[ll] = x2t;
}
// find the minimum
int ssmin = nsteps;
float x2min = x2[nsteps];
for(int ll = 0; ll <= 2*nsteps; ll++) {
......@@ -2450,9 +2541,11 @@ Float_t PSAFilterGridSearch::FitT0FromCore(pointFull &pS, int tsamp)
}
}
// if at the edge
if(ssmin == 0 || ssmin == 2*nsteps)
return ssmin+0.f-nsteps;
// parabolic interpolation around the minimum
float ym = x2[ssmin-1];
float ys = x2[ssmin ];
float yp = x2[ssmin+1];
......
......@@ -158,7 +158,7 @@ public:
void AddToSolution (pointFull &pS, pointPsa &bestPoint, int numsamp, float scaleFact);
void AddToSolution (pointFull &pS, pointPsa &bestPoint1, pointPsa &bestPoint2, float bestFac, int numsamp, float scaleFact);
void SaveTotalTrace(pointFull &pS, PsaData *pD);
Float_t FitT0FromCore (pointFull &pS, int tsamp);
Float_t FitT0AfterPSA(pointFull &pS, int tsamp);
private:
void Extend(float *dd, int n1, int n2) {
......@@ -173,6 +173,9 @@ private:
float xInterp() {return (rand()%200 - 100)*0.01f; }
float yInterp() {return (rand()%200 - 100)*0.01f; }
float zInterp() {return (rand()%200 - 125)*0.01f; }
//float xInterp() {return (gCoarseOnly) ? 0 : (rand()%200 - 100)*0.01f; }
//float yInterp() {return (gCoarseOnly) ? 0 : (rand()%200 - 100)*0.01f; }
//float zInterp() {return (gCoarseOnly) ? 0 : (rand()%200 - 125)*0.01f; }
private:
......@@ -432,9 +435,12 @@ INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const float *pReal, const
float chi2 = 0;
for(int nn = 0; nn < loopSamples; nn++) {
float diff = pReal[nn] - pBase[nn];
#ifdef SPECMETRICS // be sure that specMetrics is defined and that you are NOT USING THREADS
#ifdef SPECMETRICS
// This block is to check that the size of "metrics", is big enough to address all
// possible values of "diff". If enabled, it slows very much down the execution.
// As nDhist is not thread safe, be sure that you are NOT USING THREADS.
int idiff = int(diff + RMETRIC);
specMetrics->incr(0, idiff);
if(specMetrics) specMetrics->incr(0, idiff);
if(idiff<=0 || idiff>=2*RMETRIC-1) {
int nerrr = 1; // to catch the outliers with the debugger
}
......@@ -455,7 +461,9 @@ INLINE_ALWAYS float PSAFilterGridSearch::Chi2InnerLoop(const int *pReal, const i
float chi2 = 0;
for(int nn = 0; nn < loopSamples; nn++) {
int diff = pReal[nn] - pBase[nn];
#ifdef DORESCALE
diff >>= RESCALE;
#endif
#ifdef SPECMETRICS // be sure that specMetrics is defined and that you are NOT USING THREADS
int idiff = int(diff + RMETRIC);
specMetrics->incr(0, idiff);
......
......@@ -387,6 +387,8 @@ public:
int netSeg; // the specific net-charge segment during grid-search processing
int bestPt1; // the best matching first base point
int bestPt2; // the best matching second base point
int randPt1; // a random point of the fine grid associated to the coarse of bestPt1
int randPt2; // a random point of the fine grid associated to the coarse of bestPt2
float bestFac; // relative amplitude of the first interaction point
float chi2min; // not very meaningful outside the grid-search
......@@ -411,6 +413,8 @@ private:
netSeg = 0;
bestPt1 = 0;
bestPt2 = 0;
randPt1 = 0;
randPt2 = 0;
bestFac = 0;
chi2min = 0;
}
......
......@@ -39,6 +39,7 @@ const int NMETRIC = 2*RMETRIC; // positive and negative
const int MAXNORM = 1<<21; // normalization of the traces for the search loops (i.e. precision of the metrics vector)
const int MAXMORE = (1<<24)-MAXNORM-2; // max value of the rescaled experimental trace
const int RESCALE = 8; // divide by 2^RSCALE to match RMETRIC
#define DORESCALE
#elif GS_TYPE == GS_TYPE_INT16
typedef short gs_type; // averPt of pointPsa can have overflow problems with his type
const int MAXNORM = 4096; // normalization of the traces for the search loops (i.e. precision of the metrics vector)
......@@ -48,18 +49,19 @@ const int NMETRIC = 2*RMETRIC; // positive and negative
# error GS_TYPE invalid value
#endif
#define USE_NETCHARGES // to use the net charge signals in the grid search
//#define USE_CORETRACE // to use also the core signal in the grid search
//#define SIMPLIFY_SMALL_ENERGIES // simplified search for segments with low energy if multiplicity > 2
#define USE_NETCHARGES // to use the net charge signals in the grid search
#define USE_CORETRACE // to use also the core signal in the grid search
#define SIMPLIFY_SMALL_ENERGIES // simplified search for segments with low energy if multiplicity > 2
#ifdef USE_NETCHARGES
const int DIFFLAG = 0; // use the net-charge signals as they are
const int DIFFLAG = 0; // use the net-charge signals as they are
#else
const int DIFFLAG = 8; // use the net-charge segments with a delayed-differentiation (in units of signal samples)
const int DIFFLAG = 8; // use the net-charge segments with a delayed-differentiation (in units of signal samples)
#endif
// INITATSEGCENTER costs almost nothing in terms of cpu and is good enough as pre-positioning;
// PRECOARSESEARCH costs a lot and gives no further gaing
// PRECOARSESEARCH costs a lot and gives no further gain
#define INITATSEGCENTER // initialize the search by putting the interaction points at the center of their segment
//#define PRECOARSESEARCH // a coarse-only preliminary search to remove the effect of the other net-charge segments
......@@ -77,7 +79,7 @@ const int TIME0 = 0; // the "T=0" sample of the calculated trac
const int SAMP_STEP = 10; // (ns) time step of the experimental data
const int NSAMP = defTraceLengthPSA; // length of experimental data (taken as float data) (defTraceLengthPSA should be 60)
const int SAMP0 = (defTriggerSample/4)*4; // pretrigger of the experimental traces (defTriggerSample should be 10)
const int SAMP0 = (defTriggerSample/4)*4; // pretrigger of the experimental traces ( = 8 ) (defTriggerSample should be 10)
#define SFIXEDLEN 40 // number of samples used to calculate the Figure Of Merit in the Grid Search
// If SFIXEDLEN==0, the number of samples is determined for each event as T95 (95%) of the core amplitude
......
......@@ -53,7 +53,7 @@ const float tScale = 10.f; // tScale=10 ==> 1ns/ch
const int tGain = 10; // 1; // to increase the scale and length of ...TT1 and ...TT2 spectra in a (more or less) consistent way
const int speTTlen = int(1000*tScale);
#if defined(_DEBUG) && 1 // set to 1 to enable writing the intermediate waves while debugging
#if defined(_DEBUG) && 0 // set to 1 to enable writing the intermediate waves while debugging
# define WRITEWORKINGWAVE1(p1,p2) WriteWorkingWave(p1, p2)
#ifdef TIMING_CFD
# define WRITEWORKINGWAVE2(p1,p2) WriteWorkingWave(p1, p2)
......@@ -1280,7 +1280,7 @@ bool PreprocessingFilterPSA::ProcessEvent(int &smult, float &eSumSG1, float &e
float fHandShift = CC.pCC[iTR].tmove; // CC.pCC[iTR].tmove is interpreted as global shift and should be given in units of samples !!
float fshiftCC = (defTriggerSample+defTriggerSamplePlus) - tRefCC + fHandShift;
for(int nc = 0; nc < nCC; nc++) {
ShiftMoveTrace(fTracesCC[nc], nNsamp, fshiftCC + 1.f, tmp); // no extra shift for the cores
ShiftMoveTrace(fTracesCC[nc], nNsamp, fshiftCC, tmp); // no extra shift for the cores
memcpy(fTracesCC[nc], tmp, sizeof(tmp));
}
for(int ns = 0; ns < nSG; ns++) {
......
......@@ -15,7 +15,7 @@
// Until this is fixed, it is better to re-time the signals using the CFD, which is less sensitive to the problem.
// By enabling the following the re-timing of the signals inside the captured traces is done by a CFD (see MWD for details of this "interpolated" version)
//#define TIMING_CFD
#define TIMING_CFD
#ifdef TIMING_CFD
# include "mwdlib.h"
......
......@@ -1189,8 +1189,8 @@ bool TrackingFilter::AnalysisOfCrystals()
if(id2 < 0 || id2 >= NDETS) // this should not happen and should be handled as an error ??
continue;
dts = (int)(pCr2->tstamp - pCr1->tstamp);
dtf = pCr2->CCT0 - pCr1->CCT0 ;
dtp = pCr2->CCDT - pCr1->CCDT ;
dtf = pCr2->CCT0 - pCr1->CCT0;
dtp = pCr2->CCDT - pCr1->CCDT;
dtt = tScale*(dts + dtf) ;
OftSpec_TS->Incr(id1, id2, dts+specLenT1/2); // time-time from tstamp
OftSpec_TT->Incr(id1, id2, int(tr0+dtt) ); // time-time from tstamp+CFD
......
......@@ -82,15 +82,15 @@ protected:
UInt_t id;
UInt_t nhits; // number of hits
UInt_t hit1st; // first hit in pEXYZ
UInt_t idSeg;
UInt_t evnum;
ULong64_t tstamp;
Float_t CCE0;
Float_t CCE1;
Float_t CCT0;
Float_t CCT1;
Float_t CCDT;
Float_t Esum;
UInt_t idSeg; // ID of segment with the maximum energy release
UInt_t evnum; // event number
ULong64_t tstamp; // time stamp (samples)
Float_t CCE0; // energy of triggering core (keV)
Float_t CCE1; // energy of the other core (keV)
Float_t CCT0; // trigger time (samples) of the core relative to the timestamp, possibly corrected in the PSA
Float_t CCT1; // "" of the other core (usually the same)
Float_t CCDT; // average time of the event derived from the segments (ns)
Float_t Esum; // sum of segment energies (value not forced to CoreE[0])
} crystdata;
// General tracking structures
......
......@@ -45,6 +45,10 @@ fTrigger("data:crystal")
fStopErrorCode = 111;
fNumFilesIn = 1;
fMinSegAmpli = 500.f; // ~50 (150) keV for high(low) gain range if shaping = 1000
fSpekGain = 1.f; // good for trapezoid risetime 1000
fInputDataFile.clear();
fVerbose = false;
......@@ -221,7 +225,9 @@ void CrystalProducer::GetParameters( UInt_t *error_code )
}
else if( stringEq(keyw, "WriteBaseLines") ) {
ok = 1 == sscanf(data.c_str(), "%d", &fWriteBaselines);
ok = true;
}
else if( stringEq(keyw, "AmplitudeGain") ) {
ok = 1 == sscanf(data.c_str(), "%f", &fSpekGain);
}
else if( stringEq(keyw, "TimeStep") ) {
ok = 1 == sscanf(data.c_str(), "%f", &fTimeStep);
......
......@@ -74,6 +74,8 @@ protected:
UInt_t fStopErrorCode; // how to stop analysis on EOF (> 100 stop everything) default 111
UInt_t fWriteBaselines; // if > 0 produces the distribution of baselines, decimating events by fWriteBaselines
Int_t fNumFilesIn; // run on a sequence of input files
Float_t fMinSegAmpli; // minimum segment amplitude to detect segment multiplicity
Float_t fSpekGain; // rescaling of amplitude spectra to ease energy calibration for short shaping times
std::string fInputDataFile; // filename of raw data, used to simplify the decoding of CrystalProducerATCA.conf
Bool_t fVerbose;
#ifdef CP_LOCALSPECTRA
......
......@@ -52,7 +52,7 @@ tStampErr(0), tStampNew(0)
fProdBaseCount = 0;
#ifdef CP_LOCALSPECTRA
ProdSpec_Ener = NULL;
ProdSpec_Ampl = NULL;
ProdSpec_Base = NULL;
ProdSpec_Rate = NULL;
#endif //CP_LOCALSPECTRA
......@@ -154,12 +154,12 @@ Int_t CrystalProducerATCA::AlgoSpecificInitialise()
}