Fix Coverity defects
[u/mrichter/AliRoot.git] / TRD / AliTRDSimParam.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 // Class containing constant simulation parameters                        //
21 //                                                                        //
22 // Request an instance with AliTRDSimParam::Instance()                    //
23 // Then request the needed values                                         //
24 //                                                                        //
25 ////////////////////////////////////////////////////////////////////////////
26
27 #include <TMath.h>
28
29 #include "AliRun.h"
30
31 #include "AliTRDSimParam.h"
32 #include "AliTRDCommonParam.h"
33
34 ClassImp(AliTRDSimParam)
35
36 AliTRDSimParam *AliTRDSimParam::fgInstance   = 0;
37 Bool_t          AliTRDSimParam::fgTerminated = kFALSE;
38
39 //_ singleton implementation __________________________________________________
40 AliTRDSimParam* AliTRDSimParam::Instance()
41 {
42   //
43   // Singleton implementation
44   // Returns an instance of this class, it is created if neccessary
45   // 
46   
47   if (fgTerminated != kFALSE) {
48     return 0;
49   }
50
51   if (fgInstance == 0) {
52     fgInstance = new AliTRDSimParam();
53   }  
54
55   return fgInstance;
56
57 }
58
59 //_ singleton implementation __________________________________________________
60 void AliTRDSimParam::Terminate()
61 {
62   //
63   // Singleton implementation
64   // Deletes the instance of this class and sets the terminated flag,
65   // instances cannot be requested anymore
66   // This function can be called several times.
67   //
68   
69   fgTerminated = kTRUE;
70   
71   if (fgInstance != 0) {
72     delete fgInstance;
73     fgInstance = 0;
74   }
75
76 }
77
78 //_____________________________________________________________________________
79 AliTRDSimParam::AliTRDSimParam()
80   :TObject()
81   ,fGasGain(0.0)
82   ,fNoise(0.0)
83   ,fChipGain(0.0)
84   ,fADCoutRange(0.0)
85   ,fADCinRange(0.0)
86   ,fADCbaseline(0)
87   ,fDiffusionOn(kFALSE)
88   ,fElAttachOn(kFALSE)
89   ,fElAttachProp(0.0)
90   ,fTRFOn(kFALSE)
91   ,fTRFsmp(0)
92   ,fTRFbin(0)
93   ,fTRFlo(0.0)
94   ,fTRFhi(0.0)
95   ,fTRFwid(0.0)
96   ,fCTOn(kFALSE)
97   ,fCTsmp(0)
98   ,fPadCoupling(0.0)
99   ,fTimeCoupling(0.0)
100   ,fTimeStructOn(kFALSE)
101   ,fPRFOn(kFALSE)
102   ,fNTimeBins(0)
103 {
104   //
105   // Default constructor
106   //
107   
108   Init();
109
110 }
111
112 //_____________________________________________________________________________
113 void AliTRDSimParam::Init()
114 {
115   // 
116   // Default initializiation
117   //
118   
119   // The default parameter for the digitization
120   fGasGain           = 4000.0;
121   fChipGain          =   12.4;
122   fNoise             = 1250.0;
123   fADCoutRange       = 1023.0;          // 10-bit ADC
124   fADCinRange        = 2000.0;          // 2V input range
125   fADCbaseline       =   10;
126
127   // Diffusion on
128   fDiffusionOn       = kTRUE;
129   
130   // Propability for electron attachment
131   fElAttachOn        = kFALSE;
132   fElAttachProp      = 0.0;
133
134   // The time response function
135   fTRFOn             = kTRUE;
136
137   // The cross talk
138   fCTOn              = kTRUE;
139
140   // The pad coupling factor
141   // Use 0.46, instead of the theroetical value 0.3, since it reproduces better 
142   // the test beam data, even tough it is not understood why.
143   fPadCoupling       = 0.46;
144
145   // The time coupling factor (same number as for the TPC)
146   fTimeCoupling      = 0.4;
147
148   // Use drift time maps
149   fTimeStructOn      = kTRUE;
150   
151   // The pad response function
152   fPRFOn             = kTRUE;
153
154   // The number of time bins
155   fNTimeBins         = 30;
156
157   ReInit();
158
159 }
160
161 //_____________________________________________________________________________
162 AliTRDSimParam::~AliTRDSimParam() 
163 {
164   //
165   // Destructor
166   //
167   
168   if (fTRFsmp) {
169     delete [] fTRFsmp;
170     fTRFsmp = 0;
171   }
172
173   if (fCTsmp) {
174     delete [] fCTsmp;
175     fCTsmp  = 0;
176   }
177
178 }
179
180 //_____________________________________________________________________________
181 AliTRDSimParam::AliTRDSimParam(const AliTRDSimParam &p)
182   :TObject(p)
183   ,fGasGain(p.fGasGain)
184   ,fNoise(p.fNoise)
185   ,fChipGain(p.fChipGain)
186   ,fADCoutRange(p.fADCoutRange)
187   ,fADCinRange(p.fADCinRange)
188   ,fADCbaseline(p.fADCbaseline)
189   ,fDiffusionOn(p.fDiffusionOn)
190   ,fElAttachOn(p.fElAttachOn)
191   ,fElAttachProp(p.fElAttachProp)
192   ,fTRFOn(p.fTRFOn)
193   ,fTRFsmp(0)
194   ,fTRFbin(p.fTRFbin)
195   ,fTRFlo(p.fTRFlo)
196   ,fTRFhi(p.fTRFhi)
197   ,fTRFwid(p.fTRFwid)
198   ,fCTOn(p.fCTOn)
199   ,fCTsmp(0)
200   ,fPadCoupling(p.fPadCoupling)
201   ,fTimeCoupling(p.fTimeCoupling)
202   ,fTimeStructOn(p.fTimeStructOn)
203   ,fPRFOn(p.fPRFOn)
204   ,fNTimeBins(p.fNTimeBins)
205 {
206   //
207   // Copy constructor
208   //
209
210   Int_t iBin = 0;
211
212   fTRFsmp = new Float_t[fTRFbin];
213   for (iBin = 0; iBin < fTRFbin; iBin++) {
214     fTRFsmp[iBin] = ((AliTRDSimParam &) p).fTRFsmp[iBin];
215   }                                                                             
216
217   fCTsmp  = new Float_t[fTRFbin];
218   for (iBin = 0; iBin < fTRFbin; iBin++) {
219     fCTsmp[iBin]  = ((AliTRDSimParam &) p).fCTsmp[iBin];
220   }                                                                             
221
222 }
223
224 //_____________________________________________________________________________
225 AliTRDSimParam &AliTRDSimParam::operator=(const AliTRDSimParam &p)
226 {
227   //
228   // Assignment operator
229   //
230
231   if (this != &p) {
232     ((AliTRDSimParam &) p).Copy(*this);
233   }
234
235   return *this;
236
237 }
238
239 //_____________________________________________________________________________
240 void AliTRDSimParam::Copy(TObject &p) const
241 {
242   //
243   // Copy function
244   //
245   
246   AliTRDSimParam *target = dynamic_cast<AliTRDSimParam *> (&p);
247   if (!target) {
248     return;
249   }
250
251   target->fGasGain            = fGasGain;
252   target->fNoise              = fNoise;
253   target->fChipGain           = fChipGain;  
254   target->fADCoutRange        = fADCoutRange;
255   target->fADCinRange         = fADCinRange;
256   target->fADCbaseline        = fADCbaseline; 
257   target->fDiffusionOn        = fDiffusionOn; 
258   target->fElAttachOn         = fElAttachOn;
259   target->fElAttachProp       = fElAttachProp;
260   target->fTRFOn              = fTRFOn;
261   target->fTRFbin             = fTRFbin;
262   target->fTRFlo              = fTRFlo;
263   target->fTRFhi              = fTRFhi;
264   target->fTRFwid             = fTRFwid;
265   target->fCTOn               = fCTOn;
266   target->fPadCoupling        = fPadCoupling;
267   target->fTimeCoupling       = fTimeCoupling;
268   target->fPRFOn              = fPRFOn;
269   target->fNTimeBins          = fNTimeBins;
270
271   if (target->fTRFsmp) {
272     delete[] target->fTRFsmp;
273   }
274   target->fTRFsmp = new Float_t[fTRFbin];
275   for (Int_t iBin = 0; iBin < fTRFbin; iBin++) {
276     target->fTRFsmp[iBin] = fTRFsmp[iBin];
277   }
278
279   if (target->fCTsmp) {
280     delete[] target->fCTsmp;
281   }
282   target->fCTsmp  = new Float_t[fTRFbin];
283   for (Int_t iBin = 0; iBin < fTRFbin; iBin++) {
284     target->fCTsmp[iBin]  = fCTsmp[iBin];
285   }
286   
287 }
288
289 //_____________________________________________________________________________
290 void AliTRDSimParam::ReInit()
291 {
292   //
293   // Reinitializes the parameter class after a change
294   //
295
296   if      (AliTRDCommonParam::Instance()->IsXenon()) {
297     // The range and the binwidth for the sampled TRF 
298     fTRFbin = 200;
299     // Start 0.2 mus before the signal
300     fTRFlo  = -0.4;
301     // End the maximum drift time after the signal 
302     fTRFhi  =  3.58;
303     // Standard gas gain
304     fGasGain = 4000.0;
305   }
306   else if (AliTRDCommonParam::Instance()->IsArgon()) {
307     // The range and the binwidth for the sampled TRF 
308     fTRFbin  =  50;
309     // Start 0.2 mus before the signal
310     fTRFlo   =  0.02;
311     // End the maximum drift time after the signal 
312     fTRFhi   =  1.98;
313     // Higher gas gain
314     fGasGain = 8000.0;
315   }
316   else {
317     AliFatal("Not a valid gas mixture!");
318     exit(1);
319   }
320   fTRFwid = (fTRFhi - fTRFlo) / ((Float_t) fTRFbin);
321
322   // Create the sampled TRF
323   SampleTRF();
324
325 }
326
327 //_____________________________________________________________________________
328 void AliTRDSimParam::SampleTRF()
329 {
330   //
331   // Samples the new time response function.
332   //
333
334   Int_t ipasa = 0;
335
336   // Xenon
337   // From Antons measurements with Fe55 source, adjusted by C. Lippmann.
338   // time bins are -0.4, -0.38, -0.36, ...., 3.54, 3.56, 3.58 microseconds
339   const Int_t kNpasa     = 200;  // kNpasa should be equal to fTRFbin!
340   Float_t xtalk[kNpasa];
341   Float_t signal[kNpasa]     = { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000
342                                , 0.0002, 0.0007, 0.0026, 0.0089, 0.0253, 0.0612, 0.1319
343                                , 0.2416, 0.3913, 0.5609, 0.7295, 0.8662, 0.9581, 1.0000
344                                , 0.9990, 0.9611, 0.8995, 0.8269, 0.7495, 0.6714, 0.5987
345                                , 0.5334, 0.4756, 0.4249, 0.3811, 0.3433, 0.3110, 0.2837
346                                , 0.2607, 0.2409, 0.2243, 0.2099, 0.1974, 0.1868, 0.1776
347                                , 0.1695, 0.1627, 0.1566, 0.1509, 0.1457, 0.1407, 0.1362
348                                , 0.1317, 0.1274, 0.1233, 0.1196, 0.1162, 0.1131, 0.1102
349                                , 0.1075, 0.1051, 0.1026, 0.1004, 0.0979, 0.0956, 0.0934
350                                , 0.0912, 0.0892, 0.0875, 0.0858, 0.0843, 0.0829, 0.0815
351                                , 0.0799, 0.0786, 0.0772, 0.0757, 0.0741, 0.0729, 0.0718
352                                , 0.0706, 0.0692, 0.0680, 0.0669, 0.0655, 0.0643, 0.0630
353                                , 0.0618, 0.0607, 0.0596, 0.0587, 0.0576, 0.0568, 0.0558
354                                , 0.0550, 0.0541, 0.0531, 0.0522, 0.0513, 0.0505, 0.0497
355                                , 0.0490, 0.0484, 0.0474, 0.0465, 0.0457, 0.0449, 0.0441
356                                , 0.0433, 0.0425, 0.0417, 0.0410, 0.0402, 0.0395, 0.0388
357                                , 0.0381, 0.0374, 0.0368, 0.0361, 0.0354, 0.0348, 0.0342
358                                , 0.0336, 0.0330, 0.0324, 0.0318, 0.0312, 0.0306, 0.0301
359                                , 0.0296, 0.0290, 0.0285, 0.0280, 0.0275, 0.0270, 0.0265
360                                , 0.0260, 0.0256, 0.0251, 0.0246, 0.0242, 0.0238, 0.0233
361                                , 0.0229, 0.0225, 0.0221, 0.0217, 0.0213, 0.0209, 0.0206
362                                , 0.0202, 0.0198, 0.0195, 0.0191, 0.0188, 0.0184, 0.0181
363                                , 0.0178, 0.0175, 0.0171, 0.0168, 0.0165, 0.0162, 0.0159
364                                , 0.0157, 0.0154, 0.0151, 0.0148, 0.0146, 0.0143, 0.0140
365                                , 0.0138, 0.0135, 0.0133, 0.0131, 0.0128, 0.0126, 0.0124
366                                , 0.0121, 0.0119, 0.0120, 0.0115, 0.0113, 0.0111, 0.0109
367                                , 0.0107, 0.0105, 0.0103, 0.0101, 0.0100, 0.0098, 0.0096
368                                , 0.0094, 0.0092, 0.0091, 0.0089, 0.0088, 0.0086, 0.0084
369                                , 0.0083, 0.0081, 0.0080, 0.0078 };
370   signal[0] = 0.0;
371   signal[1] = 0.0;
372   signal[2] = 0.0;
373   // With undershoot, positive peak corresponds to ~3% of the main signal:
374   for (ipasa = 3; ipasa < kNpasa; ipasa++) {
375     xtalk[ipasa] = 0.2 * (signal[ipasa-2] - signal[ipasa-3]);
376   }
377   xtalk[0]  = 0.0;   
378   xtalk[1]  = 0.0;  
379   xtalk[2]  = 0.0;  
380
381   // Argon
382   // Ar measurement with Fe55 source by Anton
383   // time bins are 0.02, 0.06, 0.10, ...., 1.90, 1.94, 1.98 microseconds
384   const Int_t kNpasaAr = 50;
385   Float_t xtalkAr[kNpasaAr];
386   Float_t signalAr[kNpasaAr] = { -0.01,  0.01,  0.00,  0.00,  0.01
387                                , -0.01,  0.01,  2.15, 22.28, 55.53
388                                , 68.52, 58.21, 40.92, 27.12, 18.49
389                                , 13.42, 10.48,  8.67,  7.49,  6.55
390                                ,  5.71,  5.12,  4.63,  4.22,  3.81
391                                ,  3.48,  3.20,  2.94,  2.77,  2.63
392                                ,  2.50,  2.37,  2.23,  2.13,  2.03
393                                ,  1.91,  1.83,  1.75,  1.68,  1.63
394                                ,  1.56,  1.49,  1.50,  1.49,  1.29
395                                ,  1.19,  1.21,  1.21,  1.20,  1.10 };
396   // Normalization to maximum
397   for (ipasa = 0; ipasa < kNpasaAr; ipasa++) {
398     signalAr[ipasa] /= 68.52;
399   }
400   signalAr[0] = 0.0;
401   signalAr[1] = 0.0;
402   signalAr[2] = 0.0;
403   // With undershoot, positive peak corresponds to ~3% of the main signal:
404   for (ipasa = 3; ipasa < kNpasaAr; ipasa++) {
405     xtalkAr[ipasa] = 0.2 * (signalAr[ipasa-2] - signalAr[ipasa-3]);
406   }
407   xtalkAr[0]  = 0.0;   
408   xtalkAr[1]  = 0.0;  
409   xtalkAr[2]  = 0.0;  
410
411   if (fTRFsmp) {
412     delete [] fTRFsmp;
413   }
414   fTRFsmp = new Float_t[fTRFbin];
415
416   if (fCTsmp)  {
417     delete [] fCTsmp;
418   }
419   fCTsmp  = new Float_t[fTRFbin];
420
421   if      (AliTRDCommonParam::Instance()->IsXenon()) {
422     if (fTRFbin != kNpasa) {
423       AliError("Array mismatch (xenon)\n\n");
424     }
425   }
426   else if (AliTRDCommonParam::Instance()->IsArgon()) {
427     if (fTRFbin != kNpasaAr) {
428       AliError("Array mismatch (argon)\n\n");
429     }
430   }
431
432   for (Int_t iBin = 0; iBin < fTRFbin; iBin++) {
433     if      (AliTRDCommonParam::Instance()->IsXenon()) {
434       fTRFsmp[iBin] = signal[iBin];
435       fCTsmp[iBin]  = xtalk[iBin];
436     }
437     else if (AliTRDCommonParam::Instance()->IsArgon()) {
438       fTRFsmp[iBin] = signalAr[iBin];
439       fCTsmp[iBin]  = xtalkAr[iBin];
440     }
441   }
442
443 }
444
445 //_____________________________________________________________________________
446 Double_t AliTRDSimParam::TimeResponse(Double_t time) const
447 {
448   //
449   // Applies the preamp shaper time response
450   // (We assume a signal rise time of 0.2us = fTRFlo/2.
451   //
452
453   Double_t rt   = (time - .5*fTRFlo) / fTRFwid;
454   Int_t    iBin = (Int_t) rt; 
455   Double_t dt   = rt - iBin; 
456   if ((iBin >= 0) && (iBin+1 < fTRFbin)) {
457     return fTRFsmp[iBin] + (fTRFsmp[iBin+1] - fTRFsmp[iBin])*dt;
458   } 
459   else {
460     return 0.0;
461   }
462
463 }
464
465 //_____________________________________________________________________________
466 Double_t AliTRDSimParam::CrossTalk(Double_t time) const
467 {
468   //
469   // Applies the pad-pad capacitive cross talk
470   //
471
472   Double_t rt   = (time - fTRFlo) / fTRFwid;
473   Int_t    iBin = (Int_t) rt; 
474   Double_t dt   = rt - iBin; 
475   if ((iBin >= 0) && (iBin+1 < fTRFbin)) {
476     return fCTsmp[iBin] + (fCTsmp[iBin+1] - fCTsmp[iBin])*dt;
477   } 
478   else {
479     return 0.0;
480   }
481
482 }