HarpoSimElectronics.cxx

Go to the documentation of this file.

#include "HarpoSimElectronics.h"
#include "HarpoSimStrips.h"
#include "HarpoDccMap.h"
#include "HarpoConfig.h"
#include "HarpoDebug.h"

#include "TMath.h"
#include "TRandom.h"
#include <iostream>

ClassImp(HarpoSimElectronics)

HarpoSimElectronics::HarpoSimElectronics()
{

  fNoise = 10.;
  // fNoiseAGET = 0.5;
  fZeroSuppr = 5;

  fGammaPower = 3;
  fTau = 0.120;       // Shaping time [us]
  fSamplingFreq = 50; // Sampling frequency [MHz]
  fTimeRange = 511;
  fTimeOffset = 0.;      // trigger time in us

  fThrSat = 4096;
  //  fThrSatCh = 5e4;
  fThrSatCh = 26600;
  fAmpSatCh = 28900;
  fXtalkSatCh = 0.;

  Double_t noise;
  if ( ! gHConfig->Lookup("Sim.Noise",noise) )
    Info("Constructor","Use default Noise %g",fNoise);
  else
    fNoise = noise;

  Long64_t zeroSuppr;
  if ( ! gHConfig->Lookup("Sim.ZeroSuppr",zeroSuppr) )
    Info("Constructor","Use default ZeroSuppr %i",fZeroSuppr);
  else
    fZeroSuppr = zeroSuppr;

  Double_t tau;
  if ( ! gHConfig->Lookup("Sim.Tau",tau) )
    Info("Constructor","Use default Tau %g",fTau);
  else
    fTau = tau;

  Long64_t gammaPower;
  if ( ! gHConfig->Lookup("Sim.GammaPower",gammaPower) )
    Info("Constructor","Use default GammaPower %d",fGammaPower);
  else
    fGammaPower = gammaPower;

  Double_t samplingFreq;
  if ( ! gHConfig->Lookup("Sim.SamplingFreq",samplingFreq) )
    Info("Constructor","Use default SamplingFreq %g",fSamplingFreq);
  else
    fSamplingFreq = samplingFreq;

  Long64_t timeRange;
  if ( ! gHConfig->Lookup("Sim.TimeRange",timeRange) )
    Info("Constructor","Use default TimeRange %i",fTimeRange);
  else
    fTimeRange = timeRange;

  Double_t timeOffset;
  if ( ! gHConfig->Lookup("Sim.TimeOffset",timeOffset) )
    Info("Constructor","Use default TimeOffset %g",fTimeOffset);
  else
    fTimeOffset = timeOffset;

  Double_t thrSat;
  if ( ! gHConfig->Lookup("Sim.ThrSat",thrSat) )
    Info("Constructor","Use default ThrSat %g",fThrSat);
  else
    fThrSat = thrSat;

  Double_t thrSatCh;
  if ( ! gHConfig->Lookup("Sim.ThrSatCh",thrSatCh) )
    Info("Constructor","Use default ThrSatCh %g",fThrSatCh);
  else
    fThrSatCh = thrSatCh;

  Double_t ampSatCh;
  if ( ! gHConfig->Lookup("Sim.AmpSatCh",ampSatCh) )
    Info("Constructor","Use default AmpSatCh %g",fAmpSatCh);
  else
    fAmpSatCh = ampSatCh;

  Double_t xtalkSatCh;
  if ( ! gHConfig->Lookup("Sim.XtalkSatCh",xtalkSatCh) )
    Info("Constructor","Use default XtalkSatCh %g",fXtalkSatCh);
  else
    fXtalkSatCh = xtalkSatCh;

  fPedestals[XDCC] = new HarpoPedestal();
  fPedestals[XDCC]->SetMean(0.);
  fPedestals[XDCC]->SetSigma(fNoise);
  fPedestals[YDCC] = new HarpoPedestal();
  fPedestals[YDCC]->SetMean(0.);
  fPedestals[YDCC]->SetSigma(fNoise);


  fDccMapX = new HarpoDccMap();
  fDccMapY = new HarpoDccMap();
  //  std::cout << "Pedestals: " << fPedestals[XDCC] << std::endl;
  fDccMapX->SetPedestals(fPedestals[XDCC]);
  fDccMapX->SetPedestals(fPedestals[YDCC]);

}

