Added possibility to set condition on group of ge-detectors in GlobalFilter and TrackingFilter.

Added Smoke Plot (by Francesco R.) in CrystalProducerATCA to analyze differential/integral non linearity of digitizers
Fixed femul to allow more than one LOOP...ENDLOOP block in the topology file


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

WinCtest/femul.cpp

@@ -523,6 +523,9 @@ bool topologyRead(string topo)
         cout << "Format is: Loop Name Value1 Value2 ..." << endl;
         return false;
       }
+      // reset vectors 
+      loopValues.clear();
+      loopActors.clear();
       // tokenize values
       while(stringSplit(values, value1, value2)) {
         loopValues.push_back(value1);
@@ -777,17 +780,15 @@ void listKeysAndExit()
 
   PostPSAFilter             *pf3  = new PostPSAFilter;
   pf3->GetParameters(true);
-  //pf3->ListParameters();
   delete pf3;
 
   EventBuilder              *pe1  = new EventBuilder("EventBuilder");
   pe1->GetParameters(true);
   delete pe1;
 
-  //// not really developed and used
-  //GlobalFilter              *pf4 = new GlobalFilter;
-  //pf4->GetParameters(true);
-  //delete pf4;
+  GlobalFilter              *pf4 = new GlobalFilter;
+  pf4->GetParameters(true);
+  delete pf4;
 
   TrackingFilter            *pf5  = new TrackingFilter;
   pf5->GetParameters(true);

WinCtest/femul.vcproj

@@ -308,7 +308,7 @@
 				UseLibraryDependencyInputs="false"
 				AdditionalDependencies="libCore.lib libTree.lib libRIO.lib"
 				LinkIncremental="1"
-				AdditionalLibraryDirectories="C:\Boost\lib64-VC90;C:\root\lib"
+				AdditionalLibraryDirectories=""C:\Boost\boost_1_40\lib64-VC90";C:\root\lib"
 				GenerateDebugInformation="true"
 				SubSystem="1"
 				OptimizeReferences="2"

common/commonDefs.h

@@ -95,6 +95,9 @@ const float  RAD2DEG = float(57.29577951308232286464772187173366546630859375);
 const float  DEG2RAD = 1.f/RAD2DEG;
 //const float  DEG2RAD = float(0.0174532925199432954743716805978692718781530857086181640625);
 
+const float  Rad2Deg = RAD2DEG;
+const float  Deg2Rad = DEG2RAD;
+
 // default/maximum length with the present ATCA electronics
 const int     defTraceLengthRaw = 160;
 // default trace length produced by the PreprocessingFilter

filters/Global/GlobalFilter.h

@@ -15,6 +15,9 @@
 #include "Trigger.h"
 #include "AgataCompositeFrame.h"
 
+#include <cmath>
+#include <vector>
+
 #include <ctime>
 
 #include "commonDefs.h"
@@ -71,20 +74,21 @@ protected:
     float X;
     float Y;
     float Z;
-    float T;        // time of hit (samples) referred to the common reference timestampEvt
-    UInt_t Det;     // detector ID
-    UInt_t Seg;     // segment number
-    UInt_t Ind;     // sequential crystal number (from which to take all needed info)
+    float T;            // time of hit (samples) referred to the common reference timestampEvt ???
+    UInt_t Det;         // detector ID
+    UInt_t Seg;         // segment number
+    UInt_t Ind;         // sequential crystal number (from which to take all needed info)
   } exyzHit;
 
   // Output gammas from tracking
-  // The only information needed is actually the identifier of the original PSA hit (I1 and I2)
+  // The only information needed is actually the identifier of the original PSA hit iH[nH]
   // the rest (xyz) is redundant and could be removed
   typedef struct {
     float E;
     float Ecorr;
-    int   nH;         // number of hits of this gamma
-    int   type;       // 0-photo 1-compt 2-pair
+    int   nH;           // number of hits of this gamma
+
+    int   type;         // 0-photo 1-compt 2-pair
     int   I1;         // original PSA hit of the 1st interaction point
     float X1, Y1, Z1; // info recoverable from the original PSA hit I1
     int   I2;         // original PSA hit of the 2nd interaction point
@@ -93,10 +97,8 @@ protected:
   } trGamma;
 
   typedef struct {
-    float rXX, rXY, rXZ;
-    float rYX, rYY, rYZ;
-    float rZX, rZY, rZZ;
-    float TrX, TrY, TrZ;
+    float rot[3][3];
+    float tra[3];
   } rotTr3D; 
 
   rotTr3D  *rotTr;
