Just refactoring

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@1087 170316e4-aea8-4b27-aad4-0380ec0519c9

WinCtest/femul.cpp

@@ -10,7 +10,7 @@
 //  - in the driver program (#define WCT_THREADED in ChainLocker.h).
 //
 // All libraries used so far in Agata Demonstrator experiments and analyses are available, using the same source code;
-//   indeed a consistent part of that code has been developed and/or tested by various contributors using this emulator.
+//   indeed, a large part of that code has been developed and/or tested by various contributors using this emulator.
 //
 // A C++ implementation of the EventBuilder replaces that of Narval (which is written in Ada).
 //

WinCtest/femul.vcproj

@@ -202,7 +202,7 @@
 				Name="VCCLCompilerTool"
 				Optimization="2"
 				EnableIntrinsicFunctions="true"
-				AdditionalIncludeDirectories="..\PRISMA\src\lib_prisma\include;"C:\Program Files (x86)\boost\boost_1_40";..\myADF0.2;..\..\adf\standalone;..\common;..\producers\Crystal;..\producers\Crystal\includeATCA;..\producers\AncillaryTCP;..\filters\Preprocessing;..\filters\Preprocessing\includePrePSA;..\filters\Ancillary;..\filters\Ancillary\includeVME;..\filters\PSA;..\filters\PSA\includePSA;..\filters\Tracking;..\filters\Tracking\includeOFT;..\filters\Tracking\includeMGT;..\builders;C:\root\include"
+				AdditionalIncludeDirectories="..\PRISMA\src\lib_prisma\include;"C:\Program Files (x86)\boost\boost_1_40";..\..\adf;..\WinCtest;..\common;..\producers\Crystal;..\producers\Crystal\includeATCA;..\producers\AncillaryTCP;..\filters\Preprocessing;..\filters\Preprocessing\includePrePSA;..\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"
 				PreprocessorDefinitions="WIN32;NDEBUG;_CONSOLE;_CRT_SECURE_NO_WARNINGS;NRV_TYPE=NRV_OFFLINE;TF_ROOTTREE"
 				RuntimeLibrary="2"
 				EnableFunctionLevelLinking="true"

builders/EventBuilder.h

@@ -2,6 +2,7 @@
 #define EVENTBUILDER_H_INCLUDED
 
 #include "ADFConfig.h"
+#include "commonDefs.h"
 #include "EventQueue.h"
 #include "ChainLocker.h"
 

common/ChainLocker.h

@@ -10,8 +10,8 @@
 //      the DISPATCHER writes to it only if the buffer is empty and then unlocks the mutex 
 //      the BUILDER reads the data in one shot or event-by-event and unlocks when finished (after resetting the size of stored data)
 // In case of non-threaded analyses:
-//      Lock return false if it is marked as already locked
-//      The other methods (Unlock included) are applied only if they are called in a locked state and with the ID of the locker 
+//      Lock() returns false if it is marked as already locked
+//      The other methods (Unlock() included) are applied only if they are called in a locked state and with the ID of the locker 
 
 // enable using threads
 #ifdef HAS_AGAPROCONFIG

common/commonDefs.h

@@ -113,7 +113,7 @@ const int     defTriggerSample  =  10;  // 60-10=50 samples passed to the PSA
 //////// CrystalProducer ////////
 /////////////////////////////////
 #define CRP_MULTIHIST     // enable MultiHist spectra
-//#define CPR_FromGRU       // enable GRU/ViGRU
+//#define CRP_FromGRU       // enable GRU/ViGRU
 
 /////////////////////////////////////
 //////// PreprocessingFilter ////////

filters/PSA/PSAFilterGridSearch.h

@@ -12,7 +12,7 @@
 
 //#define SPECMETRICS        // to produce specMetrics in SearchFullGrid and SearchAdaptive ==> will be very slow and is not thread-safe
 
