source: trunk/MagicSoft/Simulation/Detector/StarResponse/MStarLight.cxx@ 20115

#include "MStarLight.hxx" 

using namespace std;

MStarLight::MStarLight(Float_t fadc_slices_per_ns, Int_t response_slices_fadc) {
  //
  //   default constructor
  //

  fBrightness  = 0.; 
  fTimeRange   = TIMERANGE; 

  fFadcSlicesPerNanosec = fadc_slices_per_ns; // "real" number of FADC slices per ns

  fResponseSlicesFadc = response_slices_fadc; 
  // total number of bins in the histogram containing the response of the FADC to a
  // single photoelectron. The bins are narrower than the true FADC slices by a factor
  // equal to SUBBINS (see MFadcDefine.h)


  fBinsTrig = (Int_t)(TRIG_SLICES_PER_NSEC*fTimeRange);   // Default value 4*10000=40000
  fTrigShape   = 0; 
  fAmplTrig    = 0.; 
  fFwhmTrig    = 0.; 

  fBinsFadc    = (Int_t)(fFadcSlicesPerNanosec*fTimeRange);  // Default value 0.3*10000=3000
  fFadcShape   = 0; 
  fIntegFadc   = 0.; 
  fFwhmFadc    = 0.; 

  fGainFluctuations = 1;

  fTrig     = new Float_t[fBinsTrig];
  fTrigResp = new Float_t[RESPONSE_SLICES_TRIG];
  fFadc     = new Float_t[fBinsFadc];
  fFadcResp = new Float_t[fResponseSlicesFadc]; 
  
  for (Int_t i= 0; i< fBinsTrig ; i++) 
    fTrig[i] = 0. ; 

  for (Int_t i= 0; i < fBinsFadc; i++) 
    fFadc[i] = 0.; 
  
  for (Int_t i = 0; i < RESPONSE_SLICES_TRIG; i++) 
    fTrigResp[i] = 0.; 

  for (Int_t i = 0; i < fResponseSlicesFadc; i++) 
    fFadcResp[i] = 0.; 
}

void MStarLight::Reset() {

  fBrightness  = 0. ;

  fTrigShape   = 0; 
  fAmplTrig    = 0.; 
  fFwhmTrig    = 0.; 

  fFadcShape   = 0; 
  fIntegFadc   = 0.; 
  fFwhmFadc    = 0.; 

  fGainFluctuations = 1;
  
  for (Int_t i= 0; i < fBinsTrig ; i++) 
    fTrig[i] = 0.; 

  for (Int_t i= 0; i < fBinsFadc; i++) 
    fFadc[i] = 0.; 
  
  for (Int_t i = 0; i < RESPONSE_SLICES_TRIG; i++) 
    fTrigResp[i] = 0.; 

  for (Int_t i = 0; i < fResponseSlicesFadc; i++) 
    fFadcResp[i] = 0.; 
}

Float_t MStarLight::GetBrightness () 
{
  return fBrightness; 
}

void MStarLight::SetBrightness (Float_t in ) 
{
  fBrightness = in; 
}

Float_t MStarLight::GetAmplTrig () 
{
  return fAmplTrig ; 
}

void MStarLight::SetAmplTrig (Float_t in ) 
{
  fAmplTrig = in; 
}

Float_t  MStarLight::GetFwhmTrig () 
{
  return fFwhmTrig; 
}

void MStarLight::SetFwhmTrig (Float_t in ) 
{
  fFwhmTrig = in; 
}

void MStarLight::SetFadcSlicesPerNanosec (Float_t in)
{
  fFadcSlicesPerNanosec = in;
  fBinsFadc    = (Int_t)(fFadcSlicesPerNanosec*fTimeRange);

  if (fFadc)
    delete [] fFadc;

  fFadc = new Float_t[fBinsFadc];

  for (Int_t i= 0; i < fBinsFadc; i++) 
    fFadc[i] = 0.;  
}

Float_t MStarLight::GetIntegFadc () 
{
  return fIntegFadc; 
}

void MStarLight::SetIntegFadc (Float_t in ) 
{
  fIntegFadc = in; 
}

Float_t  MStarLight::GetFwhmFadc () 
{
  return fFwhmFadc; 
}

void MStarLight::SetFwhmFadc (Float_t in ) 
{
  fFwhmFadc = in; 
}


void MStarLight::SetTrigResponse( Float_t *in )
{
  for (Int_t i = 0; i < RESPONSE_SLICES_TRIG; i++) 
    fTrigResp[i] = in[i]; 
}

void MStarLight::SetFadcResponse( Float_t *in )
{
  for (Int_t i = 0; i < fResponseSlicesFadc; i++)
    fFadcResp[i] = in[i]; 
}

