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