PHOS/AliPHOSPulseGenerator.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007, ALICE Experiment at CERN, All rights reserved.      *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // The class which simulates the pulse shape from the PHOS FEE shaper,
  19 // make sampled amplitudes, digitize them.
  20 // Use case:
  21 //   AliPHOSPulseGenerator *pulse = new AliPHOSPulseGenerator(energy,time);
  22 //   Int_t *adcHG = new Int_t[pulse->GetRawFormatTimeBins()];
  23 //   Int_t *adcLG= new Int_t[pulse->GetRawFormatTimeBins()];
  24 //   pulse->AddNoise(1.);
  25 //   pulse->MakeSamples();
  26 //   pulse->GetSamples(adcHG, adcHG) ;
  27 //   pulse->Print();
  28 //   pulse->Draw();
  29 //
  30 // Author: Yuri Kharlov, after Yves Schutz and Per Thomas Hille
  31
  32 // --- ROOT system ---
  33 #include "TMath.h"
  34 #include "TF1.h"
  35 #include "TGraph.h"
  36 #include "TCanvas.h"
  37 #include "TRandom.h"
  38
  39 // --- AliRoot header files ---
  40 #include "AliLog.h"
  41 #include "AliPHOSPulseGenerator.h"
  42
  43 // --- Standard library ---
  44 #include <cmath>
  45 #include <iostream>
  46
  47 using std::cout;
  48 using std::endl;
  49
  50 ClassImp(AliPHOSPulseGenerator)
  51
  52 Double_t AliPHOSPulseGenerator::fgCapa        = 1.;        // 1pF
  53 Int_t    AliPHOSPulseGenerator::fgOrder       = 2 ;        // order of the Gamma function
  54 Double_t AliPHOSPulseGenerator::fgTimeMax     = 2.56E-5 ;  // each sample is over 100 ns fTimeMax/fTimeBins
  55 Double_t AliPHOSPulseGenerator::fgTimePeak    = 4.1E-6 ;   // 4 micro seconds
  56 Double_t AliPHOSPulseGenerator::fgTimeTrigger = 100E-9 ;   // 100ns, just for a reference
  57 Double_t AliPHOSPulseGenerator::fgHighCharge  = 8.2;       // adjusted for a high gain range of 5.12 GeV (10 bits)
  58 Double_t AliPHOSPulseGenerator::fgHighGain    = 6.64;
  59 Double_t AliPHOSPulseGenerator::fgHighLowGainFactor = 16.; // adjusted for a low gain range of 82 GeV (10 bits)
  60
  61 //-----------------------------------------------------------------------------
  62 AliPHOSPulseGenerator::AliPHOSPulseGenerator(Double_t a, Double_t t0)
  63   : TObject(), fAmplitude(a), fTZero(t0), fDataHG(0), fDataLG(0), fDigitize(kTRUE)
  64 {
  65   // Contruct a pulsegenrator object and initializes all necessary parameters
  66   // @param a digit amplitude in GeV
  67   // @param t0 time delay in nanoseconds of signal relative the first sample.
  68   // This value should be between 0 and Ts, where Ts is the sample interval
  69
  70   fDataHG = new Double_t[fkTimeBins];
  71   fDataLG = new Double_t[fkTimeBins];
  72   Reset();
  73 }
  74
  75 //-----------------------------------------------------------------------------
  76 AliPHOSPulseGenerator::AliPHOSPulseGenerator(const AliPHOSPulseGenerator & pulse)
  77   : TObject(), fAmplitude(pulse.fAmplitude), fTZero(pulse.fTZero), fDataHG(0), fDataLG(0), fDigitize(kTRUE)
  78 {
  79   fDataHG = new Double_t[pulse.fkTimeBins];
  80   fDataLG = new Double_t[pulse.fkTimeBins];
  81   for (Int_t i=0; i<pulse.fkTimeBins; i++) {
  82     fDataHG[i] = pulse.fDataHG[i];
  83     fDataLG[i] = pulse.fDataHG[i];
  84   }
  85 }
  86
  87 //-----------------------------------------------------------------------------
  88 AliPHOSPulseGenerator::~AliPHOSPulseGenerator()
  89 {
  90   // Destructor: delete arrays of samples
  91
  92   delete [] fDataHG;
  93   fDataHG=0;
  94   delete [] fDataLG;
  95   fDataLG=0;
  96 }
  97
  98 //-----------------------------------------------------------------------------
  99 void AliPHOSPulseGenerator::Reset()
 100 {
 101   // Reset all sample amplitudes to 0
 102
 103   for (Int_t i=0; i<fkTimeBins; i++) {
 104     fDataHG[i] = 0.;
 105     fDataLG[i] = 0.;
 106   }
 107 }
 108
 109 //-----------------------------------------------------------------------------
 110 void AliPHOSPulseGenerator::AddBaseline(Double_t baselineLevel)
 111 {
 112   // Adds a baseline offset to the signal
 113   // @param baselineLevel The basline level to add
 114   for (Int_t i=0; i<fkTimeBins; i++) {
 115     fDataHG[i] += baselineLevel;
 116     fDataLG[i] += baselineLevel;
 117   }
 118   // Digitize floating point amplitudes to integers
 119   if (fDigitize) Digitize();
 120 }
 121
 122 //-----------------------------------------------------------------------------
 123 void AliPHOSPulseGenerator::AddNoise(Double_t sigma)
 124 {
 125   // Adds Gaussian uncorrelated to the sample array
 126   // @param sigma the noise amplitude in entities of ADC levels
 127
 128   for (Int_t i=0; i<fkTimeBins; i++) {
 129     fDataHG[i] = gRandom->Gaus(0., sigma) ;
 130     fDataLG[i] = gRandom->Gaus(0., sigma) ;
 131   }
 132 }
 133
 134 //-----------------------------------------------------------------------------
 135 void AliPHOSPulseGenerator::AddNoise(Double_t * /* sigma */, Double_t /* cutoff */)
 136 {
 137   //Adds correlated Gaussian noise with cutof frequency "cutoff"
 138   // @param sigma noise amplitude in entities of ADC levels
 139   // @param -30DB cutoff frequency of the noise in entities of sampling frequency
 140
 141   AliError("not implemented yet");
 142 }
 143
 144 //-----------------------------------------------------------------------------
 145 void AliPHOSPulseGenerator::AddPretriggerSamples(Int_t nPresamples)
 146 {
 147   // Adds pretrigger samples to the sample arrays and replace them
 148   // with concatinated and truncated arrays
 149
 150   Double_t *tmpDataHG = new Double_t[fkTimeBins];
 151   Double_t *tmpDataLG = new Double_t[fkTimeBins];
 152   Int_t i;
 153   for (i=0; i<fkTimeBins; i++) {
 154     tmpDataHG[i] = fDataHG[i];
 155     tmpDataLG[i] = fDataLG[i];
 156   }
 157   for (i=0; i<fkTimeBins; i++) {
 158     if (i<nPresamples) {
 159       fDataHG[i] = 0.;
 160       fDataLG[i] = 0.;
 161     }
 162     else {
 163       fDataHG[i] = tmpDataHG[i-nPresamples];
 164       fDataLG[i] = tmpDataLG[i-nPresamples];
 165     }
 166   }
 167   delete [] tmpDataHG;
 168   delete [] tmpDataLG;
 169 }
 170
 171 //-----------------------------------------------------------------------------
 172 void AliPHOSPulseGenerator::Digitize()
 173 {
 174   // Emulates ADC: rounds down to nearest integer value all amplitudes
 175   for (Int_t i=0; i<fkTimeBins; i++) {
 176     fDataHG[i] = (Double_t) ((Int_t)(fDataHG[i] + 0.5));
 177     fDataLG[i] = (Double_t) ((Int_t)(fDataLG[i] + 0.5));
 178   }
 179 }
 180
 181 //-----------------------------------------------------------------------------
 182 Double_t AliPHOSPulseGenerator::RawResponseFunction(Double_t *x, Double_t *par)
 183 {
 184   // Shape of the electronics raw reponse:
 185   // It is a semi-gaussian, 2nd order Gamma function of the general form
 186   // v(t) = n**n * Q * A**n / C *(t/tp)**n * exp(-n * t/tp) with
 187   // tp : peaking time par[0]
 188   // n  : order of the function
 189   // C  : integrating capacitor in the preamplifier
 190   // A  : open loop gain of the preamplifier
 191   // Q  : the total APD charge to be measured Q = C * energy
 192
 193   Double_t signal ;
 194   Double_t xx = x[0] - ( fgTimeTrigger + par[3] ) ;
 195
 196   if (xx < 0 || xx > fgTimeMax)
 197     signal = 0. ;
 198   else {
 199     Double_t fac = par[0] * TMath::Power(fgOrder, fgOrder) * TMath::Power(par[1], fgOrder)/fgCapa ;
 200     signal = fac * par[2] * TMath::Power(xx/fgTimePeak, fgOrder) * TMath::Exp(-fgOrder*(xx/fgTimePeak)) ;
 201   }
 202   return signal ;
 203 }
 204
 205 //__________________________________________________________________
 206 Double_t AliPHOSPulseGenerator::RawResponseFunctionMax(Double_t charge, Double_t gain)
 207 {
 208   // Maximum value of the shaper response function
 209   return ( charge * TMath::Power(fgOrder, fgOrder) * TMath::Power(gain, fgOrder)
 210            / ( fgCapa * TMath::Exp(fgOrder) ) );
 211 }
 212
 213 //__________________________________________________________________
 214 Bool_t AliPHOSPulseGenerator::MakeSamples()
 215 {
 216   // for a start time fTZero and an amplitude fAmplitude given by digit,
 217   // calculates the raw sampled response AliPHOSPulseGenerator::RawResponseFunction
 218
 219   const Int_t kRawSignalOverflow = 0x3FF ; // decimal 1023
 220   Bool_t lowGain = kFALSE ;
 221
 222   TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4);
 223
 224   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
 225     signalF.SetParameter(0, fgHighCharge) ;
 226     signalF.SetParameter(1, fgHighGain) ;
 227     signalF.SetParameter(2, fAmplitude) ;
 228     signalF.SetParameter(3, fTZero) ;
 229     Double_t time = iTime * GetRawFormatTimeMax() / GetRawFormatTimeBins() ;
 230     Double_t signal = signalF.Eval(time) ;
 231     fDataHG[iTime] += signal;
 232     if ( static_cast<Int_t>(fDataHG[iTime]+0.5) > kRawSignalOverflow ){  // larger than 10 bits
 233       fDataHG[iTime] = kRawSignalOverflow ;
 234       lowGain = kTRUE ;
 235     }
 236
 237     signalF.SetParameter(0, GetRawFormatLowCharge() ) ;
 238     signalF.SetParameter(1, GetRawFormatLowGain() ) ;
 239     signal = signalF.Eval(time) ;
 240     fDataLG[iTime] += signal;
 241     if ( static_cast<Int_t>(fDataLG[iTime]+0.5) > kRawSignalOverflow)  // larger than 10 bits
 242       fDataLG[iTime] = kRawSignalOverflow ;
 243
 244   }
 245   // Digitize floating point amplitudes to integers
 246   if (fDigitize) Digitize();
 247   return lowGain ;
 248 }
 249
 250 //__________________________________________________________________
 251 void AliPHOSPulseGenerator::GetSamples(Int_t *adcHG, Int_t *adcLG) const
 252 {
 253   // Return integer sample arrays adcHG and adcLG
 254   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
 255     adcHG[iTime] = static_cast<Int_t>(fDataHG[iTime]) ;
 256     adcLG[iTime] = static_cast<Int_t>(fDataLG[iTime]) ;
 257   }
 258 }
 259
 260 //__________________________________________________________________
 261 void AliPHOSPulseGenerator::Print(Option_t*) const
 262 {
 263   // Prints sampled amplitudes to stdout
 264   Int_t i;
 265   cout << "High gain: ";
 266   for (i=0; i<fkTimeBins; i++)
 267     cout << (Int_t)fDataHG[i] << " ";
 268   cout << endl;
 269
 270   cout << "Low  gain: ";
 271   for (i=0; i<fkTimeBins; i++)
 272     cout << (Int_t)fDataLG[i] << " ";
 273   cout << endl;
 274 }
 275
 276 //__________________________________________________________________
 277 void AliPHOSPulseGenerator::Draw(Option_t*)
 278 {
 279   // Draw graphs with high and low gain samples
 280
 281   Double_t *time = new Double_t[fkTimeBins];
 282   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
 283     time[iTime] = iTime * GetRawFormatTimeMax() / GetRawFormatTimeBins() ;
 284   }
 285   Int_t nPoints = fkTimeBins;
 286   TGraph *graphHG = new TGraph(nPoints,time,fDataHG);
 287   TGraph *graphLG = new TGraph(nPoints,time,fDataLG);
 288   graphHG->SetMarkerStyle(20);
 289   graphLG->SetMarkerStyle(20);
 290   graphHG->SetMarkerSize(0.4);
 291   graphLG->SetMarkerSize(0.4);
 292   graphHG->SetTitle("High gain samples");
 293   graphLG->SetTitle("Low gain samples");
 294
 295   TCanvas *c1 = new TCanvas("c1","Raw ALTRO samples",10,10,700,500);
 296   c1->SetFillColor(0);
 297   c1->Divide(2,1);
 298   c1->cd(1);
 299   gPad->SetLeftMargin(0.15);
 300   graphHG->Draw("AP");
 301   graphHG->GetHistogram()->SetTitleOffset(1.6,"Y");
 302   graphHG->GetHistogram()->SetXTitle("time, #musec");
 303   graphHG->GetHistogram()->SetYTitle("Amplitude, ADC counts");
 304   c1->cd(2);
 305   gPad->SetLeftMargin(0.15);
 306   graphLG->Draw("AP");
 307   graphLG->GetHistogram()->SetTitleOffset(1.6,"Y");
 308   graphLG->GetHistogram()->SetXTitle("time, #musec");
 309   graphLG->GetHistogram()->SetYTitle("Amplitude, ADC counts");
 310   c1->Update();
 311 }