/* Tracking-actor base class, J. Ljungvall 2008,
* based on implementation by Olivier Stezowski.
* 
* Modified and maintained by
*   D.Mengoni
*   D.Bazzacco
*   F.Recchia
*/

#include "TrackingFilter.h"
#include "AgataKeyFactory.h"
#include "AgataFrameFactory.h"

#include "DFAgent.h"

#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <cmath>
#include <cstdlib>

#include "misc.h"

std::string TrackingFilter::gMotherClass = "TrackingFilter";
std::string TrackingFilter::gActualClass;

using namespace std;
using namespace ADF;

//#define CHECKRAWBUF     // uncomment to check repetitions in rawbuf

const int   specLenE   = 16*1024;     // size of energy spectra

const int   tGain      = 1/*00*/;           // to increase tScale and the length of time spectra in a consistent way
const float tScale     = tGain*10.f;  // tScale=10 ==> 1ns/ch
const int   specLenT1  = tGain*1000;

const int   specLenT2  = 400;

const float fEnergyGainET = 0.05f;
const int   intmax     =  800;

TrackingFilter::TrackingFilter() :
  fTriggerEventDataPSA("event:data:psa"),
  fTriggerEventData_1X("event:data &event:data:psa |data:rancN"),
  fTriggerEventData_01("event:data !event:data:psa &data:rancN"),
  fTriggerDataPSA     ("data:psa"),
  fTriggerDataRawAnc  ("data:rancN"),
  fFramePSA(NULL),
  fFrameEvb(NULL),
  fFrameRaw(NULL)
{
  fExcludeTracking = false;

  fOdirPrefix.clear();

  trigType = NULL;
  trigNums = trigUnknown = 0;

  fWhichRotoTranslations  = "CrystalPositionLookUpTable";

  fAncRawFrameNum = 1;
  fAncRawFrameStr = "data:ranc1";

  fAncillary = false;
  fRecoiling = false;
  DopplerCorrect = NULL;

  fSource[0] = fSource[1] = fSource[2] = 0;
  fRecoil[0] = fRecoil[1] = 0; fRecoil[2] = 1;    // along beam (??)
  fRecoilVc  = 0;                                 // not moving
  fVerbose       = false;
  fDebug         = false;
  fDiscardEmpty  = true;
  fSetWriteModel = false;

  fWriteModel = ConfAgent::kStrict;

  //fBlockIn.SetModeIO(ConfAgent::kRead);
  //fBlockOut.SetModeIO(ConfAgent::kWrite);

  pEXYZ = NULL;
  rotTr = NULL;
  pGams = NULL;
  number_of_hits   = 0;
  number_of_gammas = 0;
  number_of_tracked_gammas = 0;

  fCoreSumEnergy = 0;

  rawBufSize  = 0;
  rawBuf      = NULL;
  rawBufOld   = NULL;

  numRawItems = 0;
  len_rawBuf  = 0;
  has_rawBuf  = 0;

  crystal_id     = 0;
  crystal_status = 0;
  CoreE[0] = CoreE[1] = 0;
  CoreT[0] = CoreT[1] = 0;
  timestampEvt = 0;
  evnumberEvt  = 0; 
  timestampCry = 0;
  evnumberCry  = 0; 
  evnumberAnc  = 0; 
  timestampAnc = 0;

  stop_called = false;

  hasAncillary = false;
  hasRecoil    = false;

  fTsGeAnc_min = 0x80000000;
  fTsGeAnc_max = 0x7FFFFFFF;

  //fForceSegmentsToCore = false;

  fWriteMgtData      = false;
  fWriteMgtAncillary = false;
  fWriteMgtTiming    = false;
  fWriteMgtExtended  = false;

  fWriteRootTree  = false;
  fWriteRootHits  = false;

  fWriteGsortData = false;

  fWriteTracked   = false;

  fDetEnerMin = 20.f;     // should be much wider 
  fDetEnerMax = 20000.f;
  fSumEnerMin = 20.f;     // should be much wider 
  fSumEnerMax = NUMGEDETS*fDetEnerMax;
  fDetFoldMin = 0;        // i.e. keep events with no gamma detectors
  fDetFoldMax = NUMGEDETS;

  fEnergyGain   = 1.f;

  fSizeMatr_gg1 = 0;       // nomatrix
  fSizeMatr_gg2 = 0;       // nomatrix
  fGainMatr_gg1 = 1;
  fGainMatr_gg2 = 1;

  //fRecal = false;
  //fRecalFile.clear();
  //for(int ii = 0; ii < NUMGEDETS; ii++) {
  //  fRecalCoeffs[ii][0] = 0;
  //  fRecalCoeffs[ii][1] = 1.f;
  //}

  rVx = rVy = rVz = rVc = 0;
  aPx = aPy = aPz = 0;

  fAcceptanceValue = 0; // 0 is recognized by the daughter class as default
  fGeometrySummary = "235.008 329.202 15 3 6 6 6 6 6 6";

  fp_conf    = NULL;
  fp_tracked = NULL;
  fp_hits    = NULL;

  fUseMultiHist = true;
#ifdef TRF_MULTIHIST
  TrackSpec_Fold = NULL;
  TrackSpec_EC   = NULL;
  TrackSpec_ES   = NULL;
  TrackSpec_TC   = NULL;
  TrackSpec_EE   = NULL;
  TrackSpec_TS   = NULL;
  TrackSpec_TT   = NULL;
  TrackSpec_TG   = NULL;
  TrackMatr_IEET = NULL;
  TrackSpec_EA   = NULL;
  TrackSpec_TA   = NULL;
  TrackSpec_TOF  = NULL;
  TrackSpec_Anc  = NULL;
  TrackMatr_TOF  = NULL;
  TrackMatr_AG   = NULL;
  TrackSpec_TTA  = NULL;
  TrackMatr_ETA  = NULL;
  TrackMatr_gg1  = NULL;
  TrackMatr_gg2  = NULL;
  TrackMatr_ETC  = NULL;
  TrackMatr_SGSG = NULL;
  TrackMatr_Fold = NULL;
#endif  //TRF_MULTIHIST

  RecoilVc = 0;
  RecoilGm = 1;
  RecoilCx = 0;
  RecoilCy = 0;
  RecoilCz = 1.;

#ifdef TRF_ROOTTREE
  rootfile    = NULL;
  t           = NULL;
  gammaE      = NULL;
  gammaT      = NULL;
  gammaTstamp = NULL;
  gammaTrType = NULL;
  gammaX1     = NULL;
  gammaY1     = NULL;
  gammaZ1     = NULL;
  gammaX2     = NULL;
  gammaY2     = NULL;
  gammaZ2     = NULL;
  hitE        = NULL;
  hitX        = NULL;
  hitY        = NULL;
  hitZ        = NULL;
  hitT        = NULL;
  hitTstamp   = NULL;
  hitId       = NULL;
  hitSg       = NULL;
#endif  // #ifdef TRF_ROOTTREE

#ifdef TRF_FromGRU
  Int_t port=9095;
  SpectraDB = new GSpectra();
  NetworkRoot = new GNetServerRoot (SpectraDB);
  NetworkRoot->SetPort(port);
  SpectraDB->SetfDefaultHostPort(port);
  SpectraDB->SetfDefaultHostName((char*)gSystem->HostName());
  NetworkRoot->SetVerbose(0);

  HitDistrib          = new TH1I("HitDistrib","Number of hits", 100, 0, 99);
  TrackedEnergies     = new TH1F("Trackedenergies",    "energy",    32768,0,32767);
  TrackedEnergiesCorr = new TH1F("TrackedenergiesCorr","energyCorr",32768,0,32767);
  CE_Tracking         = new TH1F("SumHitEnergy","SumHitEnergy",32768,0,32767);
  CC3D_BefAftRT       = new TH3F("Hist_position_befaft_transform","Hist_position_befaft_transform",256,-200,200,256,-200,200,256,-300,50);

  SpectraDB->AddSpectrum (HitDistrib);
  SpectraDB->AddSpectrum (TrackedEnergies);
  SpectraDB->AddSpectrum (TrackedEnergiesCorr);
  SpectraDB->AddSpectrum (CE_Tracking);
  SpectraDB->AddSpectrum (CC3D_BefAftRT);  

  NetworkRoot->StartGNetServer(false);
#endif  // TRF_FromGRU

#ifdef TRF_WRITE_GSORT
  gaspEvent  = NULL;
  gaspBuffer = NULL;
#endif // TRF_WRITE_GSORT

}

TrackingFilter::~TrackingFilter()
{
  cServer.Finish();

#ifdef TRF_FromGRU  
  delete HitDistrib;
  delete TrackedEnergies;
  delete TrackedEnergiesCorr;
  delete CE_Tracking;
  delete CC3D_BefAftRT;
#endif

#ifdef TRF_ROOTTREE
  //if(fWriteRootTree)
  //  closeTTree();
  //if(rootfile)    delete rootfile;
  //if(t)           delete t;
  if(gammaE)      delete [] gammaE;
  if(gammaT)      delete [] gammaT;
  if(gammaTstamp) delete [] gammaTstamp;
  if(gammaTrType) delete [] gammaTrType;
  if(gammaX1)     delete [] gammaX1;
  if(gammaY1)     delete [] gammaY1;
  if(gammaZ1)     delete [] gammaZ1;
  if(gammaX2)     delete [] gammaX2;
  if(gammaY2)     delete [] gammaY2;
  if(gammaZ2)     delete [] gammaZ2;
  if(hitE)        delete [] hitE;
  if(hitX)        delete [] hitX;
  if(hitY)        delete [] hitY;
  if(hitZ)        delete [] hitZ;
  if(hitT)        delete [] hitT;
  if(hitTstamp)   delete [] hitTstamp;
  if(hitId)       delete [] hitId;
  if(hitSg)       delete [] hitSg;
#endif  // #ifdef TRF_ROOTTREE

#ifdef TRF_WRITE_GSORT
  if(gaspBuffer) 
    delete gaspBuffer;
  if(gaspEvent) 
    delete gaspEvent;
#endif // TRF_WRITE_GSORT
} 

void TrackingFilter::process_config(const Char_t *directory_path, UInt_t *error_code)
{
  *error_code = 0;
  cout << gMotherClass + "::process_config() called with directory_path = " << directory_path << endl;

  // first init narval stuff
  NarvalInterface::process_config(directory_path, error_code);
  if( *error_code )
    return;

  // Get name of daughter class from from configuration directory_path/gMotherClass.conf 
  int rv = getKeyFromFile(GetGlobalConfPath() + gMotherClass + ".conf", "ActualClass", gActualClass);
  if(rv == 2)
    *error_code = 102; // Fatal error because the configuration file MUST be present
}