void MStarLight::FillResponse( Float_t ampl, Float_t time )
{
  //  fill the trigger response
  
  Int_t startbin = (Int_t) (time * ((Float_t)fBinsTrig/fTimeRange));

  Int_t icount = 0;
  Int_t idata; 
  
  for ( Int_t i = startbin ; i < startbin+RESPONSE_SLICES_TRIG ; i++)
    {
      if ( i < fBinsTrig ) 
        idata = i; 
      else if ( i >= fBinsTrig ) 
        idata = i - fBinsTrig; 

      fTrig[idata] = fTrig[idata] + ampl * fTrigResp[icount];
      icount++ ; 
    }

  //
  //   fill the FADC content
  //

  startbin = (Int_t) ( time * ((Float_t)(fBinsFadc*SUBBINS)/fTimeRange));

  Int_t ichanfadc = 0 ; 
  
  //
  //   putting the response slices in the right sig slices.
  //   Be careful, because both slices have different widths. 
  //

  //
  // Changed, Jan 2004, A. Moralejo: now the FADC does not integrate all
  // the signal within each FADC slice, but measures just the signal height
  // at one point, like the real FADC does. By default, each FADC slice 
  // contains SUBBINS bins of the response histogram
  // (fFadcResp). Warning: do not change this unless you do the corresponding 
  // modifications also in MFadc.cxx, or the signal and the noise photoelectrons 
  // will be treated differently!!
  //

  for ( Int_t i = 0 ; i < fResponseSlicesFadc; i += SUBBINS )
    {
      ichanfadc = (Int_t) ((startbin+i)/(Float_t)SUBBINS); 

      if ( ichanfadc < fBinsFadc ) 
        idata = ichanfadc; 
      else if ( ichanfadc >= fBinsFadc ) 
        idata = ichanfadc-fBinsFadc; 

      fFadc[idata] += (ampl * fFadcResp[i]); 
    }
}

void MStarLight::ElecNoise ( Float_t noiseTrig, Float_t noiseFadc ) {
  //
  //     putting some noise to the baseline
  // 

  TRandom2 wuerfel( (UInt_t) (noiseTrig*100) ) ; 
 
  for (Int_t i=0; i< fBinsTrig ; i++ )
    fTrig[i] += wuerfel.Gaus(0., noiseTrig );

  for (Int_t i=0; i< fBinsFadc ; i++ )
    fFadc[i] += wuerfel.Gaus(0., noiseFadc );
  
}

Float_t MStarLight::GetTrig( Int_t i)
{
  //------------------------------------------------------------------
  //
  //  It gets the value of the simulated trigger in the i bin
  //
  return fTrig[i];
}

Float_t MStarLight::GetFadc( Int_t i)
{
  //------------------------------------------------------------------
  //
  //  It gets the value of the simulated FADC signal in the i bin
  //
  return fFadc[i];
}

void MStarLight::StoreHisto( char *filename)
{

  Float_t baseline = 0.; 

  //    first the histograms for trigger
  //
  //    the calculated trigger signal before baseline
  // 

  TH1F  trigresp ("trigresp", "Trigger Response", fBinsTrig, 0., fTimeRange);

  for (Int_t i=0; i< fBinsTrig ; i++ )
    {
      trigresp.SetBinContent(i+1, fTrig[i]);
      baseline += fTrig[i];
    }

  baseline /= fBinsTrig;
  
  TH1F  trigbase ("trigbase", "Response after Baseline shift", 
                  fBinsTrig, 0., fTimeRange) ; 

  for (Int_t i = 0; i < fBinsTrig ; i++)
    trigbase.SetBinContent(i+1, fTrig[i]-baseline );
  
  TH1F  trigdist ("trigdist", "Noise on the baseline", 
                  1000, -25., 25. );

  for (Int_t i = 0; i < fBinsTrig ; i++)
    trigdist.Fill( (Float_t) trigbase.GetBinContent(i+1) );


  // 
  // Now the histograms for the fadc
  //

  TH1F  fadcresp ("fadcresp", "Fadc Response", fBinsFadc, 0., fTimeRange); 

  baseline = 0.;
  for (Int_t i=0; i < fBinsFadc; i++)
    {
      fadcresp.SetBinContent(i+1, fFadc[i]); 
      baseline += fFadc[i]; 
    }

  baseline /= fBinsFadc; 
  
  TH1F  fadcbase ("fadcbase", "Fadc after Baseline shift", 
                  fBinsFadc, 0., fTimeRange) ; 

  for (Int_t i=0; i< fBinsFadc ; i++ )
    fadcbase.SetBinContent(i+1, fFadc[i]-baseline ); 

  TH1F  fadcdist ("fadcdist", "Noise on fadc", 
                  1000, -100., 100. ) ; 

  for (Int_t i=0; i< fBinsFadc ; i++ )
    fadcdist.Fill( (Float_t) fadcbase.GetBinContent(i+1) ); 

  TFile outfile( filename, "RECREATE");
  
  trigresp.Write(); 
  trigbase.Write(); 
  trigdist.Write(); 

  fadcresp.Write(); 
  fadcbase.Write(); 
  fadcdist.Write(); 
  
  
  outfile.Close(); 
}

