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
  34 #include <TCanvas.h>
  35 #include <TF1.h>
  36 #include <TGraph.h>
  37 #include <TH1F.h>
  38 #include <TMath.h>
  39 #include <TRandom.h>
  40
  41 // --- AliRoot header files ---
  42 #include "AliLog.h"
  43 #include "AliPHOSPulseGenerator.h"
  44
  45 // --- Standard library ---
  46 #include <cmath>
  47 #include <iostream>
  48
  49 using std::cout;
  50 using std::endl;
  51
  52 ClassImp(AliPHOSPulseGenerator)
  53
  54 Double_t AliPHOSPulseGenerator::fgCapa        = 1.1;       // 1pF
  55 Int_t    AliPHOSPulseGenerator::fgOrder       = 2 ;        // order of the Gamma function
  56 Double_t AliPHOSPulseGenerator::fgTimePeak    = 2.1E-6 ;   // tau=2.1 micro seconds
  57 Double_t AliPHOSPulseGenerator::fgTimeTrigger = 100E-9 ;   // one tick 100 ns
  58 Double_t AliPHOSPulseGenerator::fgHighCharge  = 8.8;       // adjusted for a high gain range of 5 GeV (10 bits)
  59 Double_t AliPHOSPulseGenerator::fgHighGain    = 6.85;
  60 Double_t AliPHOSPulseGenerator::fgHighLowGainFactor = 16.; // adjusted for a low gain range of 80 GeV (10 bits)
  61
  62 //-----------------------------------------------------------------------------
  63 AliPHOSPulseGenerator::AliPHOSPulseGenerator(Double_t a, Double_t t0)
  64   : TObject(), fAmplitude(a), fTZero(t0), fDataHG(0), fDataLG(0), fDigitize(kTRUE)
  65 {
  66   // Contruct a pulsegenrator object and initializes all necessary parameters
  67   // @param a digit amplitude in GeV
  68   // @param t0 time delay in nanoseconds of signal relative the first sample.
  69   // This value should be between 0 and Ts, where Ts is the sample interval
  70
  71   fDataHG = new Double_t[fkTimeBins];
  72   fDataLG = new Double_t[fkTimeBins];
  73   Reset();
  74 }
  75
  76 //-----------------------------------------------------------------------------
  77 AliPHOSPulseGenerator::AliPHOSPulseGenerator(const AliPHOSPulseGenerator & pulse)
  78   : TObject(), fAmplitude(pulse.fAmplitude), fTZero(pulse.fTZero), fDataHG(0), fDataLG(0), fDigitize(kTRUE)
  79 {
  80   fDataHG = new Double_t[pulse.fkTimeBins];
  81   fDataLG = new Double_t[pulse.fkTimeBins];
  82   for (Int_t i=0; i<pulse.fkTimeBins; i++) {
  83     fDataHG[i] = pulse.fDataHG[i];
  84     fDataLG[i] = pulse.fDataHG[i];
  85   }
  86 }
  87
  88 //-----------------------------------------------------------------------------
  89 AliPHOSPulseGenerator::~AliPHOSPulseGenerator()
  90 {
  91   // Destructor: delete arrays of samples
  92
  93   delete [] fDataHG;
  94   fDataHG=0;
  95   delete [] fDataLG;
  96   fDataLG=0;
  97 }
  98
  99 //-----------------------------------------------------------------------------
 100 void AliPHOSPulseGenerator::Reset()
 101 {
 102   // Reset all sample amplitudes to 0
 103
 104   for (Int_t i=0; i<fkTimeBins; i++) {
 105     fDataHG[i] = 0.;
 106     fDataLG[i] = 0.;
 107   }
 108 }
 109
 110 //-----------------------------------------------------------------------------
 111 void AliPHOSPulseGenerator::AddBaseline(Double_t baselineLevel)
 112 {
 113   // Adds a baseline offset to the signal
 114   // @param baselineLevel The basline level to add
 115   for (Int_t i=0; i<fkTimeBins; i++) {
 116     fDataHG[i] += baselineLevel;
 117     fDataLG[i] += baselineLevel;
 118   }
 119   // Digitize floating point amplitudes to integers
 120   if (fDigitize) Digitize();
 121 }
 122
 123 //-----------------------------------------------------------------------------
 124 void AliPHOSPulseGenerator::AddNoise(Double_t sigma)
 125 {
 126   // Adds Gaussian uncorrelated to the sample array
 127   // @param sigma the noise amplitude in entities of ADC levels
 128
 129   for (Int_t i=0; i<fkTimeBins; i++) {
 130     fDataHG[i] = gRandom->Gaus(0., sigma) ;
 131     fDataLG[i] = gRandom->Gaus(0., sigma) ;
 132   }
 133 }
 134
 135 //-----------------------------------------------------------------------------
 136 void AliPHOSPulseGenerator::AddNoise(Double_t * /* sigma */, Double_t /* cutoff */)
 137 {
 138   //Adds correlated Gaussian noise with cutof frequency "cutoff"
 139   // @param sigma noise amplitude in entities of ADC levels
 140   // @param -30DB cutoff frequency of the noise in entities of sampling frequency
 141
 142   AliError("not implemented yet");
 143 }
 144
 145 //-----------------------------------------------------------------------------
 146 void AliPHOSPulseGenerator::AddPretriggerSamples(Int_t nPresamples)
 147 {
 148   // Adds pretrigger samples to the sample arrays and replace them
 149   // with concatinated and truncated arrays
 150
 151   Double_t *tmpDataHG = new Double_t[fkTimeBins];
 152   Double_t *tmpDataLG = new Double_t[fkTimeBins];
 153   Int_t i;
 154   for (i=0; i<fkTimeBins; i++) {
 155     tmpDataHG[i] = fDataHG[i];
 156     tmpDataLG[i] = fDataLG[i];
 157   }
 158   for (i=0; i<fkTimeBins; i++) {
 159     if (i<nPresamples) {
 160       fDataHG[i] = 0.;
 161       fDataLG[i] = 0.;
 162     }
 163     else {
 164       fDataHG[i] = tmpDataHG[i-nPresamples];
 165       fDataLG[i] = tmpDataLG[i-nPresamples];
 166     }
 167   }
 168   delete [] tmpDataHG;
 169   delete [] tmpDataLG;
 170 }
 171
 172 //-----------------------------------------------------------------------------
 173 void AliPHOSPulseGenerator::Digitize()
 174 {
 175   // Emulates ADC: rounds down to nearest integer value all amplitudes
 176   for (Int_t i=0; i<fkTimeBins; i++) {
 177     fDataHG[i] = (Double_t) ((Int_t)fDataHG[i]);
 178     fDataLG[i] = (Double_t) ((Int_t)fDataLG[i]);
 179   }
 180 }
 181
 182 //-----------------------------------------------------------------------------
 183 Double_t AliPHOSPulseGenerator::RawResponseFunction(Double_t *x, Double_t *par)
 184 {
 185   // Shape of the electronics raw reponse:
 186   // It is a semi-gaussian, 2nd order Gamma function of the general form
 187   // v(t) = n**n * Q * A**n / C *(t/tp)**n * exp(-n * t/tp) with
 188   // tp : peaking time par[0]
 189   // n  : order of the function
 190   // C  : integrating capacitor in the preamplifier
 191   // A  : open loop gain of the preamplifier
 192   // Q  : the total APD charge to be measured Q = C * energy
 193
 194   Double_t signal ;
 195   Double_t xx = x[0] - ( fgTimeTrigger + par[3] ) ;
 196
 197   if (xx < 0 || xx > GetRawFormatTimeMax())
 198     signal = 0. ;
 199   else {
 200     Double_t fac = par[0] * TMath::Power(fgOrder, fgOrder) * TMath::Power(par[1], fgOrder)/fgCapa ;
 201     signal = fac * par[2] * TMath::Power(xx/fgTimePeak, fgOrder) * TMath::Exp(-fgOrder*(xx/fgTimePeak)) ;
 202   }
 203   return signal ;
 204 }
 205
 206 //__________________________________________________________________
 207 Double_t AliPHOSPulseGenerator::RawResponseFunctionMax(Double_t charge, Double_t gain)
 208 {
 209   // Maximum value of the shaper response function
 210   return ( charge * TMath::Power(fgOrder, fgOrder) * TMath::Power(gain, fgOrder)
 211            / ( fgCapa * TMath::Exp(fgOrder) ) );
 212 }
 213
 214 //__________________________________________________________________
 215 Bool_t AliPHOSPulseGenerator::MakeSamples()
 216 {
 217   // for a start time fTZero and an amplitude fAmplitude given by digit,
 218   // calculates the raw sampled response AliPHOSPulseGenerator::RawResponseFunction
 219
 220   const Int_t kRawSignalOverflow = 0x3FF ; // decimal 1023
 221   Bool_t lowGain = kFALSE ;
 222
 223   TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4);
 224
 225   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
 226     signalF.SetParameter(0, fgHighCharge) ;
 227     signalF.SetParameter(1, fgHighGain) ;
 228     signalF.SetParameter(2, fAmplitude) ;
 229     signalF.SetParameter(3, fTZero) ;
 230     Double_t time = iTime * GetRawFormatTimeMax() / GetRawFormatTimeBins() ;
 231     Double_t signal = signalF.Eval(time) ;
 232     fDataHG[iTime] += signal;
 233     if ( static_cast<Int_t>(fDataHG[iTime]+0.5) > kRawSignalOverflow ){  // larger than 10 bits
 234       fDataHG[iTime] = kRawSignalOverflow ;
 235       lowGain = kTRUE ;
 236     }
 237
 238     signalF.SetParameter(0, GetRawFormatLowCharge() ) ;
 239     signalF.SetParameter(1, GetRawFormatLowGain() ) ;
 240     signal = signalF.Eval(time) ;
 241     fDataLG[iTime] += signal;
 242     if ( static_cast<Int_t>(fDataLG[iTime]+0.5) > kRawSignalOverflow)  // larger than 10 bits
 243       fDataLG[iTime] = kRawSignalOverflow ;
 244
 245   }
 246   // Digitize floating point amplitudes to integers
 247   if (fDigitize) Digitize();
 248   return lowGain ;
 249 }
 250
 251 //__________________________________________________________________
 252 void AliPHOSPulseGenerator::GetSamples(Int_t *adcHG, Int_t *adcLG) const
 253 {
 254   // Return integer sample arrays adcHG and adcLG
 255   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
 256     adcHG[iTime] = static_cast<Int_t>(fDataHG[iTime]) ;
 257     adcLG[iTime] = static_cast<Int_t>(fDataLG[iTime]) ;
 258   }
 259 }
 260
 261 //__________________________________________________________________
 262 void AliPHOSPulseGenerator::Print(Option_t*) const
 263 {
 264   // Prints sampled amplitudes to stdout
 265   Int_t i;
 266   cout << "High gain: ";
 267   for (i=0; i<fkTimeBins; i++)
 268     cout << (Int_t)fDataHG[i] << " ";
 269   cout << endl;
 270
 271   cout << "Low  gain: ";
 272   for (i=0; i<fkTimeBins; i++)
 273     cout << (Int_t)fDataLG[i] << " ";
 274   cout << endl;
 275 }
 276
 277 //__________________________________________________________________
 278 void AliPHOSPulseGenerator::Draw(Option_t* opt)
 279 {
 280   // Draw graphs with high and low gain samples
 281   // Option_t* opt="all": draw both HG and LG in one canvas
 282   //               "HG" : draw HG only
 283   //               "LG" : draw LG only
 284
 285   Double_t *time = new Double_t[fkTimeBins];
 286   for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) {
 287     time[iTime] = iTime * GetRawFormatTimeMax() / GetRawFormatTimeBins() ;
 288   }
 289   Int_t nPoints = fkTimeBins;
 290   TGraph *graphHG = new TGraph(nPoints,time,fDataHG);
 291   TGraph *graphLG = new TGraph(nPoints,time,fDataLG);
 292   graphHG->SetMarkerStyle(20);
 293   graphLG->SetMarkerStyle(20);
 294   graphHG->SetMarkerSize(0.4);
 295   graphLG->SetMarkerSize(0.4);
 296   graphHG->SetTitle("High gain samples");
 297   graphLG->SetTitle("Low gain samples");
 298
 299   TCanvas *c1 = new TCanvas("c1","Raw ALTRO samples",10,10,700,500);
 300   c1->SetFillColor(0);
 301
 302   if (strstr(opt,"all")){
 303     c1->Divide(2,1);
 304     c1->cd(1);
 305     gPad->SetLeftMargin(0.15);
 306   }
 307   if (strstr(opt,"LG") == 0){
 308     graphHG->Draw("AP");
 309     graphHG->GetHistogram()->SetTitleOffset(1.6,"Y");
 310     graphHG->GetHistogram()->SetXTitle("time, #musec");
 311     graphHG->GetHistogram()->SetYTitle("Amplitude, ADC counts");
 312   }
 313   if (strstr(opt,"all")){
 314     c1->cd(2);
 315     gPad->SetLeftMargin(0.15);
 316   }
 317   if (strstr(opt,"HG") == 0){
 318     graphLG->Draw("AP");
 319     graphLG->GetHistogram()->SetTitleOffset(1.6,"Y");
 320     graphLG->GetHistogram()->SetXTitle("time, #musec");
 321     graphLG->GetHistogram()->SetYTitle("Amplitude, ADC counts");
 322   }
 323   c1->Update();
 324 }
 325
 326 //__________________________________________________________________
 327 Double_t AliPHOSPulseGenerator::GeV2ADC()
 328 {
 329   //Return GeV to ADC counts conversion factor.
 330   //adc_counts = energy[GeV]*AliPHOSPulseGenerator::GeV2ADC().
 331
 332   return RawResponseFunctionMax(fgHighCharge,fgHighGain);
 333 }