void TrackingFilter::process_initialise (UInt_t *error_code)
{
  Log.ClearMessage();
  Log.GetProcessName() = gMotherClass; 
  Log.SetProcessMethod("process_initialise");
  Log.SetPID(GetPID());

  GetParameters(error_code);
  if(*error_code) {
    Log << dolog;
    return;
  }

  cout << endl;
  DefineTriggers(error_code);
  if(*error_code) {
    Log << dolog;
    return;
  }

  if(fSetWriteModel)
    GetConfAgent()->GetDFAgent()->SetModel(fWriteModel);

  // state set to kIdle so that the data can be treated 
  fFrameIO.SetStatus(BaseFrameIO::kIdle);

  if(fAncillary && fRecoiling)
    DopplerCorrect = &TrackingFilter::DopplerCorrectAnc;
  else
    DopplerCorrect = &TrackingFilter::DopplerCorrectFix;

  // algo-related initializations 
  InitGenStructures();  // parts managed by the base class

#ifdef TRF_ROOTTREE
  if(fWriteRootTree)
    openTTree();
#else
  if(fWriteRootTree)
    cout << "\nWriteRootTree --> IGNORED (Program built without TRF_ROOTTREE)\n" << endl;
#endif

#ifdef TRF_WRITE_TRACKED
  if(fWriteTracked) {
  // opening file where tracked energies will be written out
  tracked_file = fOdirPrefix + "Oft_TrackedEnergies";
  fp_tracked = fopen(tracked_file.c_str(), "w");
  if(fp_tracked == NULL) {
    cout << "Error opening " << tracked_file << endl;
    *error_code = 100;
    return;
  }
  cout << "Tracked energies will be written to " << tracked_file << endl;
  }
#else
  if(fWriteTracked)
    cout << "\nWriteTracked --> IGNORED (Program built without TRF_WRITE_TRACKED)\n" << endl;
#endif

#ifdef TFR_WRITE_HITS_MGT
  if(fWriteMgtData) {
    if(!OpenInputHits()) {
      *error_code = 100;
      Log << dolog;
      return;
    }
  }
#else
  if(fWriteMgtData)
    cout << "\nWriteInputHits --> IGNORED (Program built without TFR_WRITE_HITS_MGT)\n" << endl;
#endif

#ifdef TRF_WRITE_GSORT
  if(fWriteGsortData) {
    string prismaLUTfile = GetConfPath() + "lutPRISMA.txt";
    if(!ReadGaspLUT( prismaLUTfile ) ) {
      cout << "Error reading PRISMA LUT " << prismaLUTfile << endl;
      *error_code = 100;
      return;
    }
    gaspEvent  = new unsigned short [8192];
    gaspBuffer = new CGaspBuffer();
    gaspBuffer->SetPrefix(fOdirPrefix+"Gsort_");
    gaspBuffer->BeginOfRun( gaspLUT[599]-1 );  // run number
  }
#else
  if(fWriteGsortData)
    cout << "\nWriteGsortData --> IGNORED (Program built without TRF_WRITE_GSORT)\n" << endl;
#endif // TRF_WRITE_GSORT
  

  // configuration file to get transformation from relative to world coordinates
  conf_file = GetConfPath() + fWhichRotoTranslations;
  cout << "Opening configuration file " << conf_file << endl;
  fp_conf = fopen(conf_file.c_str(), "r");
  if(!fp_conf) {
    cout << "Could not open configuration file " << conf_file << endl;
    *error_code = 100;
    Log << dolog;
    return;
  }

  int nr = GetRotoTranslations(fp_conf);
  cout << "Found " << nr << " coefficients" << endl;
  if(nr < 3) {
    cout << "Error reading configuration file " << conf_file << endl;
    *error_code = 100;
    Log << dolog;
    return;
  }
  fclose(fp_conf);

  //if(fRecalFile.size() > 0) {
  //  bool ok = GetRecalCoeffs();
  //  if(!ok) {
  //    cout << "Error reading recalibration coefficients from " << fRecalFile << endl;
  //    *error_code = 100;
  //    Log << dolog;
  //    return;
  //  }
  //}

  number_of_tracked_gammas   = 0;

  // version-specific initializations
  *error_code = AlgoSpecificInitialise();

  cServer.SetCommandFile(GetConfPath() + gMotherClass + ".live");
#ifdef TRF_MULTIHIST
  cServer.SetHistGroup(&hGroup);
#endif  //TRF_MULTIHIST

  cServer.Start(gMotherClass);
}

  // Definition of the frames (through triggers) needed by the tracking algorithm

  // At present, only the following data frames can reach TrackingFilter
  // event:data      composite frame built by the EventMerger, containing
  //               - event:data:psa (composite frame containing data:psa frames)
  //               - data:rancN frames produced by AncillaryFilter (N=0,1,2)
  //                 one of the two frames can be absent
  // event:data:psa  composite frame built by the EventBuilder, containing
  //               - data:psa frames produced by the PSA filter(s)
  // data:psa        produced by PSA filters
  // data:rancN      produced by AncillaryFilters
  //
  // The trigger on event:data for the case of missing event:data:psa is implemented
  // independently so as to avoid producing a data:tracked frame on the output 
  // All triggers are implemented in the following but only (the two) event:data
  // and event:data:psa (if not using the Merger) are really in practice.
  // We keep also the other because this actor is (still) used to check
  // the integrity of data
  //
  // To be able to (re)read the final adf file (e.g. to perform a new tracking) one
  // should also implement the data:tracked frame produced by the tracking


void TrackingFilter::DefineTriggers(UInt_t *error_code)
{
  const std::list<std::pair <FactoryItem,FactoryItem> > & kf = GetConfAgent()->GetDFAgent()->GetListOfKnownFrames();
	std::list< std::pair<FactoryItem,FactoryItem> >::const_iterator lit;
  int ii = 0;
  cout << " **** List of known frames **** " << endl;
  for(lit = kf.begin(); lit != kf.end(); lit++, ii++) {
    cout << setw(3) << ii << "  " <<  lit->first << "   "  << lit->second << endl;
  }

  SharedFP *frame;
  UInt_t rerr = 0;
  trigNums = 0;
  string errMsg = "TrackingFilter::process_initialise() --> Error defining trigger ";
  string trigName;

  char chtmp[100];
  fAncRawFrameNum = fAncRawFrameNum > 9 ? 9 : fAncRawFrameNum;
  fAncRawFrameNum = fAncRawFrameNum < 0 ? 0 : fAncRawFrameNum;
  sprintf(chtmp, "data:ranc%d", fAncRawFrameNum); 
  fAncRawFrameStr = string(chtmp);
  cout << endl;

#if 1
  // Trigger on event:data  &event:data:psa |data:rancN inside  (normal output of EventMerger, with or without ancillary)
  rerr = 0;
  trigName = "event:data &event:data:psa |" + fAncRawFrameStr;
  frame = fTriggerEventData_1X.Add("Agata", "event:data",     0, false);        // keep input data
  if(!frame)
    rerr |= 8;
  frame = fTriggerEventData_1X.Add("Agata", "event:data:psa", 0, false);        // keep input data
  if(!frame)
    rerr |= 16;
  frame = fTriggerEventData_1X.Add("Agata", fAncRawFrameStr.c_str(),  0, false, "|");   // keep input data
  if(!frame)
    rerr |= 32;
  frame = fTriggerEventData_1X.SetOutputFrame("Agata", "data:tracked");
  if(!frame)
    rerr |= 4;
  if( !fFrameIO.Register(&fTriggerEventData_1X) )
    rerr |= 1;
  if(rerr) {
    cout << errMsg + trigName + " --> " << rerr << endl;
    *error_code = 100 + rerr;
    return;
  }  
  cout << "Trigger " << setw(2) << trigNums << "   " << "event:data &event:data:psa &" << fAncRawFrameStr << endl;
  nTriggerEventData_11 = trigNums++;    // event:data:psa  WITH   data:rancN
  cout << "Trigger " << setw(2) << trigNums << "   " << "event:data &event:data:psa !" << fAncRawFrameStr << endl;
  nTriggerEventData_10 = trigNums++;    // event:data:psa WITHOUT data:rancN
#endif

#if 1
  // Trigger on event:data !event:data:psa &data:rancN
  // Output of Merger without Ge component. Kept for diagnostics. No tracking, of course
  rerr = 0;
  trigName = "event:data !event:data:psa &" + fAncRawFrameStr;
  frame = fTriggerEventData_01.Add("Agata", "event:data",     0, false);   // keep input data
  if(!frame)
    rerr |= 8;
  frame = fTriggerEventData_01.Add("Agata", "event:data:psa", 0, false, "!");   // keep input data
  if(!frame)
    rerr |= 16;
  frame = fTriggerEventData_01.Add("Agata", fAncRawFrameStr.c_str(),  0, false);   // keep input data
  if(!frame)
    rerr |= 32;
  if( !fFrameIO.Register(&fTriggerEventData_01) )
    rerr |= 1;
  if(rerr) {
    cout << errMsg + trigName + " --> " << rerr << endl;
    *error_code = 100 + rerr;
    return;
  }  
  cout << "Trigger " << setw(2) << trigNums << "   " << trigName << endl;
  nTriggerEventData_01 = trigNums++;
#endif

#if 1
  // Trigger on event:data:psa (normal output of EventBuilder)
  rerr = 0;
  trigName = "event:data:psa";
  frame = fTriggerEventDataPSA.Add("Agata", "event:data:psa", 0, false);   // keep input data
  if(!frame)
    rerr |= 8;
  frame = fTriggerEventDataPSA.SetOutputFrame("Agata", "data:tracked");
  if(!frame)
    rerr |= 4;
  if( !fFrameIO.Register(&fTriggerEventDataPSA) )
    rerr |= 1;
  if(rerr) {
    cout << errMsg + trigName + " --> " << rerr << endl;
    *error_code = 100 + rerr;
    return;
  }
  cout << "Trigger " << setw(2) << trigNums << "   " << trigName << endl;
  nTriggerEventDataPSA = trigNums++;
#endif

#if 1
  // Trigger on data:psa frames (directly from PSA) , needed also to decode composite frames
  // used also to instantiate fFramePSA, which is needed to decode composite frames
  rerr = 0;
  trigName = "data:psa";
  frame = fTriggerDataPSA.Add("Agata", "data:psa", 0, true);   // discard the input data
  if(!frame)
    rerr |= 8;
  frame = fTriggerDataPSA.SetOutputFrame("Agata", "data:tracked");
  if(!frame)
    rerr |= 4;
  // UNCOMMENT if you want/need to trigger on this data
  //if( !fFrameIO.Register(&fTriggerDataPSA) )
  //  rerr |= 1;
  if(rerr) {
    cout << errMsg + trigName + " --> " << rerr << endl;
    *error_code = 100 + rerr;
    return;
  }
  cout << "Trigger " << setw(2) << trigNums << "   " << trigName << endl;
  nTriggerDataPSA = trigNums++;
#endif

#if 1
  // Trigger on raw ancillary frames (directly from AncillaryFilter)
  // used also to instantiate fFrameRaw, which is needed to decode composite frames
  rerr = 0;
  trigName = fAncRawFrameStr;
  frame = fTriggerDataRawAnc.Add("Agata", fAncRawFrameStr.c_str(), 0, true);   // discard input data
  if(!frame)
    rerr |= 8;
  // UNCOMMENT if you want/need to trigger on this data
  //if( !fFrameIO.Register(&fTriggerDataRawAnc) )
  //  rerr |= 1;
  if(rerr) {
    cout << errMsg + trigName + " --> " << rerr << endl;
    *error_code = 100 + rerr;
    return;
  }
  cout << "Trigger " << setw(2) << trigNums << "   " << trigName << endl;
  nTriggerDataRawAnc = trigNums++;
#endif

  if(trigNums > 0) {
    trigType = new UInt_t[trigNums];
    memset(trigType, 0, sizeof(UInt_t)*trigNums);
  }
  else {
    cout << "WARNING: no ADF triggers defined in TrackingFilter::DefineTriggers()" << endl;
  }

  //// The data frames fFramePSA and fFrameRaw needed to decompose composite frames
  //// They are taken from the respective trigger definitions
  fFramePSA = fTriggerDataPSA.GetInputSharedFP();
  if( fFramePSA ) {
    // link named items with local variables.
    GObject *glob =  GetDataPointer<PSAInterface>(fFramePSA)->Global();
    glob->LinkItem("CrystalID",     &crystal_id); 
    glob->LinkItem("CrystalStatus", &crystal_status); 
    glob->LinkItem("CoreE0",        &CoreE[0]); 
    glob->LinkItem("CoreE1",        &CoreE[1]);
    glob->LinkItem("CoreT0",        &CoreT[0]);
    glob->LinkItem("CoreT1",        &CoreT[1]);
  }
  else {
    rerr |= 2;
    *error_code = 100 + rerr;
    return;
  }
  fFrameRaw = fTriggerDataRawAnc.GetInputSharedFP();
  if ( fFrameRaw ) {
  }
  else {
    rerr |= 2;
    *error_code = 100 + rerr;
    return;
  }

  //// Alternatively (and certainly if one does not define those triggers) they
  //// can be generated like in the following example. Of course, in this case we
  //// have to manage the case of changed versions due to dynamical reconfiguration.
  //FactoryItem keyItem("Agata", "event:data:psa", Version(4,0));
  //FactoryItem fraItem("Agata", "event:data:psa", Version(65000,1));
  //fFrameEvb = MainFrameFactory::theMainFactory().NewSharedFrame(keyItem, fraItem);

  *error_code = 0;
}

UInt_t TrackingFilter::ProcessBlock(ADF::FrameBlock &inBlock, ADF::FrameBlock &outBlock)
{
  // attach the input/output buffer to the FrameIO system
  fFrameIO.Attach(&inBlock, &outBlock);

  // start the processing
  UInt_t error_code = 0;
  UInt_t eventsDone = cServer.GetCounts();

  if(trigNums)
    memset(trigType, 0, sizeof(UInt_t)*trigNums);
  trigUnknown = 0;

  Int_t errStatus;

  while( fFrameIO.Notify() ) { 

    // fill local structures with data from the input
    errStatus = SetInput();
    if( errStatus ) {
      if(errStatus < 0)
        continue;
      error_code = 1;
      LOCK_COUT;
      Log.SetProcessMethod("ProcessBlock");
      Log << error  << " During : SetInput()" << dolog;
      break;
    }

    // analysis of event before tracking
    errStatus = PreProcessEvent();

    if(!errStatus) {
      // tracking is done only if PreProcessEvent OK;
      errStatus = fExcludeTracking ?
                  TrackingFilter::Process() :   // if tracking excluded, force using the fake TrackingFilter::Process() 
                  Process();                    // daughter::Process()
      if( errStatus ) {
        error_code = 2;
        LOCK_COUT;
        Log.SetProcessMethod("ProcessBlock");
        Log << error  << " During : Process()" << dolog;
        break;
      }
    }
    else {
      // Event rejected by PreProcessEvent (e.g. because there are no PSA hits or out of gates)
#ifdef TRF_MULTIHIST
      if(TrackSpec_Anc && hasAncillary) { // spectrum of ancillaries with no gamma-hits
        for(unsigned int nn = 0; nn < numRawItems; nn++) {
          TrackSpec_Anc->Incr(1, nn, (int)rawBuf[nn]);
        }
      }
#endif
      // 
      //continue;   // UNCOMMENT THIS LINE TO AVOID PASSING THESE EVENTS TO THE POST-PROCESSING AND TO THE OUTPUT
    }

    // Analysis after tracking (writing the root tree and various histograms)
    errStatus = PostProcessEvent();

    cServer.Exec(timestampEvt);

    if(number_of_gammas<1 && fDiscardEmpty)
        continue;

    // move event to the output frame
    errStatus = SetOutput();
    if( errStatus ) {
      if(errStatus < 0)
        continue;
      error_code = 3;
      LOCK_COUT;
      Log.SetProcessMethod("ProcessBlock");
      Log << error  << " During : SetOutput()" << dolog;
      break;
    }

    // send the produced frame to the ouput buffer
    if( !fFrameIO.Record() ) {
      error_code = 4;
      LOCK_COUT;
      Log.SetProcessMethod("ProcessBlock");
      Log << error  << " During : Record()" << dolog;
      break;
    }
  }   //  while( fFrameIO.Notify() )

  fFrameIO.Detach(&inBlock, &outBlock);

  eventsDone = cServer.GetCounts() - eventsDone;

  LOCK_COUT;
  cServer.Prompt(-1, eventsDone, UInt_t(outBlock.GetSize()), false);
  cout << "  {";
  for(UInt_t ttt = 0; ttt < trigNums; ttt++)
    cout << " " << setw(5) << trigType[ttt];
  cout << " }";
  if( trigUnknown )
    cout << " ** " << trigUnknown << " **" ;
  cout << endl;

  // see note in TrackingFilter::process_stop
#ifdef TRF_ROOTTREE
  if(fWriteRootTree && stop_called) {
    closeTTree();
    fWriteRootTree = false;   // at least this protection against further access
  }
#endif
#ifdef TRF_WRITE_GSORT
  if(fWriteGsortData && stop_called) {
    gaspBuffer->EndOfRun();
    fWriteGsortData = false;
  }
#endif

  return error_code;
}