void HarpoSimElectronics::Reset()
{
  //  std::cout << "fDccMapX " << fDccMapX << std::endl;

  fDccMapX = new HarpoDccMap();
  fDccMapX->SetPedestals(fPedestals[XDCC]);
  fDccMapY = new HarpoDccMap();
  fDccMapY->SetPedestals(fPedestals[YDCC]);

}

Int_t HarpoSimElectronics::AddCharge(Int_t plane, Int_t channel, Double_t time, Double_t charge)
{  

  if(plane == XPLANE) return AddCharge(fDccMapX,channel,time,charge);
  if(plane == YPLANE) return AddCharge(fDccMapY,channel,time,charge);
  return 0 ; // ????
}

Int_t HarpoSimElectronics::AddCharge(HarpoDccMap* dccMap, Int_t channel, Double_t time, Double_t charge)
{  

  if(gHarpoDebug>2) Info("AddCharge","###Pedestals: %p",dccMap->GetPedestals());
  if(channel<0) return -1;
  if(channel>=288) return -1;

  if(gHarpoDebug>2) Info("AddCharge","ch = %d, t = %.1g, Q = %.1g",channel,time,charge);
  Double_t minCharge = 0.01;

  Double_t tautmp = fTau*fSamplingFreq;

  Double_t chargeSum = 0;
  
  time += fTimeOffset;


  Double_t data = 0;
  Bool_t signal = kFALSE;
  Int_t tmin = Int_t(time*fSamplingFreq);
  if(tmin<0) tmin = 0;
  for(Int_t t = tmin; t<fTimeRange; t++){
    Double_t cellCharge = 0;
    Double_t t1 = (t-time*fSamplingFreq)/tautmp;
    if(t1<0) continue;   
    switch(fGammaPower){
    case 4: cellCharge = charge/(3.0/128.0*tautmp)*TMath::Power(t1,4)*TMath::Exp(-4*t1); break;
    case 3: cellCharge = charge/(2.0/27.0*tautmp)*TMath::Power(t1,3)*TMath::Exp(-3*t1); break;
    case -3: cellCharge = charge/(2.0/27.0*tautmp)*TMath::Power(t1,3)*TMath::Exp(-3*t1)*TMath::Sin(t1)/0.7776; break;
    case 2: cellCharge = charge/(1.0/4.0*tautmp)*TMath::Power(t1,2)*TMath::Exp(-2*t1); break;
    default: cellCharge = charge/(3.0/128.0*tautmp)*TMath::Power(t1,4)*TMath::Exp(-4*t1);
    }
    data = dccMap->GetRawData(channel,t);
    if(data <= 0) data = cellCharge;
    else       data += cellCharge;
    if(data<minCharge && signal) break;
    if(data<minCharge && !signal) continue;
    signal = kTRUE;
    
    if(gHarpoDebug>4) Info("AddCharge","map(%d,%d) = %g",channel,t,data);
    //    dccMap->map(channel,t) = data;
    dccMap->SetRawData(channel,t,data);
    chargeSum += cellCharge;
  }

  if(gHarpoDebug>3) Info("AddCharge","*** In: %.3g, OUT: %.3g   =>   %.3g",charge, chargeSum, chargeSum/charge);

  return 1;
}




void HarpoSimElectronics::Noise()
{
  Noise(fDccMapX);
  Noise(fDccMapY);
  Saturation(fDccMapX);
  Saturation(fDccMapY);
}

void HarpoSimElectronics::Noise(HarpoDccMap* dccMap)
{
  for(Int_t ch = 0; ch<NCHAN;ch++){
    for(Int_t t = 0; t<NADC;t++){
      //      Double_t data = dccMap->map(ch,t);
      Double_t data = dccMap->GetRawData(ch,t);
      data += gRandom->Gaus(0,fNoise);
      //      dccMap->map(ch,t) = Int_t(data);
      dccMap->SetRawData(ch,t, Int_t(data));
    }
  }
}

