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