Int_t TrackingFilter::SetInput()
{
  timestampEvt = timestampCry = timestampAnc = 0xFFFFFFFFFFFFFFFFULL;   // preset to maximum value
  evnumberEvt  = evnumberCry  = evnumberAnc  = 0;                       // preset to 0

  number_of_crystals = 0;
  number_of_hits     = 0;
  number_of_gammas   = 0;
  fCoreSumEnergy     = 0;
  numRawItems = 0;
  len_rawBuf  = 0;
  has_rawBuf  = 0;

  hasAncillary = false;
  hasRecoil    = false;
  rVx = rVy = rVz = rVc = 0;
  aPx = aPy = aPz = 0;

  // Check which (of the alternative) trigger has fired and read the data into our local structures

  if( fTriggerEventData_1X.IsFired() ) {
    if(fDebug) {
      cout << "*** Trigger fired on event:data WITH event:data:psa and, optionally, " + fAncRawFrameStr << endl;
      PrintAllFrames(fTriggerEventData_1X);  // to understand how the composite frame is done
    }
    fAgataComposite = dynamic_cast< AgataCompositeFrame *> ( fTriggerEventData_1X.GetInputFrame(1) ); 
    fAgataComposite->Scan();
    /////////////////////////////////////////////////////////////
    int nbSubFrames = fAgataComposite->GetNbSubFrame();
    for(Int_t nsf = 0; nsf < nbSubFrames; nsf++) {
      if( fAgataComposite->LinkSubFrame(nsf, fFramePSA->GetFrame()) )
        GetDataPSA();   // extract the needed data from this PSA frame
      else
        break;          // unexpected frame; should this be treated as an error ?? 
    }
    if(fTriggerEventData_1X.IsIndividualFired(2)) {               
      trigType[nTriggerEventData_11]++;     // event:data:psa AND data:rancN
      /////////////////////////////////////////////////////////////
      fFrameRaw = fTriggerEventData_1X.GetInputSharedFP(2);
      if( fFrameRaw )
        GetDataRaw();         // extract the needed data from this PSA frame
      else
        return 18;            // ??
    }
    else {
      trigType[nTriggerEventData_10]++;     // event:data:psa WITHOUT data:rancN
    }
  }

  else if( fTriggerEventData_01.IsFired() ) {
    if(fDebug) {
      cout << "*** Trigger fired on event:data WITHOUT event:data:psa and WITH " + fAncRawFrameStr << endl;
      PrintAllFrames(fTriggerEventData_01);  // to understand how the composite frame is done
    }
    trigType[nTriggerEventData_01]++;
    if(fTriggerEventData_01.IsIndividualFired(2)) {
      /////////////////////////////////////////////////////////////
      fFrameRaw = fTriggerEventData_01.GetInputSharedFP(2);
      if( fFrameRaw )
        GetDataRaw();         // extract the needed data from this PSA frame
      else
        return 18;            // ??
    }
  }

  else if( fTriggerEventDataPSA.IsFired() ) {
    if(fDebug) {
      cout << "*** Trigger fired on event:data:psa" << endl;
      PrintAllFrames(fTriggerEventDataPSA);    // to understand how the composite frame is done
    }
    trigType[nTriggerEventDataPSA]++;
    fAgataComposite = dynamic_cast< AgataCompositeFrame *> ( fTriggerEventDataPSA.GetInputFrame(0) );
    fAgataComposite->Scan();
    /////////////////////////////////////////////////////////////
    int nbSubFrames = fAgataComposite->GetNbSubFrame();
    for(Int_t nsf = 0; nsf < nbSubFrames; nsf++) {
      if( fAgataComposite->LinkSubFrame(nsf, fFramePSA->GetFrame()) )
        GetDataPSA(); // extract the needed data from this PSA frame
      else
        break;        // unexpected frame; should this be treated as an error ?? 
    }
  }

  else if( fTriggerDataPSA.IsFired() ) {
    if(fDebug) {
      cout << "*** Trigger fired on data:psa" << endl;
      PrintAllFrames(fTriggerDataPSA);    // to understand how the composite frame is done
    }
    trigType[nTriggerDataPSA]++;
    fFramePSA = fTriggerDataPSA.GetInputSharedFP(0); 
    if( fFramePSA )
      GetDataPSA();   // extract the needed data from this PSA frame
    else
      return 17;      // ??
  }

  else if( fTriggerDataRawAnc.IsFired() ) {
    if(fDebug) {
      cout << "*** Trigger fired on " + fAncRawFrameStr << endl;
      PrintAllFrames(fTriggerDataRawAnc);    // to understand how the composite frame is done
    }
    trigType[nTriggerDataRawAnc]++;
    fFrameRaw = fTriggerDataRawAnc.GetInputSharedFP(0); 
    if( fFrameRaw )
      GetDataRaw();         // extract the needed data from this PSA frame
    else
      return 18;            // ??
  }

  else {
    cout << "an unknown trigger has fired" << endl;
    trigUnknown++;
  }

  // select common references
  // (should take the oldest timestamp not simply that of Ancillary)
  if(hasAncillary) {
    timestampEvt = timestampAnc;
    evnumberEvt  = evnumberAnc;
  }
  else {
    timestampEvt = timestampCry;
    evnumberEvt  = evnumberCry;
  }

  // correct the time information of the hits to refer to timestampEvt
  exyzHit *pt = pEXYZ;
  for(UShort_t nh = 0; nh < number_of_hits; nh++, pt++) {
    pt->T += cryst[pt->Ind].tstamp - timestampEvt;
  }
  // should we do something also for the ancillary ??

  return 0;
}

#if 1

// to understand how the composite frame is done
void TrackingFilter::PrintAllFrames(ADF::AgataFrameTrigger & fTR)
{
  UInt_t nfr0 = fTR.GetNbInputFrame();
  for(UInt_t ifr0 = 0; ifr0 < nfr0; ifr0++) {
    ADF::Frame *agataFrame = fTR.GetInputFrame(ifr0);
    if(!agataFrame) continue;
    ADF::AgataKey * akey0 = dynamic_cast<AgataKey *>(agataFrame->GetKey());
    if(!akey0) continue;
    cout << setw(2) << ifr0 << "  " <<  "    ";
    akey0->Print(); cout << "  " << AgataKeyFactory::theFactory()->GetMessage(akey0->GetMessage()) << endl;
    if( agataFrame->IsComposite() ) {
      ADF::AgataCompositeFrame *agataComposite = dynamic_cast< AgataCompositeFrame *>( agataFrame );
      if(!agataComposite) continue;
      agataComposite->Scan();
      UInt_t nfr1 = agataComposite->GetNbSubFrame();  
      for (UInt_t ifr1 = 0 ; ifr1 < nfr1 ; ifr1++) {
        const ADF::AgataKey * akey1 = dynamic_cast<const AgataKey *>(agataComposite->GetSubKey(ifr1));
        if(!akey1) continue;
        cout << setw(2) << ifr0 << "  " << setw(2) << ifr1 << "  ";
        akey1->Print(); cout << "  " << AgataKeyFactory::theFactory()->GetMessage(akey1->GetMessage()) << endl;
        //ADF::Frame *agataFrame = 
        //if( agataFrame->IsComposite() ) {
        //  ADF::AgataCompositeFrame *agataComposite = dynamic_cast< AgataCompositeFrame *>( agataFrame );
        //  agataComposite->Scan();
        //  UInt_t nfr2 = agataComposite->GetNbSubFrame();  
        //  for (UInt_t ifr2 = 0 ; ifr2 < nfr2 ; ifr2++) {
        //    const ADF::AgataKey * akey2 = dynamic_cast<const AgataKey *>(agataComposite->GetSubKey(ifr2));
        //    cout << setw(2) << ifr0 << "  " << setw(2) << ifr2 << "  ";
        //    akey2->Print(); cout << "  " << AgataKeyFactory::theFactory()->GetMessage(akey2->GetMessage()) << endl;
        //  }
        //}
      }
    }
  }
}

#else

void TrackingFilter::PrintAllFrames(ADF::AgataFrameTrigger & fTR)
{
  UInt_t nfr0 = fTR.GetNbInputFrame();
  for(UInt_t ifr0 = 0; ifr0 < nfr0; ifr0++) {
    ADF::Frame *agataFrame0 = fTR.GetInputFrame(ifr0);
    ADF::AgataKey   *agataKey0  = dynamic_cast<AgataKey *>(agataFrame0->GetKey());
    cout << setw(2) << ifr0 << "  " <<  "    " <<  "    ";
    //agataKey0->Print(); cout << "  " << agataKey0->GetSignature().GetItemName() << endl;
    agataKey0->Print(); cout << "  " << AgataKeyFactory::theFactory()->GetMessage(agataKey0->GetMessage()) << endl;
    if( agataFrame0->IsComposite() ) {
      ADF::AgataCompositeFrame *agataComposite1 = dynamic_cast< AgataCompositeFrame *>( agataFrame0 );
      agataComposite1->Scan(); UInt_t nfr1 = agataComposite1->GetNbSubFrame();  
      for (UInt_t ifr1 = 0 ; ifr1 < nfr1 ; ifr1++) {
        const ADF::AgataKey * agataKey1 = dynamic_cast<const AgataKey *>(agataComposite1->GetSubKey(ifr1));
        cout << setw(2) << ifr0 << "  " << setw(2) << ifr1  << "  " <<  "    ";
        //agataKey1->Print(); cout << "  " << agataKey1->GetSignature().GetItemName() << endl;
        agataKey1->Print(); cout << "  " << AgataKeyFactory::theFactory()->GetMessage(agataKey1->GetMessage()) << endl;
        if(!agataComposite1->LinkSubFrame(ifr1, fFrameEvb->GetFrame()))
          continue;
        if( fFrameEvb->GetFrame()->IsComposite() ) {
          ADF::AgataCompositeFrame *agataComposite2 = dynamic_cast< AgataCompositeFrame *>( fFrameEvb->GetFrame());
          agataComposite2->Scan(); UInt_t nfr2 = agataComposite2->GetNbSubFrame();  
          for (UInt_t ifr2 = 0 ; ifr2 < nfr2 ; ifr2++) {
            const ADF::AgataKey * agataKey2 = dynamic_cast<const AgataKey *>(agataComposite2->GetSubKey(ifr2));
            cout << setw(2) << ifr0 << "  " << setw(2) << ifr1 << "  " << setw(2) << ifr2 << "  ";
            //agataKey2->Print(); cout << "  " << agataKey2->GetSignature().GetItemName() << endl;
            agataKey2->Print(); cout << "  " << AgataKeyFactory::theFactory()->GetMessage(agataKey2->GetMessage()) << endl;
          }
        }
      }
    }
  }
}

#endif