@@ -109,10 +111,13 @@ protected:
 
   int numCrystals;
   int numCrystalsGood;
-  crystdata cryst[NUMGEDETS];
+  crystdata *cryst;
 
   Float_t fCoreSumEnergy;
 
+  Int_t   fNumGeDets;         // number of Ge crystals (defaults to NUMGEDETS)
+  Int_t   fIDtoCryst[180];    // contains (index+1) in crys
+
   UInt_t *trigType;           // statistics of fired triggers
   UInt_t trigNums;            // number of triggers defined
   UInt_t trigUnknown;         // triggers fired and not known
@@ -151,7 +156,15 @@ protected:
   MultiHist<unsigned short> *GlobMatr_gg2;    // gamma-gamma of tracked and doppler corrected data data
   MultiHist<unsigned short> *GlobMatr_ETC;    // gamma-trigger
   MultiHist<unsigned short> *GlobMatr_SGSG;   // SG-SG counts for all pairs
-#endif  //GLF_MULTIHIST
+  MultiHist<unsigned int  > *GlobMatr_Fold;   // fold distributions
+  MultiHist<unsigned short> *GlobMatr_ZYX;    // 3D representation of points and lines joinning them
+  MultiHist<unsigned int  > *GlobSpec_XYZ;    // 3D representation of placed points
+
+  int SpecMap180[180];  // in principle there will not be higher ID's than this but it should be done better
+  int SM180(int ii) {return SpecMap180[ii];}
+#else
+  int SM180(int ii) {return ii;}
+#endif  //TRF_MULTIHIST
   cycleServer cServer;  // manager of cyclic operations
 
   // fixed values to correct when there is no Ancillary defining the recoil vector
@@ -163,8 +176,8 @@ protected:
 
 protected:
   //! input and ouput data frames
-  ADF::SharedFP *fFramePSA[NUMGEDETS];
-  ADF::SharedFP *fFrameRaw;
+  ADF::SharedFP **fFramePSA;
+  ADF::SharedFP  *fFrameRaw;
 
   ULong64_t timestampEvt;       // common reference for gamma timing (smallest tstamp among crystals or tstamp of ancillary)
   ULong64_t timestampCry;       // this is the smallest timestamp among the crystals; equals timestampEvt if no ancyllary
@@ -188,6 +201,18 @@ protected:
   std::string fOdirPrefix;
 
   std::string fWhichRotoTranslations;
+  typedef struct {
+    int id;
+    int ax;
+    float ang;
+  } extraRot;
+  std::vector<extraRot> eRot;
+  typedef struct {
+    int id;
+    float dx, dy, dz;
+  } extraShift;
+  std::vector<extraShift> eShift;
+
   std::string fGeometrySummary;
 
   bool    fDiscardEmpty;
@@ -214,6 +239,17 @@ protected:
   Float_t fSumEnerMax;
   Int_t   fDetFoldMin;
   Int_t   fDetFoldMax;
+  Int_t   fHitFoldMin;
+  Int_t   fHitFoldMax;
+
+  Float_t fGroupSumEnerMin;
+  Float_t fGroupSumEnerMax;
+  Int_t   fGroupDetFoldMin;
+  Int_t   fGroupDetFoldMax;
+  Int_t   fGroupHitFoldMin;
+  Int_t   fGroupHitFoldMax;
+
+  std::vector< std::vector<int> > GeGroups;
 
   Float_t fEnergyGain;          // gain of energy spectra
 
@@ -222,6 +258,12 @@ protected:
   Int_t   fSizeMatr_gg2;        // <=0 ==> nomatrix        
   Float_t fGainMatr_gg2;
 
+  Bool_t  fMatrZYX;
+  Float_t fMatrZYX_zOffset;     // mm  to subtract to the z coordinate for the positions inside the crystals
+  Float_t fMatrZYX_zScale1;     // mm->cm for the Z positions inside the crystals (fixed to 0.1)
+  Int_t   fMatrZYX_zRange;      // range of z values to intercept the lines joining the hits
+  Float_t fMatrZYX_zScale2;     // mm->?? for fMatrZYX_zRange
+
   //Bool_t      fRecal;
   //std::string fRecalFile;
   //float       fRecalCoeffs[NUMGEDETS][2];  // 0=offset, 1=gain
