PreprocessingFilterPSA.cpp 61.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
#include "PreprocessingFilterPSA.h"

#include <iostream>
#include <fstream>
#include <string>
#include <cmath>

#include "misc.h"

using namespace std;
using namespace ADF;

13 14 15
// The parts related to the timing tests done by Francesco using triangular
// pulses, which were enabled by #define PULSER_ON_CC_A1_A2, have been removed
// In case of need, look to revisions < 1064 
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

const int   specLenE = 16*1024;       // size of energy spectra
const int   specOffE =  0;            // 2*1024;  // offset to see the negative amplitudes in the segments
const float fEnergyGainET = 0.05f;    // normalize output to 20 keV/ch

const int matEEdim = 2000;  const float matEEgain = 1.0f;  // at  1 keV/ch
const int matPPdim = 1000;  const float matPPgain = 0.1f;  // at 10 keV/ch; on the main diagonal if tauCC==TauCC_used_in_Preprocessing
const int matETdim = 1536;
const int matSSdim = 1/*00*/;
const int matTTdim = 500;

const float calcTT_eminCC =  30.f;    // energy threshold when making the tCC-tSG analysis

const float xFactor = 500.f;          // scaling factor for the differential cross-talk correction

// The parabolic fit of the risetime has been removed because the extra degree of freedom leads to strange results (e.g. no zero-crossing) 

const int   trigChansCC  = 5;     // normally 4
const int   trigChansSG  = 5;     // normally 4
const float trigFraction = .05f;  // start fitting strightline at this fraction of max amplitude (but not lower than trigLevel1)
//const float trigLevel1   = 20.f;  // 25 keV ~5 sigma of noise ==> should be calculated as a running average of the second threshold
const float trigLevel1   = 10.f;  // lowered to 10 keV because of the gaussian smoothing performend before fitting the stright line

const float tStep    = 10.f;        // ns/sample
const float tScale   = tStep;       // 1ns/ch
const int   tGain    = 10;  // 1;   // to increase the scale and length of ...TT1 and ...TT2 spectra in a (more or less) consistent way
const int   speTTlen = int(1000*tScale);

44
#if 0 && defined(_DEBUG)      // set to 1 to enable writing the intermediate waves while debugging
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138
# define WRITEWORKINGWAVE1(p1,p2)    WriteWorkingWave(p1, p2)
#ifdef TIMING_CFD
# define WRITEWORKINGWAVE2(p1,p2)    WriteWorkingWave(p1, p2)
#else
# define WRITEWORKINGWAVE2(p1,p2)
#endif
#else
# define WRITEWORKINGWAVE1(p1,p2)
# define WRITEWORKINGWAVE2(p1,p2)
#endif

PreprocessingFilterPSA::PreprocessingFilterPSA() :
#ifdef PPF_MULTIHIST
  PreSpecEner(NULL), PreSpecEsum(NULL), PreMatrEcEsum(NULL),
  PreMatrEeEtrCC(NULL), PreMatrEeEtr(NULL), PreMatrEsEs(NULL), PreMatrEcTc(NULL), PreMatrEsTs(NULL), PreMatrEeTr(NULL),
  PreMatrTT1(NULL), PreMatrTT2(NULL), PreMatrTT3(NULL), PreMatrXT(NULL),
  PreSpecTT1(NULL), PreSpecTT2(NULL),
#endif  //PPF_MULTIHIST
  fpPrepTraces(NULL)
{
  // will be initialized from members in the mother class
  eMinCC  = 0;
  eMaxCC  = 0;
  eMinSG  = 0;

  segF1   = 0;
  segF2   = 0;

  fCC  = fSG  = NULL;
  fTC  = fTS  = fTT = fWW = NULL;
  fTFA = fCFD = NULL;

  firstEventPP = true;

#ifdef TIMING_CFD
  mwd  = NULL;
#endif
  ResetTraceLength();

  nEvDump = -1;

#ifdef PPF_FromGRU  
  CCE_raw  = new TH1F("CCE_raw","CCE_raw",32768,0,32767);
  CCE1_raw = new TH1F("CCE1_raw","CCE1_raw",32768,0,32767);
  CCE      = new TH1F("CCE","CCE",32768,0,32767);
  CCE1     = new TH1F("CCE1","CCE1",32768,0,32767);
  seg_raw  = new TH1F("seg","all 36 seg",32768*36,0,32*1024*36-1);

  // segment traces
    string SGname = "SGTrace_00"; 
  for (int i=0; i<CrystalInterface::kNbSegments; i++){
    cout << "SGname : " << SGname << endl;
    SGTraces[i]= new TH1F(SGname.c_str(),SGname.c_str(),ADF::CrystalInterface::kDefaultLength,0,ADF::CrystalInterface::kDefaultLength);
  }

  // core traces
  string CCname = "CCTrace_0"; 
  for (int i=0; i<CrystalInterface::kNbCores; i++) {
    cout << "CCname : " << CCname << endl;
    CCTraces[i]= new TH1F(CCname.c_str(),CCname.c_str(),ADF::CrystalInterface::kDefaultLength,0,ADF::CrystalInterface::kDefaultLength); 
  }

  hs = new THStack("hs","test stacked histograms");

  PPSpectraDB = new GSpectra();
  PPNetworkRoot = new GNetServerRoot (PPSpectraDB);

  PPSpectraDB->AddSpectrum (CCE_raw);
  PPSpectraDB->AddSpectrum (CCE1_raw);
  PPSpectraDB->AddSpectrum (seg_raw);
  PPSpectraDB->AddSpectrum (CCE);
  PPSpectraDB->AddSpectrum (CCE1);

  for (int i=0; i<CrystalInterface::kNbSegments; i++){
    PPSpectraDB->AddSpectrum (SGTraces[i]);
  }
  for (int i=0; i<CrystalInterface::kNbCores; i++){
    PPSpectraDB->AddSpectrum (CCTraces[i]);
  }

  // GRU-v8
  Int_t port =9091;
  PPNetworkRoot->SetPort(port);
  PPSpectraDB->SetfDefaultHostPort(port);
  PPSpectraDB->SetfDefaultHostName((char*)gSystem->HostName());
  // 
  // //////HERE TO RESTART THE SERVER//////
  PPNetworkRoot->StartGNetServer(false);
#endif // PPF_FromGRU   

}

PreprocessingFilterPSA::~PreprocessingFilterPSA()
{
139 140 141 142 143 144 145
  for(int nn = 0; nn < CrystalInterface::kNbSegments; nn++) {
    delete [] sTracesSG[nn];
  }
  for(int nn = 0; nn < CrystalInterface::kNbCores; nn++) {
    delete [] sTracesCC[nn];
  }

146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171
  cServer.Finish();

#ifdef PPF_FromGRU  
  delete CCE_raw;
  delete CCE1_raw;
  delete seg_raw;
  delete CCE;
  delete CCE1;
  for(int nn = 0; nn < CrystalInterface::kNbSegments; nn++) {
    delete [] SGTraces[nn];
  }
  for(int nn = 0; nn < CrystalInterface::kNbCores; nn++) {
    delete [] CCTraces[nn];
  }
#endif // PPF_FromGRU
}

Int_t PreprocessingFilterPSA::AlgoSpecificInitialise()
{
  // read the configuration file
  setupFile = GetConfPath() + gActualClass + ".conf";
  cout << "\n" << gActualClass + "::AlgoSpecificInitialise() --> reading from " << setupFile << endl;

  // calibrations and other parameters
  CC.SetTriggerLevel(fCoreEnerMin);
  CC.SetSegMinEnergy(fPrepSegEmin);
172
  if( !CC.ReadCalibCoeffs(setupFile, false) ) {
173 174 175 176 177 178 179 180 181 182 183 184 185 186
    cout << "Error decoding file " << setupFile << endl;
    return 130;
  }

  // cross talk
  if(fPrepXtalkFile.length() > 0) {
    string fnameF2 = GetConfPath() + fPrepXtalkFile;
    if( !CC.ReadXtalkCoeffs(fnameF2, fVerbose) ) {
      cout << "Error decoding file " << fnameF2 << endl;
      return 140;
    }
    CC.CalcDifferentialXtalk(xFactor);
  }

187 188 189 190 191 192 193
  for(int nn = 0; nn < CrystalInterface::kNbSegments; nn++) {
    sTracesSG[nn] = new Float_t [defTraceLengthRaw];
  }
  for(int nn = 0; nn < CrystalInterface::kNbCores; nn++) {
    sTracesCC[nn] = new Float_t [defTraceLengthRaw];
  }

194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527
  ResetTraceLength();

  // these are initialized to the maximum length
  if(fCC ) delete [] fCC;  fCC  = new float [defTraceLengthRaw]; memset(fCC,  0, sizeof(float)*defTraceLengthRaw);
  if(fSG ) delete [] fSG;  fSG  = new float [defTraceLengthRaw]; memset(fSG,  0, sizeof(float)*defTraceLengthRaw);
  if(fTC ) delete [] fTC;  fTC  = new float [defTraceLengthRaw]; memset(fTC,  0, sizeof(float)*defTraceLengthRaw);
  if(fTS ) delete [] fTS;  fTS  = new float [defTraceLengthRaw]; memset(fTS,  0, sizeof(float)*defTraceLengthRaw);
  if(fTT ) delete [] fTT;  fTT  = new float [defTraceLengthRaw]; memset(fTT,  0, sizeof(float)*defTraceLengthRaw);
  if(fWW ) delete [] fWW;  fWW  = new float [defTraceLengthRaw]; memset(fWW,  0, sizeof(float)*defTraceLengthRaw);
  if(fTFA) delete [] fTFA; fTFA = new float [defTraceLengthRaw]; memset(fTFA, 0, sizeof(float)*defTraceLengthRaw);
  if(fCFD) delete [] fCFD; fCFD = new float [defTraceLengthRaw]; memset(fCFD, 0, sizeof(float)*defTraceLengthRaw);

  initPresort();    // local spectra

  return 0;
}