void TrackingFilter::GetDataPSA()
{
  if(fFramePSA->GetFrame()->Read() == 0)
    return;

  // only if not flagged as bad in previous filters
  if(crystal_status)
    return;

  ULong64_t timestamp = ((AgataKey *)fFramePSA->GetFrame()->GetKey())->GetTimeStamp(); 
  UInt_t    evnumber  = ((AgataKey *)fFramePSA->GetFrame()->GetKey())->GetEventNumber();

  timestampCry = min(timestampCry, timestamp);
  evnumberCry  = evnumber;

  PSAInterface *data = GetDataPointer<PSAInterface>(fFramePSA);

  int nbhits = data->GetNbHits();
#ifdef TRF_FromGRU
  HitDistrib->Fill (data->GetNbHits());
#endif

  if(CoreE[0] < fDetEnerMin || CoreE[0] > fDetEnerMax)
    return;

  fCoreSumEnergy += CoreE[0];

  // loop over number of interaction points in fFramePSA
  float esumhits = 0;
  float averageT = 0;
  int   hit1st   = number_of_hits;
  exyzHit *pt = pEXYZ + number_of_hits;
  for(UShort_t nh = 0; nh < nbhits; nh++, pt++) {
    const PSAHit *ahit = dynamic_cast<PSAHit*>(data->GetHit(nh));
    //pt->Det = crystal_id;         // needed in previous versions
    pt->Det = ahit->GetID(1);
    pt->Seg = ahit->GetID(0);
    pt->Ind = number_of_crystals;   // used to recover crystal-level infos (tstamp, ...)
    pt->E  = (float)ahit->GetE(); // passed in units of keV
    pt->X = (float)ahit->GetX();
    pt->Y = (float)ahit->GetY();
    pt->Z = (float)ahit->GetZ();
    //pt->T = (float)ahit->GetT();  // not clear what this is --> replaced with common time
    pt->T = CoreT[0];               // same time for all hits; this should be in units of samples
    averageT += pt->T* pt->E;       // nonsense as long as all times are the same
    esumhits += pt->E;              // build the sum-of-segments energy

#ifdef TRF_FromGRU      
    CC3D_BefAftRT ->Fill( pt->X, pt->Y, pt->Z);
#endif

    Transform(pt, crystal_id);

#ifdef TRF_FromGRU      
    CC3D_BefAftRT ->Fill(pt->X, pt->Y, pt->Z);
#endif

  }   // loop over interaction points

  // The keyword RecalFile allows to re-scale the energy of the core and of the segments.
  // This can be used to track gain changes and to compensate for the increase of energy
  // due to the present procedure to correct the charge trapping due to neutron damage.

  // The keyword ForceSegmentsToCore forces the segment hits to sum-up to the energy of the core.
  // In this way the n-damage-corrected energy resolution is recovered to the value of CC.
  // In principle this should be done also to esumhits (i.e. esumhits=CoreE[0]), but this
  // is left unchanged to check and define the recalibration coefficient for the segments.
  // All this works well at least for low g-ray multiplicity (e.g. 60C0).
  // - A BENEFIT OF THIS TREATMENT IS A STRONG REDUCTION IN THE TRACKED SPECTRUM OF
  // THE STEP ON THE LEFT OF THE PEAKS (clearly seen in Oft__2-10-16384-UI__EE#02.spec)
  // CAUSED BY HAVING MISSED THE "VERY-LOW-AMPLITUDE" SEGMENTS
  // - A PROBLEM OF THIS TREATMENT IS THE APPEARANCE OF LEFT AND RIGHT TAILS IN THE PEAKS
  // WHEN THE HITS OF ONE CRYSTALS END UP IN DIFFERENT GAMMAS 

  //if(fRecal) {
  //  // cores
  //  float yval0 = fRecalCoeffs[crystal_id][0]; // offset
  //  float yval1 = fRecalCoeffs[crystal_id][1]; // gain
  //  CoreE[0] = yval0 + yval1*CoreE[0];
  //  CoreE[1] = yval0 + yval1*CoreE[1];
  //}

  //if(fForceSegmentsToCore) {
  //  float yval = CoreE[0]/esumhits;
  //  exyzHit *pt = pEXYZ + number_of_hits;
  //  for(UShort_t nh = 0; nh < nbhits; nh++, pt++) {
  //    pt->E *= yval;
  //  }
  //}

  number_of_hits += nbhits;
  if(esumhits > 0)
    averageT /= esumhits;
  else
    averageT  = 0;

#ifdef TRF_FromGRU
  CE_Tracking->Fill(esumhits*10);
#endif

  // save data relative to crystal
  cryst[number_of_crystals].id      = crystal_id;
  cryst[number_of_crystals].nhits   = nbhits;
  cryst[number_of_crystals].hit1st  = hit1st;
  cryst[number_of_crystals].evnum   = evnumber;
  cryst[number_of_crystals].tstamp  = timestamp;
  cryst[number_of_crystals].CCE0    = CoreE[0];
  cryst[number_of_crystals].CCE1    = CoreE[1];
  cryst[number_of_crystals].CCT0    = CoreT[0];
  cryst[number_of_crystals].CCT1    = CoreT[1];
  cryst[number_of_crystals].Esum    = (Float_t)esumhits;  // value not forced to CoreE[0]

  number_of_crystals++;
}

void TrackingFilter::GetDataRaw()
{
  ADF::RawFrame *raw_frame = dynamic_cast<ADF::RawFrame *> (fFrameRaw->GetFrame()); 

  ULong64_t timestamp = ((AgataKey *)raw_frame->GetKey())->GetTimeStamp();  
  UInt_t    evnumber  = ((AgataKey *)raw_frame->GetKey())->GetEventNumber();

  timestampAnc = min(timestampAnc, timestamp);
  evnumberAnc  = evnumber;

  if(fAncRawFrameNum == 1) {
    // load the values from the frame into the ADFObject attached to the frame
    UInt_t lenChar = raw_frame->Read();  // length of the data buffer in bytes
    //data:ranc1 is composed by a UInt_t with the number of Float_t data-items that follow
    raw_frame->RawBuffer().Export((Char_t*)&numRawItems, sizeof(UInt_t));
    if(numRawItems > 10000) {
      cout << "TrackingFilter::GetDataRaw accepts max 10000 channels" << endl;
      exit(-1); // <==== should be replaced by a "fatal" error to stop the analysis
    }
    //// realloc removed as rawBufSize is now fixed to 10000 for TF_ROOT
    //if(numRawItems > rawBufSize) {
    //  if(rawBuf) delete rawBuf;
    //  rawBufSize = numRawItems;
    //  rawBuf     = new Float_t[rawBufSize];
    //}
    len_rawBuf = numRawItems; // used for the root tree
    has_rawBuf = 1;           // used for the root tree
    raw_frame->RawBuffer().Export((Char_t*)rawBuf, numRawItems*sizeof(Float_t));
  }
  else {
    // decoding of other formats
    numRawItems = 6;  // the 6 extra values present in data:ranc1
    len_rawBuf  = 0;  // used for the root tree
    has_rawBuf  = 0;  // used for the root tree
  }

  rawBufTstamp[0] = timestamp;   // used for the root tree
  hasAncillary = true;
}

Int_t TrackingFilter::SetOutput()
{
  Frame                 *trackframe = NULL;
  GammaTrackedInterface *trackdata  = 0x0;

  if( fTriggerEventData_1X.IsFired() ) {
    trackframe = fTriggerEventData_1X.GetOutputFrame();
    trackdata  = GetDataPointer<GammaTrackedInterface>(trackframe);
  }
  else if( fTriggerEventData_01.IsFired() ) {
    trackframe = fTriggerEventData_01.GetOutputFrame(); 
    return 0; // nothing to write
  }
  else  if( fTriggerEventDataPSA.IsFired() ) {
    trackframe = fTriggerEventDataPSA.GetOutputFrame();
    trackdata  = GetDataPointer<GammaTrackedInterface>(trackframe);
  }
  else if( fTriggerDataPSA.IsFired() ) {
    trackframe = fTriggerDataPSA.GetOutputFrame();      
		trackdata  = GetDataPointer<GammaTrackedInterface>(trackframe);
  }
  else if( fTriggerDataRawAnc.IsFired() ) {
    trackframe = fTriggerDataRawAnc.GetOutputFrame();      
    return 0; // nothing to write
  }
  else {
    cout << "TrackingFilter::SetOutput not a recognized trigger" << endl;
    return -1;   // continue without writing to the output
  }
	if ( trackdata == 0x0 )
		return 1;
	
  trackframe->Reset(); // MANDATORY to start filling a new Frame to set # of gammas to 0

  ///////////////////////////////////////////////////////////////
  // evnum should be the same for all items (if valid)
  // But, which tstamp should be written out ??
  // probably all of them ==> should keep the PSA frame(s)
  // At the moment save the last obtained from the PSA frame(s) --> since 26/11/2011 use the common reference
  ///////////////////////////////////////////////////////////////

  ((AgataKey *)trackframe->GetKey())->SetEventNumber(evnumberEvt);
  ((AgataKey *)trackframe->GetKey())->SetTimeStamp(timestampEvt);

  number_of_tracked_gammas += number_of_gammas;

  // loop on the tracked gammas
  trGamma *pG = pGams;
  for(Int_t ii = 0; ii < number_of_gammas; ii++, pG++ ) {
    TrackedHit *ahit = trackdata->NewGamma();
    ahit->SetStatus(0);
    ahit->SetE(pG->E);                      // energy, passed over in units of keV 
    ahit->SetXYZ(pG->X1, pG->Y1, pG->Z1);   // 3D position of the tracked gamma is the first hit
    exyzHit *pP = pEXYZ + pG->I1;           // original PSA hit of the the first interaction point of the tracked gamma
    for(Int_t nh = 0; nh < pG->nH; nh++) {  // Decision on how many hits to write taken later in  if(nh==0) and if(pG-type!=2)
      Hit *pH = ahit->NewHit();             // BEST OF ALL WOULD BE TO WRITE ALL HITS IN THE PROPER SEQUENCE!! (OR TO USE A FLAG)
      pH->SetXYZ(pP->X, pP->Y, pP->Z);      // info taken from the original PSA hit
      pH->SetE (pP->E);
      pH->SetT (pP->T);                     // info taken from the original PSA hit
      pH->SetID(pP->Seg,0);
      pH->SetID(pP->Det,1);
      if(nh==0) {
        ahit->SetT(pP->T);                  // time of the tracked gamma taken from its first hit
        //ahit->SetID(pP->Seg, 0);    // there is no such field in TrackedHit
        //ahit->SetID(pP->Ind, 1);    //  ""
        if(pG->type!=2)                     // Photoelectric and Pair Production write only on hit. When reading back the
          break;                            // adf file they differ in Egamma_photo = Ehit vs Egamma_PP ~ Ehit + 1022
        pP = pEXYZ + pG->I2;
      }
    }
  }

  Int_t bytes_out = trackframe->Write();

  if(bytes_out == 0)
    return 1;   // error
  else
    return 0;   // normal
}

int TrackingFilter::PreProcessEvent()
{ 
  if(number_of_hits == 0)
    return 1;

  if(number_of_crystals < fDetFoldMin || number_of_crystals > fDetFoldMax)
    return 1;

  if(fCoreSumEnergy < fSumEnerMin ||  fCoreSumEnergy > fSumEnerMax)
    return 1;

  if( !AnalysisOfCrystals() )
    return 1;

  if(fAncillary && hasAncillary) {
    if(fRecoiling) {
      rVx = rawBuf[numRawItems-3 + 0]; 
      rVy = rawBuf[numRawItems-3 + 1]; 
      rVz = rawBuf[numRawItems-3 + 2];
      rVc = rVx*rVx+rVy*rVy+rVz*rVz;
      if(rVc) {
        rVc  = (float)sqrt(rVc);
        rVx /= rVc; 
        rVy /= rVc; 
        rVz /= rVc;
      }
    }
    aPx = rawBuf[numRawItems-6 + 0];
    aPy = rawBuf[numRawItems-6 + 1];
    aPz = rawBuf[numRawItems-6 + 2];
    if( !AnalysisGeAncillary() ) 
      return 1;
  }

#ifdef TFR_WRITE_HITS_MGT
  if(fWriteMgtData)
    WriteInputHits();
#endif

#ifdef TRF_ROOTTREE
  if(fWriteRootTree)
    CopyHits2Tree();
#endif

  if(fAncillary && !hasAncillary)
    return 1;

  return 0;
}