-//! Implementation of a simple grid search method. 
+//! Implementation of a two-step grid-search method. 
 /**
    Coded by Roberto Venturelli
    Ported first to Narval by Joa Ljungvall
@@ -70,7 +70,6 @@ private:
   SignalBasis fBasis;
   float       mappedMetric[NMETRIC];
   char        hmask[NSEGS][NSEGS+4];  // Some padding to have a "good" length
-  float       chi2NetSegs[NSEGS];     // the chi2 of the net-charge segments 
   float       chi2AllSegs[NCHAN];     // the chi2 of all segments and CC 
 
   bool        gUseAdaptive;
@@ -116,15 +115,16 @@ private:
   int   ProcessTwoHits (pointFull &PF);
   void  InitialSolution(pointFull &PF, bool keep = false);
   void  PreSearchCoarse(pointFull &PF);
-  int   SearchFullGrid (pointFull &PF, int netChSeg);
+  int   SearchFullGrid (pointFull &PF, int netChSeg, float netChEner);
   int   SearchAdaptive1(pointFull &PF, int netChSeg, bool bCoarseOnly = false, bool bFirstOnly = false); // 1 hit
   int   SearchAdaptive2(pointFull &PF, int netChSeg, bool bCoarseOnly = false, bool bFirstOnly = false); // 2 hits
+  float EnergyPartition(pointFull &PF, int netChSeg, int pt1, int pt2, float & chi2Opt);
   void  MakeSearchWave (pointFull &PF, float fact, float maxVal);
   void  MakeLocalMask  (pointFull &PF);
   float CalcChi2Residue(pointFull &PF);
   float CalcChi2GridPts(int segNum, int pdNum, pointFull &PF);
   bool  WriteBaseAverPt(float scale, std::string ext);
-  bool  WriteBasePoint (float scale, pointPsa &pt, FILE *FP);
+  bool  WriteBasePoint (pointPsa &pt, float scale, FILE *FP);
 
   bool IsItWorthwhile(pointFull &PF) {return true;}  // just to start
   int  TryToFit(pointFull &PF);

filters/PSA/includePSA/GridSearchClasses.h

@@ -21,30 +21,40 @@
 // Because of this, the experimental traces cannot be longer than NTIME samples
 // For each signal we keep  TCHAN  = NSEGS+3 single traces:
 // - the first NSEGS traces are for the segment.
-// - the 3nd last trace     [at TCHAN-3==INDCC==NSEGS] is the (possibly differentiated) core
+// - the 3rd last trace     [at TCHAN-3==INDCC==NSEGS] is the (possibly differentiated) core
 // - the 2nd last trace     [at TCHAN-2] is a non-differentiated copy of the net charge segment.
 // - the     last trace     [at TCHAN-1] is a non-differentiated copy of the core.
-// The basis traces are normalized to MAXNORM (normally 1000.f) and the net charge segment is already differentiated.
+// The basis traces are normalized to MAXNORM (normally 1000.f). 
+// If DIFFLAG > 0, the net charge segment is already differentiated.
 // To use the SSE extensions, the structure is aligned to a 16-byte boundary and NTIME must be a multiple of 4.
 
+enum ENextXYZ {kNextMx=0, kNextPx=1, kNextMy=2, kNextPy=3, kNextMz=4, kNextPz=5};
+//#define kNextMx 0
+//#define kNextPx 1
+//#define kNextMy 2
+//#define kNextPy 3
+//#define kNextMz 4
+//#define kNextPz 5
+
+#define numNeigh 27
+
 class pointPsa  
 {
 public:
   float   x, y, z;                    // coordinates of the interaction (mm)
   int     netChargeSegment;           // its net-charge segment
   MEMALIGN_16                         // compiler-specific alignment directive
-  float   amplitude[TCHAN][NTIME];    // 36 segments + CC + non differentiated copies of net-charge segment and the core
+    float   amplitude[TCHAN][NTIME];    // 36 segments + CC + non differentiated copies of net-charge segment and the core
   int     addr_first;                 // valid/useful start of the amplitude
   int     addr_last;                  // valid/useful last  of the amplitude (margin included)
   float   sumTransients;              // sum of all transient signals
   float   sumDiffNetSeg;              // area of the differentiated net-charge segment
-  int     ixP, ixM;                   // index of x neighbours in the same segment (-1 if not existing)
-  int     iyP, iyM;                   // index of y neighbours in the same segment
-  int     izP, izM;                   // index of z neighbours in the same segment
+  int     nextXYZ[6];                 // [ixM, ixP, iyM, iyP, izM, izP] cartesian axis neighbours in the same segment (-1 if not existing)
+  //int     nextAll[numNeigh];          // all neigbours of a point; first is the point itself
 
   bool    isFullCharge;               // the net charge signal has reached full amplitude (the 2 bool moved here to get 40 bytes)
   bool    isInnerPoint;               // the point is not on the boundary of the segment
-  
+
   int     numFinePoints;              // if > 0 this is part of the coarse grid and this is the number of associated fine points
 
   float   ePath;                      // total distance travelled by the electrons
@@ -70,16 +80,14 @@ public:
     addr_last  = NTIME-1;
     sumTransients = 0;
     sumDiffNetSeg = 0;
-    ixP = ixM = -1;
-    iyP = iyM = -1;
-    izP = izM = -1;
+    for(int nn = 0; nn < 6; nn++) nextXYZ[nn] = -1;
+    //for(int nn = 0; nn < numNeigh; nn++) nextAll[nn] = -1;
     isFullCharge = false;
     isInnerPoint = false;
     numFinePoints = 0;
     ePath = hPath= 0;
     eS_CC = hS_CC = eS_SG = hS_SG = 0;
     eE_CC = hE_CC = eE_SG = hE_SG = 0;
-    //trapCC = trapSG = 1.f;
   }
 
   void FullCopy(pointPsa & PF) {
@@ -90,9 +98,8 @@ public:
     addr_last  = PF.addr_last;
     sumTransients = PF.sumTransients;
     sumDiffNetSeg = PF.sumDiffNetSeg;
-    ixP = PF.ixP; ixM = PF.ixM;
-    iyP = PF.iyP; iyM = PF.iyM;
-    izP = PF.izP; izM = PF.izM;
+    for(int nn = 0; nn < 6; nn++) nextXYZ[nn] = PF.nextXYZ[nn];
+    //for(int nn = 0; nn < numNeigh; nn++) nextAll[nn] = PF.nextAll[nn];
     isFullCharge  = PF.isFullCharge;
     isInnerPoint  = PF.isInnerPoint;
     numFinePoints = PF.numFinePoints;
@@ -106,14 +113,13 @@ public:
     hE_CC = PF.hE_CC;
     eE_SG = PF.eE_SG;
     hE_SG = PF.hE_SG;
-    //trapCC = PF.trapCC; trapSG = PF.trapSG;
   }
 
   void clearAmplitude() {
     memset(amplitude, 0, sizeof(amplitude));
   }
 
-  void SumOfSignals() {
+  void sumOfSignals() {
     sumTransients = 0;
     sumDiffNetSeg = 0;
     for(int ns = 0; ns < NSEGS; ns++) {
@@ -128,7 +134,7 @@ public:
     }
   }
 
-  float GetSumSeg(int ns) {
+  float getSumSeg(int ns) {
     if(ns < 0 || ns > NCHAN)
       return 0;
     float sum = 0;
@@ -221,7 +227,7 @@ public:
     memcpy(amplitude[INDCC],            amplitude[INDQC], sizeof(float)*NTIME);
   }
 
-  // moving average filter
+  // moving-average filter
   void smooth(int nn) {
     if(nn<2 || nn>=NTIME)
       return;
@@ -270,7 +276,7 @@ public:
   }
 
   // exp -> step
-  void deconvExpToStep(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++) {
@@ -317,8 +323,8 @@ public:
       float v0 = 0;
       psg = amplitude[sg];
       for(int ii = 1; ii < NTIME; ii++) {
-       v0 = v1; v1 = v2; v2 = psg[ii];
-       psg[ii-1] = 0.1f * (v0 + 8.f*v1 + v2); 
+        v0 = v1; v1 = v2; v2 = psg[ii];
+        psg[ii-1] = 0.1f * (v0 + 8.f*v1 + v2); 
       }
       v0 = v1; v1 = v2;
       psg[NTIME-1] = 0.1f * (v0 + 8.f*v1 + v2);
@@ -401,18 +407,18 @@ public:
   int     netChargeSegnum[NSEGS];     // normally given in the order of energy release
   float   netChargeEnergy[NSEGS];
   MEMALIGN_16
-  float   tAmplitude[TCHAN][NTIME];   // the total original experimental data
+    float   tAmplitude[TCHAN][NTIME];   // the total original experimental data
   MEMALIGN_16
-  float   wAmplitude[TCHAN][NTIME];   // the experimental data to be used and modified by the grid search
+    float   wAmplitude[TCHAN][NTIME];   // the experimental data to be used and modified by the grid search
   MEMALIGN_16
-  float   sAmplitude[TCHAN][NTIME];   // experimental trace manipulated for being used by the "chi2" loops
+    float   sAmplitude[TCHAN][NTIME];   // experimental trace manipulated for being used by the "chi2" loops
   MEMALIGN_16
-  float   rAmplitude[TCHAN][NTIME];   // the "fitted" trace
+    float   rAmplitude[TCHAN][NTIME];   // the "fitted" trace
 
   int     resPt1[NSEGS];              // the best 1st solution point for the numNetCharges segments
   int     resPt2[NSEGS];              // the best 2nd solution point                 "
   float   resFac[NSEGS];              // relative amplitude of 1st point             "
- 
+
   // working temporaries
   int     netSeg;                     // the specific net-charge segment during grid-search processing
   int     bestPt1;                    // the best matching first  base point
@@ -421,12 +427,12 @@ public:
   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
- 
+
   pointExp() {
     initIt();
   }
 private:
-  void initIt() {
+  void InitIt() {
     idGe    = 0;
     enerGe  = 0;
     sumSegm = 0;
@@ -477,8 +483,8 @@ public:
     chi2min = PF.chi2min;
   }
 
-  // Produces sAmplitude from wAmplitude:  to get the right dispersion when
-  // accessing the mapped-metrix value, sAmplitude is further amplified by RESCALE 
+  // Produces sAmplitude from wAmplitude. To get the right dispersion when
+  // accessing the mapped-metric value, sAmplitude is further amplified by RESCALE 
   // Differentiates netChargeSegment and Core, if needed.
   void MakeSearchWave(float fact, char *mask) {
     memset(sAmplitude, 0, sizeof(sAmplitude));
@@ -526,6 +532,7 @@ public:
 
   }
 
+  // Modify wAmplitude (float) adding the two base points passed in pPsa1 and pPsa2 
   void AddBaseTrace(pointPsa &pPsa1, pointPsa &pPsa2, const float bestFac, const float fact, const int samp_first, const int time_first) {
 
     if(bestFac == 1.f) {

filters/PSA/includePSA/SignalBasis.h

@@ -65,8 +65,9 @@ public:
   bool        bCoarseFine;    // the coarse-fine structures are valid
   bool        bCoarseOnly;    // only the coarse part of the coarse-fine structures is valid
   int         numPts[NSEGS];  // number of points in the segments
-  pointPsa      *Pts[NSEGS];  // the points of the segment, generated in Read()
+  pointPsa   *segPts[NSEGS];  // the points of the segment, generated in Read()
   pointPsa    averPt[NCHAN];  // the "center" of the segment (in INDCC the overall center)
+  float      *ddPtPt[NSEGS];  // precalculated point-point distance inside a segment
 
   float xTalkProp[NSEGS][NSEGS];
   float xTalkDiff[NSEGS][NSEGS];
@@ -108,14 +109,15 @@ public:
   SignalBasis() {
     bCenters = bFullGrid = bCoarseFine = bCoarseOnly = xTalkCorr = false;
     memset(numPts,    0,     sizeof(numPts));
-    memset(Pts,       0,     sizeof(Pts));
+    memset(segPts,    0,     sizeof(segPts));
+    memset(ddPtPt,    0,     sizeof(ddPtPt));
     memset(cflist,    0,     sizeof(cflist));
     memset(segBroken, false, sizeof(segBroken));
   }
   virtual ~SignalBasis() {
     for(int ns=0; ns<NSEGS; ns++) {
-      if(Pts[ns]) {
-        delete[] Pts[ns]; Pts[ns] = NULL;
+      if(segPts[ns]) {
+        delete[] segPts[ns]; segPts[ns] = NULL;
       }
       if(bCoarseFine) {
         delete [] cflist[ns].indFine;   cflist[ns].indFine   = NULL;
@@ -123,32 +125,39 @@ public:
       }
     }
   }
-  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  MakeCoarseFineList(bool both, bool verbose);
-  bool CoarseFineForOneSegment(int ns, bool both, bool verbose);
-  bool AssignIt(float px, float py, float pz, float cx, float cy, float cz);
-  bool IsCoarse(int jj, cgrid *cgr, int numC);
-  int  NearestPoint(float px, float py, float pz, int ns, float& dist);
-  int  NearestCoarse(float px, float py, float pz, cgrid *cgr, int numC, float& dist);
-  void SmoothBasisSignals(int nsmoo);
-  void ExperimentalResponse();
-  void ExponentialResponse(const float *TauSGCC);
-  void CalculateSegmentCenters();
-  void FindNeighbours();
-  int  StoreNetCharges();
-  void CalculateSumOfSignals();
-  int  DifferentiateNetChargeSignals(int lag);
-  void ReorderSamples();
-  bool ReadXtalkCoeffs(std::string cname, float xfactor, bool verbose);
-  bool CalcXtalkCoeffs(float *xTprop, float xfactor);
-  bool XtalkCorrection(int badSegment);
+  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);
+  void  CalcPtPtDistance(int ns);
+  float PtPtDistance(int ns, int n1, int n2);
+  int   MakeCoarseFineList(bool both, bool verbose);
+  bool  CoarseFineForOneSegment(int ns, bool both, bool verbose);
+  bool  AssignIt(float px, float py, float pz, float cx, float cy, float cz);
+  bool  IsCoarse(int jj, cgrid *cgr, int numC);
+  int   NearestPoint(float px, float py, float pz, int ns, float& dist);
+  int   NearestCoarse(float px, float py, float pz, cgrid *cgr, int numC, float& dist);
+  void  SmoothBasisSignals(int nsmoo);
+  void  ExperimentalResponse();
+  void  ExponentialResponse(const float *TauSGCC);
+  void  CalculateSegmentCenters();
+  void  FindNeighbours();
+  int   StoreNetCharges();
+  void  CalculateSumOfSignals();
+  int   DifferentiateNetChargeSignals(int lag);
+  void  ReorderSamples();
+  bool  ReadXtalkCoeffs(std::string cname, float xfactor, bool verbose);
+  bool  CalcXtalkCoeffs(float *xTprop, float xfactor);
+  bool  XtalkCorrection(int badSegment);
   float GetStepCoarse() {return float(cstep);}
   float GetStepFine()   {return float(fstep);}
 };
 
+inline float SignalBasis::PtPtDistance(int ns, int n1, int n2)
+{
+  return ddPtPt[ns][n1*numPts[ns] + n2];
+}
+
 inline bool SignalBasis::AssignIt(float px, float py, float pz, float cx, float cy, float cz)
 {
   //const float surface = 0;              // to exclude the surface
@@ -184,7 +193,7 @@ inline bool SignalBasis::IsCoarse(int jj, cgrid *cgr, int numC)
 inline int SignalBasis::NearestPoint(float px, float py, float pz, int ns, float& dist)
 {
   int       npts = numPts[ns];  // the number of points in this segment
-  pointPsa *ppts = Pts[ns];     // the points of this segment
+  pointPsa *ppts = segPts[ns];     // the points of this segment
 
   int   cnear = -1;
   float cdist = (float)1.e20;