@@ -237,14 +279,44 @@ protected:
 
   void InitGenStructures();
   int  GetRotoTranslations(FILE *fp_conf);
-  void Transform(exyzHit *pt, int nb_det) {    // transformation to global frame
+  void TransformToGlobalFrame(exyzHit *pt, int nb_det) {
     float Xr = pt->X; 
     float Yr = pt->Y;
     float Zr = pt->Z;
     rotTr3D *prt = rotTr + nb_det;
-    pt->X = prt->rXX*Xr + prt->rXY*Yr + prt->rXZ*Zr + prt->TrX;
-    pt->Y = prt->rYX*Xr + prt->rYY*Yr + prt->rYZ*Zr + prt->TrY;
-    pt->Z = prt->rZX*Xr + prt->rZY*Yr + prt->rZZ*Zr + prt->TrZ; 
+    pt->X = prt->rot[0][0]*Xr + prt->rot[0][1]*Yr + prt->rot[0][2]*Zr + prt->tra[0];
+    pt->Y = prt->rot[1][0]*Xr + prt->rot[1][1]*Yr + prt->rot[1][2]*Zr + prt->tra[1];
+    pt->Z = prt->rot[2][0]*Xr + prt->rot[2][1]*Yr + prt->rot[2][2]*Zr + prt->tra[2]; 
+  }
+  void rotTr3D_rotX (int which, float adeg) {return rotTr3D_rotAx(which, adeg, 0);}
+  void rotTr3D_rotY (int which, float adeg) {return rotTr3D_rotAx(which, adeg, 1);}
+  void rotTr3D_rotZ (int which, float adeg) {return rotTr3D_rotAx(which, adeg, 2);}
+  void rotTr3D_shift(int which, float dx, float dy, float dz) {
+    rotTr[which].tra[0] += dx;
+    rotTr[which].tra[1] += dy;
+    rotTr[which].tra[2] += dz;
+  }
+private:
+  void rotTr3D_rotAx(int which, float adeg, int ax) {
+    rotTr3D *prt = rotTr + which;
+    rotTr3D  old; memcpy(&old, prt, sizeof(rotTr3D));
+    float cs = cos(adeg*Deg2Rad);
+    float sn = sin(adeg*Deg2Rad);
+    float rm[3][3] = {{0}};
+    rm[(ax+0)%3][(ax+0)%3] = 1; rm[(ax+0)%3][(ax+1)%3] =  0; rm[(ax+0)%3][(ax+2)%3] =   0;
+    rm[(ax+1)%3][(ax+0)%3] = 0; rm[(ax+1)%3][(ax+1)%3] = cs; rm[(ax+1)%3][(ax+2)%3] = -sn;
+    rm[(ax+2)%3][(ax+0)%3] = 0; rm[(ax+2)%3][(ax+1)%3] = sn; rm[(ax+2)%3][(ax+2)%3] =  cs;
+    float tmp;
+    for(int i = 0; i < 3; ++i) {
+      for(int j = 0; j < 3; ++j) {
+        tmp = 0;
+        for(int k = 0; k < 3; ++k) tmp += rm[i][k]*old.rot[k][j];
+        prt->rot[i][j] = tmp;
+      }
+      tmp = 0;
+      for(int k = 0; k < 3; ++k) tmp += rm[i][k]*old.tra[k];
+      prt->tra[i] = tmp;
+    }
   }
 
   //bool GetRecalCoeffs();

filters/Tracking/TrackingFilter.cpp

@@ -139,8 +139,6 @@ TrackingFilter::TrackingFilter() :
 
   fNumGeDets = NUMGEDETS;
 
-  cryst = new crystdata[fNumGeDets];
-
   fDetEnerMin = 20.f;       // should be much wider 
   fDetEnerMax = 20000.f;
   fSumEnerMin = 20.f;       // should be much wider 
@@ -150,6 +148,13 @@ TrackingFilter::TrackingFilter() :
   fHitFoldMin = 0;          // i.e. keep events with no gamma detectors
   fHitFoldMax = 1000000;    // big enough
 
+  fGroupSumEnerMin = fSumEnerMin;
+  fGroupSumEnerMax = fSumEnerMax;
+  fGroupDetFoldMin = fDetFoldMin;
+  fGroupDetFoldMax = fDetFoldMax;
+  fGroupHitFoldMin = fHitFoldMin;
+  fGroupHitFoldMax = fHitFoldMax;
+
   fEnergyGain   = 1.f;
 
   fSizeMatr_gg1 = 0;       // nomatrix
@@ -157,6 +162,12 @@ TrackingFilter::TrackingFilter() :
   fGainMatr_gg1 = 1;
   fGainMatr_gg2 = 1;
 
+  fMatrZYX         = false;   // mm  to subtract to the z coordinate for the positions inside the crystals
+  fMatrZYX_zOffset = 50.f;    // mm  to subtract to the z coordinate for the positions inside the crystals
+  fMatrZYX_zScale1 = 0.1f;    // mm->cm for the Z positions inside the crystals
+  fMatrZYX_zRange  = 10;      // range of z values to intercept the lines joining the hits
+  fMatrZYX_zScale2 = 1.0f;    // units for fMatrZYX_zRange
+
   //fRecal = false;
   //fRecalFile.clear();
   //for(int ii = 0; ii < NUMGEDETS; ii++) {
@@ -199,6 +210,8 @@ TrackingFilter::TrackingFilter() :
   TrackMatr_ETC  = NULL;
   TrackMatr_SGSG = NULL;
   TrackMatr_Fold = NULL;
+  TrackMatr_ZYX  = NULL;
+  TrackSpec_XYZ  = NULL;
 
   for(int ii = 0; ii < 180; ii++)
     SpecMap180[ii] = ii;
@@ -343,6 +356,8 @@ void TrackingFilter::process_initialise (UInt_t *error_code)
     return;
   }
 