Int_t PreprocessingFilterPSA::Process()
{
  float sumE1 = 0;    // sum of net charge segments
  float sumE2 = 0;    // sum of net charge segments after crosstalk correction
  int   sMult = 0;    // how many segments with net charge
  
#ifdef PPF_FromGRU
   CCE_raw  -> Fill((int)CoreE[0]);
   CCE1_raw -> Fill((int)CoreE[1]);
   for (int i=0; i<CrystalInterface::kNbSegments; i++){
     seg_raw->Fill(SegE[i]+i*specLenE);
   }
   for (int i=0; i<CrystalInterface::kNbSegments; i++){
    SGTraces[i]->Reset();
    for (int k=0; k<CrystalInterface::kDefaultLength; k++){
      SGTraces[i]->Fill(k,(UShort_t)fTracesSG[i][k]);//Weigthed Histo
//       std::cout << "trace_raw" <<"\t"<< k <<"\t"<< k <<"\t"<< SegmentTraces[i][k]<< "\n";
    }
  }
  for (int i=0; i<CrystalInterface::kNbCores; i++){
    CCTraces[i]->Reset();
    for (int k=0; k<CrystalInterface::kDefaultLength; k++){
      CCTraces[i]->Fill(k,(UShort_t)fTracesCC[i][k]);//Weigthed Histo
//       std::cout << "trace_raw" <<"\t"<< k <<"\t"<< k <<"\t"<< CoreTraces[i][k]<< "\n";
      }
   }
#endif  // PPF_FromGRU 

  // Normally the traces have always the same length but we must check it because the value stored in
  // the event can differ from the user-defined one (which may certainly happen for the first event)
  if(nNsamp != (int)fTraceLengthRaw)
    ResetTraceLength();

  bool evok = ProcessEvent(sMult, sumE1, sumE2);

  if(!evok) {
    crystal_status = 1;
    return 0;
  }
  crystal_status = 0;

#ifdef PPF_MULTIHIST
  // preprocessed projections
  if(PreSpecEner) {
    for(int nn = 0; nn < CC.nSG; nn++) {
      if(SegE[nn] > 0)
        PreSpecEner->Incr(0, nn, (int)(SegE[nn]*fEnergyGain+specOffE));
    }
    for(int nn = 0; nn < CC.nCC; nn++) {
      PreSpecEner->Incr(0, CC.nSG + nn,        (int)(CoreE[nn]*fEnergyGain+specOffE));  // with gain
      PreSpecEner->Incr(0, CC.nSG + nn+CC.nCC, (int)(CoreE[nn]            +specOffE));  // at 1 keV/ch
    }
  }

  if(PreSpecEsum) {
    int iee1 = (int)(sumE1*fEnergyGain+specOffE);
    PreSpecEsum->Incr(0,       0, iee1);
    PreSpecEsum->Incr(0, sMult, iee1);
    int iee2 = (int)(sumE2*fEnergyGain+specOffE);
    PreSpecEsum->Incr(1,       0, iee2);
    PreSpecEsum->Incr(1, sMult, iee2);
  }
  if(PreMatrEcEsum) {
    PreMatrEcEsum->Incr(0, (int)CoreE[0], (int)sumE1);
    PreMatrEcEsum->Incr(1, (int)CoreE[0], (int)sumE2);
  }
#endif  //PPF_MULTIHIST

  if(fWriteNumTraces > 0) {
    int rv = WriteTraces();
    if(rv)
      return rv;
  }

  // As the PSA does not know on which core we are triggering, we must pass it out as first
  if(CC.itrig == 1) {
    swap(CoreE[0], CoreE[1]);
    swap(CoreT[0], CoreT[1]);
    memcpy(fTC,          fTracesCC[0], sizeof(float)*nNsamp);
    memcpy(fTracesCC[0], fTracesCC[1], sizeof(float)*nNsamp);
    memcpy(fTracesCC[1], fTC,          sizeof(float)*nNsamp);
  }

  cServer.Exec(timestamp);

  return 0;
}

// traces and energies written in keV (if properly normalized...)
int PreprocessingFilterPSA::WriteTraces()
{
  // Opening the file must be done here because in AlgoSpecificInitialize fTraceLengthRaw
  // is not yet set to the actual value as obtained form the data itself
  if(fWriteNumTraces>0 && fpPrepTraces==NULL) {
    std::ostringstream filename;
    filename << "Prep__" << fWriteNumTraces << "-44-" << fTraceLengthRaw << "-S__Traces.samp";
    fnPrepTraces = fOdirPrefix+filename.str();
    cout << "Opening Preprocessing traces file " << fnPrepTraces << endl;
    fpPrepTraces = fopen(fnPrepTraces.c_str(), "wb");
    if(!fpPrepTraces) {
      cout << "Error: could not open " <<  fnPrepTraces << endl;
      return 1;
    }
  }

  if(!fpPrepTraces)
    return 0;

  UShort_t ss[defTraceLengthRaw];   // short int buffer to write
  UShort_t ee[defTraceLengthRaw];   // energies
  memset(ss, 0, sizeof(ss));
  memset(ee, 0, sizeof(ee));
  size_t nwritten = 0;
  // segments
  for(UInt_t nn = 0; nn < CrystalInterface::kNbSegments; nn++) {
    Float_t  *pft = fTracesSG[nn];
    for(UInt_t ii = 0; ii < fTraceLengthRaw; ii++) {
      //Short_t stmp = Short_t(pft[ii]*fEnergyGain);
      Short_t stmp = Short_t(pft[ii]);
      ss[ii] = UShort_t(stmp);
    }
    nwritten = fwrite(ss, sizeof(Short_t), fTraceLengthRaw, fpPrepTraces);
    if(nwritten != fTraceLengthRaw)
      return 1;
    //ee[nn] = (short)(SegE[nn]*fEnergyGain);
    ee[nn] = Short_t(SegE[nn]);
  }
  int eeind = CrystalInterface::kNbSegments;
  // the 2 cores
  for(UInt_t nn = 0; nn < CrystalInterface::kNbCores; nn++) {
    UShort_t *ps = ss;
    Float_t  *pf = fTracesCC[nn];
    for(UInt_t ii = 0; ii < fTraceLengthRaw; ii++) {
      //Short_t stmp = Short_t((*pf++)*fEnergyGain);
      Short_t stmp = Short_t(*pf++);
      *ps++ = UShort_t(stmp);
    }
    nwritten = fwrite(ss, sizeof(short), fTraceLengthRaw, fpPrepTraces);
    //ee[eeind++] = (short)(CoreE[nn]*fEnergyGain);
    ee[eeind++] = Short_t(CoreE[nn]);
  }
  eeind++;
  ee[eeind++] = crystal_id;
  ee[eeind++] = 1000*crystal_status;
  eeind++;
  eeind += CopyEventTags((unsigned short *)(ee+eeind));
  eeind++;
  float ttrig = CoreT[CC.itrig]*tStep;                // ns
  ee[eeind++] = UShort_t(ttrig);
  ee[eeind++] = UShort_t((ttrig-int(ttrig))*1000.f);  // fractional part in ps

  // extra stuff to have 42 channels

  // 38 - triggering core unshifted
  for(UInt_t ii = 0; ii < fTraceLengthRaw; ii++) {
    //ss[ii] = (UShort_t)(fTC[ii]*fEnergyGain);
    ss[ii] = UShort_t(fTC[ii]);
  }
  nwritten = fwrite(ss, sizeof(Short_t), fTraceLengthRaw, fpPrepTraces);
  if(nwritten != fTraceLengthRaw)
    return 1;

  // 39 - sum of net charge segments unshifted
  for(UInt_t ii = 0; ii < fTraceLengthRaw; ii++) {
    //ss[ii] = (UShort_t)(fTS[ii]*fEnergyGain);
    ss[ii] = UShort_t(fTS[ii]);
  }
  nwritten = fwrite(ss, sizeof(Short_t), fTraceLengthRaw, fpPrepTraces);
  if(nwritten != fTraceLengthRaw)
    return 1;

  // 40 - sum of net charge segments and core unshifted
  for(UInt_t ii = 0; ii < fTraceLengthRaw; ii++) {
    //ss[ii] = (UShort_t)(fTT[ii]*fEnergyGain);
    ss[ii] = UShort_t(fTT[ii]);
  }
  nwritten = fwrite(ss, sizeof(Short_t), fTraceLengthRaw, fpPrepTraces);
  if(nwritten != fTraceLengthRaw)
    return 1;

  // 41 - energies, id, status, eventnumber and timestamp
  nwritten = fwrite(ee, sizeof(Short_t), fTraceLengthRaw, fpPrepTraces);
  if(nwritten != fTraceLengthRaw)
    return 1;

  // 42 - TFA used to re-trigger the wave (empty if stright-line)
  for(UInt_t ii = 0; ii < fTraceLengthRaw; ii++) {
    ss[ii] = UShort_t(fTFA[ii]);
  }
  nwritten = fwrite(ss, sizeof(Short_t), fTraceLengthRaw, fpPrepTraces);
  if(nwritten != fTraceLengthRaw)
    return 1;

  // 43 - CFD used to re-trigger the wave (actual wave stright-line)
  for(UInt_t ii = 0; ii < fTraceLengthRaw; ii++) {
    ss[ii] = UShort_t(fCFD[ii]);
  }
  nwritten = fwrite(ss, sizeof(Short_t), fTraceLengthRaw, fpPrepTraces);
  if(nwritten != fTraceLengthRaw)
    return 1;

  fWriteNumTraces--;

  if(fWriteNumTraces < 1) {
    cout << "Closing Preprocessing traces file " << fnPrepTraces << endl;
    fclose(fpPrepTraces);
    fpPrepTraces = NULL;
  }

  return 0;
}

// this should be used for debugging purposes as it is slow and not thread safe
bool PreprocessingFilterPSA::WriteWorkingWave(float *pW, bool append)
{
  string fnWaves = "debug__100-100-100-F__prepr.fsamp";

  FILE  *fpWaves = NULL;      

  if(firstEventPP) {
    firstEventPP = false;
    append = false;
  }

  if(append)
    fpWaves = fopen(fnWaves.c_str(), "ab");
  else
    fpWaves = fopen(fnWaves.c_str(), "wb");

  if(!fpWaves) {
    cout << "Error: could not open " <<  fnWaves << endl;
    return false;
  }

  const size_t wlen = fTraceLengthRaw;
  size_t nw = fwrite(pW, sizeof(float), wlen, fpWaves);
  fclose(fpWaves);

  if(nw != wlen) {
    cout << "Error writing " <<  fnWaves << endl;
    return false;
  }
  return true;
}


