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