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