Add digitization of amplitudes, baseline and noise
[u/mrichter/AliRoot.git] / 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 }