bool PreprocessingFilterPSA::initPresort()
{
  srand(1);

  eMinCC    = float(fCoreEnerMin);
  eMaxCC    = float(fCoreEnerMax);
  if(eMinCC > eMaxCC) swap(eMinCC, eMaxCC);

  eMinSG   = fPrepSegEmin;

  segF1    = fPrepFoldMin;
  segF2    = fPrepFoldMax;
  if(segF1 > segF2) swap(segF1, segF2);

#ifdef PPF_MULTIHIST
  if(fUseMultiHist) {
    PreMatrEeEtrCC   = new MultiHist<unsigned short>(2, matPPdim, matPPdim);
    PreMatrEeEtrCC->setFileName(fOdirPrefix+"Prep?EeEtrCC.matr");
    PreMatrEeEtrCC->setComment("energy vs energy-from-traces for the triggering CC");
    hGroup.add(PreMatrEeEtrCC);

    //PreMatrEeEtr   = new MultiHist<unsigned short>(CC.nSG+CC.nCC, matEEdim, matEEdim);
    //PreMatrEeEtr->setFileName(fOdirPrefix+"Prep?EeEtr.matr");
    //PreMatrEeEtr->setComment("energy vs energy-from-traces for segments and cores");
    //hGroup.add(PreMatrEeEtr);

    PreMatrEsEs = new MultiHist<unsigned short>(matSSdim*CC.nSG, matSSdim*CC.nSG);
    PreMatrEsEs->setFileName(fOdirPrefix+"Prep?EsEs.matr");
    PreMatrEsEs->setComment("segment-segment energy for all pairs");
    hGroup.add(PreMatrEsEs);

    PreMatrEcTc   = new MultiHist<unsigned short>(2, matETdim, matETdim);
    PreMatrEcTc->setFileName(fOdirPrefix+"Prep?EcTc.matr");
    PreMatrEcTc->setComment("Energy(keV)-Time(ns) for the triggering core");
    hGroup.add(PreMatrEcTc);

    //PreMatrEeTr  = new MultiHist<unsigned short>(CC.nSG+CC.nCC+2, matETdim, 100);
    //PreMatrEeTr->setFileName(fOdirPrefix+"Prep?EeTr.matr");
    //PreMatrEeTr->setComment("Energy(keV)-T10-90(samp)");
    //hGroup.add(PreMatrEeTr);

    //PreMatrEsTs  = new MultiHist<unsigned short>(CC.nSG, 100, matTTdim);   // energy in steps of 20 keV, time in ns
    //PreMatrEsTs->setFileName(fOdirPrefix+"Prep?EsTs.matr");
    //PreMatrEsTs->setComment("Segments: Energy(keV)-Time(ns, relative to core)");
    //hGroup.add(PreMatrEsTs);

    //PreMatrTT1  = new MultiHist<unsigned short>(matTTdim*CC.nSG/6, matTTdim*CC.nSG/6);
    //PreMatrTT1->setFileName(fOdirPrefix+"Prep?TT1.matr");
    //PreMatrTT1->setComment("relative and absolute Core_to_Segment time, before shifting traces");
    //hGroup.add(PreMatrTT1);

    //PreMatrTT2  = new MultiHist<unsigned short>(matTTdim*CC.nSG/6, matTTdim*CC.nSG/6);
    //PreMatrTT2->setFileName(fOdirPrefix+"Prep?TT2.matr");
    //PreMatrTT2->setComment("relative and absolute Core_to_Segment time, after shifting the traces");
    //hGroup.add(PreMatrTT2);

    //PreMatrTT3  = new MultiHist<unsigned short>(5, 1000, 1000);
    //PreMatrTT3->setFileName(fOdirPrefix+"Prep?TT3.matr");
    //PreMatrTT3->setComment("timing checks");
    //hGroup.add(PreMatrTT3);

    //PreMatrXT = new MultiHist<unsigned short>(100*CC.nSG*CC.nSG, 1536);
    //PreMatrXT->setFileName(fOdirPrefix+"Prep?XT.matr");
    //PreMatrXT->setComment("CrossTalk matrix at Fold 1");
    //hGroup.add(PreMatrXT);

    PreSpecTT1 = new MultiHist<unsigned int>(2, 40, speTTlen);
    PreSpecTT1->setFileName(fOdirPrefix+"Prep?TT1.spec");
    PreSpecTT1->setComment("time spectra, before shifting the traces)");
    hGroup.add(PreSpecTT1);

528 529 530 531
    PreSpecTT2 = new MultiHist<unsigned int>(2, 40, speTTlen);
    PreSpecTT2->setFileName(fOdirPrefix+"Prep?TT2.spec");
    PreSpecTT2->setComment("time spectra after shifting the traces)");
    hGroup.add(PreSpecTT2);
532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714

    PreSpecEner = new MultiHist<unsigned int>(2, CC.nSG + 2*CC.nCC, specLenE);
    PreSpecEner->setFileName(fOdirPrefix+"Prep?Ener.spec");
    PreSpecEner->setComment("calibrated energy spectrum of segments and cores");
    hGroup.add(PreSpecEner);

    PreSpecEsum = new MultiHist<unsigned int>(2, 10, specLenE);   // up to multiplicity 9
    PreSpecEsum->setFileName(fOdirPrefix+"Prep?Esum.spec");
    PreSpecEsum->setComment("sum energy of segments before and after xTalk correction");
    hGroup.add(PreSpecEsum);

    //PreMatrEcEsum   = new MultiHist<unsigned short>(2, matEEdim, matEEdim);
    //PreMatrEcEsum->setFileName(fOdirPrefix+"Prep?EcEsum.matr");
    //PreMatrEcEsum->setComment("energy of core vs sum-of-segments before and after Xtalk");
    //hGroup.add(PreMatrEcEsum);
  }
#endif  //PPF_MULTIHIST

  return true;
}

// The new version with triggering on the sum of core + net-charge segments
// Done in FindTriggerTime by fitting a stright line to the first samples of the trace

// This method is too complicated and should be split into a series of smaller blocks 

