FarProcessor.cxx

Go to the documentation of this file.

#include <iostream>
#include <math.h>
#include "ReadoutSim/FarProcessor.h"

using namespace rsim;

FarProcessor::FarProcessor()
{
  MIPSignal=20.0; // 20 pe's/MIP
  DarkCurrent=0.250; //Prob of elec in sample period
  ASIC_Noise=1.50; //150e- divided by gain of 100
  ExcessNoise=2.5; // APD amplification excess noise factor
  
  maxADC=4096;

  T_F = 10.0;
  T_R = 2.0;

  pRandom          = new TRandom(100);
}

FarProcessor::~FarProcessor()
{
  delete pRandom;
}

void FarProcessor::Configure(const cfg::Config& c)
{
  c("MIPSignal").Get(MIPSignal);
  c("DarkCurrent").Get(DarkCurrent);
  c("ASIC_Noise").Get(ASIC_Noise);
  c("ExcessNoise").Get(ExcessNoise);
  c("maxADC").Get(maxADC);
  c("T_F").Get(T_F);
  c("T_R").Get(T_R);
  c("MessageLevel").Get(msglevel);
}

void FarProcessor::ProcessChannel(std::vector<unsigned short>& ADC, std::vector<unsigned short>& TDC,
                                  std::vector<sim::HitProfile> & MCInfo)
{
  //All times in units of sample time.  Assumed to be 500ns for FD.
  
  MaxHeight=TMath::Power(T_F/T_R,-T_R/(T_F-T_R))-
    TMath::Power(T_F/T_R,-T_F/(T_F-T_R));  //peak of a pulse with risetime T_R and falltime T_F
  //cout << T_R << " " << T_F  << " " <<MaxHeight << endl;
  
  //  TCanvas *test = new TCanvas("test","test",800,800);
  const Int_t maxArray=1000;
  static Double_t Input[maxArray];
  static Double_t Output[maxArray];
  static Double_t Interp_Out[8*maxArray];
  for (Int_t Interp_index = 0 ; Interp_index < 8*maxArray ; Interp_index++)
    Interp_Out[Interp_index]=0.0;
  static Double_t DCS[maxArray];
  static Double_t Convolution[maxArray];
  static Double_t Index[maxArray];
  static Double_t PerfectSignal[8][20];  //Perfect signal in 8 steps of phase, 20 steps of digitization
  
  static Int_t PSmade=0;
  
  //Clear input array
  for (Int_t i=0;i<maxArray;i++) Input[i]=0.0;
  
  // For MC information. First prepare the vector of pair<index, ADC> for each array component
  sim::HitProfile tempprofile;
  const unsigned int knoise = tempprofile.NoiseId();
  static std::vector< std::pair<unsigned int, Double_t> > InputComponent[maxArray];
  static std::vector< std::pair<unsigned int, Double_t> > OutputComponent[maxArray];
  static std::vector< std::pair<unsigned int, Double_t> > Interp_OutComponent[8*maxArray];
  static std::vector< std::pair<unsigned int, Double_t> > ConvolutionComponent[maxArray];
  
  //Clear MC array
  for (Int_t i=0;i<maxArray;i++) {
    InputComponent[i].clear();
    OutputComponent[i].clear();
  }
  for (Int_t Interp_index = 0 ; Interp_index < 8*maxArray ; Interp_index++)
    Interp_OutComponent[Interp_index].clear();
  
  //Add signal
  // Input[20]=1.0*MIPSignal;
  //  Input[50]=1.0*MIPSignal;
  //  Input[60]=1.0*MIPSignal;
  //  Input[150]=0.50*MIPSignal;
  //  Input[300]=1.0*MIPSignal;
  //  Input[250]=1.0*MIPSignal;
  //  Input[350]=5.0*MIPSignal;
  
  //really add signals
  
  static int startsignal=50;
  
  std::pair<unsigned int, Double_t> tempindex;
  
  for(unsigned int i=0;i<ADC.size();i++)
  {
    //TDC is recorded in ns... each bin in 500 ns
    //first TDC should be in at least bin 50 to allow for ramp up of any noise
    int tbin = (TDC[i]-TDC[0])/500+startsignal;
    if (tbin>maxArray) 
    {
      //printf("error.... length of signal exceeds array size\n");
      continue;
    }
    
    Input[tbin]+=ADC[i];
    tempindex.first  = i;
    tempindex.second = ADC[i];
    InputComponent[tbin].push_back(tempindex);
    //printf("%d %f - ",i, TDC[i]);
  }//printf("\n");
  
  TDC.clear();
  ADC.clear();
  
  std::vector<sim::HitProfile> copyMCInfo = MCInfo;
  MCInfo.clear();
  
  
  //printf("adding dark current\n");
  
  for (Int_t i=0;i<maxArray;i++) {
    //add random dark current
    Double_t darkcurrent = pRandom->Poisson(DarkCurrent);
    Input[i]+=darkcurrent;
    
    tempindex.first = knoise;
    tempindex.second = darkcurrent;
    InputComponent[i].push_back(tempindex);
    
    if (Input[i]>0.0) {
      //smear by apd excess noise response
      Double_t smeared = pRandom->Gaus(Input[i],TMath::Sqrt(Input[i]*(ExcessNoise-1))); 
      
      for(unsigned int j=0; j<InputComponent[i].size(); j++) { 
        InputComponent[i][j].second *= smeared/Input[i];
      }
      
      Input[i] = smeared;
    }
    Output[i]=0.0;
    Index[i]=(Double_t)i;
  }
  
  // Create a perfect signal in 8 steps of phase to match 8x interpolation.
  
  //printf("perfect signal\n");
  
  if(PSmade<1)
    for (Int_t phase=0;phase<8;phase++){
      PerfectSignal[phase][0]=0.0;
      // PerfSig->SetBinContent(0,phase,0.0);
      for (Int_t step=1;step<20;step++){
        PerfectSignal[phase][step]=exp(-(step-1+(float)phase/8.0)/T_F)*(1-exp(-(step-1+(float)phase/8.0)/T_R));
        //  PerfSig->SetBinContent(step+1,phase+1,PerfectSignal[phase][step]);   
      }
    }
  PSmade=1;
  
  //printf("output and dcs\n");
  
  Double_t ConvFudge = 0.2; //Fudge factor to get convolution on same scale.
  
  for (Int_t a=0 ;  a < 5*T_F; a++){
    Double_t fun=exp(-(a)/T_F)*(1-exp(-(a)/T_R))/MaxHeight;
    for (Int_t i=0; i<maxArray ; i++){
      if(i-a>=0) {
        Output[i]+=Input[i-a]*fun;              
        
        if(fun==0) continue;

        for(unsigned int j=0; j<InputComponent[i-a].size(); j++) {

          tempindex.first  = InputComponent[i-a][j].first;
          tempindex.second = fun * InputComponent[i-a][j].second;

          bool duplicated = false;
          for(unsigned int k=0; k<OutputComponent[i].size(); k++) {
            if(OutputComponent[i][k].first == tempindex.first) {
              OutputComponent[i][k].second += tempindex.second;
              duplicated = true;
            }
          }
          if (!duplicated) OutputComponent[i].push_back(tempindex);
        }

      }
    }
  }
  
  tempindex.first = knoise;
  
  Double_t asicnoise;
  
  asicnoise = pRandom->Gaus(0,ASIC_Noise);
  Output[0] += asicnoise;
  tempindex.second = asicnoise;
  OutputComponent[1].push_back(tempindex);
  
  asicnoise=pRandom->Gaus(0,ASIC_Noise);
  Output[1] += asicnoise;
  tempindex.second = asicnoise;
  OutputComponent[2].push_back(tempindex);
  
  for (Int_t i=2; i<maxArray ; i++){
    asicnoise = pRandom->Gaus(0,ASIC_Noise);
    Output[i] += asicnoise;
    tempindex.second = asicnoise;
    OutputComponent[i].push_back(tempindex);
    
    DCS[i-1]=Output[i]-Output[i-2];
  }
  
  // clock_t tho;
  //tho=clock();
  
  //printf("int out\n");
  
  Double_t x_interp[640];
  
  //this should be put in a seperate static class
  for (Int_t Interp_index = -320;  Interp_index < 320 ; Interp_index++)
  {
    x_interp[Interp_index+320]=(Double_t)Interp_index*TMath::Pi()/8.0;
    
    if( fabs(x_interp[Interp_index+320]) > 1.0E-3)
      x_interp[Interp_index+320]=1.0*TMath::Sin(x_interp[Interp_index+320])/x_interp[Interp_index+320];
    else
      x_interp[Interp_index+320]=1;
    
  }
  //Interpolate Output to add 8X the points with sin(x)/x
  //Truncate after 40 periods of sinx oscillation
  for (Int_t Interp_index = -320;  Interp_index < 320 ; Interp_index++)
  {
    //Double_t x_interp=(Double_t)Interp_index*TMath::Pi()/8.0;
    
    //if( fabs(x_interp) > 1.0E-3)
    //  {
    //          Double_t fun=1.0*TMath::Sin(x_interp)/x_interp;
    
    //for (Int_t Sample_index = 40 ; Sample_index < maxArray && Sample_index*8+Interp_index >0 && Sample_index*8+Interp_index < maxArray ; Sample_index++)
    for (Int_t Sample_index = -Interp_index/8 ; Sample_index < (maxArray -Interp_index)/8  ; Sample_index++)
    {     
      Interp_Out[Sample_index*8+Interp_index] +=
        Output[Sample_index]*x_interp[Interp_index+320];
      /*
      for(unsigned int j=0; j<OutputComponent[Sample_index].size(); j++) {
        tempindex.first  = OutputComponent[Sample_index][j].first;
        tempindex.second = x_interp[Interp_index+320] * OutputComponent[Sample_index][j].second;
        Interp_OutComponent[Sample_index*8+Interp_index].push_back(tempindex);
      }
      */
    }
    
    //} else {
    
    //  for (Int_t Sample_index = -Interp_index/8 ; Sample_index < (maxArray -Interp_index)/8  ; Sample_index++)
    //  Interp_Out[Sample_index*8+Interp_index] += Output[Sample_index];
    
    //}
  }
  
  //tho-=clock(); //printf ("It took you %d ticks \n", -tho );
  
  //printf("conv\n");
  for (Int_t i=0; i<maxArray ; i++){
    //Convolution with perfect signal (Still needs to use all phases.)
    Convolution[i]=0.0;
    ConvolutionComponent[i].clear();
    if (i>20) {
      for (Int_t i_t=i;i-i_t<20;i_t--) {
        
        Convolution[i] += ConvFudge*PerfectSignal[0][20-(i-i_t)]*Output[i_t];

        for(unsigned int j=0; j<OutputComponent[i_t].size(); j++) {
          
          tempindex.first  = OutputComponent[i_t][j].first;
          tempindex.second = ConvFudge * PerfectSignal[0][20-(i-i_t)] * OutputComponent[i_t][j].second;
          
          bool duplicated = false;
          for(unsigned int k=0; k<ConvolutionComponent[i].size(); k++) {
            if(ConvolutionComponent[i][k].first == tempindex.first) {
              ConvolutionComponent[i][k].second += tempindex.second;
              duplicated = true;
            }
          }
          if (!duplicated) {
            ConvolutionComponent[i].push_back(tempindex);
          }
        }

      }  
    }
  }    

  //printf("making output\n");
  
  Int_t maxbin=-1;
  Float_t maxval=-1;
  Float_t thresh=20;

  //std::vector<float> ADC(0);
  //std::vector<float> TDC(0);
  
  for (Int_t i=0; i<maxArray ; i++){
    if(Convolution[i]>thresh)
    {
      if(Convolution[i]>maxval)
      {
        maxbin=i;
        maxval=Convolution[i];
      }
    }else{
      if(maxval>0)
      {
        ADC.push_back(maxval>maxADC ? maxADC : (short)maxval);
        TDC.push_back((short)(maxbin-19-startsignal));  

        tempprofile.Clear();
        
        for(unsigned int j=0; j<ConvolutionComponent[maxbin].size(); j++) {
          
          unsigned int hitindex = ConvolutionComponent[maxbin][j].first;
          
          if(hitindex==knoise) {
            tempprofile.Push(tempprofile.NoiseId(),ConvolutionComponent[maxbin][j].second);
          } else {
            for(unsigned int k=0; k<copyMCInfo[hitindex].NComponents(); k++) {
              tempprofile.Push(copyMCInfo[hitindex].Component(k,(float)ConvolutionComponent[maxbin][j].second));
            }
          }       
        }
        
        tempprofile.Combine();
        tempprofile.NormalizeFraction();

        MCInfo.push_back(tempprofile);
        
      }
      maxbin=-1;
      maxval=-1;
    }
  }
  
  if(maxval>0)
  {
    ADC.push_back(maxval>maxADC ? maxADC : (short)maxval);
    TDC.push_back((short)(maxbin-19-startsignal));      
    
    tempprofile.Clear();
    
    for(unsigned int j=0; j<ConvolutionComponent[maxbin].size(); j++) {
      
      unsigned int hitindex = ConvolutionComponent[maxbin][j].first;
      
      if(hitindex==knoise) {
        tempprofile.Push(tempprofile.NoiseId(),ConvolutionComponent[maxbin][j].second);
      } else {
        for(unsigned int k=0; k<copyMCInfo[hitindex].NComponents(); k++)
          tempprofile.Push(copyMCInfo[hitindex].Component(k,(float)ConvolutionComponent[maxbin][j].second));
      }   
    }
    
    tempprofile.Combine();
    tempprofile.NormalizeFraction();

    MCInfo.push_back(tempprofile);
  }

  maxbin=-1;
  maxval=-1;
  
  if(ADC.size()<1)
  {
    ADC.push_back(0);
    TDC.push_back(0);
  } else if (ADC.size()>1) {
    //printf("possible error... multiple peaks... should dump waveform\n");
  } 
  
  /*
    for(unsigned int ii=0;ii<ADC.size();ii++)
    {
    printf("peak at %d time %d\n",ADC[ii], TDC[ii]);
    }
  */    

}

---