Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

Commit 07d378a0 authored by dino's avatar dino
Browse files

Exposing more parameters to control the PSAFilterGridSearch with 2 hits in a...

Exposing more parameters to control the PSAFilterGridSearch with 2 hits in a segment. Several minor fixes.

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@1170 170316e4-aea8-4b27-aad4-0380ec0519c9
parent f1c4b66f
......@@ -28,7 +28,7 @@ using namespace ADF;
// would freeze the analysis, must be modified in favor of a dynamic mechanism. That means that if the
// DISPATCHER cannot send the data because the exchange buffer is owned by the MERGER, the chain sets a
// flag and returns (allowing the other chains to be run in turn). At the next run of this chain the
// flag is tested and the Dispatcher is called again and this mechanism is repeated until successfull
// flag is tested and the Dispatcher is called again; this mechanism is repeated until successful
// exchange and only then the chain executes the main body.
void ChainOfActors::Run()
......@@ -52,7 +52,7 @@ void ChainOfActors::Run()
// check EOD condition
if(hasEnded) {
DispatchEOD(); // if succesful, sets isClosed
DispatchEOD(); // if successful, sets isClosed
return; // anyway, nothing more to do
}
......
......@@ -82,7 +82,7 @@ DISPATCHER EventBuilder pseudo actor to connect the output of the ch
#include "PSAFilter.h" // 31 310
#include "PSAFilterGridSearch.h" // 311
#include "AncillaryFilter.h" // 33 330 must find out why this cannot be after TrackingFilterMGT.h
#include "AncillaryFilter.h" // 33 330 must find out why this cannot be included after TrackingFilterMGT.h
#include "AncillaryFilterVME.h" // 331
#include "AncillaryFilterATCA.h" // 332
......@@ -102,17 +102,17 @@ DISPATCHER EventBuilder pseudo actor to connect the output of the ch
// BuilderName // 50 destination_chain
#ifdef WCT_THREADED
boost::thread **chThread; // to record the pointers of the threads (needed to be able to join them)
int detThreads = 0; // counter of PRODUCER chains running in parallel
bool stopControl = false; // to force the control thread to exit
bool runThreaded = true; // set it to false (from the command line) to run the unthreaded main loop
vector <boost::thread *> chThread; // pointers of the threads, needed to be able to join them
boost::thread *ctrlThread; // pointer of the control thread for cyclic actions
bool stopControl = false; // to force the control thread to exit
bool runThreaded = true; // set it to false (from the command line) to run the unthreaded main loop
#else
bool runThreaded = false; // will of course run unthreaded
bool runThreaded = false; // will of course run unthreaded
#endif
#if defined(COUT_LOCKED)
boost::mutex ioMutex; // mutex for LOCK_COUT
boost::mutex ioMutex; // mutex for LOCK_COUT
#endif
// Replacement of Xavier's variable to ask producers to return from process_block even
......@@ -131,8 +131,8 @@ std::string confPrefix("Conf/");
std::string StopFile("StopFile"); // presence of this file in cwd ends the analysis in an ordered way
std::string KillFile("KillFile"); // presence of this file in cwd ends the analysis immediately
int maxTurn = 0; // if using threads this limit works only approximately in eventloop2()
int theTurn = 0; // counter of turns
int maxTurn = 0; // if using threads, the analysis will run for approximately maxTurn seconds
int theTurn = 0; // counter of turns
int nSIGINT = 0; // count number of SIGINT (CTRL_C)
bool stopSig = false; // used to check StopFile
bool killSig = false; // used to check KillFile
......@@ -1757,20 +1757,16 @@ int eventLoop2()
int TOTCHAINS = (int)theChains.size();
chThread = new boost::thread* [TOTCHAINS];
detThreads = 0;
// launch the threads of all worker chains
for(int ich = 0; ich < TOTCHAINS ; ich++) {
chThread[ich] = new boost::thread(&threadWorker, ich);
chThread.push_back( new boost::thread(&threadWorker, ich) );
cout << "STARTED CHAIN " << ich << endl;
detThreads++;
}
cout << endl;
// launch another thread for cyclic actions like checking the presence of StopFile or KillFile
// launch a thread for cyclic actions, like checking the presence of StopFile or KillFile
stopControl = false;
boost::thread *ctrlThread = new boost::thread(&threadControl);
ctrlThread = new boost::thread(&threadControl);
// join() should be replaced with timed_join(...) to avoid blocking forever in case of
// never-exiting threads. The control thread could take up the job of checking this case
......@@ -1779,7 +1775,6 @@ int eventLoop2()
for(int ich = 0; ich < TOTCHAINS ; ich++) {
chThread[ich]->join();
cout << "JOINED CHAIN " << ich << endl;
detThreads--;
}
stopControl = true; // force also the control thread to exit if that has not already happened
......@@ -1833,7 +1828,7 @@ void threadControl()
else {
// In Unix the ^C is not captured by catch_int() as defined in main() <== removed
// This functionality is replaced by checking the presence of a file named
// StopFile (for ordered close-up) or KillFile for immediate stop
// StopFile (for ordered close-up) or KillFile (for immediate stop)
if(CheckFile(StopFile, true)) {
stopSig = true;
}
......
......@@ -45,7 +45,7 @@
Name="VCCLCompilerTool"
Optimization="0"
EnableIntrinsicFunctions="true"
AdditionalIncludeDirectories="..\WinCtest;..\common;..\producers\Crystal;..\producers\Crystal\includeATCA;..\producers\AncillaryTCP;..\filters\Preprocessing;..\filters\Ancillary;..\filters\Ancillary\includeVME;..\filters\PSA;..\filters\PSA\includePSA;..\filters\Tracking;..\filters\Tracking\includeOFT;..\filters\Tracking\includeMGT;..\filters\PostPSA;..\builders;..\filters\Global;..\..\include;..\PRISMA\src\lib_prisma\include;&quot;C:\Program Files (x86)\boost\boost_1_40&quot;;C:\root\include"
AdditionalIncludeDirectories="..\WinCtest;..\common;..\producers\Crystal;..\producers\Crystal\includeATCA;..\producers\AncillaryTCP;..\filters\Preprocessing;..\filters\Ancillary;..\filters\Ancillary\includeVME;..\filters\PSA;..\filters\PSA\includePSA;..\filters\Tracking;..\filters\Tracking\includeOFT;..\filters\Tracking\includeMGT;..\filters\PostPSA;..\builders;..\filters\Global;..\..\include;..\PRISMA\src\lib_prisma\include;C:\Boost\boost_1_55_0;C:\root\include"
PreprocessorDefinitions="WIN32;_DEBUG;_CONSOLE;_CRT_SECURE_NO_WARNINGS;NRV_TYPE=NRV_OFFLINE"
MinimalRebuild="true"
BasicRuntimeChecks="3"
......@@ -70,7 +70,7 @@
Name="VCLinkerTool"
UseLibraryDependencyInputs="true"
LinkIncremental="2"
AdditionalLibraryDirectories="&quot;C:\Program Files (x86)\boost\boost_1_40\lib&quot;"
AdditionalLibraryDirectories="&quot;C:\Boost\boost_1_55_0\lib32-msvc-9.0&quot;"
GenerateDebugInformation="true"
SubSystem="1"
TargetMachine="1"
......@@ -282,7 +282,7 @@
Optimization="3"
EnableIntrinsicFunctions="true"
FavorSizeOrSpeed="1"
AdditionalIncludeDirectories="..\PRISMA\src\lib_prisma\include;&quot;C:\Program Files (x86)\boost\boost_1_40&quot;;..\common;..\producers\Crystal;..\producers\Crystal\includeATCA;..\producers\AncillaryTCP;..\filters\Preprocessing;..\filters\Preprocessing\includePreprocessing;..\filters\Ancillary;..\filters\Ancillary\includeVME;..\filters\PSA;..\filters\PSA\includePSA;..\filters\Tracking;..\filters\Tracking\includeOFT;..\filters\Tracking\includeMGT;..\builders;C:\root\include;..\filters\PostPSA;..\filters\Global;..\..\include"
AdditionalIncludeDirectories="..\PRISMA\src\lib_prisma\include;C:\Boost\boost_1_55_0;..\common;..\producers\Crystal;..\producers\Crystal\includeATCA;..\producers\AncillaryTCP;..\filters\Preprocessing;..\filters\Preprocessing\includePreprocessing;..\filters\Ancillary;..\filters\Ancillary\includeVME;..\filters\PSA;..\filters\PSA\includePSA;..\filters\Tracking;..\filters\Tracking\includeOFT;..\filters\Tracking\includeMGT;..\builders;C:\root\include;..\filters\PostPSA;..\filters\Global;..\..\include"
PreprocessorDefinitions="WIN64;NDEBUG;_CONSOLE;_CRT_SECURE_NO_WARNINGS;NRV_TYPE=NRV_OFFLINE"
ExceptionHandling="1"
RuntimeLibrary="2"
......@@ -308,7 +308,7 @@
UseLibraryDependencyInputs="false"
AdditionalDependencies="libCore.lib libTree.lib libRIO.lib"
LinkIncremental="1"
AdditionalLibraryDirectories="&quot;C:\Boost\boost_1_40\lib64-VC90&quot;;C:\root\lib"
AdditionalLibraryDirectories="&quot;C:\Boost\boost_1_55_0\lib64-msvc-9.0&quot;;C:\root\lib"
GenerateDebugInformation="true"
SubSystem="1"
OptimizeReferences="2"
......
......@@ -52,6 +52,7 @@ bool adetParams::ReadCalibCoeffs(std::string setupFile, bool verbose)
"segm",
"trig", // files before end-january 2011
"tntf" // files after end-january 2011
// no extra key == tntf with trace_gain=energy_gain
};
string line, name, data, data_tntf;
......@@ -87,34 +88,38 @@ bool adetParams::ReadCalibCoeffs(std::string setupFile, bool verbose)
}
int fileVersion = 0;
if(nn_tntf==0) {
// the old format has no keywork tntf
if (nn_trig && nn_tntf) {
cout << "keywords trig and tntf are mutually exclusive" << endl;
return false;
}
else if (nn_trig) {
// the old format has trig and no tntf
fileVersion = 1;
if(nn_trig!=1) {
cout << "there must at most 1 trig statement" << endl;
cout << "there can be at most 1 trig statement" << endl;
return false;
}
fileVersion = 1;
}
else {
// the new format, with the keywork tntf
if(nn_tntf!=1) {
fileVersion = 2;
if(nn_tntf > 1) {
cout << "there can be at most 1 tntf statement" << endl;
return false;
}
if(nn_trig && nn_trig!=1) {
cout << "there can be at most 1 trig statement" << endl;
return false;
}
fileVersion = 2;
istringstream dd(data_tntf);
dd >> tntFactor;
if(!dd.good() && (dd.rdstate()&ios_base::badbit)) {
cout << "Error reading the tntf factor from { " << line << " }" << endl;
return false;
if(nn_tntf) {
istringstream dd(data_tntf);
dd >> tntFactor;
if(!dd.good() && (dd.rdstate()&ios_base::badbit)) {
cout << "Error reading the tntf factor from { " << line << " }" << endl;
return false;
}
//if(tntFactor <=0) {
// cout << "The scaling factor for the traces must be positive --> " << tntFactor << endl;
// return false;
//}
}
if(tntFactor <=0) {
cout << "The scaling factor for the traces must be positive --> " << tntFactor << endl;
return false;
else {
tntFactor = 0; // traces have the same gain as energy (new AGATA/GALILEO electronics)
}
}
......@@ -273,7 +278,7 @@ bool adetParams::ReadCalibV02(std::ifstream &cfg, std::string *keys, bool verbos
pCC[corenum].tfall = t1;
pCC[corenum].trise = t2;
pCC[corenum].egain = gE;
pCC[corenum].sgain = gE*t1*t2/tntFactor;
pCC[corenum].sgain = (tntFactor>0) ? gE*t1*t2/tntFactor : gE;
pCC[corenum].tmove = gT;
//pCC[corenum].tfast = tf;
//pCC[corenum].lambdaH = lH;
......@@ -302,7 +307,7 @@ bool adetParams::ReadCalibV02(std::ifstream &cfg, std::string *keys, bool verbos
pSG[segnum].tfall = t1;
pSG[segnum].trise = t2;
pSG[segnum].egain = gE;
pSG[segnum].sgain = gE*t1*t2/tntFactor;
pCC[segnum].sgain = (tntFactor>0) ? gE*t1*t2/tntFactor : gE;
pSG[segnum].emink = gThr;
pSG[segnum].tmove = gT;
//pCC[segnum].tfast = tf;
......
......@@ -69,7 +69,7 @@ public:
bool bSegEmin;
float vSegEmin;
float tntFactor; // usually 2^21 = 2097152
float tntFactor; // usually 2^21 = 2097152; if <=0 traces get the same gain as energies
void SetTriggerLevel(float value) {vtrig = value;}
void SetSegMinEnergy(float value) {vSegEmin = value; bSegEmin = true;}
......
......@@ -327,7 +327,7 @@ bool StrightLine(float *data, int nchan, float &A, float&B, float ref)
yx1 += y*x;
}
// the [x] components and the determinant censtant and could/should be precalculated
// the [x] components and the determinant are constant and could/should be precalculated
// N => 1 2 3 4 5 6 7 8 9 10
// x => 0 1 2 3 4 5 6 7 8 9
// sx0 => 1 2 3 4 5 6 7 8 9 10 == N
......
......@@ -24,14 +24,14 @@ const int specLenT = 1000;
const int matLen = 100; // 100 mm
const int matOff = matLen/2; // origin in the middle
const int matPack = 2; // pixels are 2 mm^3
const int numXYZ = 10; // sub PsaMatrXYZR for up to numXYZ-1 hits
const int numXYZ = 10; // sub PsaMatrXYZR for up to numXYZ-1 hits
const size_t nPsaHitsValues = 16; // number of parameters written in Psa__0-nn-F__Hits.fdat
const int defLoopsTZ = 4; // default number of GridSearch--FitT0AfterPSA loops (0 == no T0 fit)
const int maxLoopsTZ = 10; // maximum number of GridSearch--FitT0AfterPSA loops
const float shiftDTmin = 0.15f; // samples; stop looping if dT0 smaller than this
const float shiftT0max = 5.00f; // samples; stop looping if TotalShift larger than this
const int defLoopsTZ = 4; // default number of GridSearch--FitT0AfterPSA loops (0 == no T0 fit)
const int maxLoopsTZ = 10; // maximum number of GridSearch--FitT0AfterPSA loops
const float shiftDTmin = 0.15f; // samples; stop looping if dT0 smaller than this
const float shiftT0max = 5.00f; // samples; stop looping if TotalShift larger than this
const float tauDefault = 35.f; // default time constant of the preamp response function for segments and core (in nanoseconds)
......@@ -89,6 +89,13 @@ PSAFilterGridSearch::PSAFilterGridSearch()
fHitSegThreshold = minSegEnergy; // keV
h2Chi2Start = 2.00f; // Start higher than the 1-hit
h2Chi2Final = 0.95f; // Minimum final improvement
h2EminRel = 0.05f; // minimum relative energy of the light component
h2EminAbs = 50.0f; // minimum absolute energy of the light component (keV)
h2DminCoarse = 13.0f; // minimum distance of coarse grid points for the two hits solution
h2DminFine = 7.0f; // minimum distance of fine grid points for the two hits solution
memset(hmask, ' ', sizeof(hmask));
for(int nn = 0; nn < NSEGS; nn++) { // segments and core
hmask[nn][INDCC+1] = 0; // +1 to close each line as a string
......@@ -327,13 +334,14 @@ Int_t PSAFilterGridSearch::AlgoInitialise()
cout << "Grid Search performed using a ONLY the COARSE grid with a step of " << cstep << " mm" << endl;
}
else if(gUseAdaptive) {
cout << "Grid Search performed using a COARSE-FINE grid with steps of " << cstep << " and " << fstep << " mm";
cout << "Grid Search performed using a COARSE-FINE grid with steps of " << cstep << " and " << fstep << " mm" << endl;
if(gTryTwoHits) {
cout << " ==> try also USING TWO HITS";
cout << " ==> also trying with TWO HITS, with limits on "
<< " chi2{" << h2Chi2Start << " " << h2Chi2Final << "}"
<< " ener{" << h2EminRel << " " << h2EminAbs << "}"
<< " dist{" << h2DminCoarse << " " << h2DminFine << "}" << endl;
}
cout << endl;
}
cout << "Using the trace of the Core = " << ((gUseTraceCC)?"true":"false") << endl;
cout << "Using the Net-Charge segments = " << ((gUseTraceNC)?"true":"false") << endl;
if((gUseTraceNC || gUseTraceCC) && gDelayDiff > 0)
......@@ -461,6 +469,18 @@ void PSAFilterGridSearch::AlgoParameters(std::string keyw, std::string data, Boo
ok = true;
GP_END;
GP_PAR("TwoHitsChi2", "%f %f", "coarse and final chi2 relative limits for the two-hits solution");
ok = 2 == sscanf(data.c_str(), "%f %f", &h2Chi2Start, &h2Chi2Final);
GP_END;
GP_PAR("TwoHitsEner", "%f %f", "minimum relative and absolute energy of the light component");
ok = 2 == sscanf(data.c_str(), "%f %f", &h2EminRel, &h2EminAbs);
GP_END;
GP_PAR("TwoHitsDist", "%f %f", "minimum distance (mm) of coarse and fine grid points for the two hits solution");
ok = 2 == sscanf(data.c_str(), "%f %f", &h2DminCoarse, &h2DminFine);
GP_END;
GP_PAR("DistanceMetric", "%f", "metrics to calculate the figure of merit");
ok = 1 == sscanf(data.c_str(), "%f", &gDistMetric);
GP_END;
......@@ -613,12 +633,30 @@ int PSAFilterGridSearch::PreampResponse(bool exponential)
if(fTauSegment[ns] > 0)
tauSGCC[ns] = fTauSegment[ns];
}
}
fBasis.ExponentialResponse(tauSGCC, fTauDecay); // ns
}
bool smooth = true;
fBasis.ExponentialResponse(tauSGCC, fTauDecay, smooth); // ns
cout << " Decay constants (ns) of the exponential filter applied to the signal basis:";
//ios_base::fmtflags oldflags = cout.flags();
streamsize oldprec = cout.precision(1);
cout << fixed;
for(int slice = 0; slice < NSLIC; slice++) {
cout << "\n slice " << setw(1) << slice << " ";
for(int sector = 0; sector < NSECT; sector++) {
cout << setw(8) << tauSGCC[sector*NSLIC+slice];
}
}
cout << "\n core " << setw(8) << tauSGCC[NSEGS] << endl;
if(smooth)
cout << " Signals basis smoothed with an 1-8-1 kernel" << endl;
cout.precision(oldprec);
cout.unsetf(ios_base::floatfield);
//cout.flags(oldflags);
}
else {
// alternatively, convolution with the experimental response derived from the pulser (kernel stored in fBasis)
fBasis.ExperimentalResponse();
cout << "Signal basis modified using a hard-wired experimental response" << endl;
}
return 0;
}
......@@ -728,13 +766,13 @@ int PSAFilterGridSearch::PrepareMetric()
Int_t PSAFilterGridSearch::Process(int firstSlot, int evCount)
{
// To allow parallel processing, all info&data for the present event are kept into the stack-object PF.
// To allow parallel processing, all info&data for the present event are kept into the object PF.
// After being filled by PrepareEvent(), PF is modified by the analysis methods (to which it is usually passed by reference).
// Finally the results are moved back from PF into the proper data slot.
pointFull PF;
if(gSegCenter) {
// a fake PSA that places the inteaction points at the center of their net-charge segment
// a fake PSA that places the interaction points at the center of their net-charge segment
for(int evn = 0; evn < evCount; evn++) { // loop on the events
PsaSlot *pSlot = theSlots + firstSlot + evn; // all info from/to the caller is contained here
pSlot->chi2Start = 0;
......@@ -798,8 +836,8 @@ Int_t PSAFilterGridSearch::Process(int firstSlot, int evCount)
if(!PF.isValid) {
pSlot->crystal_status = 1; // not to be used (anyway there are no hits)
pSlot->numNetSegs = 0;
pSlot->isSelected = false;
pSlot->numNetSegs = 0;
pSlot->isSelected = false;
break; // continue;
}
......@@ -839,15 +877,15 @@ Int_t PSAFilterGridSearch::Process(int firstSlot, int evCount)
pSlot->retValue = 0; // could be used as error message
if(extraLoop) {
//// Perform a final time adjustment ??
//// BETTER NOT as it produces a slight worsening of Tgg. Not yet verified for AdaptiveTwoHits
//// Perform a final time adjustment ?? BETTER NOT as it produces a slight
//// worsening of Tgg (this has not yet been checked for AdaptiveTwoHits)
//if(gFitTZero) {
// float DT = FitT0AfterPSA(PF);
// pSlot->shiftTZero += DT;
//}
// final total reduced chi2 from the residuals in wAmplitude
pSlot->chi2Final = PF.TotalResidualChi2();
break; // this is the only way to exit the "for(int rpt..." loop
break; // this is the only way (for valid events) to exit the "for(int rpt..." loop
}
// 4444444444444444444444444444 adjust T0
......@@ -985,7 +1023,7 @@ int PSAFilterGridSearch::ProcessTheEvent(pointFull &PF, bool doFull, bool twoHit
if(!twoHits)
return chiDone;
// attempts to fit (kept only as reminder)
// attempts to fit (kept only as reminders)
//TryToFit(PF);
//RandomSplit2(PF, 1000);
......@@ -1022,9 +1060,9 @@ int PSAFilterGridSearch::ProcessTwoHits(pointFull &PF)
int nCalled = 0;
int chiDone = 0;
// global analysis (now a fake one)
if(!IsItWorthwhile(PF))
return chiDone;
//// global analysis (now a fake one)
//if(!IsItWorthwhile(PF))
// return chiDone;
// make a full copy of the previous solution in case we decide to fall back to it
pointFull PF1;
......@@ -1160,7 +1198,7 @@ int PSAFilterGridSearch::ProcessTwoHits(pointFull &PF)
if(numSegTwoHits) {
float ch2Final = CalcChi2Residue(PF);
if(ch2Final > 0.95f*ch2Initial) {
if(ch2Final > h2Chi2Final*ch2Initial) {
PF.FullCopy(PF1); // WRITEPARTIALWAVE and WRITEWORKINGWAVE lost their meaning
return -chiDone;
}
......@@ -1872,9 +1910,8 @@ int PSAFilterGridSearch::SearchFullGrid(pointFull &PF, int netChSeg, float netCh
//#define FULL_TWO_HITS
#ifdef FULL_TWO_HITS
// Lower limit on the energy of the partial hits
// (this is also a 2*ener2Min threshold on the total energy of the segment)
const float ener2Min = 50.f;
if(netChEner < 2*ener2Min) {
// (this is also a 2*h2EminAbs threshold on the total energy of the segment)
if(netChEner < 2*h2EminAbs) {
return 0; // --> don't proceed with the fine search
}
// Constraints on the distance between the two points
......@@ -1893,10 +1930,10 @@ int PSAFilterGridSearch::SearchFullGrid(pointFull &PF, int netChSeg, float netCh
float chi2Opt = 0;
float ePartOpt = EnergyPartition(PF, netChSeg, iPts1, iPts2, chi2Opt);
// skip off range or too extreme energy partition
if(ePartOpt < 0.05f || ePartOpt > 0.95f)
if(ePartOpt < h2EminRel || ePartOpt > (1.f-h2EminRel))
continue;
// skip solutions with too small energy (i.e. more balanced partition for low-energy events)
if(netChEner*ePartOpt < ener2Min || netChEner*(1.f-ePartOpt) < ener2Min) // keV
if(netChEner*ePartOpt < h2EminAbs || netChEner*(1.f-ePartOpt) < h2EminAbs) // keV
continue;
#if FIXED_METRIC == FIXED_METRIC_SQUARE
//// if the metric is 2, chiSquare has been calculated in EnergyPartition
......@@ -2141,22 +2178,17 @@ int PSAFilterGridSearch::SearchAdaptive1(pointFull &PF, int netChSeg, bool bCoar
}
// Coarse search performed looping over all ORDERED pairs of points, with the the two hits sharing the energy in the optimum way.
// This solution is then refined using the fine grids around the coarse solution
// This solution is then refined using the fine-grid points around the coarse solution
int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoarseOnly, bool bFirstOnly)
{
// Lower limit on the energy of the partial hits
// (this is also a 2*ener2Min threshold on the total energy of the segment)
const float ener2Min = 50.f;
if(PF.enerSeg < 2*ener2Min) {
// (this is also a 2*h2EminAbs threshold on the total energy of the segment)
if(PF.enerSeg < 2*h2EminAbs) {
return 0; // --> don't proceed with the fine search
}
// Constraints on the distance between the two points
const float ddMin = 13.f; // mm
const float ffMin = 7.f; // mm
const float nnMin = ddMin+ffMin;
const float chi2FactorStart = 2.00f; // Start higher than the 1-hit
const float chi2FactorFinal = 0.95f; // Minimum final improvement
//const float ddMin = 13.f; // mm
//const float ffMin = 7.f; // mm
int time_first = TIME0; // start of basis
int samp_first = PF.samp_first;
......@@ -2168,10 +2200,10 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
// the coarse-fine structure of this segment
SignalBasis::cfgrid &sgrid = fBasis.cflist[netChSeg];
int h1_bestPt = PF.bestPt1; // remember previous result in case it is decided to reuse it
int h1_bestPt = PF.bestPt1; // remember previous result in case it is decided to reuse it
float h1_chi2min = PF.chi2min;
float chi2min = h1_chi2min*chi2FactorStart; // Inital limit for chi2, respect to the 1-hit solution
float chi2min = h1_chi2min*h2Chi2Start; // Inital limit for chi2, respect to the 1-hit solution
int nInner = 0;
int chiDone = 0;
......@@ -2196,21 +2228,21 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
int kPts2 = sgrid.indFine[sgrid.posCoarse[kloop2]];
// avoid points too close to each other
float dd = fBasis.PtPtDistance(netChSeg, kPts1, kPts2);
if(dd < ddMin)
if(dd < h2DminCoarse)
continue;
float chi2Opt = 0;
float ePartOpt = EnergyPartition(PF, netChSeg, kPts1, kPts2, chi2Opt);
// skip off range or too extreme energy partition
if(ePartOpt < 0.05f || ePartOpt > 0.95f)
if(ePartOpt < h2EminRel || ePartOpt > (1.f-h2EminRel))
continue;
// skip solutions with too small energy (i.e. more balanced partition for low-energy events)
if(PF.enerSeg*ePartOpt < ener2Min || PF.enerSeg*(1.f-ePartOpt) < ener2Min) // keV
if(PF.enerSeg*ePartOpt < h2EminAbs || PF.enerSeg*(1.f-ePartOpt) < h2EminAbs) // keV
continue;
#if FIXED_METRIC == FIXED_METRIC_SQUARE
//// if the metric is 2, chiSquare has been calculated in EnergyPartition
chi2 = chi2Opt;
#else
// otherwise chiSquare must be calculated explicitely
// otherwise chiSquare is calculated explicitely
float bScale1 = PF.baseScale*RESCALE*ePartOpt;
float bScale2 = PF.baseScale*RESCALE*(1.f - ePartOpt);
chiDone++;
......@@ -2287,15 +2319,15 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
for(int jloop = jloop1; jloop < jloop2; jloop++) {
int jPts = sgrid.indFine[jloop];
float dd = fBasis.PtPtDistance(netChSeg, jPts, fPts); // order not relevant
if(dd < ffMin)
if(dd < h2DminFine)
continue;
float chi2Opt = 0;
float ePartOpt = EnergyPartition(PF, netChSeg, jPts, fPts, chi2Opt);
if(ipt == 1) ePartOpt = 1.f-ePartOpt;
if(ePartOpt < 0.05f || ePartOpt > 0.95f)
if(ePartOpt < h2EminRel || ePartOpt > (1.f-h2EminRel))
continue;
// skip solutions with too small energy (implies more balanced partition for low-energy events)
if(PF.enerSeg*ePartOpt < ener2Min || PF.enerSeg*(1.f-ePartOpt) < ener2Min) // keV
if(PF.enerSeg*ePartOpt < h2EminAbs || PF.enerSeg*(1.f-ePartOpt) < h2EminAbs) // keV
continue;
#if FIXED_METRIC == FIXED_METRIC_SQUARE
//// if the metric is 2, chiSquare has been calculated in EnergyPartition
......@@ -2357,10 +2389,10 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
float chi2Opt = 0;
float ePartOpt = EnergyPartition(PF, netChSeg, jPts, fxPt, chi2Opt);
// skip off range or too extreme energy partition
if(ePartOpt < 0.05f || ePartOpt > 0.95f)
if(ePartOpt < h2EminRel || ePartOpt > (1.f-h2EminRel))
continue;
// skip solutions with too small energy (i.e. more balanced partition for low-energy events)
if(PF.enerSeg*ePartOpt < ener2Min || PF.enerSeg*(1.f-ePartOpt) < ener2Min) // keV
if(PF.enerSeg*ePartOpt < h2EminAbs || PF.enerSeg*(1.f-ePartOpt) < h2EminAbs) // keV
continue;
#if FIXED_METRIC != FIXED_METRIC_SQUARE
//// if the metric is 2, chiSquare has been calculated in EnergyPartition
......@@ -2413,10 +2445,10 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
float chi2Opt = 0;
float ePartOpt = EnergyPartition(PF, netChSeg, fxPt, jPts, chi2Opt);
// skip off range or too extreme energy partition
if(ePartOpt < 0.05f || ePartOpt > 0.95f)
if(ePartOpt < h2EminRel || ePartOpt > (1.f-h2EminRel))
continue;
// skip solutions with too small energy (i.e. more balanced partition for low-energy events)
if(PF.enerSeg*ePartOpt < ener2Min || PF.enerSeg*(1.f-ePartOpt) < ener2Min) // keV
if(PF.enerSeg*ePartOpt < h2EminAbs || PF.enerSeg*(1.f-ePartOpt) < h2EminAbs) // keV
continue;
#if FIXED_METRIC != FIXED_METRIC_SQUARE
//// if the metric is 2, chiSquare has been calculated in EnergyPartition
......@@ -2461,7 +2493,7 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
}
// back to 1 point if no real improvement over the original 1-hit solution
if(PF.chi2min >= h1_chi2min*chi2FactorFinal || PF.bestFac > 0.95 || PF.bestFac < 0.05) {
if(PF.chi2min >= h1_chi2min*h2Chi2Final || PF.bestFac < h2EminRel || PF.bestFac > (1.f-h2EminRel)) {
PF.bestPt1 = h1_bestPt;
PF.bestPt2 = -1;
PF.bestFac = 1.f;
......@@ -2471,7 +2503,7 @@ int PSAFilterGridSearch::SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoar
// back to 1 point if the 2 points are too near
float fineDD = fBasis.PtPtDistance(netChSeg, PF.bestPt1, PF.bestPt2);
if(fineDD < ffMin) {
if(fineDD < h2DminFine) {
PF.bestPt1 = h1_bestPt;
PF.bestPt2 = -1;
PF.bestFac = 1.f;
......@@ -2525,7 +2557,7 @@ float PSAFilterGridSearch::EnergyPartition(pointFull &PF, int netChSeg, int kPts
// if the metric is 2, chiSquare is given by
chi2Opt = Sumy2y2 - ePartOpt*Sum12y2;
#else
// otherwise has to be calculated explicitely
// otherwise the caller has to calculate it explicitely
chi2Opt = -1;
#endif
......
......@@ -88,6 +88,13 @@ private:
float gDistMetric; // initialized to PMETRIC but user settable; may create inconsistencies if USE_SSE_VERSION
float fDistRescale; // used to consider RESCALE when calculating chi2 for non rescaled residues
float h2Chi2Start; // max chi2(coarse two hits)/chi2(one hit)
float h2Chi2Final; // Minimum final improvement of chi2 ti accept two hits
float h2EminRel; // minimum relative energy of the light component
float h2EminAbs; // Minimum absolute energy of the light component
float h2DminCoarse; // minimum distance (mm) of coarse grid points for the two hits solution
float h2DminFine; // minimum distance (mm) of fine grid points for the two hits solution
bool gPsaMatrXYZR; // enable/disable this spectrum
float tauSGCC[NSEGS+4]; // the 36 segments and the CC (and some spare)
......
......@@ -339,7 +339,7 @@ public:
// }
// delta->exp, one value for each segment and the CC
void convDeltaToExp(const float *Tau, float tDecayLong = 0) {
void convDeltaToExp(const float *Tau, float tDecayLong, bool smooth) {
for(int sg = 0; sg < NCHAN; sg++) { // including CC
float tau = Tau[sg];
float aa = (float)exp(-TIME_STEP/tau);
......@@ -363,20 +363,20 @@ public:
amplitude[sg][nn] = vv2;
}