bool  PreprocessingFilterPSA::ProcessEvent(int &segMult, float &eSumSG1, float  &eSumSG2)
{
  int   nSG = CC.nSG;     // 36
  int   nCC = CC.nCC;     //  2
  int   iTR = CC.itrig;   //  0

  segMult = nSegFold = 0;
  eSumSG1 = 0;  // sum of net charge segments before xTalk correction
  eSumSG2 = 0;  // sum of net charge segments after  xTalk correction

  ////  calibrate CC energies
  for(int nc = 0; nc < nCC; nc++) {
    float ener  = CoreE[nc];
    ener *= CC.pCC[nc].egain;
    CoreE[nc] = ener;
  }

  //// check energy gate in triggering core
  float eCore  = CoreE[iTR];
  if(eCore < eMinCC || eCore > eMaxCC)
    return false;   // out of energy range

  //// remove offset and calibrate trace for cores
  for(int icc = 0; icc < nCC; icc++) {
    float *fptr = fTracesCC[icc];
    float vertcal = CC.pCC[icc].sgain;
    float valbase = AverageValue(fptr, nNbase);
    for(int nn = 0; nn < nNsamp; nn++) {
      fptr[nn] = (fptr[nn]-valbase)*vertcal;
    }
  }

  //// test if the core energy derived from the trace is close enough to the "good" energy of core
#ifdef PPF_CHECK_EE_ES
  // This is a kind of pileup rejection, which helps rejecting pile-up events at high singles rate
  // THESE LIMITS SHOULD BE REPLACED WITH PROPER CONFIGURATION PARAMETERS <=====================
  const float purDelta = 70.0f;
  const float purSlope = 0.98f;
  float sCore = AverageValue(fTracesCC[iTR] + (nNsamp -2*nNbase), 2*nNbase);
#ifdef PPF_MULTIHIST
  if(PreMatrEeEtrCC)
    PreMatrEeEtrCC->Incr(0, (int)(matPPgain*eCore), (int)(matPPgain*sCore));   // all events
#endif
  if(sCore+purDelta < eCore*purSlope) {  // This check should be excluded if testing the BuiltinPulser  (??)
#ifdef PPF_MULTIHIST
    if(PreMatrEeEtrCC)
      PreMatrEeEtrCC->Incr(1, (int)(matPPgain*eCore), (int)(matPPgain*sCore)); // rejected events
#endif
    return false;
  }
#endif

  //// calibrate segments and adjust their traces
  float eSumSG0 = 0;  // total sum, including subthreshold and negative cross-talk components
  for(int ns = 0; ns < nSG; ns++) {
    float ener  = SegE[ns];
    ener *=  CC.pSG[ns].egain;
    SegE[ns] = ener;
    eSumSG0 += ener;
    float *fptr = fTracesSG[ns];
    float vertcal = CC.pSG[ns].sgain;
    float valbase =  AverageValue(fptr, nNbase);
    for(int nn = 0; nn < nNsamp; nn++) {
      fptr[nn] = (fptr[nn]-valbase)*vertcal;
    }
  }
  eSumSG1 = eSumSG0; 

#ifdef PPF_MULTIHIST
    // calibrated spectra before special treatment of dead and unstable segments and before cross-talk corrections
  if(PreSpecEner) {
    for(int nn = 0; nn < CC.nSG; nn++) {
      float ener = SegE[nn];
      if(ener > CC.pSG[nn].emink)
        PreSpecEner->Incr(1, nn, (int)(ener*fEnergyGain+specOffE));
    }
    // until we introduce the Xtalk for the core, these spectra are the same as in [0]
    for(int nn = 0; nn < CC.nCC; nn++) {
      PreSpecEner->Incr(1, CC.nSG + nn,        (int)(CoreE[nn]*fEnergyGain+specOffE));  // with gain
      PreSpecEner->Incr(1, CC.nSG + nn+CC.nCC, (int)(CoreE[nn]            +specOffE));  // at 1 keV/ch
    }
  }
#endif

  //// fix problematic segments (only one of the two cases)
  if(fDeadSegment >= 0) {
    // dead segment with ghost peaks
    fOriginalEnergy = SegE[fDeadSegment];
    eSumSG1 -= fOriginalEnergy;
    SegE[fDeadSegment] = 0;                 // will be replaced by CC - Sum_other_Segments

    float cSumReal = eSumSG1/fDeadFactorS;  // corrected for xTalk
    float cMissing = eCore - cSumReal;      // first approx of missing energy
    float cCorrCC  = cMissing*fDeadFactorC; // to be added to the core
    CoreE[iTR] = eCore + cCorrCC;           // corrected core energy
    float fCorrCore = CoreE[iTR]/eCore;     // needed by DeadSegmentMake

    eCore = CoreE[iTR];                     // corrected core energy

    cMissing  = eCore - cSumReal;           // better approx of missing energy
    cMissing *= fDeadFactorS;               // scaled down for xTalk

    SegE[fDeadSegment] = cMissing;
    eSumSG1           += cMissing;

    DeadSegmentMake(fDeadFactorS, fCorrCore);     // generate its trace whether it has net charge or not
  }
  else if(fUnstableSegment >= 0) {
    // gain drifting segment
    fOriginalEnergy = SegE[fUnstableSegment];
    if(fOriginalEnergy > CC.pSG[fUnstableSegment].emink) { // don't touch if too low energy
      eSumSG1 -= fOriginalEnergy;
      // replace the energy
      float cMissing = eCore*fUnstableFactorS - eSumSG1;
      SegE[fUnstableSegment] = cMissing;
      eSumSG1 += cMissing;
      // rescale the traces (replacing with Core-sumsegs would increase its noise)
      float  fgain = SegE[fUnstableSegment]/fOriginalEnergy; 
      float *fseg  = fTracesSG[fUnstableSegment];
      for(int ss = 0; ss < nNsamp; ss++) {
        fseg[ss] *= fgain;
      }
    }
  }

  //// now that things are adjusted, calculate the real fold of segments
  nSegFold = 0;
  sumSegs  = 0;
  for(int ns = 0; ns < nSG; ns++) {
    float ener  = SegE[ns];
    //if(ener > eMinSG) {             // to use the common threshold
    if(ener > CC.pSG[ns].emink) {   // but it is better to use the individual thresholds
      sumSegs  += ener;
      netChargeSegs[nSegFold++] = ns;
    }
    //else {
    //  // Although zeroing the channels below threshold seems the right thing to do, this cannot be done
    //  // because the the "negative" values are neeed to properly treat the cross-talk correction.
    //  // In fact, if the negative amplitudes produced (in the other segments) by very high energy hits are zeroed,
    //  // the "corrected" values for the now-zero channels can be easily 20-30 keV, i.e. over threshold, producing
    //  // artificial events of very large multiplicity. What we can do to remove the peak at zero in the segment spectra
    //  // is to zero the under-threshold AFTER the cross-talk correction 
    //  SegE[ns] = 0;
    //}
  }

  //// give-up if there are no net-charge segments:
  // this can happen if the event has a low energy, which is over threshold in the core but is spread
  // among several segments, all of which are below threshold. Keep also in mind that at this point the
  // segment threshold are higher than nominal as the signals are reduced by the proportional cross-talk
  if(nSegFold < 1)
    return false;

  segMult = nSegFold;
  eSumSG2 = eSumSG1 = sumSegs;      // set to this for the case of no crosstalk correction

#ifdef PPF_MULTIHIST
715
  // THIS HAS NOT BEEN DONE BY THIS PROGRAM SINCE LONG TIME --> CHECK if still working or, better, USE SortEnergy
716 717 718 719 720
  // xTalk only for segment multiplicity = 1
  if(PreMatrXT && nSegFold == 1)
    calcXTmatrix(SegE);
#endif

721 722
  //// cross-talk correction for the segment energies
  // (cross talk of traces is applied to the signal basis in PSAFilterGridSearch)
723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764
  if(CC.doecalF2) {
    segMult =
      xTalkCorrAll(SegE, sumSegs); // the number of hit segments can change due to the increase of energy produced by the xTalk correction
    //xTalkCorrNet(SegE, sumSegs); // BEFORE USING THIS VERSION, CHECK IF IT IS UP-TO-DATE 
    eSumSG2 = sumSegs;
    if(fDeadSegment>=0 && SegE[fDeadSegment]>CC.pSG[fDeadSegment].emink ) {
      // if it comes out with an overthreshold value replace it with core-sumOthers, which has a better noise than the original estimate
      SegE[fDeadSegment] = eCore - (eSumSG2 - SegE[fDeadSegment]); // fix it from the sum of overthreshold segments
    }
    else if(fUnstableSegment>=0 && SegE[fUnstableSegment]>CC.pSG[fUnstableSegment].emink ) {
      SegE[fUnstableSegment] = eCore - (eSumSG2 - SegE[fUnstableSegment]); // fix it from the sum of overthreshold segments
    }
    // final re-calculation of folds and sumsegs
    nSegFold = 0;
    sumSegs  = 0;
    for(int ns = 0; ns < nSG; ns++) {
      float ener = SegE[ns];
      if(ener > CC.pSG[ns].emink) {
        sumSegs  += ener;
        netChargeSegs[nSegFold++] = ns;
      }
    }
    eSumSG2 = sumSegs;
  }

  //// check condition on the number of net-charge segments
  if(nSegFold < segF1 || nSegFold > segF2)
    return false;

#ifdef PPF_MULTIHIST
  // enerEnergy-enerTrace matrices
  if(PreMatrEeEtr)
    calcEEmatrix(SegE, CoreE, fTracesSG, fTracesCC);
  // enerSG-enerSG matrix only for segment multiplicity = 2
  if(PreMatrEsEs && nSegFold == 2)
    calcSSmatrix(SegE);
#endif

  ////////////////////////////////////
  ////// TIMING, done with fTT ///////
  ////////////////////////////////////

765 766 767 768 769 770 771 772 773 774 775 776 777 778 779
  // Apply time shifts specified in PreprocessingFilterPSA.conf
  // to the cores and the net-charge segments needed for timing
  // Being the reference for all time shifts the triggering core is not moved
  for(int ic = 0; ic < nCC; ic++) {
    WRITEWORKINGWAVE1(fTracesCC[ic], true);
    float xmove = (ic==iTR) ? 0 : CC.pCC[ic].tmove/tStep;
    ShiftMoveTrace(fTracesCC[ic], nNsamp, xmove, sTracesCC[ic]);
    WRITEWORKINGWAVE1(sTracesCC[ic], true);
  }
  for(int is = 0; is < nSegFold; is++) {
    int ns = netChargeSegs[is];
    WRITEWORKINGWAVE1(fTracesSG[ns], true);
    ShiftMoveTrace(fTracesSG[ns], nNsamp, CC.pSG[ns].tmove/tStep, sTracesSG[ns]);
    WRITEWORKINGWAVE1(sTracesSG[ns], true);
  }
780

781 782 783 784
  //// copy the triggering core into the "netcharge" wave
  memcpy(fTC, sTracesCC[iTR], nNsamp*sizeof(float));

  //// calculate the sum of shifted net-charge waves for the trigger
785 786 787 788 789
  // WE SHOULD PROBABLY EXCLUDE FROM THE SUM THE VERY LOW ENERGY SEGMENTS,
  // BY SETTING A THRESHOLD ON A FRACTION OF THE CORE ENERGY
  memset(fTS, 0, sizeof(float)*nNsamp);
  for(int is = 0; is < nSegFold; is++) {
    int ns = netChargeSegs[is];
790
    ShiftAccuTrace(sTracesSG[ns], nNsamp, 0, fTS);
791 792 793 794 795 796 797 798 799 800
  }

  //// sum of core and net-charge segments
  for(int ii = 0; ii < nNsamp; ii++) {
    fTT[ii] = fTC[ii] + fTS[ii];
  }
  WRITEWORKINGWAVE1(fTT, true);

  //// select which of the calculated waves to use for triggering
  float *fTU = fTrigUseOnlyCC ? fTC : fTT;    // core alone  or  core + net-charge segments
801
  float tRefCC = FindTriggerTime(fTU, trigFraction, trigChansCC);
802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827
  WRITEWORKINGWAVE1(fTU, true);
  WRITEWORKINGWAVE2(mwd->MWDptrTFA(), true);  // timing filter wave from the trigger
  WRITEWORKINGWAVE2(mwd->MWDptrCFD(), true);  // cfd wave in the trigger

  //// save a copy of the TFA and CFD shaping for WriteTraces
#ifdef TIMING_CFD
  if(fWriteNumTraces > 0) {
    memcpy(fTFA, mwd->MWDptrTFA(), nNsamp*sizeof(float));
    memcpy(fCFD, mwd->MWDptrCFD(), nNsamp*sizeof(float));
    fCFD[0] = tRefCC*100;
  }
#else
  if(fWriteNumTraces > 0) {
    memset(fTFA, 0, nNsamp*sizeof(float));
    fTFA[fitChan1] = fTFA[fitChan2] = 1000;
    memcpy(fCFD, fTU, nNsamp*sizeof(float));
    fCFD[0] = tRefCC*100;
  }
#endif

  //// what if timing too bad ??
  bBadCFD = (tRefCC < 1.f);
  if( bBadCFD )
    return false;                 // discard the event

  //// put both trigger times to the same value
828 829 830 831 832
  //   Due to the fact that the GTS tree alignment is done only up to the Preprocessing electronics an
  //   (and not down to the Digitizer) the correction CC.pCC[iTR].tmove/tStep does not adjust the
  //   Tgamma-gamma of all detector pairs to be in the same position. Therefore it does not make sense to apply it.
  //   In fact, in PostPSAFilter.conf we can give the (empyrically determined) proper value with the keyword ShiftCC
  CoreT[0] = CoreT[1] = tRefCC /*+ CC.pCC[iTR].tmove/tStep*/; // not to be used
833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853

#if 0
  //// select only events in a given time range
  if(tRefCC < 42.0f || tRefCC > 42.6f) {
    return false;
  }
#endif

#ifdef PPF_MULTIHIST
  if(PreMatrEcTc) {
    // Energy-TriggerTime matrix
    int tzero = int(tRefCC*tScale);   // in ns
    int enecc = int(CoreE[iTR]);      // in keV
    PreMatrEcTc->Incr(0, enecc, tzero);
  }
  if(PreMatrEeTr) {
    // Energy-Risetime matrix for all channels
    calcRiseTime();
  }
  if(PreSpecTT1) {
    // absolute and relative time spectra (done before shifting the waves)
854
    calcTT(tRefCC, PreSpecTT1, PreMatrTT1, PreMatrEsTs);
855 856 857 858 859 860 861 862 863
  }
#endif  // PPF_MULTIHIST

  //// time alignment of core and traces

  // Shift CC's and SG's by the amount that moves the trigger point to defTriggerSample+defTriggerSamplePlus (normally 30)
  // The final shift to (defTriggerSample/4)*4 is done in PreprocessingFilter::SetOutput() 

  // fHandShift is an ad-hoc common shift to have the traces well positioned in the PSA.
864
  float fHandShift = CC.pCC[iTR].tmove/tStep;     // is interpreted as global shift
865
  float fshiftCC   = (defTriggerSample+defTriggerSamplePlus) - tRefCC + fHandShift;
866

867
  for(int nc = 0; nc < nCC; nc++) {
868 869 870
    float xmove = (nc==iTR) ? fshiftCC : fshiftCC + CC.pCC[nc].tmove/tStep;
    ShiftMoveTrace(fTracesCC[nc], nNsamp, xmove, sTracesCC[nc]);
    memcpy(fTracesCC[nc], sTracesCC[nc], sizeof(float)*nNsamp);
871 872 873
  }
  for(int ns = 0; ns < nSG; ns++) {
    // Segments shifted respect to the core as specified (in ns) in the last column of the calibration file
874 875
    ShiftMoveTrace(fTracesSG[ns], nNsamp, fshiftCC + CC.pSG[ns].tmove/tStep, sTracesSG[ns]);
    memcpy(fTracesSG[ns], sTracesSG[ns], sizeof(float)*nNsamp);
876 877 878
  }

  // ?????????????????
879
  // After the shift we should, for consistency, modify also timestamp and the timing values reported in CoreT
880 881
  //timestamp -= int(fshiftCC + 10000.f) -10000;
  //CoreT[0] = CoreT[1] = defTriggerSample+defTriggerSamplePlus+fHandShift;  // or tRefCCs if calculated (see next block)
882
  // However, doing nothing like until now is also fine because the absolute time tstamp+CoreT[0] is (and must be) the same in both cases !!
883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146
  // (doing nothing has also the advantage to avoid rounding problems)
  // One has just to remember the inconsistency between the position of the traces and their absolute time
  // ?????????????????

#ifdef PPF_MULTIHIST
  // absolute and relative time spectra for the shifted waves
  if(PreSpecTT2) {
#if 1
    // To be sure that the traces have been shifted correctly, rebuild fTC, fTS and fTT (just fTT) with no shifts and retrigger
    // To avoid loosing the "unshifted" view in the last traces of the Prep__1000-42-100-S__Traces.samp, we put it directly in fWW
    // Remember to build this wave consistently with what done before (i.e. which fTU)
    memcpy(fWW, fTracesCC[iTR], sizeof(float)*nNsamp);
    if(fTU == fTT) {
      for(int is = 0; is < nSegFold; is++) {
        int   ns    = netChargeSegs[is];
        float *fptr = fTracesSG[ns];
        for(int ii = 0; ii < nNsamp; ii++)
          fWW[ii] += fptr[ii];
      }
    }
    float tRefCCs = FindTriggerTime(fWW, trigFraction, trigChansCC);
#else  // #if 1
    float tRefCCs = tRefCC + fshiftCC;   // trusting the old value
#endif // #if 1

    //ULong64_t newtstamp = timestamp - (int(fshiftCC + 10000.f)-10000); 
    //ULong64_t ts1 = timestamp + int(tRefCC);
    //ULong64_t ts2 = newtstamp + int(tRefCCs);

    if(PreMatrEcTc) {
      // Energy-TriggerTime matrix
      int tzero = int(tRefCCs*tScale);    // in ns
      int enecc = int(CoreE[iTR]);        // in keV
      PreMatrEcTc->Incr(1, enecc, tzero);
    }

    calcTT(tRefCCs, PreSpecTT2, PreMatrTT2, NULL);
  }
#endif  //PPF_MULTIHIST

#if 1
  //// a small gaussian smoothing of the traces
  // 25-50-25 has a bandwidth of 18 MHz, attenuation of ~ 10^8 at 50 MHz, too strong
  // 10-80-10 has a bandwidth of 33 MHz, attenuation of ~ 0.6  at 50 MHz
  const float wgg = 0.8f;           // weight of central sample
  const float wg0 = (1.f - wgg)/2;  // weight of side samples
  const float wg1 = wgg/wg0;        // weight of central sample with side samples scaled to 1  
  for(int ns = 0; ns < nSG; ns++) {
    float v2 = fTracesSG[ns][0];
    float v1 = v2;                  // extend first value
    float v0;
    float *psg = fTracesSG[ns];
    for(int ii = 1; ii < nNsamp; ii++) {
      v0 = v1; v1 = v2; v2 = psg[ii];
      psg[ii-1] = wg0 * (v0 + wg1*v1 + v2); 
    }
    v0 = v1; v1 = v2; /* v2 = 0;*/  // extend last value
    psg[nNsamp-1] = wg0 * (v0 + wg1*v1 + v2);
  }
  for(int nc = 0; nc < nCC; nc++) {
    float v2 = fTracesCC[nc][0];
    float v1 = v2;                  // extend first value
    float v0;
    float *pcc = fTracesCC[nc];
    for(int ii = 1; ii < nNsamp; ii++) {
      v0 = v1; v1 = v2; v2 = pcc[ii];
      pcc[ii-1] = wg0 * (v0 + wg1*v1 + v2); 
    }
    v0 = v1; v1 = v2; /* v2 = 0;*/  // extend last value
    pcc[nNsamp-1] = wg0 * (v0 + wg1*v1 + v2);
  }
#endif

  return true;
}