+  cryst = new crystdata[fNumGeDets];
+
   cout << endl;
   DefineTriggers(error_code);
   if(*error_code) {
@@ -442,6 +457,14 @@ void TrackingFilter::process_initialise (UInt_t *error_code)
   }
   fclose(fp_conf);
 
+  // other rotations and translations in the global frame
+  for(size_t nn = 0; nn < eRot.size(); ++nn) {
+    rotTr3D_rotAx(eRot[nn].id, eRot[nn].ang, eRot[nn].ax);
+  }
+  for(size_t nn = 0; nn < eShift.size(); ++nn) {
+    rotTr3D_shift(eShift[nn].id, eShift[nn].dx, eShift[nn].dy, eShift[nn].dz);
+  }
+
   //if(fRecalFile.size() > 0) {
   //  bool ok = GetRecalCoeffs();
   //  if(!ok) {
@@ -823,6 +846,7 @@ Int_t TrackingFilter::SetInput()
   number_of_crystals = 0;
   number_of_hits     = 0;
   number_of_gammas   = 0;
+  memset(fIDtoCryst, 0, sizeof(fIDtoCryst));
   fCoreSumEnergy     = 0;
   numRawItems = 0;
   len_rawBuf  = 0;
@@ -1063,6 +1087,8 @@ void TrackingFilter::GetDataPSA()
 
   fCoreSumEnergy += CoreE[0];
 
+  fIDtoCryst[cryst[fNumGeDets].id] = fNumGeDets+1;
+
   // loop over number of interaction points in fFramePSA
   float esumhits = 0;
   int   hit1st   = number_of_hits;