void TrackingFilter::GetParameters(UInt_t *error_code)
{
  *error_code = 0;

  string configFileName = GetConfPath() + gMotherClass + ".conf";
  ifstream configfile( configFileName.c_str() );
  if(!configfile.good()) {
    cout << endl << "Error opening " << configFileName << endl;
    *error_code = 102;
    return;
  }
  cout << endl << gMotherClass + "::GetParameters() reading from --> " << configFileName << endl;

  string line, keyw, data;
  bool ok = true;
  while(getline(configfile, line)) {
    if(!stringSplit(line, " \t", keyw, data))
      continue;       // empty or comment lines

    cout << line;
    ok = false;
    if( stringEq(keyw, "ActualClass") ) {
      ok = data.size() > 0;
      gActualClass  = data;
    }
    else if( stringEq(keyw, "SaveDataDir") || stringEq(keyw, "WriteDataDir") || stringEq(keyw, "WriteDataPrefix") ) {
      ok = data.size() > 0;
      fOdirPrefix = data;
      forceTailSlash(fOdirPrefix);
    }
    else if( stringEq(keyw, "ExcludeTracking") ) {
      fExcludeTracking = true;
      ok = true;
    }
    else if( stringEq(keyw, "Ancillary") ) {
      fAncillary = true;
      ok = true;
    }
    else if( stringEq(keyw, "AncillaryRawFrame") ) {
      ok = 1 == sscanf(data.c_str(), "%d", &fAncRawFrameNum);
    }
    else if( stringEq(keyw, "Recoiling") ) {
      fRecoiling = true;
      fAncillary = true;
      ok = true;
    }
    else if( stringEq(keyw, "RotoTranslations") ) {
      ok = data.size() > 0;
      fWhichRotoTranslations  = data;
    }
    else if( stringEq(keyw, "GeometrySummary") ) {
      ok = data.size() > 0;
      fGeometrySummary = data;
    }
    else if( stringEq(keyw, "RecoilBeta") ) {
      ok = 1 == sscanf(data.c_str(), "%f", &fRecoilVc);
    }
    else if( stringEq(keyw, "RecoilDirection") ) {
      ok = 3 == sscanf(data.c_str(), "%f %f %f", &fRecoil[0], &fRecoil[1], &fRecoil[2]);
    }
    else if( stringEq(keyw, "SourcePosition") ) {
      ok = 3 == sscanf(data.c_str(), "%f %f %f", &fSource[0], &fSource[1], &fSource[2]);
    }
    else if( stringEq(keyw, "KeepEmpty") ) {
      fDiscardEmpty = false;
      ok = true;
    }
    else if( stringEq(keyw, "DiscardEmpty") ) {
      fDiscardEmpty = true;
      ok = true;
    }
    else if( stringEq(keyw, "OutputModel") ) {
      if( stringEq(data, "kSafe") ) {
        fWriteModel    = ConfAgent::kSafe;
        fSetWriteModel = true;
        ok = true;
      }
      else if( stringEq(data, "kStrict") ) {
        fWriteModel    = ConfAgent::kStrict;
        fSetWriteModel = true;
        ok = true;
      }
      else if( stringEq(data, "kGrowing") ) {
        fWriteModel    = ConfAgent::kGrowing;
        fSetWriteModel = true;
        ok = true;
      }
      else
        ok = false;
    }
    else if( stringEq(keyw, "AcceptanceValue") ) {
      ok = 1 == sscanf(data.c_str(), "%f", &fAcceptanceValue);
    }
    else if( stringEq(keyw, "EnergyGain") ) {
      ok = 1 == sscanf(data.c_str(), "%f", &fEnergyGain);
    }
    //else if( stringEq(keyw, "RecalFile") ) {
    //  ok = data.size() > 0;
    //  fRecalFile = data;
    //}
    //else if( stringEq(keyw, "ForceSegmentsToCore") ) {
    //  fForceSegmentsToCore = true;
    //  ok = true;
    //}
    else if( stringEq(keyw, "CoreEnergyGate") ) {
      ok = 2 == sscanf(data.c_str(), "%f %f", &fDetEnerMin, &fDetEnerMax);
    }
    else if( stringEq(keyw, "SumEnergyGate") ) {
      ok = 2 == sscanf(data.c_str(), "%f %f", &fSumEnerMin, &fSumEnerMax);
    }
    else if( stringEq(keyw, "DetectorFoldGate") ) {
      ok = 2 == sscanf(data.c_str(), "%d %d", &fDetFoldMin, &fDetFoldMax);
    }
    else if(stringEq(keyw, "WriteMgtData") || stringEq(keyw, "WriteMgtHits") ||
            stringEq(keyw, "WriteOftHits") || stringEq(keyw, "WriteInputHits") ) {
#ifdef TFR_WRITE_HITS_MGT
      fWriteMgtData       = true;
      fWriteMgtAncillary  = false;
      fWriteMgtTiming     = false;
      fWriteMgtExtended   = false;
      if( stringHas(data, "Ancillary") )
        fWriteMgtAncillary = true;
      if( stringHas(data, "Timing") )
        fWriteMgtTiming = true;
      if( stringHas(data, "Extended") )
        fWriteMgtExtended = true;
#else
      cout << " --> TFR_WRITE_HITS_MGT NOT DEFINED --> ignored" << endl;
#endif
      ok = true;
    }
    else if( stringEq(keyw, "WriteRootTree") ) {
#ifdef TRF_ROOTTREE
      fWriteRootTree = true;
      fWriteRootHits = false;
      if( stringHas(data, "InputHits") ) {
        fWriteRootHits = true;
      }
#else
      cout << " --> TRF_ROOTTREE NOT DEFINED --> ignored" << endl;
#endif
      ok = true;
    }
    else if( stringEq(keyw, "WriteGsortData") ) {
#ifdef TRF_WRITE_GSORT
      fWriteGsortData = true;
#else
      cout << " --> TRF_WRITE_GSORT NOT DEFINED --> ignored" << endl;
#endif
      ok = true;
    }
    else if( stringEq(keyw, "WriteTracked") ) {
#ifdef TRF_WRITE_TRACKED
      fWriteTracked = true;
#else
      cout << " --> TRF_WRITE_TRACKED NOT DEFINED --> ignored" << endl;
#endif
      ok = true;
    }
    else if( stringEq(keyw, "TimeWindowGeAnc") || stringEq(keyw, "TimeWindowGeAncillary") ) {
      ok = 2 == sscanf(data.c_str(), "%d %d", &fTsGeAnc_min, &fTsGeAnc_max);
    }
    else if( stringEq(keyw, "Matrixgg1") ) {
      ok = 2 == sscanf(data.c_str(), "%d %f", &fSizeMatr_gg1, &fGainMatr_gg1);
    }
    else if( stringEq(keyw, "Matrixgg2") ) {
      ok = 2 == sscanf(data.c_str(), "%d %f", &fSizeMatr_gg2, &fGainMatr_gg2);
    }
    else if( stringEq(keyw, "NoMultiHist") || stringEq(keyw, "NoLocalSpectra") ) {
      fUseMultiHist = false;
      ok = true;
    }
    else if( stringEq(keyw, "Verbose") ) {
      fVerbose = true;
      ok = true;
    }
    else if( stringEq(keyw, "Debug") ) {
      fDebug = true;
      ok = true;
    }
    else {
      cout << "   --> ignored";
      ok = true;
    }

    if(!ok) {
      cout << "   --> missing argument(s)" << endl;
      *error_code = 103;
      return;
    }
    cout << endl;
  }

}

void TrackingFilter::InitGenStructures()
{
  cout << "\nTrackingFilter::InitGenStructures()" << endl;

  rotTr = (rotTr3D *)calloc(   180, sizeof(rotTr3D));
  pEXYZ = (exyzHit *)calloc(intmax, sizeof(exyzHit));
  pGams = (trGamma *)calloc(intmax, sizeof(trGamma));

  number_of_hits   = 0;
  number_of_gammas = 0;
  number_of_tracked_gammas = 0;

  rawBufSize = 10000;
  rawBuf     = new Float_t[rawBufSize]; //maximum number of channels in ancillaries
  memset(rawBuf,    0, sizeof(Float_t)*rawBufSize);

  len_rawBuf = 0;
  has_rawBuf = 0;

#ifdef CHECKRAWBUF
  rawBufOld = new Float_t[rawBufSize];  //maximum number of channels in ancillaries
  memset(rawBufOld, 0, sizeof(Float_t)*rawBufSize);
#endif

#ifdef TRF_ROOTTREE
  gammaE       = new Float_t[intmax];
  gammaT       = new Float_t[intmax];
  gammaTstamp  = new ULong64_t[intmax];
  gammaTrType  = new Int_t[intmax];
  gammaX1      = new Float_t[intmax];
  gammaY1      = new Float_t[intmax];
  gammaZ1      = new Float_t[intmax];
  gammaX2      = new Float_t[intmax];
  gammaY2      = new Float_t[intmax];
  gammaZ2      = new Float_t[intmax];
  hitE         = new Float_t[intmax];
  hitX         = new Float_t[intmax];
  hitY         = new Float_t[intmax];
  hitZ         = new Float_t[intmax];
  hitT         = new Float_t[intmax];
  hitTstamp    = new ULong64_t[intmax];
  hitId        = new Int_t[intmax];
  hitSg        = new Int_t[intmax];
#endif // TRF_ROOTTREE

#ifdef TRF_MULTIHIST
  if(fUseMultiHist) {
    TrackSpec_Fold = new MultiHist<unsigned int>(4, NUMGEDETS);
    TrackSpec_Fold->setFileName(fOdirPrefix+"Track?Fold.spec");
    TrackSpec_Fold->setComment("multiplicity of crystals");
    hGroup.add(TrackSpec_Fold);

    TrackMatr_Fold = new MultiHist<unsigned int>(NUMGEDETS, NUMGEDETS);
    TrackMatr_Fold->setFileName(fOdirPrefix+"Track?Fold.matr");
    TrackMatr_Fold->setComment("number of gammas vs number of crystals");
    hGroup.add(TrackMatr_Fold);

    TrackSpec_EC = new MultiHist<unsigned int>(2, NUMGEDETS, specLenE);
    TrackSpec_EC->setFileName(fOdirPrefix+"Track?EC.spec");
    TrackSpec_EC->setComment("core and sumSegs energy sorted by detector");
    hGroup.add(TrackSpec_EC);

    TrackSpec_ES = new MultiHist<unsigned int>(1, NUMGEDETS, specLenE);
    TrackSpec_ES->setFileName(fOdirPrefix+"Track?ES.spec");
    TrackSpec_ES->setComment("sum-energy spectra (gain 1 keV/chan)");
    hGroup.add(TrackSpec_ES);

    TrackSpec_EE = new MultiHist<unsigned int>(2, 10, specLenE);
    TrackSpec_EE->setFileName(fOdirPrefix+"Track?EE.spec");
    TrackSpec_EE->setComment("input and tracked spectra");
    hGroup.add(TrackSpec_EE);

    // no meaning if the integer part of the trigger-point has been moved to timestamp
    TrackSpec_TC = new MultiHist<unsigned int>(NUMGEDETS, specLenT1);
    TrackSpec_TC->setFileName(fOdirPrefix+"Track?TC.spec");
    TrackSpec_TC->setComment("CFD trigger relative to beginning of trace (ns)");
    hGroup.add(TrackSpec_TC);

    //TrackMatr_ETC = new MultiHist<unsigned short>(1000, NUMGEDETS*1024);
    //TrackMatr_ETC->setFileName(fOdirPrefix+"Track?ETC.matr");
    //TrackMatr_ETC->setComment("gamma-Trigger");
    //hGroup.add(TrackMatr_ETC);

    TrackSpec_TS = new MultiHist<unsigned int>(NUMGEDETS, NUMGEDETS, specLenT1);
    TrackSpec_TS->setFileName(fOdirPrefix+"Track?TS.spec");
    TrackSpec_TS->setComment("T_gamma-gamma from timestamp (samples)");
    hGroup.add(TrackSpec_TS);

    TrackSpec_TT = new MultiHist<unsigned int>(NUMGEDETS, NUMGEDETS, specLenT1);
    TrackSpec_TT->setFileName(fOdirPrefix+"Track?TT.spec");
    TrackSpec_TT->setComment("T_gamma-gamma from CFD (ns)");
    hGroup.add(TrackSpec_TT);

    TrackSpec_TG = new MultiHist<unsigned int>(2, 4*specLenT1);
    TrackSpec_TG->setFileName(fOdirPrefix+"Track?TG.spec");
    TrackSpec_TG->setComment("Tgamma and Tgg of tracked gammas");
    hGroup.add(TrackSpec_TG);

    //TrackMatr_IEET = new MultiHist<unsigned char>(NUMGEDETS, 100, 100, specLenT2);
    //TrackMatr_IEET->setFileName(fOdirPrefix+"Track?IEET.matr");
    //TrackMatr_IEET->setComment("det_index-energy-energy-time");
    //hGroup.add(TrackMatr_IEET);

    if(fSizeMatr_gg1 > 0) {
      TrackMatr_gg1 = new MultiHist<unsigned short>(fSizeMatr_gg1, fSizeMatr_gg1);
      TrackMatr_gg1->setFileName(fOdirPrefix+"Track?gg1.matr");
      TrackMatr_gg1->setComment("gamma-gamma of untracked data");
      hGroup.add(TrackMatr_gg1);
    }

    if(fSizeMatr_gg2 > 0) {
      TrackMatr_gg2 = new MultiHist<unsigned short>(fSizeMatr_gg2, fSizeMatr_gg2);
      TrackMatr_gg2->setFileName(fOdirPrefix+"Track?gg2.matr");
      TrackMatr_gg2->setComment("gamma-gamma of tracked and doppler-corrected data");
      hGroup.add(TrackMatr_gg2);
    }

    TrackMatr_SGSG = new MultiHist<unsigned short>(36*NUMGEDETS, 36*NUMGEDETS);
    TrackMatr_SGSG->setFileName(fOdirPrefix+"Track?SGSG.matr");
    TrackMatr_SGSG->setComment("segment-segment coincidences for all pairs");
    hGroup.add(TrackMatr_SGSG);

    if(fAncillary) {
      TrackSpec_EA = new MultiHist<unsigned int>(2, NUMGEDETS, specLenE); // 0 alone, 1 Ancillary
      TrackSpec_EA->setFileName(fOdirPrefix+"Track?EA.spec");
      TrackSpec_EA->setComment("core energy sorted by detector, in coincidence with ancillary");
      hGroup.add(TrackSpec_EA);

      TrackSpec_TA = new MultiHist<unsigned int>(2, NUMGEDETS, specLenT1);
      TrackSpec_TA->setFileName(fOdirPrefix+"Track?TA.spec");
      TrackSpec_TA->setComment("Difference of tstamps between crystals and ancillary");
      hGroup.add(TrackSpec_TA);

      TrackSpec_Anc = new MultiHist<unsigned int>(2, 256, 4096);
      TrackSpec_Anc->setFileName(fOdirPrefix+"Track?Anc.spec"); 
      TrackSpec_Anc->setComment("Projection of all Ancillary Parameters");
      hGroup.add(TrackSpec_Anc);

      //TrackSpec_TOF = new MultiHist<unsigned int>(NUMGEDETS, 32, specLenT1);
      //TrackSpec_TOF->setFileName(fOdirPrefix+"Track?TOF.spec");
      //TrackSpec_TOF->setComment("Timestamp-corrected time spectra crystals-ancillary");
      //hGroup.add(TrackSpec_TOF);

      //TrackMatr_TOF = new MultiHist<unsigned short>(NUMGEDETS, 32, 400, 200);
      //TrackMatr_TOF->setFileName(fOdirPrefix+"Track?TOF.matr");
      //TrackMatr_TOF->setComment("Timestamp-corrected Egamma-TOF");
      //hGroup.add(TrackMatr_TOF);

      //TrackMatr_AG = new MultiHist<unsigned short>(96, 256, 2048);
      //TrackMatr_AG->setFileName(fOdirPrefix+"Track?AG.matr");
      //TrackMatr_AG->setComment("num-Eanc-Egamma");
      //hGroup.add(TrackMatr_AG);

      //TrackSpec_TTA = new MultiHist<unsigned short>(96, 1024, specLenT1);
      //TrackSpec_TTA->setFileName(fOdirPrefix+"Track?TTA.matr");
      //TrackSpec_TTA->setComment("T_gamma-anc from CFD (ns)");
      //hGroup.add(TrackSpec_TTA);

      //TrackMatr_ETA = new MultiHist<unsigned short>(4096, 1024);
      //TrackMatr_ETA->setFileName(fOdirPrefix+"Track?EgTa-g.matr");
      //TrackMatr_ETA->setComment("E(gamma)--T(gamma_tstamp-ancillary_tstamp)");
      //hGroup.add(TrackMatr_ETA);
    }
  }
#endif //#ifdef TRF_MULTIHIST

  // for Doppler correction
  if(fRecoilVc > 0 && fRecoilVc < 1.) {
    RecoilVc = fRecoilVc;
    RecoilGm = float( 1. / sqrt(1.f - RecoilVc*RecoilVc) );
    float dd = (fRecoil[0]*fRecoil[0] + fRecoil[1]*fRecoil[1] + fRecoil[2]*fRecoil[2]);
    if(dd < 1.e-6) {
      cout << "Size of direction cosines vector is too small ==> recoil direction set along z" << endl;
      RecoilCx = RecoilCy = 0;
      RecoilCz = 1.;
    }
    else {
      dd = sqrt(dd);
      RecoilCx = fRecoil[0]/dd;
      RecoilCy = fRecoil[1]/dd;
      RecoilCz = fRecoil[2]/dd;
    }
  }
  else {
    RecoilVc = 0;
    RecoilGm = 1.;
    RecoilCx = 0;
    RecoilCy = 0;
    RecoilCz = 1.;
  }
}