// call it only if really needed (i.e. at initialization and for the first event) 
void PreprocessingFilterPSA::ResetTraceLength()
{
    nNsamp = fTraceLengthRaw;
    nNbase = int(0.15*nNsamp + 0.5);
    if(nNbase < 1) nNbase = 1;
#ifdef TIMING_CFD
    // reset also mwd (if used)
    if(mwd) delete mwd;
    mwd = new MWD(nNsamp);
    // all this should be settable by the caller
    mwd->nDcwin     = nNsamp/5; // the first 20% --> check if this is valid
    mwd->nSmooth    = fCFD_CC_int;    // 0;
    mwd->nWidthT    = fCFD_CC_diff;   // 10; // 6;
    mwd->nDelayCFD  = fCFD_CC_delay;  // 5;
    mwd->fFractCDF  = fCFD_CC_fract;  // 0.25
    mwd->fThreshTFA = fCFD_CC_thresh; // 20;
    mwd->fThreshCFA = 0; 
    mwd->nMinWidthP = 2;
    mwd->nMinWidthN = 6;
#endif
}

#ifdef TIMING_CFD
// The parameters of the call are those of the stright-line fit
float PreprocessingFilterPSA::FindTriggerTime (float *pT,  float trigFract, int trigChans)
{
  unsigned int tsamp = 0;
  float        fract = 0;
  mwd->MWtimingCFD(pT, tsamp, fract);
  if(tsamp < 1) {
    bBadCFD = true;
    mwd->MWDptrCFD()[0] = 0;
    return 0;
  }
  float ttrig = tsamp+fract;
  mwd->MWDptrCFD()[0] = ttrig;
  return ttrig;
}

#else // #ifdef TIMING_CFD

// This is the "default"
// As the data is supposed to arrive calibrated, the threshold is always given in keV
// Trigger point returned in units of samples