@@ -1295,13 +1321,34 @@ int TrackingFilter::PreProcessEvent()
 { 
   if(number_of_crystals < fDetFoldMin || number_of_crystals > fDetFoldMax)
     return 1;
-
   if(number_of_hits < fHitFoldMin || number_of_hits > fHitFoldMax)
     return 1;
-
   if(fCoreSumEnergy < fSumEnerMin || fCoreSumEnergy > fSumEnerMax)
     return 1;
 
+  // determine fold and sum-energy of detector groups and apply gates
+  for(size_t gr = 0; gr < GeGroups.size(); gr++) {
+    float sumEnergy = 0;
+    float countDets = 0;
+    float countSegs = 0;
+    for(size_t nn = 0; nn < GeGroups[gr].size(); nn++) {
+      int gid = GeGroups[gr][nn];
+      int ind = fIDtoCryst[gid];
+      if(ind) {
+        crystdata *pc = cryst + (ind -1);
+        sumEnergy += pc->CCE0;
+        countSegs += pc->nhits;
+        ++countDets;
+      }
+    }
+    if(sumEnergy < fGroupSumEnerMin || sumEnergy > fGroupSumEnerMax)
+      return 1;
+    if(countDets < fGroupDetFoldMin || countDets > fGroupDetFoldMax)
+      return 1;
+    if(countSegs < fGroupDetFoldMin || countSegs > fGroupDetFoldMax)
+      return 1;
+  }
+
   if( !AnalysisOfCrystals() )
     return 1;
 
@@ -1407,6 +1454,33 @@ void TrackingFilter::ConfParameters(std::string keyw, std::string data, Bool_t &
     fWhichRotoTranslations  = data;
   GP_END;
 
+  GP_PAR("GeRotX", "%d %f", "rotate the given crystal (ID) by the given angle (deg) around the global X axis ");
+    extraRot rtmp;
+    ok = 2 == sscanf(data.c_str(), "%d %f", &rtmp.id, &rtmp.ang);
+    rtmp.ax = 0;
+    eRot.push_back(rtmp);
+  GP_END;
+
+  GP_PAR("GeRotY", "%d %f", "rotate the given crystal (ID) by the given angle (deg) around the global Y axis ");
+    extraRot rtmp;
+    ok = 2 == sscanf(data.c_str(), "%d %f", &rtmp.id, &rtmp.ang);
+    rtmp.ax = 1;
+    eRot.push_back(rtmp);
+  GP_END;
+
+  GP_PAR("GeRotZ", "%d %f", "rotate the given crystal (ID) by the given angle (deg) around the global Z axis ");
+    extraRot rtmp;
+    ok = 2 == sscanf(data.c_str(), "%d %f", &rtmp.id, &rtmp.ang);
+    rtmp.ax = 2;
+    eRot.push_back(rtmp);
+  GP_END;
+
+  GP_PAR("GeShift", "%d %f %f %f", "shift the given crystal (ID) by the given dx, dy, dz (mm) in the global frame ");
+    extraShift stmp;
+    ok = 4 == sscanf(data.c_str(), "%d %f %f %f", &stmp.id, &stmp.dx, &stmp.dy, &stmp.dz);
+    eShift.push_back(stmp);
+  GP_END;
+
   GP_PAR("GeometrySummary", "%s", "for the header of the Mgt_Hits.txt file of input hits");
     ok = data.size() > 0;
     fGeometrySummary = data;
@@ -1475,6 +1549,38 @@ void TrackingFilter::ConfParameters(std::string keyw, std::string data, Bool_t &
     ok = 2 == sscanf(data.c_str(), "%d %d", &fHitFoldMin, &fHitFoldMax);
   GP_END;
 
+  GP_PAR("GeGroup", "%d ...", "define a group of germanium detectors containing the listed IDs");
+    vector<int> gid;
+    int num = stringCount(data);
+    string tmp1 = data;
+    string tmp2, tmp3;
+    ok = true;
+    for(int nn = 0; ok && nn < num; nn++) {
+      stringSplit(tmp1, tmp2, tmp3);
+      int id;
+      ok |= 1 == sscanf(tmp2.c_str(), "%d", &id);
+      if(ok) {
+        gid.push_back(id);
+        tmp1 = tmp3;
+      }
+    }
+    if(ok) {
+      GeGroups.push_back(gid);
+    }
+  GP_END;
+  
+  GP_PAR("GroupSumEnergyGate", "%f %f", "gate on the sum-energy of the crystals in the defined groups");
+    ok = 2 == sscanf(data.c_str(), "%f %f", &fGroupSumEnerMin, &fGroupSumEnerMax);
+  GP_END;
+
+  GP_PAR("GroupNumDetsGate", "%d %d", "event selection based on number of fired detectors in the defined groups");
+    ok = 2 == sscanf(data.c_str(), "%d %d", &fGroupDetFoldMin, &fGroupDetFoldMax);
+  GP_END;
+
+  GP_PAR("GroupNumHitsGate", "%d %d", "event selection based on total number of hits in the defined groups");
+    ok = 2 == sscanf(data.c_str(), "%d %d", &fGroupHitFoldMin, &fGroupHitFoldMax);
+  GP_END;
+
   // ???????
   GP_PAR("WriteMgtData", "%s ...", "Ancillary &| Timing &| Extended  write file of input hits for analysis with mgt");
 #ifdef TFR_WRITE_HITS_MGT
@@ -1556,6 +1662,11 @@ void TrackingFilter::ConfParameters(std::string keyw, std::string data, Bool_t &
     ok = 2 == sscanf(data.c_str(), "%d %f", &fSizeMatr_gg2, &fGainMatr_gg2);
   GP_END;
 
+  GP_PAR("MatrixZYX", "%f %f %d %f", "hits and intercepts: zOffset and zScale for hits; range and scale of z-planes");
+    ok = 4 == sscanf(data.c_str(), "%f %f %d %f", &fMatrZYX_zOffset, &fMatrZYX_zScale1, &fMatrZYX_zRange, &fMatrZYX_zScale2);
+    fMatrZYX = ok;
+  GP_END;
+
   GP_PAR("Details", "", "list input events in symbolic form (very slow and incomplete)");
     fDetails = true;
     ok = true;
@@ -1617,7 +1728,9 @@ void TrackingFilter::InitGenStructures()
   }
 #endif
 
-  rotTr = (rotTr3D *)calloc(   180, sizeof(rotTr3D));
+  rotTr = (rotTr3D *)calloc(   180, sizeof(rotTr3D)); 
+  for(int nn = 0; nn < 180; nn++) rotTr[nn].rot[0][0] = rotTr[nn].rot[1][1] = rotTr[nn].rot[2][2] = 1;
+
   pEXYZ = (exyzHit *)calloc(intmax, sizeof(exyzHit));
   pGams = (trGamma *)calloc(intmax, sizeof(trGamma));
 
@@ -1682,6 +1795,11 @@ void TrackingFilter::InitGenStructures()
     TrackSpec_ES->setComment("sum-energy spectra (gain 1 keV/chan)");
     hGroup.add(TrackSpec_ES);
 
+    //TrackSpec_XYZ = new MultiHist<unsigned int>(fNumGeDets, 3, 4000);
+    //TrackSpec_XYZ->setFileName(fOdirPrefix+"Track?XYZ.spec");
+    //TrackSpec_XYZ->setComment("3D position of interaction points");
+    //hGroup.add(TrackSpec_XYZ);
+
     TrackSpec_EE = new MultiHist<unsigned int>(2, 15, specLenE);
     TrackSpec_EE->setFileName(fOdirPrefix+"Track?EE.spec");
     TrackSpec_EE->setComment("input and tracked spectra");
@@ -1737,6 +1855,13 @@ void TrackingFilter::InitGenStructures()
     TrackMatr_SGSG->setComment("segment-segment coincidences for all pairs");
     hGroup.add(TrackMatr_SGSG);
 
+    if(fMatrZYX) {
+      TrackMatr_ZYX = new MultiHist<unsigned short>(200, 200, 200);
+      TrackMatr_ZYX->setFileName(fOdirPrefix+"Track?ZYX.matr");
+      TrackMatr_ZYX->setComment("3D position of interaction points");
+      hGroup.add(TrackMatr_ZYX);
+    }
+
     if(fAncillary) {
       TrackSpec_EA = new MultiHist<unsigned int>(2, fNumGeDets, specLenE); // 0 alone, 1 Ancillary
       TrackSpec_EA->setFileName(fOdirPrefix+"Track?EA.spec");
@@ -1823,9 +1948,9 @@ int TrackingFilter::GetRotoTranslations(FILE *fp)
 
     rotTr3D *prt = rotTr + nd;
 
-    prt->TrX = x; // in mm
-    prt->TrY = y;
-    prt->TrZ = z;
+    prt->tra[0] = x; // in mm
+    prt->tra[1] = y;
+    prt->tra[2] = z;
 
     // matrix elements of rotation
     if(fscanf(fp, "%d  %f %f %f", &k, &x, &y, &z) != 4)
@@ -1833,25 +1958,25 @@ int TrackingFilter::GetRotoTranslations(FILE *fp)
     if(k != 1)
       return -1;
 
-    prt->rXX = x;
-    prt->rXY = y;
-    prt->rXZ = z;
+    prt->rot[0][0] = x;
+    prt->rot[0][1] = y;
+    prt->rot[0][2] = z;
 
     if(fscanf(fp, "%d  %f %f %f", &k, &x, &y, &z) != 4)
       return -1;
     if(k != 2)
       return -1;
-    prt->rYX = x;
-    prt->rYY = y;
-    prt->rYZ = z;
+    prt->rot[1][0] = x;
+    prt->rot[1][1] = y;
+    prt->rot[1][2] = z;
 
     if(fscanf(fp, "%d  %f %f %f", &k, &x, &y, &z) != 4)
       return -1;
     if(k != 3)
       return -1;
-    prt->rZX = x;
-    prt->rZY = y;
-    prt->rZZ = z;
+    prt->rot[2][0] = x;
+    prt->rot[2][1] = y;
+    prt->rot[2][2] = z;
 
     nt++;
   }
@@ -2236,6 +2361,45 @@ bool TrackingFilter::AnalysisOfCrystals()
       }
     }
   }