int TrackingFilter::GetRotoTranslations(FILE *fp)
{
  int nd, k;
  int nt = 0;
  float x, y, z;

  while(true) {
    // detector number and translation
    if(fscanf(fp, "%d %d %f %f %f", &nd, &k, &x, &y, &z) != 5)
      break;

    if(nd < 0 || nd >= 180 || k != 0)
      return -1;

    rotTr3D *prt = rotTr + nd;

    prt->TrX = x; // in mm
    prt->TrY = y;
    prt->TrZ = z;

    // matrix elements of rotation
    if(fscanf(fp, "%d  %f %f %f", &k, &x, &y, &z) != 4)
      return -1;
    if(k != 1)
      return -1;

    prt->rXX = x;
    prt->rXY = y;
    prt->rXZ = 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;

    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;

    nt++;
  }
  return nt;
}

//bool TrackingFilter::GetRecalCoeffs()
//{
//  std::string conf_file = GetConfPath() + fRecalFile;
//  cout << "Opening energy recalibration coefficients file " << conf_file << endl;
//  FILE *fp_conf = fopen(conf_file.c_str(), "r");
//  if(!fp_conf) {
//    cout << "Could not open file " << conf_file << endl;
//    return false;
//  }
//
//  for(int ii = 0; ii < NUMGEDETS; ii++) {
//    fRecalCoeffs[ii][0] = 0;
//    fRecalCoeffs[ii][1] = 1.f;
//  }
//
//  int nt = 0;
//
//  while(true) {
//    int   nd;
//    float xx, yy;
//    // detector number and (linear) recalibration coefficients
//    if(fscanf(fp_conf, "%d %f %f", &nd, &xx, &yy) != 3)
//      break;
//
//    if(nd < 0 || nd >= NUMGEDETS)
//      return false;
//
//    fRecalCoeffs[[nd]0] = xx; // offset
//    fRecalCoeffs[nd][1] = yy; // gain
//    nt++;
//  }
//  cout << "Found coefficients for " << nt << " detector(s)" << endl;
//  fRecal = true;
//  fclose(fp_conf);
//  return true;
//}

// recoil along a fixed direction (beam)

#ifdef TFR_WRITE_HITS_MGT

bool TrackingFilter::OpenInputHits()
{
  // opening file to write the hits at the input of OFT in mgt format
  hits_file = fOdirPrefix + "Mgt_Hits.txt";

  fp_hits = fopen(hits_file.c_str(), "wb");
  if(!fp_hits) {
    cout << "Error opening " << hits_file << endl;
    return false;
  }
  cout << "Input hits will be written to " << hits_file << endl;
  fprintf(fp_hits, "AGATA 6.6.6\n");

  // Enrico's definition of OUTPUT_MASK 0--> disabled; 1--> enabled
  // 0  nDet            (leftmost character)
  // 1  Energy
  // 2  AbsolutePosition(x y z)
  // 3  RelativePosition(x' y' z') 
  // 4  RelativePosition for mgs (x'' y'' z'')
  // 5  nSeg
  // 6  time
  // 7  interaction    (rightmost character)  (for us this is timeStamp+eventNumber)
  fprintf(fp_hits, "OUTPUT_MASK 111001");
  if(fWriteMgtTiming)
    fprintf(fp_hits, "1");
  else
    fprintf(fp_hits, "0");
  if(fWriteMgtExtended)
    fprintf(fp_hits, "1");
  else
    fprintf(fp_hits, "0");
  fprintf(fp_hits, "\n");

  fprintf(fp_hits, "DISTFACTOR 1\n");
  fprintf(fp_hits, "ENERFACTOR 1\n");
  fprintf(fp_hits, "GEOMETRY\n");
  fprintf(fp_hits, "AGATA\n");
  fprintf(fp_hits, "SUMMARY %s\n", fGeometrySummary.c_str());
  fprintf(fp_hits, "ENDGEOMETRY\n");
  fprintf(fp_hits, "$\n");
  //size_t lposi = ftell(fp_hits);
  return true;
}

void TrackingFilter::WriteInputHits()
{
  if(!fp_hits)
    return;  // file was not opened

  if(fAncillary) {
    if(!hasAncillary)
      return;
    if(!hasRecoil)
      return;
    if(fRecoiling && rVc <=0.)
      return;
  }

  float etot = 0;
  for(int ii = 0; ii < number_of_hits; ii++)
    etot += pEXYZ[ii].E;

  // first come the Ancillary, if present
  if(fAncillary) {
    fprintf(fp_hits, "-101 ");  // for mgt the line tagged -101 defines the recoil vector
    if(fWriteMgtAncillary) {
      // write the complete ancillary buffer with of without zero suppression
#if 1      // zero suppression 
      for (int i = 0; i < 96; i++) {    // why 96 ???
        if (rawBuf[i] > 0.)
          fprintf(fp_hits, "%d %.5f ", i, rawBuf[i]);
      }
#else
      for (int i = 0; i < 96; i++) fprintf(fp_hits, "%.5f ", rawBuf[i]);
#endif
    }
    else {
      // write the recoil vector and the xyz coordinates in the frame of Dante to be able to re-do the kinematics (?? what if no Dante ??)
      fprintf(fp_hits, "%.5f %.5f %.5f %.5f %.2f %.2f %.2f", rVc, rVx, rVy, rVz, aPx, aPy, aPz);
    }
    fprintf(fp_hits, "\n");
  }

  fprintf(fp_hits, "-1 %.2f 0 0 0 %d\n", etot, evnumberCry);

  exyzHit * pt = pEXYZ;

  for(int nn = 0; nn < number_of_hits; nn++, pt++) {
    UInt_t icry = pt->Ind;
    fprintf(fp_hits, "%d %.2f %.2f %.2f %.2f %d", pt->Det, pt->E, pt->X, pt->Y, pt->Z, pt->Seg);
    if(fWriteMgtTiming) {
      float timeseg = cryst[icry].CCT0 + float(cryst[icry].tstamp - timestampEvt);
      //for (int icry = 0; icry < number_of_crystals; icry++) {
      //  if (cryst[icry].id == pt->Id) {
      //    timeseg += (Float_t)cryst[icry].CCT0;
      //    timeseg += cryst[icry].tstamp - timestampAnc;
      //    timeseg *= 10;
      //  }
      //}
      fprintf(fp_hits, " %.1f", timeseg*10.f);
    }
    if(fWriteMgtExtended)
      fprintf(fp_hits, " %lld %d", cryst[pt->Ind].tstamp, cryst[pt->Ind].evnum);
    fprintf(fp_hits, "\n");
  }
}

#endif // #ifdef TFR_WRITE_HITS_MGT

int TrackingFilter::PostProcessEvent()
{
#ifdef CHECKRAWBUF
  // a test to debug repetitions
  if(rawBufOld) {
    int nsame = 0;
    for(unsigned int nn = 0; nn < numRawItems; nn++) {
      if(rawBuf[nn]>100.f && rawBufOld[nn] == rawBuf[nn])
        nsame++;
    }
    memcpy(rawBufOld, rawBuf, numRawItems*sizeof(Float_t));
  }
#endif

#ifdef TRF_ROOTTREE
  if(fWriteRootTree) {
    CopyGamma2Tree(pGams, number_of_gammas);
    t->Fill();
  }
#endif
#ifdef TRF_WRITE_GSORT
  if(fWriteGsortData) {
    PrepareGaspEvent();      
    gaspBuffer->AddToBuffer( gaspEvent, gaspEvLen+1 );
  }
#endif

#ifdef TRF_MULTIHIST
  if(TrackSpec_Fold)
    TrackSpec_Fold->Incr(2, number_of_gammas);
  if(TrackMatr_Fold)
    TrackMatr_Fold->Incr(number_of_crystals, number_of_gammas);
  if(TrackSpec_Anc && fAncillary && numRawItems) {
    for(unsigned int nn = 0; nn < numRawItems; nn++) {
      TrackSpec_Anc->Incr(0, nn, (int)rawBuf[nn]);
    }
  }
#endif

  if(!number_of_gammas)
    return 0;   // nothing more to do if there are non gammas

  // Questo dovrebbe migliorare la qualitą della ricostruzione va fatto meglio, verificando
  // che gli hit di un cristallo non appartengano a pił di un gamma.
  // Bisogna comunque valutare il costo in efficienza

  if(number_of_crystals && number_of_gammas>number_of_crystals)
    return 0; 

  trGamma *pg = pGams;

  // loop on tracked gammas in the event
  for(Int_t i = 0; i < number_of_gammas; i++, pg++ ) {
#ifdef TRF_WRITE_TRACKED
    if(fWriteTracked)  {
      // write tracked gammas to file (in keV mm)
      fprintf(fp_tracked, "%.2f %.1f %.1f %.1f %.1f %.1f %.1f\n",
        pg->E, pg->X1, pg->Y1, pg->Z1, pg->X2, pg->Y2, pg->Z2);
    }
#endif
#ifdef TRF_MULTIHIST
    if(TrackSpec_EE) {
      TrackSpec_EE->Incr(0, 2,                   int(fEnergyGain*pg->E));   // 0-2 Tracked
      TrackSpec_EE->Incr(0, 2 + pg->type,        int(fEnergyGain*pg->E));   // 0-3 Photo=1   0-4 Compt=2 0-5 pair=3
      if(hasRecoil) {
        TrackSpec_EE->Incr(1, 2,                 int(fEnergyGain*pg->E));   // 1-2 Tracked in coinc with ancillary
        TrackSpec_EE->Incr(1, 2 + pg->type,      int(fEnergyGain*pg->E));   // 0-3 Photo=1   0-4 Compt=2 0-5 pair=3
      }
    }
#endif  //TRF_MULTIHIST
#ifdef TRF_FromGRU
    if(pg->nH > 1)
      TrackedEnergies->Fill(fEnergyGain*pg->E);
#endif

    // Doppler correction

    //float Ecorr = fAncillary ? DopplerCorrectAnc(pg) : DopplerCorrectFix(pg); // not using function pointer
    float Ecorr = pg->Ecorr = (this->*DopplerCorrect)(pg);    // using the function pointer

    if(Ecorr > 0) {
#ifdef TRF_MULTIHIST
      if(TrackSpec_EE) {
        TrackSpec_EE->Incr  (0, 6, int(fEnergyGain*Ecorr)); // 0-6 Doppler corrected
        if(hasRecoil)
          TrackSpec_EE->Incr(1, 6, int(fEnergyGain*Ecorr)); // 1-6 Doppler corrected in coincidence with ancillary
      }
#endif  //TRF_MULTIHIST
#ifdef TRF_FromGRU
      TrackedEnergiesCorr->Fill(fEnergyGain*Ecorr);
#endif
    }

  }  // loop over the gammas

#ifdef TRF_MULTIHIST
  // time of the individual gammas and time between gammas using only info from the output
  if(TrackSpec_TG) {
    for(int i1 = 0; i1 < number_of_gammas-1; i1++) {
      Double_t gamTime1 = pGams[i1].T;
      TrackSpec_TG->Incr(0, int(gamTime1));   // time of the individual gammas (re global time reference ??) 
      for(int i2 = i1+1; i2 < number_of_gammas; i2++) {
        Double_t gamTime2 = pGams[i2].T;
        TrackSpec_TG->Incr(1, int(gamTime2-gamTime1 + 2*specLenT1));  // time between gammas
      }
    }
  }
  // gamma-gamma symmetrized matrix of tracked and doppler-corrected data
  if(TrackMatr_gg2 && number_of_gammas > 1) {
    for(int i1 = 0; i1 < number_of_gammas-1; i1++) {
      int ie1 = int(fGainMatr_gg2*pGams[i1].Ecorr);
      for(int i2 = i1+1; i2 < number_of_gammas; i2++) {
        int ie2 = int(fGainMatr_gg2*pGams[i2].Ecorr);
        TrackMatr_gg2->Incr(ie1, ie2);
        TrackMatr_gg2->Incr(ie2, ie1);  // matrix is symmetrized
      }
    }
  }
#endif

  return 0;
}