// linear fit of the risetime
float PreprocessingFilterPSA::FindTriggerTime (float *pT,  float trigFract, int trigChans)
{
  // the threshold used to calculate the stright line is
  // max(trigLevel1, valAmpli*trigFract)

  // move to working buffer (if needed)
  if(pT != fWW)
    memcpy(fWW, pT, sizeof(float)*nNsamp);

#define NGSMOO 1
#if NGSMOO == 1
  // gaussian smoothing once
  float oldval = fWW[0];
  for(int nn = 0; nn < nNsamp-1; nn++) {
    float newval = 0.25f*(oldval + 2.f*fWW[nn] + fWW[nn+1]); 
    oldval = fWW[nn];
    fWW[nn] = newval;
  }
  fWW[nNsamp-1] = 0.25f*(oldval + 3.f*fWW[nNsamp-1]);
#elif NGSMOO == 2
  // gaussian smoothing twice
  float oldval2 = fWW[0];
  float oldval1 = fWW[0];
  float newval;
  int nn;
  for(nn = 0; nn < nNsamp-2; nn++) {
    newval = 0.0625f*(oldval2 + 4.f*oldval1 + 6.f*fWW[nn] + 4.f*fWW[nn+1] + fWW[nn+2]); 
    oldval2 = oldval1;
    oldval1 = fWW[nn];
    fWW[nn] = newval;
  }
  nn = nNsamp-2;
  newval = 0.0625f*(oldval2 + 4.f*oldval1 + 6.f*fWW[nn] + 5.f*fWW[nn+1]);
  oldval2 = oldval1; oldval1 = fWW[nn]; fWW[nn] = newval;
  nn = nNsamp-1;
  newval = 0.0625f*(oldval2 + 4.f*oldval1 + 11.f*fWW[nn]);
  oldval2 = oldval1; oldval1 = fWW[nn]; fWW[nn] = newval;
#endif

  // determine amplitude (valStart should be already ~0 from the calibration )
  float valStart = AverageValue(fWW, nNbase);                   // signal average at the beginning
  float valEnd   = AverageValue(fWW+nNsamp-2*nNbase, 2*nNbase); // signal average at the end
  float valAmpli = abs(valEnd - valStart);                      // signal amplitude
  if(valAmpli < trigLevel1)                                     // signal below threshold
    return 0;
  float valFract = valAmpli*trigFract;                          // trigger level as a fraction of amplitude
  // Setting level=trigLevel1 for all signals gives worst results, probably because for large
  // signals we get into the region were the risetime is limited by the preamp bandwidth.
  float level = max(trigLevel1, valFract); // trigLevel1;

  // starting from the top, go down until the first channel below level
  fitChan1 = 0;
  level += valStart;
  for(int nn = nNsamp-2*nNbase; nn > 0; nn--) {
    if(fWW[nn] < level) {
      fitChan1 = nn + 1;             // this is the first channel above the trigger level
      break;
    }
  }

  // special treatment for low amplitude signals  ???
  if(valAmpli < 50.f) {
    int iValFract = int(valFract);
    trigChans += 4 - iValFract;
  }
  fitChan2 = fitChan1+trigChans-1;

  // normal way
  float A = 0, B = 0;
  if( !StrightLine(fWW+fitChan1, trigChans, A, B) )
    return 0;   // ??

  if(B == 0)
    return 0;

  float tZero = fitChan1 + (valStart - A) / B;
  if(tZero < 1 || tZero >= nNsamp-1)
    return 0;   // ??

#if 0
  // THERE IS A BUG IN THIS SECTION, that causes a core dump in the Release version

  // generate an up-sampled version of the TrigChan channels
  const int UPNCHAN = 10;   // could create problems if this ends-up off the wave
  const int UPSCALE = 16;
  float upsampled[UPNCHAN*UPSCALE]; memset(upsampled, 0, sizeof(upsampled));
  float *chh = fWW+fitChan1-1; // start with the previous
  int iup = 0;
  for(int nn = 0; nn <= trigChans; nn++) {
    float dv = (chh[nn+1]-chh[nn])/UPSCALE;
    float vv = chh[nn];
    for(int jj = 0; jj < UPSCALE; jj++) {
      upsampled[iup++] = vv;
      vv += dv;
    }
  }
  // find again the first channel above the threshold
  int fitChan1up = 0;
  for(int nn = iup-1; nn > 0; nn--) {
    if(upsampled[nn] < level) {
      fitChan1up = nn + 1;             // this is the first channel above the trigger level
      break;
    }
  }

  float Aup = 0, Bup = 0;
  if( !StrightLine(upsampled+fitChan1up, trigChans*UPSCALE, Aup, Bup) )
    return 0;   // ??

  if(Bup == 0)
    return 0;

  float tZeroup = fitChan1up + (valStart - Aup) / Bup;
  tZeroup = fitChan1-1 + tZeroup/UPSCALE;
  if(tZeroup < 1 || tZeroup >= nNsamp-1)
    return 0;   // ??

  tZero = tZeroup;

#endif

  // avoid the rest (which is also incompatible with the upsampled part)
  return tZero;

  // now recalculate the baseline up to a few channels less than tZero

  int newBase = (int)tZero-5;
  //if(newBase < nNbase)
  //  return tZero;       // nothing to improve (to be reconsidered)
  if(newBase < 5)
    return 0;       // too early

  valStart = AverageValue(fWW, newBase);  // signal average at the beginning
  //valEnd  -= valStart;

  float tZero1 = fitChan1 + (valStart - A) / B;
  if(tZero1 < 1 || tZero1 >= nNsamp-1)
    return 0;

  return tZero1;

1147
  //// now try with sqrt(moving sum)   ????????????????????????
1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399
  //// given how the last baseline has been calculated, the running sum of values respect to it should be zero at chan1-2
  //float accu = 0;
  ////for(int nn = chan1-1; nn < chan1 + trigChans; nn++) {
  //for(int nn = 0; nn < nNsamp; nn++) {
  //  accu += (fWW[nn]-valBase);
  //  if(accu >= 0)
  //    fWW[nn] =  sqrt(accu);
  //  else
  //    fWW[nn] = -sqrt(-accu);
  //}
  //if( !StrightLine(fWW+chan1, trigChans, A, B) )
  //  return tZero1;   // be happy with previous

  //float tZero2 = chan1 - A / B;
  //if(tZero2 < 0 || tZero2 >= nNsamp)
  //  return tZero1;   //  keep previous

  //return tZero2;
}

bool PreprocessingFilterPSA::StrightLine(float *data, int nchan, float &A, float&B)
{
  if(nchan < 2) {
    if(nchan < 1)
      return false;
    A = data[0];
    B = 0;
    return true;
  }

  double ref = 0.5;   // passing throuth the middle of the channel
  double sx0 = 0, sx1 = 0, sx2 = 0;
  double yx0 = 0, yx1 = 0;
  for(int ii = 0; ii < nchan; ii++) {
    double x = ii + 0.5;
    double y = data[ii];
    sx0 += 1;
    sx1 += x;
    sx2 += x*x;
    yx0 += y;
    yx1 += y*x;
  }

  // the [x] components and the determinant censtant and could/should be precalculated
  //  N  =>  1   2    3    4    5    6    7     8     9    10
  //  x  =>  0   1    2    3    4    5    6     7     8     9
  // sx0 =>  1   2    3    4    5    6    7     8     9    10  ==  N
  // sx1 =>  0   1    3    6   10   15   21    28    36    45  ==  N*(N-1)/2
  // sx2 =>  0   1    5   14   30   55   91   140   204   285  == (2*N^3 - 3*N^2 + N)/6
  // det =>  0   1    6   20   50  105  196   336   540   825  == (N^4 - N^2)/12

  double deter = sx0*sx2-sx1*sx1;
  if(deter >= 0.) {
    A = float(( yx0*sx2 - yx1*sx1)/deter);
    B = float((-yx0*sx1 + yx1*sx0)/deter);
    return true;
  }
  else {
    A = B = 0;
    return false;
  }
}

bool PreprocessingFilterPSA::Parabola(float *data, int nchan, float &A, float &B, float &C)
{
  if(nchan < 3) {
    C = 0;
    return StrightLine(data, nchan, A, B);
  }

  // the [x] components and the determinant are constant and could/should be precalculated
  //  N  =>  1   2    3    4    5    6     7     8      9     10
  //  x  =>  0   1    2    3    4    5     6     7      8      9
  //  x2 =>  0   1    4    9   16   25    36    49     64     81
  //  x3 =>  0   1    8   27   64  125   216   343    512    729
  //  x4 =>  0   1   16   81  256  625  1296  2401   4096   6561
  // sx0 =>  1   2    3    4    5    6     7     8      9     10  ==  N
  // sx1 =>  0   1    3    6   10   15    21    28     36     45  ==  N*(N-1)/2
  // sx2 =>  0   1    5   14   30   55    91   140    204    285  == (2*N^3 - 3*N^2 + N)/6
  // sx3 =>  0   1    9   36  100  225   441   784   1296   2025  ==  sx1^2
  // sx4 =>  0   1   17   98  354  979  2275  4676   8772  15333  == (6*N^5 - 15*N^4 + 10*N^3 - N)/30
  // det =>  0   0    4   80  700 3290 16464 56448 166320 435600  ==  ...

  double sx0 = 0, sx1 = 0, sx2 = 0, sx3 = 0, sx4 = 0;
  double yx0 = 0, yx1 = 0, yx2 = 0;
  sx0 = nchan;
  for(int ii = 0; ii < nchan; ii++) {
    double x = ii+0.5, t = ii+0.5;
    double y = data[ii];
    sx1 += t;
    yx0 += y;
    yx1 += y*t;
    t   *= x;
    sx2 += t;
    yx2 += y*t;
    t   *= x;
    sx3 += t;
    sx4 += x*t;
  }

  double deter = sx0*sx2*sx4 + sx1*sx3*sx2 + sx2*sx1*sx3 - sx0*sx3*sx3 - sx1*sx1*sx4 - sx2*sx2*sx2;
  if(deter >= 0.) {
    A = float( ( yx0*(sx2*sx4-sx3*sx3) + yx1*(sx2*sx3-sx1*sx4) + yx2*(sx1*sx3-sx2*sx2) ) / deter );
    B = float( ( yx0*(sx3*sx2-sx1*sx4) + yx1*(sx0*sx4-sx2*sx2) + yx2*(sx2*sx1-sx0*sx3) ) / deter );
    C = float( ( yx0*(sx1*sx3-sx2*sx2) + yx1*(sx1*sx2-sx0*sx3) + yx2*(sx0*sx2-sx1*sx1) ) / deter );
    return true;
  }
  else {
    C = 0;
    return StrightLine(data, nchan, A, B);
  }
}

#endif  // #else // #ifdef TIMING_CFD

void PreprocessingFilterPSA::ShiftWave(float *data, int nshift)
{
  if(nshift == 0) {
    return;
  }
  else if(nshift > 0) {
    int ii = int(nNsamp)-1;
    for(; ii >= nshift; ii--)
      data[ii] = data[ii-nshift];
    for(ii = 1; ii < nshift; ii++)
      data[ii] = data[0];
  }
  else {
    int ii = 0;
    int jj = int(nNsamp)+nshift;
    for(; ii < jj; ii++)
      data[ii] = data[ii-nshift];
    ii = int(nNsamp)+nshift;
    jj = int(nNsamp)-1;
    for(; ii < jj; ii++)
      data[ii] = data[nNsamp-1];
  }
}

void PreprocessingFilterPSA::ShiftMoveTrace(float *data, int nsamp, float delta, float *accu)
{
  if(delta == 0) {
    for(int nn = 0; nn < nsamp; nn++) {
      accu[nn] = data[nn];
    }
    return;
  }
  if(delta > 0) {
    int   idelta = int(delta);
    float fdelta = delta-idelta;
    float gdelta = 1.f - fdelta;
    float first  = data[0]*fdelta + data[1]*gdelta;
    for(int nn = 0; nn <= idelta; nn++) {
      accu[nn] = first;
    }
    for(int jj = 0, nn = idelta+1; nn < nsamp; jj++, nn++) {
      accu[nn] = data[jj]*fdelta + data[jj+1]*gdelta;
    }
  }
  else {
    int   idelta = -int(delta);
    float gdelta = abs(delta)-idelta;
    float fdelta = 1.f - gdelta;
    int nn = 0;
    for(int jj = idelta; jj < nsamp-1; jj++, nn++) {
      accu[nn] = data[jj]*fdelta + data[jj+1]*gdelta;
    }
    float last  = data[nsamp-2]*fdelta + data[nsamp-1]*gdelta;
    for(; nn < nsamp; nn++) {
      accu[nn] = last;
    }
  }
}

void PreprocessingFilterPSA::ShiftAccuTrace(float *data, int nsamp, float delta, float *accu)
{
  if(delta == 0) {
    for(int nn = 0; nn < nsamp; nn++) {
      accu[nn] += data[nn];
    }
    return;
  }
  if(delta > 0) {
    int   idelta = int(delta);
    float fdelta = delta-idelta;
    float gdelta = 1.f - fdelta;
    float first  = data[0]*fdelta + data[1]*gdelta;
    for(int nn = 0; nn <= idelta; nn++) {
      accu[nn] += first;
    }
    for(int jj = 0, nn = idelta+1; nn < nsamp; jj++, nn++) {
      accu[nn] += data[jj]*fdelta + data[jj+1]*gdelta;
    }
  }
  else {
    int   idelta = -int(delta);
    float gdelta = abs(delta)-idelta;
    float fdelta = 1.f - gdelta;
    int nn = 0;
    for(int jj = idelta; jj < nsamp-1; jj++, nn++) {
      accu[nn] += data[jj]*fdelta + data[jj+1]*gdelta;
    }
    float last  = data[nsamp-2]*fdelta + data[nsamp-1]*gdelta;
    for(; nn < nsamp; nn++) {
      accu[nn] += last;
    }
  }
}

