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