// time and energy spectra for singles and gamma-gamma
bool TrackingFilter::AnalysisOfCrystals()
{
#ifdef TRF_MULTIHIST
  int   dts = 0;
  float dtf = 0;
  //float dtp = 0;
  float dtt = 0;
  const float tr0 = 0.5f*specLenT1;
  const float tr1 = 0.5f*specLenT2;

  if(TrackSpec_Fold)
    TrackSpec_Fold->Incr(0, number_of_crystals);
  float esumCC = 0;
  for(int nn1 = 0; nn1 < number_of_crystals; nn1++) {
    crystdata *pCr1 = cryst + nn1;
    int id1 = pCr1->id;
    if(id1 < 0 || id1 >= NUMGEDETS)   // this should not happen and should be handled as an error ??
      continue;
    esumCC += pCr1->CCE0;
    if(TrackSpec_EC) {
      TrackSpec_EC->Incr(0, id1, int(fEnergyGain*pCr1->CCE0));     // energy of the core sorted by detector
      TrackSpec_EC->Incr(1, id1, int(fEnergyGain*pCr1->Esum));     // sum energy of the segments
    }
    if(TrackSpec_EE) {
      TrackSpec_EE->Incr(0,   0, int(fEnergyGain*pCr1->CCE0));     // 0-0 energy of the core common spectrum
      TrackSpec_EE->Incr(0,   1, int(fEnergyGain*pCr1->Esum));     // 0-1 energy of the sumsegs common spectrum
    }
    if(TrackSpec_TC)
      TrackSpec_TC->Incr(   id1, int(tScale*pCr1->CCT0));         // trigger point of the input traces from PreprocessingFilterPSA
    //if(TrackMatr_ETC)
    //  TrackMatr_ETC->Incr(int(pCr1->CCE0/10), id1*1024+int(tScale*pCr1->CCT0));     // energy-trigger point from preprocessing
    int seg1 = 36*id1 + pEXYZ[pCr1->hit1st  ].Seg;
    if(TrackMatr_SGSG && pCr1->nhits==2) {
      int seg2 = 36*id1 + pEXYZ[pCr1->hit1st+1].Seg;
      TrackMatr_SGSG->Incr(seg1, seg2);
      TrackMatr_SGSG->Incr(seg2, seg1);
    }
    if(pCr1->nhits!=1)
      seg1 = -1;
    for(int nn2 = 0; nn2 < number_of_crystals; nn2++) {
      if(nn2 == nn1)
        continue;
      crystdata *pCr2 = cryst + nn2;
      int id2 = pCr2->id;
      if(id2 < 0 || id2 >= NUMGEDETS)   // this should not happen and should be handled as an error ??
        continue;
      dts = (int)(pCr2->tstamp - pCr1->tstamp);
      dtf =       pCr2->CCT0   - pCr1->CCT0;
      dtt =       tScale*(dts + dtf) ;
      if(TrackSpec_TS)
        TrackSpec_TS->Incr(id1, id2, dts+specLenT1/2);      // time-time from tstamp
      if(TrackSpec_TT)
        TrackSpec_TT->Incr(id1, id2, int(tr0+dtt) );        // time-time from tstamp+CFD
      //if(TrackMatr_IEET)
      //  TrackMatr_IEET->Incr(nn1, int(pCr1->CCE0*fEnergyGainET), int(pCr2->CCE0*fEnergyGainET), int(tr1+dtt));
      if(TrackMatr_gg1) {
        // gamma-gamma symmetrized matrix CC
        int ie1 = int(fGainMatr_gg1*pCr1->CCE0);
        int ie2 = int(fGainMatr_gg1*pCr2->CCE0);
        TrackMatr_gg1->Incr(ie1, ie2);
        TrackMatr_gg1->Incr(ie2, ie1);
      }
      if(TrackMatr_SGSG && pCr2->nhits==1 && seg1>=0) {
        int seg2 = 36*id2 + pEXYZ[pCr2->hit1st].Seg;
        TrackMatr_SGSG->Incr(seg1, seg2);
      }
    }
  }
  if(TrackSpec_ES && number_of_crystals > 0 && number_of_crystals < NUMGEDETS) {
    TrackSpec_ES->Incr(0,                  0, int(esumCC));     // total sum energy
    TrackSpec_ES->Incr(0, number_of_crystals, int(esumCC));     // sum energy [fold]
  }

#endif  //TRF_MULTIHIST
  return true;
}

// timestamp difference between ancillary and Ge
// to check the timegate  [fTsGeAnc_min -- fTsGeAnc_max]
bool TrackingFilter::AnalysisGeAncillary()
{
  if(!fAncillary)
    return false;   // should not be called if presence of ancillary has not been requested

  if(!hasAncillary)
    return false;   // no ancillary is present in this event

#ifdef TRF_MULTIHIST
  if(TrackSpec_Fold)
    TrackSpec_Fold->Incr(1, number_of_crystals);
#endif // TRF_MULTIHIST

  for(int nn1 = 0; nn1 < number_of_crystals; nn1++) {
    crystdata *pCr1 = cryst + nn1;
    int ii1 = pCr1->id;
    if(ii1 < 0 || ii1 >= NUMGEDETS) continue;           // this should not happen
    int eG = int(fEnergyGain*pCr1->CCE0);
    int eS = int(fEnergyGain*pCr1->Esum);
    int dT = specLenT1/2;   // timestamp difference calculated (and used) with zero at specLenT1/2
    // this way to have the ULL diff positive     
    if(pCr1->tstamp > timestampAnc)
      dT += int(pCr1->tstamp - timestampAnc);
    else
      dT -= int(timestampAnc - pCr1->tstamp);
    bool isInTime = (dT >= fTsGeAnc_min && dT <= fTsGeAnc_max);
    hasRecoil |= isInTime;

#ifdef TRF_MULTIHIST
    if(TrackSpec_EA)
      TrackSpec_EA->Incr(0, ii1, eG);         // energy of the core sorted by detector
    if(TrackSpec_TA)
      TrackSpec_TA->Incr(0, ii1, dT);         // tstamp difference ge-Anc

    if(TrackMatr_ETA && (pCr1->tstamp > timestampAnc) ) {
      int dT = int(pCr1->tstamp - timestampAnc);
      int eG = int(pCr1->CCE0);
      TrackMatr_ETA->Incr(eG, dT);
    }

#if 0
    //////// leftover from a Dante or Prisma experiment ??? //////////////
    // TOF with respect to the ancillary
    // The TDC is the 3rd converter
    // The discretization due to GTS/AGAVA is already corrected in AncillaryFilterVME
    if(TrackSpec_TOF) {
      //if(pCr1->CCE0 > 500.f) { // energy threshold on the germanium energy
      for(unsigned int nns = 64; nns < 96; nns++) {   // only the TDCs
        float valAnc = rawBuf[nns];
        if(valAnc > 10.f) {
          valAnc *= .1f; // ns->samples
          float dtf = pCr1->CCT0 - valAnc;
          float dtt = tScale*(dT + dtf) ;
          dtt -= 7000.f;     // ~ tScale*(specLenT1/2 + 180)  with 180=delay of Ge relative to Anc  
          TrackSpec_TOF->Incr(ii1, nns-64, int(dtt));
          if(TrackMatr_TOF) {
            dtt -= 200;      // adjust to a proper value
            float dee = pCr1->CCE0/10;   // steps of 10 keV
            TrackMatr_TOF->Incr(ii1, nns-64, int(dee), int(dtt));
          }
        }
      }
      //}
    }
#endif

#if 0
    //  TOF with respect to the ancillary
    if(TrackSpec_TTA) {
      if(pCr1->CCE0 > 500.f) {                 // energy threshold on the germanium energy
        const float t0 = 0.5f*specLenT1;
        float dtt  = tScale*(dT-t0 + pCr1->CCT0);
        dtt -= 1600.f; // to get it ~middle of the range
        unsigned int nnmx = TrackSpec_TTA->getLenSpec(0);
        nnmx = nnmx > numRawItems ? numRawItems : nnmx;
        for(unsigned int nns = 0; nns < nnmx; nns++) {
          float valAnc = rawBuf[nns];
          if(valAnc < 20.f) continue;
          TrackSpec_TTA->Incr(nns, int(valAnc), int(dtt));
        }
      }
    } 
#endif

    // ge-Ancillary in the time window
    if(isInTime) {
      if(TrackSpec_EA)
        TrackSpec_EA->Incr(1, ii1, eG);                 // coinc. energy of the core
      if(TrackSpec_TA)
        TrackSpec_TA->Incr(1, ii1, dT);                 // tstamp difference ge-Anc
      //if( rVc ) {
        if(TrackSpec_EE) {
          TrackSpec_EE->Incr(1, 0, eG);                 // 1-0 total energy spectrum of the core
          TrackSpec_EE->Incr(1, 1, eS);                 // 1-1 energy of the sumsegs common spectrum
        }
      //}
      if(TrackMatr_AG) {
        eG = int(pCr1->CCE0*0.5f);            // fixed gain 2 keV/chan
        if(eG > 0 && eG < 2048) {
          unsigned int nnmx = TrackMatr_AG->getLenSpec(0);
          nnmx = (nnmx > numRawItems) ? numRawItems : nnmx;
          for(unsigned int nns = 0; nns < nnmx; nns++) {
            if(rawBuf[nns] > 200 && rawBuf[nns] < 4000) { // adjusted with proper values
              TrackMatr_AG->Incr(nns, int(rawBuf[nns]/16.f), eG);
            }
          }
        }
      }
    }
#endif  //TRF_MULTIHIST
  }
  return true;
}

float TrackingFilter::DopplerCorrectFix(trGamma *pg)
{
  if(!RecoilVc)
    return pg->E;

  float dx = pg->X1 - fSource[0];
  float dy = pg->Y1 - fSource[1];
  float dz = pg->Z1 - fSource[2];
  float dd = dx*dx + dy*dy + dz*dz;
  if(dd == 0)
    return 0;

  float cosTheta = (dx*RecoilCx + dy*RecoilCy + dz*RecoilCz)/sqrt(dd);
  float cDoppler = RecoilGm*(1.f - RecoilVc * cosTheta);
  return cDoppler*pg->E;
}

// recoil direction taken from ancillary
float TrackingFilter::DopplerCorrectAnc(trGamma *pg)
{
  if(!rVc)
    return pg->E;

  float dx = pg->X1 - fSource[0];
  float dy = pg->Y1 - fSource[1];
  float dz = pg->Z1 - fSource[2];
  float dd = dx*dx + dy*dy + dz*dz;
  if(dd == 0)
    return 0;

  float cosTheta = (dx*rVx + dy*rVy + dz*rVz)/sqrt(dd);
  float recGamma = sqrt(1.f - rVc*rVc);
  float cDoppler = (1.f - rVc * cosTheta)/recGamma;
  return cDoppler*pg->E;
}

void TrackingFilter::process_reset (UInt_t *error_code)
{
  *error_code = 0;
  Log.ClearMessage();
  Log.SetProcessMethod("process_reset"); 
  fFrameIO.Print(Log()); 

  NarvalFilter::process_reset(error_code);

  cServer.Reset();

  Log << dolog;
}

void TrackingFilter::process_start (UInt_t *error_code)
{
  Log.GetProcessName() = "TrackingFilter";
  Log.SetProcessMethod("process_start"); 
  cServer.Start(gMotherClass);

  stop_called = false;
  *error_code = 0;
  Log << dolog;
}

void TrackingFilter::process_stop (UInt_t *error_code)
{
  cout << "TrackingFilter::process_stop called with GetPID() " << GetPID() << endl;
  *error_code = 0;


  // Remember that in narval this is an asynchronous call and that process_block will be
  // called again and again until the intermediate buffers are completely consumed.
  // So, process_stop should simply take note of having been called and nothing more.
  // Actions like closing files (e.g. closeTTree) are likely to end up in a program-crash
  // because the rest of the filter might still need to access them. 

  stop_called = true;

#if NRV_TYPE == NRV_ONLINE
  return;
#endif

  cServer.Finish();
#ifdef TRF_ROOTTREE
  if(fWriteRootTree) {
    closeTTree();
    fWriteRootTree = false;   // at least this protection against further access
  }
#endif

#ifdef TRF_WRITE_GSORT
  if(fWriteGsortData) {
    gaspBuffer->EndOfRun();
    fWriteGsortData = false;
  }
#endif
}