void PreprocessingFilterPSA::DeadSegmentMake(float fractS, float fractC)
{
  float *fDead = fTracesSG[fDeadSegment];
  memset(fDead, 0, sizeof(float)*nNsamp);

  for(int ns = 0; ns < CrystalInterface::kNbSegments; ns++) {
    if(ns == fDeadSegment)
      continue;
    float *fRest = fTracesSG[ns];
    for(int ss = 0; ss < nNsamp; ss++)
      fDead[ss] += fRest[ss];
  }

  float *fCore = fTracesCC[CC.itrig];
  for(int ss = 0; ss < nNsamp; ss++) {
    fCore[ss] *= fractC;                        // adjustment of CC
    fDead[ss]  = fractS*fCore[ss] - fDead[ss];  // value for deadSeg
  }

}

#ifdef PPF_MULTIHIST

// specTT[0][0...NSEG-1]    triggering point of segments
// specTT[0][NSEG, NSEG+1]  triggering point of the two cores cRefCC and cRefCC1
// specTT[0][NSEG+2]        triggering point of the external trigger tRefCC
// specTT[1][0...NSEG-1]    timing (in ns) segment re external trigger (tRefSG-tRefCC)*tScale+speTTlen/2 (for eCore > calcTT_eminCC)
// specTT[1][NSEG, NSEG+1]  timing of the two cores re external trigger
// specTT[0][NSEG+2]        timing between the two cores ()
// Remember that the recalculated local triggers are not (yet?) moved according to the calibrationo table, which so far affects only tRefCC
void PreprocessingFilterPSA::calcTT(float tRefCC, MultiHist<unsigned int>*specTT, MultiHist<unsigned short>*matrTT, MultiHist<unsigned short>*matrETsg)
{
  //int aa1 = int((tRefCC-int(tRefCC))*tScale*tGain);
  //int aa2 = int(tRefCC*tScale*tGain);
  //int aa3 = aa2%100;

  if(!specTT)
    return;

  specTT->Incr(0, CC.nSG+2, int(tRefCC*tScale*tGain));  // "external" trigger point (usually from the CC+Netcharges)
  specTT->Incr(0, CC.nSG+2, int((tRefCC-int(tRefCC))*tScale*tGain));  // to verify flatness of subsample distribution

  // trigger on the triggering core
1400 1401
  float cRefCC = FindTriggerTime(sTracesCC[CC.itrig], trigFraction, trigChansCC);
  WRITEWORKINGWAVE1(sTracesCC[CC.itrig], true);
1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413
  WRITEWORKINGWAVE2(mwd->MWDptrTFA(), true);       // timing filter wave from the trigger
  WRITEWORKINGWAVE2(mwd->MWDptrCFD(), true);       // cfd wave in the trigger
  if(cRefCC < 1.f)
    return;

  specTT->Incr(0, CC.nSG, int(cRefCC*tScale*tGain));        // triggering-core trigger point 
  specTT->Incr(0, CC.nSG, int((cRefCC-int(cRefCC))*tScale*tGain));  // to verify flatness of subsample distribution
  float tCCE  = (cRefCC-tRefCC)*tScale*tGain;               // time relative to the external trigger
  specTT->Incr(1, CC.nSG, int(tCCE + speTTlen/2));  // relative time placed at the center of spectrum

  // trigger on the second core
  int   iCC1 = (CC.itrig+1)%ADF::CrystalInterface::kNbCores;
1414 1415
  float cRefCC1 = FindTriggerTime(sTracesCC[iCC1], trigFraction, trigChansCC);
  WRITEWORKINGWAVE1(sTracesCC[iCC1], true);
1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427
  WRITEWORKINGWAVE2(mwd->MWDptrTFA(), true);       // timing filter wave from the trigger
  WRITEWORKINGWAVE2(mwd->MWDptrCFD(), true);       // cfd wave in the trigger

  if(cRefCC1 > 1.f) {
    specTT->Incr(0, CC.nSG+1, int(cRefCC1*tScale*tGain));       // trigger point for the second core
    specTT->Incr(0, CC.nSG+1, int((cRefCC1-int(cRefCC1))*tScale*tGain));  // to verify flatness of subsample distribution
    float tCC1  = (cRefCC1-tRefCC)*tScale*tGain;                // time relative to the external trigger
    specTT->Incr(1, CC.nSG+1, int(tCC1 + speTTlen/2));    // relative time placed at the center of spectrum
    float tCC2  = (cRefCC1-cRefCC)*tScale*tGain;                // timing between the two cores
    specTT->Incr(1, CC.nSG+2, int(tCC2 + speTTlen/2));    // relative time placed at the center of spectrum
  }

1428
//#define TT_FOR_ALL_SEGMENT_MULTIPLICITIES
1429 1430 1431 1432 1433 1434 1435 1436

#ifndef TT_FOR_ALL_SEGMENT_MULTIPLICITIES
  if(nSegFold != 1)
    return;
#endif

  for(int nnseg = 0; nnseg < nSegFold; nnseg++) {
    int theSeg = netChargeSegs[nnseg];
1437 1438
    float tRefSG = FindTriggerTime(sTracesSG[theSeg], trigFraction, trigChansSG);
    WRITEWORKINGWAVE1(sTracesSG[theSeg], true);
1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761 1762 1763 1764 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 1776 1777 1778 1779 1780
    WRITEWORKINGWAVE2(mwd->MWDptrTFA(), true);       // timing filter wave from the trigger
    WRITEWORKINGWAVE2(mwd->MWDptrCFD(), true);       // cfd wave in the trigger
    if(tRefSG < 1.f)
      continue;

    specTT->Incr(0, theSeg, int(tRefSG*tScale*tGain));        // trigger point for the segment
    specTT->Incr(0, theSeg, int((tRefSG-int(tRefSG))*tScale*tGain));  // to verify flatness of subsample distribution
    float tCCSG  = (tRefSG-tRefCC)*tScale*tGain;                // time relative to the external trigger
    if(CoreE[CC.itrig] > calcTT_eminCC) {
      specTT->Incr(1, theSeg, int(tCCSG + speTTlen/2)); // relative time placed at the center of spectrum
    }
    if(matrETsg) {
      matrETsg->Incr(theSeg, (int)(fEnergyGainET*SegE[theSeg]), int(tCCSG+matTTdim/2));  // seg-gamma-time
    }
    if(matrTT) {
      int t1 = int(cRefCC*tScale - 100);    // offset depends on electronics 
      int t2 = int(tRefSG*tScale - 100);    // and should be adjusted manually
      if(t1>0 && t1<matTTdim && t2>0 && t2<matTTdim) {
        int n1 = (theSeg/6)*matTTdim;  // sectors
        int n2 = (theSeg%6)*matTTdim;  // slices
        matrTT->Incr(t1+n1, t2+n2);
      }
    }
  }
}

void PreprocessingFilterPSA::calcTTPulser(float tRefCC, MultiHist<unsigned int>*specTT, MultiHist<unsigned short>*matrTT, MultiHist<unsigned short>*matrETsg)
{
  if(!specTT)
    return;

  specTT->Incr(0, CC.nSG+2, int(tRefCC*tScale*tGain));  // "external" trigger point (usually from the CC+Netcharges)

  // trigger on the core
  float cRefCC = FindTriggerTime(fTracesCC[CC.itrig], trigFraction, trigChansCC);
  if(cRefCC < 1.f)
    return;

#if 0
  // discard if background contains extra noise
  int   ncback = int(cRefCC) - 5;
  const float mVariance = 2.f;
  float cValue    = AverageValue   (fTracesCC[CC.itrig], ncback);
  float cVariance = AverageVariance(fTracesCC[CC.itrig], ncback, cValue);
  if(cVariance > mVariance)
    return;
  // practically no improvement
#endif

  specTT->Incr(0, CC.nSG, int(cRefCC*tScale*tGain));  // triggering-core trigger point 
  float tCCE  = (cRefCC-tRefCC)*tScale*tGain;         // time relative to the external trigger
  specTT->Incr(1, CC.nSG+2, int(tCCE + speTTlen/2));  // relative time placed at the center of spectrum

  // the second core
  int   iCC1 = (CC.itrig+1)%ADF::CrystalInterface::kNbCores;
  float cRefCC1 = FindTriggerTime(fTracesCC[iCC1], trigFraction, trigChansCC);

  if(cRefCC1 > 1.f) {
    specTT->Incr(0, CC.nSG+1, int(cRefCC1*tScale*tGain));
    float tCC1  = (cRefCC1-cRefCC)*tScale*tGain;
    specTT->Incr(1, CC.nSG+1, int(tCC1 + speTTlen/2));
  }

  float tsg0 = 0, tsg1 = 0;
  for(int iSegLast = 0; iSegLast < nSegFold; iSegLast++) {
    float tRefSG = FindTriggerTime(fTracesSG[iSegLast], trigFraction, trigChansSG);
    if(iSegLast == 0) tsg0 = tRefSG;
    if(iSegLast == 1) tsg1 = tRefSG;
    if(tRefSG > 1.f) {
      specTT->Incr(0, iSegLast, int(tRefSG*tScale*tGain));        // trigger point for the segment
      float tCCSG  = (tRefSG-cRefCC)*tScale*tGain;                // time relative to the triggering core
      if(CoreE[CC.itrig] > calcTT_eminCC)
        specTT->Incr(1, iSegLast, int(tCCSG + speTTlen/2)); // relative time placed at the center of spectrum
      if(matrETsg) {
        matrETsg->Incr(iSegLast, (int)(fEnergyGainET*CoreE[0]), int(tCCSG+matTTdim/2));  // seg-gamma-time
      }
      if(matrTT) {
        int t1 = (unsigned int)(cRefCC*tScale - 100);    // offset depends on electronics 
        int t2 = (unsigned int)(tRefSG*tScale - 100);    // and should be adjusted manually
        if(t1>0 && t1<matTTdim && t2>0 && t2<matTTdim) {
          unsigned int n1 = (iSegLast/6)*matTTdim;  // sectors
          unsigned int n2 = (iSegLast%6)*matTTdim;  // slices
          matrTT->Incr(t1+n1, t2+n2);
        }
      }
    }
  }
  if(nSegFold==2) {
    specTT->Incr(1, CC.nSG, int((tsg1-tsg0)*tScale*tGain + speTTlen/2)); // relative time placed at the center of spectrum
  }
}

