38ee115a6c8b0190667b0fcba85f97855a3c6280
[u/mrichter/AliRoot.git] / TRD / AliTRDv1.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 /*
17 $Log$
18 Revision 1.17.2.1  2000/05/08 14:59:16  cblume
19 Made inline function non-virtual. Bug fix in setting sensitive chamber
20
21 Revision 1.17  2000/02/28 19:10:26  cblume
22 Include the new TRD classes
23
24 Revision 1.16.4.1  2000/02/28 18:04:35  cblume
25 Change to new hit version, introduce geometry class, and move digitization and clustering to AliTRDdigitizer/AliTRDclusterizerV1
26
27 Revision 1.16  1999/11/05 22:50:28  fca
28 Do not use Atan, removed from ROOT too
29
30 Revision 1.15  1999/11/02 17:20:19  fca
31 initialise nbytes before using it
32
33 Revision 1.14  1999/11/02 17:15:54  fca
34 Correct ansi scoping not accepted by HP compilers
35
36 Revision 1.13  1999/11/02 17:14:51  fca
37 Correct ansi scoping not accepted by HP compilers
38
39 Revision 1.12  1999/11/02 16:35:56  fca
40 New version of TRD introduced
41
42 Revision 1.11  1999/11/01 20:41:51  fca
43 Added protections against using the wrong version of FRAME
44
45 Revision 1.10  1999/09/29 09:24:35  fca
46 Introduction of the Copyright and cvs Log
47
48 */
49
50 ///////////////////////////////////////////////////////////////////////////////
51 //                                                                           //
52 //  Transition Radiation Detector version 2 -- slow simulator                //
53 //                                                                           //
54 //Begin_Html
55 /*
56 <img src="picts/AliTRDfullClass.gif">
57 */
58 //End_Html
59 //                                                                           //
60 //                                                                           //
61 ///////////////////////////////////////////////////////////////////////////////
62
63 #include <TMath.h>
64 #include <TVector.h>
65 #include <TRandom.h>
66
67 #include "AliRun.h"
68 #include "AliMC.h"
69 #include "AliConst.h"
70
71 #include "AliTRDv1.h"
72 #include "AliTRDmatrix.h"
73 #include "AliTRDgeometry.h"
74
75 ClassImp(AliTRDv1)
76
77 //_____________________________________________________________________________
78 AliTRDv1::AliTRDv1(const char *name, const char *title) 
79          :AliTRD(name, title) 
80 {
81   //
82   // Standard constructor for Transition Radiation Detector version 1
83   //
84
85   fIdSens        =  0;
86
87   fIdChamber1    =  0;
88   fIdChamber2    =  0;
89   fIdChamber3    =  0;
90
91   fSensSelect    =  0;
92   fSensPlane     = -1;
93   fSensChamber   = -1;
94   fSensSector    = -1;
95
96   fDeltaE        = NULL;
97
98   SetBufferSize(128000);
99
100 }
101
102 //_____________________________________________________________________________
103 AliTRDv1::~AliTRDv1()
104 {
105
106   if (fDeltaE) delete fDeltaE;
107
108 }
109  
110 //_____________________________________________________________________________
111 void AliTRDv1::CreateGeometry()
112 {
113   //
114   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
115   // This version covers the full azimuth. 
116   //
117
118   // Check that FRAME is there otherwise we have no place where to put the TRD
119   AliModule* FRAME = gAlice->GetModule("FRAME");
120   if (!FRAME) return;
121
122   // Define the chambers
123   AliTRD::CreateGeometry();
124
125 }
126
127 //_____________________________________________________________________________
128 void AliTRDv1::CreateMaterials()
129 {
130   //
131   // Create materials for the Transition Radiation Detector version 1
132   //
133
134   AliTRD::CreateMaterials();
135
136 }
137
138 //_____________________________________________________________________________
139 void AliTRDv1::Init() 
140 {
141   //
142   // Initialise Transition Radiation Detector after geometry has been built.
143   //
144
145   AliTRD::Init();
146
147   printf("          Slow simulator\n\n");
148   if (fSensSelect) {
149     if (fSensPlane   >= 0)
150       printf("          Only plane %d is sensitive\n",fSensPlane);
151     if (fSensChamber >= 0)   
152       printf("          Only chamber %d is sensitive\n",fSensChamber);
153     if (fSensSector  >= 0)
154       printf("          Only sector %d is sensitive\n",fSensSector);
155   }
156   printf("\n");
157
158   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
159   const Float_t kPoti = 12.1;
160   // Maximum energy (50 keV);
161   const Float_t kEend = 50000.0;
162   // Ermilova distribution for the delta-ray spectrum
163   Float_t Poti = TMath::Log(kPoti);
164   Float_t Eend = TMath::Log(kEend);
165   fDeltaE  = new TF1("deltae",Ermilova,Poti,Eend,0);
166
167   // Identifier of the sensitive volume (drift region)
168   fIdSens     = gMC->VolId("UL05");
169
170   // Identifier of the TRD-driftchambers
171   fIdChamber1 = gMC->VolId("UCIO");
172   fIdChamber2 = gMC->VolId("UCIM");
173   fIdChamber3 = gMC->VolId("UCII");
174
175   for (Int_t i = 0; i < 80; i++) printf("*");
176   printf("\n");
177
178 }
179
180 //_____________________________________________________________________________
181 void AliTRDv1::SetSensPlane(Int_t iplane)
182 {
183   //
184   // Defines the hit-sensitive plane (0-5)
185   //
186
187   if ((iplane < 0) || (iplane > 5)) {
188     printf("Wrong input value: %d\n",iplane);
189     printf("Use standard setting\n");
190     fSensPlane  = -1;
191     fSensSelect =  0;
192     return;
193   }
194
195   fSensSelect = 1;
196   fSensPlane  = iplane;
197
198 }
199
200 //_____________________________________________________________________________
201 void AliTRDv1::SetSensChamber(Int_t ichamber)
202 {
203   //
204   // Defines the hit-sensitive chamber (0-4)
205   //
206
207   if ((ichamber < 0) || (ichamber > 4)) {
208     printf("Wrong input value: %d\n",ichamber);
209     printf("Use standard setting\n");
210     fSensChamber = -1;
211     fSensSelect  =  0;
212     return;
213   }
214
215   fSensSelect  = 1;
216   fSensChamber = ichamber;
217
218 }
219
220 //_____________________________________________________________________________
221 void AliTRDv1::SetSensSector(Int_t isector)
222 {
223   //
224   // Defines the hit-sensitive sector (0-17)
225   //
226
227   if ((isector < 0) || (isector > 17)) {
228     printf("Wrong input value: %d\n",isector);
229     printf("Use standard setting\n");
230     fSensSector = -1;
231     fSensSelect =  0;
232     return;
233   }
234
235   fSensSelect = 1;
236   fSensSector = isector;
237
238 }
239
240 //_____________________________________________________________________________
241 void AliTRDv1::StepManager()
242 {
243   //
244   // Slow simulator. Every charged track produces electron cluster as hits 
245   // along its path across the drift volume. The step size is set acording
246   // to Bethe-Bloch. The energy distribution of the delta electrons follows
247   // a spectrum taken from Ermilova et al.
248   //
249
250   Int_t    iIdSens, icSens;
251   Int_t    iIdSpace, icSpace;
252   Int_t    iIdChamber, icChamber;
253   Int_t    pla = 0;
254   Int_t    cha = 0;
255   Int_t    sec = 0;
256   Int_t    iPdg;
257
258   Float_t  hits[4];
259   Float_t  random[1];
260   Float_t  charge;
261   Float_t  aMass;
262
263   Double_t pTot;
264   Double_t qTot;
265   Double_t eDelta;
266   Double_t betaGamma, pp;
267
268   TLorentzVector pos, mom;
269   TClonesArray  &lhits = *fHits;
270
271   const Double_t kBig     = 1.0E+12;
272
273   // Ionization energy
274   const Float_t  kWion    = 22.04;
275   // Maximum energy for e+ e- g for the step-size calculation
276   const Float_t  kPTotMax = 0.002;
277   // Plateau value of the energy-loss for electron in xenon
278   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
279   //const Double_t kPlateau = 1.70;
280   // the averaged value (26/3/99)
281   const Float_t  kPlateau = 1.55;
282   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
283   const Float_t  kPrim    = 48.0;
284   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
285   const Float_t  kPoti    = 12.1;
286
287   // PDG code electron
288   const Int_t    pdgElectron = 11;
289
290   // Set the maximum step size to a very large number for all 
291   // neutral particles and those outside the driftvolume
292   gMC->SetMaxStep(kBig); 
293
294   // Use only charged tracks 
295   if (( gMC->TrackCharge()       ) &&
296       (!gMC->IsTrackStop()       ) && 
297       (!gMC->IsTrackDisappeared())) {
298
299     // Inside a sensitive volume?
300     iIdSens = gMC->CurrentVolID(icSens);
301     if (iIdSens == fIdSens) { 
302
303       iIdSpace   = gMC->CurrentVolOffID(4,icSpace  );
304       iIdChamber = gMC->CurrentVolOffID(1,icChamber);
305
306       // Calculate the energy of the delta-electrons
307       eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
308       eDelta = TMath::Max(eDelta,0.0);
309
310       // The number of secondary electrons created
311       qTot = (Double_t) ((Int_t) (eDelta / kWion) + 1);
312
313       // The hit coordinates and charge
314       gMC->TrackPosition(pos);
315       hits[0] = pos[0];
316       hits[1] = pos[1];
317       hits[2] = pos[2];
318       hits[3] = qTot;
319
320       // The sector number (0 - 17)
321       // The numbering goes clockwise and starts at y = 0
322       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
323       if (phi < 90.) 
324         phi = phi + 270.;
325       else
326         phi = phi -  90.;
327       sec = ((Int_t) (phi / 20));
328
329       // The chamber number 
330       //   0: outer left
331       //   1: middle left
332       //   2: inner
333       //   3: middle right
334       //   4: outer right
335       if      (iIdChamber == fIdChamber1)
336         cha = (hits[2] < 0 ? 0 : 4);
337       else if (iIdChamber == fIdChamber2)       
338         cha = (hits[2] < 0 ? 1 : 3);
339       else if (iIdChamber == fIdChamber3)       
340         cha = 2;
341
342       // The plane number
343       // The numbering starts at the innermost plane
344       pla = icChamber - TMath::Nint((Float_t) (icChamber / 7)) * 6 - 1;
345
346       // Check on selected volumes
347       Int_t addthishit = 1;
348       if (fSensSelect) {
349         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
350         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
351         if ((fSensSector  >= 0) && (sec != fSensSector )) addthishit = 0;
352       }
353
354       // Add this hit
355       if (addthishit) {
356
357         new(lhits[fNhits++]) AliTRDhit(fIshunt
358                                       ,gAlice->CurrentTrack()
359                                       ,fGeometry->GetDetector(pla,cha,sec)
360                                       ,hits);
361
362         // The energy loss according to Bethe Bloch
363         gMC->TrackMomentum(mom);
364         pTot = mom.Rho();
365         iPdg = TMath::Abs(gMC->TrackPid());
366         if ( (iPdg != pdgElectron) ||
367             ((iPdg == pdgElectron) && (pTot < kPTotMax))) {
368           aMass     = gMC->TrackMass();
369           betaGamma = pTot / aMass;
370           pp        = kPrim * BetheBloch(betaGamma);
371           // Take charge > 1 into account
372           charge = gMC->TrackCharge();
373           if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
374         }
375         // Electrons above 20 Mev/c are at the plateau
376         else {
377           pp = kPrim * kPlateau;
378         }
379       
380         // Calculate the maximum step size for the next tracking step
381         if (pp > 0) {
382           do 
383             gMC->Rndm(random,1);
384           while ((random[0] == 1.) || (random[0] == 0.));
385           gMC->SetMaxStep( - TMath::Log(random[0]) / pp);
386         }
387
388       }
389       else {
390         // set step size to maximal value
391         gMC->SetMaxStep(kBig); 
392       }
393
394     }
395
396   }
397
398 }
399
400 //_____________________________________________________________________________
401 Double_t AliTRDv1::BetheBloch(Double_t bg) 
402 {
403   //
404   // Parametrization of the Bethe-Bloch-curve
405   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
406   //
407
408   // This parameters have been adjusted to averaged values from GEANT
409   const Double_t kP1 = 7.17960e-02;
410   const Double_t kP2 = 8.54196;
411   const Double_t kP3 = 1.38065e-06;
412   const Double_t kP4 = 5.30972;
413   const Double_t kP5 = 2.83798;
414
415   // This parameters have been adjusted to Xe-data found in:
416   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
417   //const Double_t kP1 = 0.76176E-1;
418   //const Double_t kP2 = 10.632;
419   //const Double_t kP3 = 3.17983E-6;
420   //const Double_t kP4 = 1.8631;
421   //const Double_t kP5 = 1.9479;
422
423   if (bg > 0) {
424     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
425     Double_t aa = TMath::Power(yy,kP4);
426     Double_t bb = TMath::Power((1./bg),kP5);
427              bb = TMath::Log(kP3 + bb);
428     return ((kP2 - aa - bb)*kP1 / aa);
429   }
430   else
431     return 0;
432
433 }
434
435 //_____________________________________________________________________________
436 Double_t Ermilova(Double_t *x, Double_t *)
437 {
438   //
439   // Calculates the delta-ray energy distribution according to Ermilova.
440   // Logarithmic scale !
441   //
442
443   Double_t energy;
444   Double_t dpos;
445   Double_t dnde;
446
447   Int_t    pos1, pos2;
448
449   const Int_t nV = 31;
450
451   Float_t vxe[nV] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120  
452                     , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
453                     , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
454                     , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
455                     , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
456                     , 9.4727, 9.9035,10.3735,10.5966,10.8198
457                     ,11.5129 };
458
459   Float_t vye[nV] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
460                     , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
461                     ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
462                     ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
463                     ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
464                     ,  0.04 ,  0.023,  0.015,  0.011,  0.01
465                     ,  0.004 };
466
467   energy = x[0];
468
469   // Find the position 
470   pos1 = pos2 = 0;
471   dpos = 0;
472   do {
473     dpos = energy - vxe[pos2++];
474   } 
475   while (dpos > 0);
476   pos2--; 
477   if (pos2 > nV) pos2 = nV;
478   pos1 = pos2 - 1;
479
480   // Differentiate between the sampling points
481   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
482
483   return dnde;
484
485 }