void MStarLight::WriteBinary( char *filename)
{
  //
  //  write the information to the binary file

  FILE *datei ; 

  datei = fopen ( filename, "w" ) ; 
  
  if ( ! datei )
    {
      cout << " ERROR: Can't open the file " << filename << endl;
      exit (230);
    }

  Float_t version = VERSIONSR;

  //   write them out 

  fwrite ( &version,               sizeof(Float_t), 1, datei ); 
  fwrite ( &fBrightness,           sizeof(Float_t), 1, datei ); 
  fwrite ( &fTimeRange,            sizeof(Float_t), 1, datei ); 
  fwrite ( &fFadcSlicesPerNanosec, sizeof(Float_t), 1, datei ); 

  fwrite ( &fBinsTrig  , sizeof(Int_t),   1, datei ); 
  fwrite ( &fTrigShape , sizeof(Int_t),   1, datei ); 
  fwrite ( &fAmplTrig  , sizeof(Float_t), 1, datei ); 
  fwrite ( &fFwhmTrig  , sizeof(Float_t), 1, datei ); 

  fwrite ( &fBinsFadc  , sizeof(Int_t),   1, datei ); 
  fwrite ( &fFadcShape , sizeof(Int_t),   1, datei ); 
  fwrite ( &fIntegFadc , sizeof(Float_t), 1, datei ); 
  fwrite ( &fFwhmFadc  , sizeof(Float_t), 1, datei ); 

  fwrite ( &fGainFluctuations, sizeof(Int_t), 1, datei ); 

  fwrite (fTrig, sizeof(Float_t), fBinsTrig, datei); 

  //    We want to store the FADC signal taking into account the AC
  //    coupling
  //
  //    We calculate and substract the baseline
  // 

  Float_t baseline = 0.; 

  for (Int_t i=0; i< fBinsFadc ; i++ )
    baseline += fFadc[i]; 

  baseline /= fBinsFadc; 

  Float_t *temp = new Float_t[fBinsFadc];

  for (Int_t i=0; i < fBinsFadc; i++ )
    temp[i] = fFadc[i]-baseline;
  
  fwrite (temp, sizeof(Float_t), fBinsFadc, datei); 

  delete [] temp;

  fclose ( datei ); 
  
}

void MStarLight::ReadBinary( char *filename) {
  
  //
  //  read the things from the binary file

  FILE *datei ; 

  datei = fopen ( filename, "r" ) ; 
  
  if ( ! datei ) {
    cout << " ERROR: Can't open the file " << filename 
         << endl ; 
    cout<<  "        The database for the NSB may be too small. "
        <<  endl;
    cout<<  "        See the How to Use of the Camera simulation for more information"
        <<  endl;
    exit (230) ; 
  }

  Float_t read_version;
  Float_t current_version = VERSIONSR;

  // Check that we read the right version of the Database
  fread ( &read_version, sizeof(Float_t), 1, datei ) ; 

  if(Int_t(read_version)!=Int_t(current_version)){
    cout<<" ERROR: You are trying to read database VERSION "<<
      read_version << endl;
    cout<<"        You must generate a database for the current VERSION "<<
      current_version << endl;
    cout<<"        See the NSB database section in the Camera How to Use note."
        <<endl;
    exit (230);
  }

  fread ( &fBrightness, sizeof(Float_t), 1, datei ); 
  fread ( &fTimeRange , sizeof(Float_t), 1, datei ); 
  fread ( &fFadcSlicesPerNanosec, sizeof(Float_t), 1, datei ); 

  fread ( &fBinsTrig  , sizeof(Int_t),   1, datei ); 
  fread ( &fTrigShape , sizeof(Int_t),   1, datei ); 
  fread ( &fAmplTrig  , sizeof(Float_t), 1, datei ); 
  fread ( &fFwhmTrig  , sizeof(Float_t), 1, datei ); 

  fread ( &fBinsFadc  , sizeof(Int_t),   1, datei ); 
  fread ( &fFadcShape , sizeof(Int_t),   1, datei ); 
  fread ( &fIntegFadc , sizeof(Float_t), 1, datei ); 
  fread ( &fFwhmFadc  , sizeof(Float_t), 1, datei ); 

  fread ( &fGainFluctuations, sizeof(Int_t), 1, datei ); 

  if (fTrig)
    delete [] fTrig;

  if (fFadc)
    delete [] fFadc;

  fTrig = new Float_t[fBinsTrig];
  fFadc = new Float_t[fBinsFadc];

  fread (fTrig, sizeof(Float_t), fBinsTrig, datei); 
  fread (fFadc, sizeof(Float_t), fBinsFadc, datei); 

  fclose ( datei ) ; 
}