+
+  if(TrackSpec_XYZ) {
+    for(int nn = 0; nn < number_of_hits; nn++) {
+      exyzHit pt = pEXYZ[nn];
+      int id = SM180(pt.Det);
+      TrackSpec_XYZ->incr(id, 0, int(2000 + pt.X));
+      TrackSpec_XYZ->incr(id, 1, int(2000 + pt.Y));
+      TrackSpec_XYZ->incr(id, 2, int(2000 + pt.Z));
+    }
+  }
+
+  if(TrackMatr_ZYX) {
+    if(fMatrZYX_zScale1) {
+      for(int nn = 0; nn < number_of_hits; nn++) {
+        exyzHit pt = pEXYZ[nn];
+        float z = (pt.Z > 0) ? (pt.Z - fMatrZYX_zOffset)*fMatrZYX_zScale1 : (pt.Z + fMatrZYX_zOffset)*fMatrZYX_zScale1;
+        TrackMatr_ZYX->incr(int(z + 100.f), int(pt.Y + 100.f), int(pt.X + 100.f));
+      }
+    }
+    if(fMatrZYX_zScale2) {
+      for(int nn1 = 0; nn1 < number_of_hits-1; nn1++) {
+        exyzHit pt1 = pEXYZ[nn1];
+        for(int nn2 = nn1+1; nn2 < number_of_hits; nn2++) {
+          exyzHit pt2 = pEXYZ[nn2];
+          // stright line from pt1 to pt2
+          float vx = pt2.X - pt1.X;
+          float vy = pt2.Y - pt1.Y;
+          float vz = pt2.Z - pt1.Z;
+          for(int iz = -fMatrZYX_zRange; iz <= fMatrZYX_zRange; iz++) { // cm
+            float vv = (iz/fMatrZYX_zScale2 - pt1.Z)/vz;
+            float xx = pt1.X + vx *vv;
+            float yy = pt1.Y + vy *vv;
+            TrackMatr_ZYX->incr(iz+100, int(yy+100.f), int(xx+100.f));
+          }
+        }
+      }
+    }
+  } // TrackMatr_ZYX
+
   if(TrackSpec_EE) {
     TrackSpec_EE->Incr(0,   8, int(fEnergyGain*esumCC));                    // [0][8] input sum energy
     if(number_of_crystals==1)

filters/Tracking/TrackingFilter.h

@@ -16,6 +16,9 @@
 #include "Trigger.h"
 #include "AgataCompositeFrame.h"
 
+#include <cmath>
+#include <vector>
+
 //#include "AgataKeyFactory.h"
 //#include "AgataFrameFactory.h"
 //#include "AgataCompositeFrame.h"
@@ -123,10 +126,12 @@ protected:
   } trGamma;
 
   typedef struct {
-    float rXX, rXY, rXZ;
-    float rYX, rYY, rYZ;
-    float rZX, rZY, rZZ;
-    float TrX, TrY, TrZ;
+    float rot[3][3];
+    float tra[3];
+    //float rXX, rXY, rXZ;
+    //float rYX, rYY, rYZ;
+    //float rZX, rZY, rZZ;
+    //float TrX, TrY, TrZ;
   } rotTr3D;