void PreprocessingFilterPSA::calcRiseTime()
{
  if(!PreMatrEeTr)
    return;

  float tRise;
  for(int nnseg = 0; nnseg < nSegFold; nnseg++) {
    int isg = netChargeSegs[nnseg];
    tRise = FindRiseTime(fTracesSG[isg]);
    PreMatrEeTr->Incr(isg, (int)SegE[isg], (int)tRise);
  }

  for(int nc = 0; nc < CC.nCC; nc++) {
    tRise = FindRiseTime(fTracesCC[nc]);
    PreMatrEeTr->Incr(CC.nSG+nc, (int)CoreE[nc], (int)tRise);
  }

  tRise = FindRiseTime(fTS);
  PreMatrEeTr->Incr(CC.nSG+CC.nCC, (int)CoreE[0], (int)tRise);

  tRise = FindRiseTime(fTT);
  PreMatrEeTr->Incr(CC.nSG+CC.nCC+1, (int)CoreE[0], (int)tRise);
}

float PreprocessingFilterPSA::FindRiseTime(float *data)
{
  float valStart = AverageValue(data, nNbase);                   // signal average at the beginning
  float valEnd   = AverageValue(data+nNsamp-2*nNbase, 2*nNbase); // signal average at the end
  float valAmpli = abs(valEnd - valStart);

  float val90 = valAmpli*0.90f;  // *****
  float val10 = valAmpli*0.10f;
  
  int t90 = nNsamp-1;
  int t10 = 0;

  // starting from the top, go down until the first channel below t90
  float level = val90;
  level += valStart;
  for(int nn = nNsamp-2*nNbase; nn > 0; nn--) {
    if(data[nn] < level) {
      t90 = nn;
      break;
    }
  }
  // go down from t90 until the first channel below t10
  level = val10;
  level += valStart;
  for(int nn = t90; nn > 0; nn--) {
    if(data[nn] < level) {
      t10 = nn;
      break;
    }
  }

  return float(t90-t10);
}

void PreprocessingFilterPSA::CheckTiming(float *fWW)
{
  if(!PreMatrTT3)
    return;

#ifdef TIMING_CFD
  float tFract   = mwd->fFractCDF;
  int   tSmooth  = mwd->nSmooth;

  mwd->fFractCDF  = 0.20f;
  mwd->nSmooth    = 0;
  float cRefCC1 = FindTriggerTime(fWW, trigFraction, trigChansCC);

  mwd->fFractCDF  = 0.30f;
  mwd->nSmooth    = 5;
  float cRefCC2 = FindTriggerTime(fWW, trigFraction, trigChansCC);

  mwd->fFractCDF  = tFract;   // leave it as it was
  mwd->nSmooth    = tSmooth;
#else
  float cRefCC1 = FindTriggerTime(fWW, .10f,  5);
  float cRefCC2 = FindTriggerTime(fWW, .10f, 10);
#endif
  if(cRefCC1 > 1.f && cRefCC2 > 1.f) {
    // 1-2
    int iRefCC1 = int(tStep*cRefCC1);
    int iRefCC2 = int(tStep*cRefCC2);
    PreMatrTT3->Incr(0, iRefCC1, iRefCC2);
    // int(1)-2
    iRefCC1 = 10*int(cRefCC1);
    iRefCC2 = int(tStep*cRefCC2);
    PreMatrTT3->Incr(1, iRefCC1, iRefCC2);
    // int(2)-1
    iRefCC1 = int(tStep)*int(cRefCC2);
    iRefCC2 = int(tStep*cRefCC1);
    PreMatrTT3->Incr(2, iRefCC1, iRefCC2);
    // int(1)-1
    iRefCC1 = int(tStep)*int(cRefCC1);
    iRefCC2 = int(tStep*cRefCC1);
    PreMatrTT3->Incr(3, iRefCC1, iRefCC2);
    // int(2)-2
    iRefCC1 = int(tStep)*int(cRefCC2);
    iRefCC2 = int(tStep*cRefCC2);
    PreMatrTT3->Incr(4, iRefCC1, iRefCC2);
  }
}

void PreprocessingFilterPSA::calcEEmatrix(float *eSG, float *eCC, float **sSG, float **sCC)
{
  if(!PreMatrEeEtr)
    return;

  for(int nn = 0; nn < nSegFold; nn++) {
    int ns = netChargeSegs[nn];
    int ee = (int)(matEEgain*eSG[ns]);
    int es = (int)(matEEgain*AverageValue(sSG[ns]+nNsamp-2*nNbase, 2*nNbase));
    PreMatrEeEtr->Incr(ns, ee, es);
  }
  for(int nc = 0; nc < CC.nCC; nc++) {
    int ee = (int)(matEEgain*eCC[nc]);
    int es = (int)(matEEgain*AverageValue(sCC[nc]+nNsamp-2*nNbase, 2*nNbase));
    PreMatrEeEtr->Incr(CC.nSG+nc, ee, es);
  }
}

void PreprocessingFilterPSA::calcSSmatrix(float *eS)
{
  if(!PreMatrEsEs)
    return;

  if(matSSdim==1) {
    for(int s1 = 0; s1 < CC.nSG-1; s1++) {
      if(eS[s1] < CC.pSG[s1].emink) continue;
      for(int s2 = s1+1; s2 < CC.nSG; s2++) {
        if(eS[s2] < CC.pSG[s2].emink) continue;
        PreMatrEsEs->Incr(s1, s2);
        PreMatrEsEs->Incr(s2, s1);
      }
    }
  }
  else {
    static const float fEmax = (matSSdim-1)/fEnergyGainET;
    for(int s1 = 0; s1 < CC.nSG-1; s1++) {
      float e1 = eS[s1];
      if(e1<CC.pSG[s1].emink || e1>fEmax) continue;
      for(int s2 = s1+1; s2 < CC.nSG; s2++) {
        float e2 = eS[s2];
        if(e2<CC.pSG[s2].emink || e2>fEmax) continue;
        int ie1 = s1*matSSdim + (int)(fEnergyGainET*e1);
        int ie2 = s2*matSSdim + (int)(fEnergyGainET*e2);
        PreMatrEsEs->Incr(ie1, ie2);
        PreMatrEsEs->Incr(ie2, ie1);
      }
    }
  }
}

void PreprocessingFilterPSA::calcXTmatrix(float *pSG_Ener)
{
  if(!PreMatrXT)
    return;

  unsigned int iSegLast = netChargeSegs[nSegFold-1];
  int sgEnerY  = int(pSG_Ener[iSegLast]);

  if(sgEnerY > 10 && sgEnerY < 1530) {
    // y axis = triggering segment : units  1   keV/ch   range 10 to 1530 keV
    // x axis = other segments     : units  0.1 keV/ch   range -7 to 3    keV
    for(int ii = 0; ii < CC.nSG; ii++) {
      if(ii == iSegLast) continue;
      int sgEnerX = int(fEnergyGain*pSG_Ener[ii] + 100000.f) - 100000;  // this is like floor()
      if(sgEnerX > -70 && sgEnerX < 30) {
        int valx = (iSegLast*CC.nSG + ii)*100 + (sgEnerX + 70);
        PreMatrXT->Incr(valx, sgEnerY);
      }
    } 
  }
}
#endif  //PPF_MULTIHIST

// Second order calibration for the segments a la "Bart",
// using the xtalk matrix derived from the position of the "negative"
// energy peaks in the multiplicity-one Esegment-Esegment matrices 
// The coefficients are those of the inverted 36x36 matrix

// Calculations done for all segments.
// Sum energy only if above segment threshold, otherwise the sum of noise of all segments spoils the energy resolution)
// The energy resolution is slightly better than xTalkCorrNet() but xTalkCorrNet() is of course faster

// the various 36 should be replaced with CC.nSG

int PreprocessingFilterPSA::xTalkCorrAll(float  *pSG_Ener, float &segsum)
{
  float  xener[36];
  int    segfold = 0;

  segsum = 0;
  for(int n1 = 0; n1 < 36; n1++) {
    float  aa = 0;
    float *pc = CC.xTalkProp[n1];
    for(int n2 = 0; n2 < 36; n2++) {
      aa += pc[n2]*pSG_Ener[n2];
    }
    xener[n1] = aa;               // in order to pass over also the sub-threshold ones
    //if(aa > eMinSG) {             // common threshold set on the corrected energy
    if(aa > CC.pSG[n1].emink) {   // individual threshold set on the corrected energy
      segsum += aa;
      segfold++;
    }
    else
      xener[n1] = 0;    ////// Must check the consequences of this
  }

  memcpy(pSG_Ener, xener, sizeof(xener));
  return segfold;
}

// TO BE CHECKED
// Correcting only net charge segments from net charge segments
// all other segments are set to zero ?????????????????????????
// should remove the dependence from the external netChargeSegs ??
int PreprocessingFilterPSA::xTalkCorrNet(float  *pSG_Ener, float &segsum)
{
  // the singles are treated differently
  segsum = 0;
  if(nSegFold < 2) {
    int i1 = netChargeSegs[0];
    if(nSegFold == 1)
      segsum = pSG_Ener[i1];
    memset(pSG_Ener, 0, sizeof(float)*36);  // 
    if(nSegFold == 1)
      pSG_Ener[i1] = segsum;
    return nSegFold;
  }

  float  xener[36] = {0};

  for(int n1 = 0; n1 < nSegFold; n1++) {
    int   i1 = netChargeSegs[n1];
    float aa = 0;
    for(int n2 = 0; n2 < nSegFold; n2++) {
      int   i2 = netChargeSegs[n2];
      float ee = CC.xTalkProp[i1][i2]*pSG_Ener[i2];
      aa += ee;
    }
    xener[i1] = aa;
    segsum   += aa;
  }

  memcpy(pSG_Ener, xener, sizeof(xener));
  return nSegFold;
}