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