void HarpoSimElectronics::Saturation(HarpoDccMap* dccMap)
{

  Bool_t isSat = kTRUE;
  for(Int_t ch = 0; ch<NCHAN;ch++){
    Double_t qtot = 0, qreadold = 0; 
    Bool_t isPrevSat = isSat;
    isSat = kFALSE;
    for(Int_t t = 0; t<NADC;t++){
      //      Double_t data = dccMap->map(ch,t);
      Double_t data = dccMap->GetRawData(ch,t);
      qtot += data;
      Double_t qread = ChannelSat(qtot);
      data = qread - qreadold;
      qreadold = qread;
      if(qtot>fThrSatCh){ // Cross-talk when saturated
        isSat = kTRUE;
        if(!isPrevSat && ch>0){
          Double_t dataprev = dccMap->GetRawData(ch-1,t)+data*fXtalkSatCh;
          if(dataprev>fThrSat) dataprev = fThrSat;
          dccMap->SetRawData(ch-1,t, Int_t(dataprev));
        }
        if(ch<NCHAN-1){
          Double_t datanext = dccMap->GetRawData(ch+1,t)+data*fXtalkSatCh;
          if(datanext>fThrSat) datanext = fThrSat;
          dccMap->SetRawData(ch+1,t, Int_t(datanext));
        }
      }
      if(data>fThrSat) data = fThrSat;
      dccMap->SetRawData(ch,t, Int_t(data));
    }
  }
}

Double_t HarpoSimElectronics::ChannelSat(Double_t qreal)
{
  // if(qreal<fThrSatCh) return qreal;
  Double_t qread;
  if(fAmpSatCh>0){
    if(qreal<fThrSatCh) return qreal;
    qread = fThrSatCh + fAmpSatCh*TMath::TanH((qreal-fThrSatCh)/fAmpSatCh);
    return qread;
  }else{
    qread = fThrSatCh + (qreal-fThrSatCh)/(-fAmpSatCh);
    return qread;
  }
  
  // Double_t A = 0.0357923;
  // Double_t B = 103545;
  // Double_t C = (1-A);
  // Double_t lambda = 5.;
  // qread = qreal*C/TMath::Power(1+TMath::Power(qreal/B*C,lambda),1./lambda) + qreal*A;
  // return qread;
  Double_t thr = 40000;
  if(qreal<thr) return qreal;
  Double_t A = 0.8;
  Double_t thr2 = 60000;
  Double_t thr2real = thr+(thr2-thr)/A;
  if(qreal<thr2real) return thr + (qreal-thr)*A;
  Double_t A2 = 0.4;
  Double_t thr3 = 100000;
  Double_t thr3real = thr2real+(thr3-thr2)/A2;
  if(qreal<thr3real) return thr2 + (qreal-thr2real)*A2;
  Double_t A3 = 0.03;
  return thr3 + (qreal-thr3real)*A3;
}

void HarpoSimElectronics::ZeroSuppression()
{
  ZeroSuppression(fDccMapX);
  ZeroSuppression(fDccMapY);
}

void HarpoSimElectronics::ZeroSuppression(HarpoDccMap* dccMap)
{
  for(Int_t ch = 0; ch<304;ch++){
    for(Int_t t = 0; t<511;t++){
      //      Double_t data = dccMap->map(ch,t);
      Double_t data = dccMap->GetData(ch,t);
      if(data>4096) data = 4096; 
      if(data < fZeroSuppr*fNoise)
        dccMap->SetRawData(ch,t, -1);
      else
        dccMap->SetData(ch,t, data);
    }
  }
}

// void HarpoSimElectronics::AGETmultiplicity()
// {

//   for(Int_t t = 0; t<NADC;t++){
//     Double_t mult[fNchipAGET];
//     for(Int_t i = 0; i<fNchipAGET; i++) mult[i] = 0;
//     for(Int_t ch = 0; ch<304;ch++){
//       Double_t dataX = fDccMapX->GetData(ch,t);
//       if(dataX>fThrAGET){
//      Int_t chip = GetAGETchip(XPLANE,ch);
//      mult[chip] += 1 + gRandom->Gaus(0,fNoiseAGET);
//       }
//       Double_t dataY = fDccMapY->GetData(ch,t);
//       if(dataY>fThrAGET){
//      Int_t chip = GetAGETchip(YPLANE,ch);
//      mult[chip] += 1 + gRandom->Gaus(0,fNoiseAGET);
//       }
//     }
//     for(Int_t i = 0; i<fNaget; i++){
//       if(mult[i] == 0) continue;
//       fAGET->Add(i,t,mult[i]);
//     }
//   }
// }