#ifdef TRF_ROOTTREE
void TrackingFilter::openTTree ()
{
  string fn1;
  fn1 = fOdirPrefix + "AGATREE.root";

  rootfile = new TFile(fn1.c_str(),"RECREATE");
  rootfile->cd();

  t = new TTree("Events","Events");
  t->SetAutoSave(1000000);
  //t->SetMaxTreeSize(100000000);

  t->Branch("EventNumber", &evnumberEvt,  "evnumber/i");  // 32 bitnsigned integer; not really meaningful because only 24 bits are valid
  t->Branch("TimeStamp",   &timestampEvt, "timestamp/l"); // 64 bit unsigned integer, is the common timestamp identifying the event 

  t->Branch("len_rawBuf", &len_rawBuf, "len_rawBuf/I");
  t->Branch("rawBuf",          rawBuf, "rawBuf[len_rawBuf]/F");
  t->Branch("has_rawBuf", &has_rawBuf, "has_rawBuf/I");
  t->Branch("rawBufTstamp",    rawBufTstamp, "rawBufTstamp[has_rawBuf]/l");

  t->Branch("number_of_gammas",&number_of_gammas,"number_of_gammas/I");
  t->Branch("gammaE", gammaE, "gammaE[number_of_gammas]/F");
  t->Branch("gammaT", gammaT, "gammaT[number_of_gammas]/F");
  t->Branch("gammaTstamp", gammaTstamp, "gammaTstamp[number_of_gammas]/l"); // 26-11-2011 always zero as gammaT given relative to timestampEvt
  t->Branch("gammaTrType", gammaTrType, "gammaTrType[number_of_gammas]/I");
  t->Branch("gammaX1",gammaX1,"gammaX1[number_of_gammas]/F");
  t->Branch("gammaY1",gammaY1,"gammaY1[number_of_gammas]/F");
  t->Branch("gammaZ1",gammaZ1,"gammaZ1[number_of_gammas]/F");
  t->Branch("gammaX2",gammaX2,"gammaX2[number_of_gammas]/F");
  t->Branch("gammaY2",gammaY2,"gammaY2[number_of_gammas]/F");
  t->Branch("gammaZ2",gammaZ2,"gammaZ2[number_of_gammas]/F");

  if(fWriteRootHits) {
    t->Branch("number_of_hits",&number_of_hits,"number_of_hits/I");
    t->Branch("hitE",hitE,"hitE[number_of_hits]/F");
    t->Branch("hitX",hitX,"hitX[number_of_hits]/F");
    t->Branch("hitY",hitY,"hitY[number_of_hits]/F");
    t->Branch("hitZ",hitZ,"hitZ[number_of_hits]/F");
    t->Branch("hitT",hitT,"hitT[number_of_hits]/F");
    t->Branch("hitTstamp",hitTstamp,"hitTstamp[number_of_hits]/l"); // original timestamps from the PSA
    t->Branch("hitId",hitId,"hitId[number_of_hits]/I");
    t->Branch("hitSg",hitSg,"hitSg[number_of_hits]/I");
  }

}

void TrackingFilter::CopyHits2Tree()
{
  exyzHit * pt = pEXYZ;
  for(int nn = 0; nn < number_of_hits; nn++, pt++) {
    hitE[nn]  = (Float_t)pt->E;
    hitX[nn]  = (Float_t)pt->X;
    hitY[nn]  = (Float_t)pt->Y;
    hitZ[nn]  = (Float_t)pt->Z;
    hitT[nn]  = (Float_t)cryst[pt->Ind].CCT0;
    hitTstamp[nn] = cryst[pt->Ind].tstamp;
    hitId[nn] = (Int_t)pt->Det;
    hitSg[nn] = (Int_t)pt->Seg;
  }
}

void TrackingFilter::CopyGamma2Tree (trGamma *pg, Int_t number_of_gammas)
{
  number_of_gammas = number_of_gammas > intmax ? intmax : number_of_gammas;
  for(Int_t i = 0; i < number_of_gammas; i++, pg++ ) {
    gammaE[i]  = (Float_t)pg->E;
    gammaT[i]  = (Float_t)pg->T;
    gammaTstamp[i]  = 0; //  they are all relative to timestampEvt, which is already known;
    gammaTrType[i]  = pg->type;
    gammaX1[i] = (Float_t)pg->X1;
    gammaY1[i] = (Float_t)pg->Y1;
    gammaZ1[i] = (Float_t)pg->Z1;
    gammaX2[i] = (Float_t)pg->X2;
    gammaY2[i] = (Float_t)pg->Y2;
    gammaZ2[i] = (Float_t)pg->Z2;
  }
}

void TrackingFilter::closeTTree ()
{
  if(fWriteRootTree) {
    t->Write();
    rootfile = t->GetCurrentFile();
    rootfile->Close();
    cout << "TrackingFilter::closeTTree() done" << endl;
  }
}

#endif // #ifdef TRF_ROOTTREE

#ifdef TRF_WRITE_GSORT

bool TrackingFilter::ReadGaspLUT(std::string prismaLUTfile)  
{
  FILE*  fp_lut_prisma;
  fp_lut_prisma = fopen( prismaLUTfile.data(), "r" );
  if(fp_lut_prisma==NULL) {
    cout <<  "Could not open PRISMA LUT file " << prismaLUTfile << endl;
    return false;
  }

  for( int k=0; k<600; k++ ){
    gaspLUT[k] = -1;  // empty channels will be negative
  }
  char label[256];
  int dumm1=0, dumm2=0;
  while( fgets(label, 256, fp_lut_prisma) ){
    if( label[0]=='#' )
      continue;
    if(sscanf( label, "%d %d", &dumm1, &dumm2 ) == 2) 
      gaspLUT[dumm1] = dumm2;
  }
  //for( int k=0; k<600; k++ ){
  //  cout << "gaspLUT[" << k << "] = " << gaspLUT[k] << endl;
  //}

  fclose( fp_lut_prisma );
  cout << "PRISMA LUT file " << prismaLUTfile << " read succesfully" << endl;
  return true;
}

void TrackingFilter::PrepareGaspEvent()  
{		         
  // GAMMAS (from AGATA)
  number_of_gammas = number_of_gammas > intmax ? intmax : number_of_gammas;

  // MONITOR
  int num_monitor=0;
  for( int jj=0; jj<100; jj++ ){
    if( gaspLUT[jj]<0 )
      continue;
    if(rawBuf[ gaspLUT[jj] ] > 0 )
      num_monitor ++;
  }        
  // MCP
  int num_mcp=0;
  for( int jj=100; jj<200; jj++ ){
    if( gaspLUT[jj]<0 )
      continue;
    if(rawBuf[ gaspLUT[jj] ] > 0 )
      num_mcp = 1;
      //cout << " rawBuf[ gaspLUT[" << jj << "] ] = " << rawBuf[ gaspLUT[jj] ] << " " << rawBuf[ gaspLUT[jj] ] << endl;
  }  		         

  // PPAC
  int num_ppac=0;
  bool sectFired[10];
  for( int s=0; s<10; s++ ){
    sectFired[s]=false;
  }

  for( int sect=0; sect<10; sect++ ){
    for( int param=0; param<6; param++ ){
      if( gaspLUT[200+6*sect+param]<0 )
        continue;        
      if( rawBuf[ gaspLUT[200+6*sect+param] ] > 0 )
        sectFired[sect] = true;
    }
    if( sectFired[sect] ) 
      num_ppac++;
  }
  // IC
  int num_ic=0;
  bool icFired[10];
  for( int s=0; s<10; s++ ){
    icFired[s]=false;
  }

  for( int sect=0; sect<10; sect++ ){
    for( int param=0; param<4; param++ ){
      if( gaspLUT[300+4*sect+param]<0 )
        continue;        
      if( rawBuf[ gaspLUT[300+4*sect+param] ] > 0 )
        icFired[sect] = true;
    }
    if( icFired[sect] ) 
      num_ic++;
  }
  // ICside
  int num_ics=0;
  bool icsFired[2];
  for( int s=0; s<2; s++ ){
    icsFired[s]=false;
  }

  for( int sect=0; sect<2; sect++ ){
    for( int param=0; param<4; param++ ){
      if( gaspLUT[400+4*sect+param]<0 )
        continue;        
      if( rawBuf[ gaspLUT[400+4*sect+param] ] > 0 )
        icsFired[sect] = true;
    }
    if( icsFired[sect] ) 
      num_ics++;
  }
  // LaBr3 (DANTE chan...)
  int num_labr=0;
  bool labrFired[8];
  for( int s=0; s<8; s++ ){
    labrFired[s]=false;
  }

  for( int sect=0; sect<8; sect++ ){
    for( int param=0; param<2; param++ ){
      if( gaspLUT[500+2*sect+param]<0 )
        continue;        
      if( rawBuf[ gaspLUT[500+2*sect+param] ] > 0 )
        labrFired[sect] = true;
    }
    if( labrFired[sect] ) 
      num_labr++;
  }
  gaspEvLen  = 0;
  gaspEvLen += 7;                   // number of detectors firing for each event + gammas
  gaspEvLen += number_of_gammas*5;  // 4params + crystal id (ficticious value in this case)
  gaspEvLen += num_monitor*2;
  gaspEvLen += num_mcp*4;
  gaspEvLen += num_ppac*7;
  gaspEvLen += num_ic*5;
  gaspEvLen += num_ics*5;
  gaspEvLen += num_labr*3;
  // Filling of the event array...
  int index = 0;
  gaspEvent[index++] = gaspEvLen + 0xf000; 
  gaspEvent[index++] = number_of_gammas;
  gaspEvent[index++] = num_monitor;
  gaspEvent[index++] = num_mcp;
  gaspEvent[index++] = num_ppac;
  gaspEvent[index++] = num_ic;
  gaspEvent[index++] = num_ics;
  gaspEvent[index++] = num_labr;
  // Filling detectors ...
  Float_t theta=0., phi=0., rho=0.;
  if( number_of_gammas>0 ){
    for(Int_t i = 0; i < number_of_gammas; i++ ) {
      gaspEvent[index++] = (unsigned short)i; // fictictious...gamma number
      gaspEvent[index++] = (unsigned short)pGams[i].E;
      gaspEvent[index++] = (unsigned short)pGams[i].T;
      rho   = sqrt( pGams[i].X1*pGams[i].X1 + pGams[i].Y1*pGams[i].Y1 + pGams[i].Z1*pGams[i].Z1 );
      theta = RAD2DEG * acos ( pGams[i].Z1/rho );
      phi   = RAD2DEG * atan2( pGams[i].Y1, pGams[i].X1 );
      if( phi<0 ) phi += 360.f;
      gaspEvent[index++] = (unsigned short)(10.f*theta);   // 0.1 deg/chan
      gaspEvent[index++] = (unsigned short)(10.f*phi);     // 0.1 deg/chan
    }   
  }
  if( num_monitor>0 ){
    for( int kk=0; kk<3; kk++ ){
      if( rawBuf[ gaspLUT[kk] ]>0 ){
        gaspEvent[index++] = kk;
        gaspEvent[index++] = (unsigned short)rawBuf[ gaspLUT[kk] ];
      }
    } 
  }
  if( num_mcp>0 ){
    gaspEvent[index++] = 0;
    for( int kk=0; kk<3; kk++ ){
      if( rawBuf[ gaspLUT[100+kk] ]>0 ){
        gaspEvent[index++] = (unsigned short)rawBuf[ gaspLUT[100+kk] ];
      }
    }  
  }
  if( num_ppac>0 ){
    for( int kk=0; kk<10; kk++ ){
      if( !(sectFired[kk]) )
        continue;
      gaspEvent[index++] = kk;     
      for( int jj=0; jj<6; jj++ ){
        gaspEvent[index++] = (unsigned short)rawBuf[ gaspLUT[200+6*kk+jj] ];
      }  
    }  
  }
  if( num_ic>0 ){
    for( int kk=0; kk<10; kk++ ){
      if( !(icFired[kk]) )
        continue;
      gaspEvent[index++] = kk;     
      for( int jj=0; jj<4; jj++ ){
        gaspEvent[index++] = (unsigned short)rawBuf[ gaspLUT[300+4*kk+jj] ];
      }  
    }  
  }
  if( num_ics>0 ){
    for( int kk=0; kk<2; kk++ ){
      if( !(icsFired[kk]) )
        continue;
      gaspEvent[index++] = kk;     
      for( int jj=0; jj<4; jj++ ){
        gaspEvent[index++] = (unsigned short)rawBuf[ gaspLUT[400+4*kk+jj] ];
      }  
    }  
  }
  if( num_labr>0 ){
    for( int kk=0; kk<8; kk++ ){
      if( !(labrFired[kk]) )
        continue;
      gaspEvent[index++] = kk;     
      for( int jj=0; jj<2; jj++ ){
        gaspEvent[index++] = (unsigned short)rawBuf[ gaspLUT[500+2*kk+jj] ];
      }  
    }  
  } 

} 	
#endif // TRF_WRITE_GSORT