Commit 48929af3 authored by dino's avatar dino
Browse files

Added a loop over PSA-FitT0AfterPSA to optimize T0 in the PSA.

This lowers the speed by ~30% and we will see if it is worthwhile keepint it.

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@1032 170316e4-aea8-4b27-aad4-0380ec0519c9
parent 7dbf5f95
......@@ -955,6 +955,10 @@
RelativePath="..\common\cycleServer.h"
>
</File>
<File
RelativePath="..\common\dataPairFile.h"
>
</File>
<File
RelativePath="..\common\misc.cpp"
>
......
......@@ -373,6 +373,8 @@ Int_t PSAFilter::SetInput(int slot)
pD->evnumber = ((AgataKey *)fFrameCrystal->GetKey())->GetEventNumber();
pD->timestamp = ((AgataKey *)fFrameCrystal->GetKey())->GetTimeStamp();
pD->numNetSegs = 0;
pD->shiftT0 = 0;
pD->tCycles = 0;
return 0;
}
......
......@@ -93,6 +93,7 @@ struct PsaData
UInt_t evnumber;
ULong64_t timestamp;
Float_t shiftT0; // samples
Int_t tCycles;
// input data (length of traces managed by the owner of the object)
Float_t *fTracesSG[ADF::CrystalInterface::kNbSegments];
......
......@@ -49,7 +49,7 @@ const float lowSegEnergy = 50.f; // segment-energies smaller than this are de
// Therefore it makes sense only in Debug mode by setting break points at the instruction immediately after the call
// to PrepareEvent in Process().
// For the same reasons it is incompatible with threads and multiple crystal analyses.
#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 WRITEWORKINGWAVE(p1,p2) WriteWorkingWave(p1, p2)
# define WRITEPARTIALWAVE(p1,p2,p3) WritePartialWave(p1, p2, p3)
#else
......@@ -64,6 +64,10 @@ PSAFilterGridSearch::PSAFilterGridSearch()
fpPsaHits = NULL;
fPsaIndex = 0;
#ifdef TSHIFT_CORRECTION
tstmpFileO = NULL;
#endif
// Chi2InnerLoop dafault values
SetLoopSamples(NTIME);
SetLoopFactors(1.0f, 0.0f);
......@@ -341,9 +345,23 @@ Int_t PSAFilterGridSearch::AlgoSpecificInitialise()
fMetrics[i] = fv;
if(fv > largest) largest = fv;
}
//const float max_xx_max = 50.f*50.f;
//for(int i=0; i<NMETRIC; 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;
// if(fv > max_xx_max) fMetrics[i] = max_xx_max;
//}
largest++; // this line is just to be able to set a breakpoint
}
#ifdef TSHIFT_CORRECTION
if(doFit_T0) {
tstmpFileO = new dataPairFile<ULong64_t, float>(fOdirPrefix+"Tadjust.fdat", false, true, 1000);
}
#endif
#ifdef PSA_LOCALSPECTRA
specEner = new nDhist<unsigned int>(2, 40, specLenE);
specEner->setFileName(fOdirPrefix+"Psa?Ener.spec");
......@@ -608,8 +626,8 @@ Int_t PSAFilterGridSearch::Process(int slot, int nslots)
// a fake PSA that places the inteaction points at the center of their net-charge segment
for(int ntot = 0; ntot < nslots; ntot++) {
int indexDD = slot + ntot;
PsaData *pD = DD + indexDD; // all info from/to the caller is contained in this
PrepareEvent(pD, &pS, false); // reset pS, identify the hits and prepare the traces
PsaData *pD = DD + indexDD; // all info from/to the caller is contained in this
PrepareEvent(pD, &pS, true, false); // reset pS, identify the hits and prepare the traces
if(pS.isValid)
SetToSegCenter(pD, &pS); // hits are centered in their segments
else {
......@@ -622,74 +640,90 @@ Int_t PSAFilterGridSearch::Process(int slot, int nslots)
return 0; // nothing more to do
}
const int nrpt = 4;
const float dT0min = 0.2f;
const float shiftT0max = 4.00f;
for(int ntot = 0; ntot < nslots; ntot++) {
// 1111111111111111111111111111 prepare one event
int indexDD = slot + ntot;
PsaData *pD = DD + indexDD; // all info from/to the caller is contained in this
PrepareEvent(pD, &pS); // reset pS, identify the hits and prepare the traces
WRITEWORKINGWAVE((float *)pS.tAmplitude, true); // false --> also rewind the file
for(int rpt = 0; rpt < nrpt; rpt++) {
if(!pS.isValid) {
pD->crystal_status = 1; // not to be used (anyway there are no hits)
pD->numNetSegs = 0;
pD->isSelected = false;
continue;
}
// 1111111111111111111111111111 prepare one event
int indexDD = slot + ntot;
PsaData *pD = DD + indexDD; // all info from/to the caller is contained in this
PrepareEvent(pD, &pS, (rpt==0), true); // reset pS, identify the hits and prepare the traces
WRITEWORKINGWAVE((float *)pS.tAmplitude, true); // false --> also rewind the file
if(!pS.isValid) {
pD->crystal_status = 1; // not to be used (anyway there are no hits)
pD->numNetSegs = 0;
pD->isSelected = false;
continue;
}
int sMult = pS.numNetCharges;
pS.samp_first = pS.addr_first;
pS.samp_last = pS.addr_last;
pD->tCycles++;
int sMult = pS.numNetCharges;
pS.samp_first = pS.addr_first;
pS.samp_last = pS.addr_last;
#ifdef USE_SSE_VERSION
pS.samp_last = pS.samp_first + SFIXEDLEN - 1;
#else
if(SFIXEDLEN > 0) {
// search using a fixed number of samples
pS.samp_last = pS.samp_first + SFIXEDLEN - 1;
}
else {
// search up to a fraction of signal amplitude (of the core)
pS.samp_last = pS.addr_last;
pS.SearchLimits(0.95f, fHitSegThreshold, pS.samp_first, pS.samp_last);
}
#else
if(SFIXEDLEN > 0) {
// search using a fixed number of samples
pS.samp_last = pS.samp_first + SFIXEDLEN - 1;
}
else {
// search up to a fraction of signal amplitude (of the core)
pS.samp_last = pS.addr_last;
pS.SearchLimits(0.95f, fHitSegThreshold, pS.samp_first, pS.samp_last);
}
#endif
//if(sMult > 1) { // initial solution not calculated for M=1
//if(sMult > 1) { // initial solution not calculated for M=1
#ifdef INITATSEGCENTER
InitialSolution(pS); // everything placed at the center of the segment; solution points recorded in pS.resPt1
WRITEWORKINGWAVE((float *)pS.wAmplitude, true);
// The first time, everything is placed at the center of the segment;
// After shifting the input trace, we start from the previous solution
// Points recorded in pS.resPt1
InitialSolution(pS, (rpt!=0) );
WRITEWORKINGWAVE((float *)pS.wAmplitude, true);
#endif
#ifdef PRECOARSESEARCH
PreSearchCoarse(pS); // make also a Coarse-Only search
WRITEWORKINGWAVE((float *)pS.wAmplitude, true);
PreSearchCoarse(pS); // make also a Coarse-Only search
WRITEWORKINGWAVE((float *)pS.wAmplitude, true);
#endif
//}
//}
// Clear the "fitted" trace
memset(pS.rAmplitude, 0, sizeof(pS.rAmplitude));
// Clear the "fitted" trace
memset(pS.rAmplitude, 0, sizeof(pS.rAmplitude));
// 2222222222222222222222222222 call the grid-search algo for one event
int ncalls = ProcessTheEvent(pS);
WRITEWORKINGWAVE((float *)pS.rAmplitude, true);
// 2222222222222222222222222222 call the grid-search algo for one event
int ncalls = ProcessTheEvent(pS);
WRITEWORKINGWAVE((float *)pS.rAmplitude, true);
//TryToFit(pS);
//RandomSplit2(pS, 1000);
//TryToFit(pS);
//RandomSplit2(pS, 1000);
// 3333333333333333333333333333 move the search-results from pS to pD
pD->numNetSegs = pS.numNetCharges;
pD->isSelected = pS.isSelected;
for(int nn = 0; nn < pS.numNetCharges; nn++) {
memcpy(pD->PsaOut + nn, pS.PsaOut + nn, sizeof(PsaOut_t));
}
if(doFit_T0) {
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
if(fWriteNumTraces > 0) {
SaveTotalTrace(pS, pD);
// 3333333333333333333333333333 move the search-results from pS to pD
pD->numNetSegs = pS.numNetCharges;
pD->isSelected = pS.isSelected;
for(int nn = 0; nn < pS.numNetCharges; nn++) {
memcpy(pD->PsaOut + nn, pS.PsaOut + nn, sizeof(PsaOut_t));
}
float dT0 = 0;
if(doFit_T0) {
dT0 = FitT0AfterPSA(pS, pS.samp_first + pS.samp_used);
pD->shiftT0 += dT0; // Units of samples. Accumulated if different runs
}
// While writing the traces, save the results in the proper slot
if(fWriteNumTraces > 0) {
SaveTotalTrace(pS, pD);
}
pD->retValue = 0; // could be used as error message
if(abs(dT0)<dT0min || abs(pD->shiftT0)>shiftT0max)
break;
}
pD->retValue = 0; // could be used as error message
} // for(int ntot = 0; ntot < nslots; ntot++)
......@@ -946,15 +980,16 @@ int PSAFilterGridSearch::ProcessTwoHits(pointFull &pS)
return ncalls;
}
void PSAFilterGridSearch::PrepareEvent(PsaData *pD, pointFull *pS, bool traces)
void PSAFilterGridSearch::PrepareEvent(PsaData *pD, pointFull *pS, bool reset, bool traces)
{
if(pD->crystal_status) {
pS->isValid = false; // not a good event (is this still possible ?)
return;
}
pS->Reset(); // clear everything
pS->enerGe = pD->CoreE[0]; // copy Core energy
if(reset)
pS->Reset(); // clear everything
pS->enerGe = pD->CoreE[0]; // copy Core energy
#ifdef SELECTEVENTS
pS->isSelected = false;
......@@ -1006,20 +1041,75 @@ void PSAFilterGridSearch::PrepareEvent(PsaData *pD, pointFull *pS, bool traces)
// Now fill tAmplitude with up to NTIME samples (including the pretrigger)
float *shifted = NULL;
float shiftT0 = -pD->shiftT0;
if(shiftT0)
shifted = new float[NSAMP];
int ntoskip = defTriggerSample - SAMP0; // alignment to a multiple of 16 bytes ( SAMP0=(defTriggerSample/4)*4 )
int ntocopy = (NSAMP-ntoskip) > NTIME ? NTIME : NSAMP-ntoskip;
pS->addr_first = SAMP0; // this is the nominal position of T0 used in the PSA
pS->addr_last = ntocopy-1;
for(Int_t atseg=0; atseg<NSEGS; atseg++) {
memcpy(pS->tAmplitude[atseg], &pD->fTracesSG[atseg][ntoskip], ntocopy*sizeof(float));
if(shiftT0) {
ShiftMoveTrace(pD->fTracesSG[atseg], NSAMP, shiftT0, shifted);
memcpy(pS->tAmplitude[atseg], &shifted[ntoskip], ntocopy*sizeof(float));
}
else {
memcpy(pS->tAmplitude[atseg], &pD->fTracesSG[atseg][ntoskip], ntocopy*sizeof(float));
}
if(NTIME>ntocopy) Extend(pS->tAmplitude[atseg], ntocopy, NTIME);
}
// copy also the CC
memcpy(pS->tAmplitude[NSEGS], &pD->fTracesCC[0][ntoskip], ntocopy*sizeof(float));
if(shiftT0) {
ShiftMoveTrace(pD->fTracesCC[0], NSAMP, shiftT0, shifted);
memcpy(pS->tAmplitude[NSEGS], &shifted[ntoskip], ntocopy*sizeof(float));
}
else {
memcpy(pS->tAmplitude[NSEGS], &pD->fTracesCC[0][ntoskip], ntocopy*sizeof(float));
}
if(NTIME>ntocopy) Extend(pS->tAmplitude[NSEGS], ntocopy, NTIME);
// save tAmplitude into wAmplitude, where is can be dynamically modified by the search algorithm
memcpy(pS->wAmplitude, pS->tAmplitude, sizeof(pS->wAmplitude));
if(shiftT0)
delete [] shifted;
}
void PSAFilterGridSearch::ShiftMoveTrace(float *data, int nsamp, float delta, float *accu)
{
if(delta == 0) {
for(int nn = 0; nn < nsamp; nn++) {
accu[nn] = data[nn];
}
return;
}
if(delta > 0) {
int idelta = int(delta);
float fdelta = delta-idelta;
float gdelta = 1.f - fdelta;
float first = data[0]*fdelta + data[1]*gdelta;
for(int nn = 0; nn <= idelta; nn++) {
accu[nn] = first;
}
for(int jj = 0, nn = idelta+1; nn < nsamp; jj++, nn++) {
accu[nn] = data[jj]*fdelta + data[jj+1]*gdelta;
}
}
else {
int idelta = -int(delta);
float gdelta = abs(delta)-idelta;
float fdelta = 1.f - gdelta;
int nn = 0;
for(int jj = idelta; jj < nsamp-1; jj++, nn++) {
accu[nn] = data[jj]*fdelta + data[jj+1]*gdelta;
}
float last = data[nsamp-2]*fdelta + data[nsamp-1]*gdelta;
for(; nn < nsamp; nn++) {
accu[nn] = last;
}
}
}
void PSAFilterGridSearch::SetToSegCenter(PsaData *pD, pointFull *pS)
......@@ -1900,15 +1990,16 @@ void PSAFilterGridSearch::MakeLocalMask(pointFull &pS)
localMask[pS.netSeg] = '1'; // (ri)enable the segment-of-interest
}
// place everything at the center
void PSAFilterGridSearch::InitialSolution(pointFull &pS)
// place everything at the center or keep what already found in the previous minimization
void PSAFilterGridSearch::InitialSolution(pointFull &pS, bool keep)
{
for(int snum = 0; snum < pS.numNetCharges; snum++) {
int netChSeg = pS.netChargeSegnum[snum];
float netChEner = pS.netChargeEnergy[snum];
float scaleFact = netChEner / MAXNORM;
int bestPt = fBasis.cflist[netChSeg].indFine[0];
int bestPt = (keep && pS.isInitialized) ? pS.resPt1[snum] : fBasis.cflist[netChSeg].indFine[0];
pS.resPt1[snum] = bestPt;
pS.resPt2[snum] = -1;
pS.resFac[snum] = 1.f;
......@@ -2039,6 +2130,11 @@ Int_t PSAFilterGridSearch::PostProcess(int slot)
}
}
#ifdef TSHIFT_CORRECTION
if(doFit_T0 && tstmpFileO)
tstmpFileO->put(pD->timestamp, pD->shiftT0);
#endif
#ifdef PSA_LOCALSPECTRA
if(pD->isSelected) {
PsaOut_t * pOut = pD->PsaOut;
......@@ -2233,14 +2329,15 @@ 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
int bestT1 = int((pD->CoreT[0]+pD->shiftT0)*SAMP_STEP);
specTzero->incr(NSEGS+1, bestT1); // corrected distribution
int bestT2 = int((pD->CoreT[0]-pD->shiftT0)*SAMP_STEP);
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
int bestT0 = int( pD->CoreT[0] *SAMP_STEP); // original distribution
specTzero->incr(NSEGS, bestT0);
int bestT1 = int((pD->CoreT[0]+pD->shiftT0)*SAMP_STEP); // corrected distribution
specTzero->incr(NSEGS+1, bestT1);
int bestT2 = int((pD->CoreT[0]-pD->shiftT0)*SAMP_STEP); // corrected distribution (with reversed sign)
specTzero->incr(NSEGS+2, bestT2);
int bestDT = int(pD->shiftT0*SAMP_STEP + specLenT/2); // the correction
specTzero->incr(NSEGS+3, bestDT);
specTzero->incr(NSEGS+3, pD->tCycles); // number of search-tShift cycles
}
if(specStat)
specStat->Incr(0, numNetSegs); // distribution of number of hit segments
......@@ -2516,16 +2613,15 @@ Float_t PSAFilterGridSearch::FitT0AfterPSA(pointFull &pS, int tsamp)
}
}
data = tData;
rfit = tRfit;
rfit = tRfit + SAMP0/2;
const int nsteps = 3;
float x2[2*nsteps+1]; memset(x2, 0, sizeof(x2));
for(int ll = 0; ll <= 2*nsteps; ll++) {
float *tfit = rfit + nsteps - ll;
float *data = tData + SAMP0/2 - nsteps + ll;
float x2t = 0;
// 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]);
// calculate chi2 using SAMP0/2 of pretrigger and 20 of risetime
for(int ii = nsteps; ii < SAMP0/2+20; ii++) {
float diff = float(data[ii] - rfit[ii]);
x2t += diff*diff;
}
x2[ll] = x2t;
......
......@@ -8,6 +8,12 @@
//
//#include "commonDefs.h"
//#define TSHIFT_CORRECTION
#ifdef TSHIFT_CORRECTION
#include "dataPairFile.h"
#endif
#include "PSAFilter.h"
#include "GridSearchParams.h"
#include "SignalBasis.h"
......@@ -98,6 +104,10 @@ private:
nDhist<unsigned short> *matrRZE;
#endif //PSA_LOCALSPECTRA
#ifdef TSHIFT_CORRECTION
dataPairFile<ULong64_t, float> *tstmpFileO;
#endif
std::string fnPsaTraces;
FILE *fpPsaTraces;
......@@ -117,11 +127,12 @@ private:
void SegmentMapMake (int neighbours); // based on Manhattan distance
void SegmentMapTest (); // exclude elements related to broken segments
void PrepareEvent (PsaData *pD, pointFull *pS, bool traces = true);
void PrepareEvent (PsaData *pD, pointFull *pS, bool reset, bool traces);
void ShiftMoveTrace (float *data, int nsamp, float delta, float *accu);
void SetToSegCenter (PsaData *pD, pointFull *pS);
int ProcessTheEvent(pointFull &pS);
int ProcessTwoHits (pointFull &pS);
void InitialSolution(pointFull &pS);
void InitialSolution(pointFull &pS, bool keep = false);
void PreSearchCoarse(pointFull &pS);
int SearchFullGrid (pointFull *pS, int netChSeg);
int SearchAdaptive1(pointFull *pS, int netChSeg, bool bCoarseOnly = false, bool bFirstOnly = false); // 1 hit
......
......@@ -87,6 +87,13 @@ PreprocessingFilterPSA::PreprocessingFilterPSA() :
fTC = fTS = fTT = fWW = NULL;
fTFA = fCFD = NULL;
#ifdef TSHIFT_CORRECTION
tstmpFileI = NULL;
fTimestampOld = 0;
fValueOld = 0;
#endif
#ifdef PPF_FromGRU_
CCE_raw = new TH1F("CCE_raw","CCE_raw",32768,0,32767);
......@@ -208,6 +215,11 @@ Int_t PreprocessingFilterPSA::AlgoSpecificInitialise()
mwd->nMinWidthN = 6;
#endif
#ifdef TSHIFT_CORRECTION
tstmpFileI = new dataPairFile<ULong64_t, float>(GetConfPath()+"Tadjust.fdat", false, false, 1000);
#endif
// the preset value should come from the configuration file
tAverage.Initialize(30.f);
......@@ -1279,6 +1291,8 @@ bool PreprocessingFilterPSA::ProcessEvent(int &smult, float &eSumSG1, float &e
// fHandShift is an ad-hoc common shift to have the traces well positioned in the PSA.
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;
//float extraShift = getTstampFileValue(timestamp);
//fshiftCC += extraShift;
for(int nc = 0; nc < nCC; nc++) {
ShiftMoveTrace(fTracesCC[nc], nNsamp, fshiftCC, tmp); // no extra shift for the cores
memcpy(fTracesCC[nc], tmp, sizeof(tmp));
......
......@@ -21,6 +21,12 @@
# include "mwdlib.h"
#endif
//#define TSHIFT_CORRECTION
#ifdef TSHIFT_CORRECTION
#include "dataPairFile.h"
#endif
#ifdef PPF_FromGRU_
// #include <TAttFill.h>
#include <TROOT.h>
......@@ -172,6 +178,12 @@ public:
nDhist<unsigned int > *specTT2;
#endif //PPF_LOCALSPECTRA
#ifdef TSHIFT_CORRECTION
dataPairFile<ULong64_t, float> *tstmpFileI;
ULong64_t fTimestampOld;
float fValueOld;
#endif
// various working buffers
float *fCC; // for the core
float *fSG; // for the segments
......@@ -203,6 +215,7 @@ public:
int xTalkCorrNet (float *pSG_Ener, float &segsum);
int xTalkCorrAll (float *pSG_Ener, float &segsum);
int xTalkTraces ();
float getTstampFileValue(ULong64_t &tst);
#ifdef PPF_LOCALSPECTRA
void calcEEmatrix (float *eSg, float *eCc, float **sSg, float **sCc);
......@@ -238,4 +251,23 @@ inline float PreprocessingFilterPSA::AverageVariance(float *data, int nchan, flo
return vv/nchan;
}
#ifdef TSHIFT_CORRECTION
inline float PreprocessingFilterPSA::getTstampFileValue(ULong64_t &tst)
{
if(!tstmpFileI || !tstmpFileI->status)
return 0;
if(tst == fTimestampOld)
return fValueOld;
while(tst > fTimestampOld) {
tstmpFileI->get(fTimestampOld, fValueOld);
}
if(tst == fTimestampOld)
return fValueOld;
else
return 0;
}
#endif
#endif // PREPROCESSINGFILTERPSA_H_INCLUDED
......@@ -46,8 +46,9 @@ fTrigger("data:crystal")
fNumFilesIn = 1;
fMinSegAmpli = 500.f; // ~50 (150) keV for high(low) gain range if shaping = 1000
fProjeM1 = false;
fSpekGain = 1.f; // good for trapezoid risetime 1000
fMinSegAmpli = 500.f; // ~50 (150) keV for high(low) gain range if shaping = 1000
fInputDataFile.clear();
......@@ -226,7 +227,14 @@ void CrystalProducer::GetParameters( UInt_t *error_code )
else if( stringEq(keyw, "WriteBaseLines") ) {
ok = 1 == sscanf(data.c_str(), "%d", &fWriteBaselines);
}
else if( stringEq(keyw, "AmplitudeGain") ) {
else if( stringEq(keyw, "ProjeM1") ) {
ok = 2 == sscanf(data.c_str(), "%f %f", &fMinSegAmpli, &fSpekGain);
}
else if( stringEq(keyw, "ProjeThreshold") ) {
ok = 1 == sscanf(data.c_str(), "%f", &fMinSegAmpli);
fProjeM1 = true;
}
else if( stringEq(keyw, "ProjeGain") ) {
ok = 1 == sscanf(data.c_str(), "%f", &fSpekGain);
}
else if( stringEq(keyw, "TimeStep") ) {
......
......@@ -74,6 +74,7 @@ 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
Bool_t fProjeM1; // extend the total projection file to subspectra with segMult=1 for the segments and the core
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
......
......@@ -154,7 +154,7 @@ Int_t CrystalProducerATCA::AlgoSpecificInitialise()
}
#ifdef CP_LOCALSPECTRA
ProdSpec_Ampl = new nDhist<unsigned int>(3, CrystalInterface::kNbSegments + CrystalInterface::kNbCores, specLenE);
ProdSpec_Ampl = new nDhist<unsigned int>((fProjeM1)?3:1, CrystalInterface::kNbSegments + CrystalInterface::kNbCores, specLenE);
ProdSpec_Ampl->setFileName(fOdirPrefix+"Prod?Ampli.spec");
ProdSpec_Ampl->setComment("amplitude spectra of segments and cores");
hGroup.add(ProdSpec_Ampl);
......@@ -771,14 +771,23 @@ Int_t CrystalProducerATCA::Process()
if(ProdSpec_Ampl) {
int segMult = 0;
int segLast = 0;
for(int nn = 0; nn < CrystalInterface::kNbSegments; nn++) {
float fAmpli = SegE[nn]*fSpekGain;
if(fAmpli > fMinSegAmpli) {
segLast = nn;
segMult++;
if(!fProjeM1) {
for(int nn = 0; nn < CrystalInterface::kNbSegments; nn++) {
float fAmpli = SegE[nn]*fSpekGain;
int iampli = int(fAmpli + specOffE);
ProdSpec_Ampl->Incr(0, nn, iampli);
}
}
else {
for(int nn = 0; nn < CrystalInterface::kNbSegments; nn++) {
float fAmpli = SegE[nn]*fSpekGain;
if