383c8e20778a6a47aabe93024a0c7a00ca4fc4f1
[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.40  2003/04/08 08:14:53  morsch
19 AddTrackReference moved to AliModule.
20
21 Revision 1.39  2003/02/11 16:54:07  hristov
22 Updated AliTrackReference class (S.Radomski)
23
24 Revision 1.38  2003/01/31 11:41:06  cblume
25 Fix bug in StepManager in treating geometry with holes
26
27 Revision 1.37  2003/01/28 14:38:18  cblume
28 Add track length to track references
29
30 Revision 1.36  2002/11/21 22:38:47  alibrary
31 Removing AliMC and AliMCProcess
32
33 Revision 1.35  2002/10/14 14:57:44  hristov
34 Merging the VirtualMC branch to the main development branch (HEAD)
35
36 Revision 1.33.6.2  2002/07/24 10:09:31  alibrary
37 Updating VirtualMC
38
39 Revision 1.34  2002/06/13 08:11:56  cblume
40 Add the track references
41
42 Revision 1.33  2002/02/20 14:01:40  hristov
43 Compare a TString with a string, otherwise the conversion cannot be done on Sun
44
45 Revision 1.32  2002/02/13 16:58:37  cblume
46 Bug fix reported by Jiri. Make atoi input zero terminated in StepManager()
47
48 Revision 1.31  2002/02/11 14:25:27  cblume
49 Geometry update, compressed hit structure
50
51 Revision 1.30  2001/05/21 16:45:47  hristov
52 Last minute changes (C.Blume)
53
54 Revision 1.29  2001/05/16 14:57:28  alibrary
55 New files for folders and Stack
56
57 Revision 1.28  2001/05/07 08:03:22  cblume
58 Generate also hits in the amplification region
59
60 Revision 1.27  2001/03/30 14:40:15  cblume
61 Update of the digitization parameter
62
63 Revision 1.26  2000/11/30 17:38:08  cblume
64 Changes to get in line with new STEER and EVGEN
65
66 Revision 1.25  2000/11/15 14:30:16  cblume
67 Fixed bug in calculating detector no. of extra hit
68
69 Revision 1.24  2000/11/10 14:58:36  cblume
70 Introduce additional hit with amplitude 0 at the chamber borders
71
72 Revision 1.23  2000/11/01 14:53:21  cblume
73 Merge with TRD-develop
74
75 Revision 1.17.2.5  2000/10/15 23:40:01  cblume
76 Remove AliTRDconst
77
78 Revision 1.17.2.4  2000/10/06 16:49:46  cblume
79 Made Getters const
80
81 Revision 1.17.2.3  2000/10/04 16:34:58  cblume
82 Replace include files by forward declarations
83
84 Revision 1.17.2.2  2000/09/18 13:50:17  cblume
85 Include TR photon generation and adapt to new AliTRDhit
86
87 Revision 1.22  2000/06/27 13:08:50  cblume
88 Changed to Copy(TObject &A) to appease the HP-compiler
89
90 Revision 1.21  2000/06/09 11:10:07  cblume
91 Compiler warnings and coding conventions, next round
92
93 Revision 1.20  2000/06/08 18:32:58  cblume
94 Make code compliant to coding conventions
95
96 Revision 1.19  2000/06/07 16:27:32  cblume
97 Try to remove compiler warnings on Sun and HP
98
99 Revision 1.18  2000/05/08 16:17:27  cblume
100 Merge TRD-develop
101
102 Revision 1.17.2.1  2000/05/08 14:59:16  cblume
103 Made inline function non-virtual. Bug fix in setting sensitive chamber
104
105 Revision 1.17  2000/02/28 19:10:26  cblume
106 Include the new TRD classes
107
108 Revision 1.16.4.1  2000/02/28 18:04:35  cblume
109 Change to new hit version, introduce geometry class, and move digitization and clustering to AliTRDdigitizer/AliTRDclusterizerV1
110
111 Revision 1.16  1999/11/05 22:50:28  fca
112 Do not use Atan, removed from ROOT too
113
114 Revision 1.15  1999/11/02 17:20:19  fca
115 initialise nbytes before using it
116
117 Revision 1.14  1999/11/02 17:15:54  fca
118 Correct ansi scoping not accepted by HP compilers
119
120 Revision 1.13  1999/11/02 17:14:51  fca
121 Correct ansi scoping not accepted by HP compilers
122
123 Revision 1.12  1999/11/02 16:35:56  fca
124 New version of TRD introduced
125
126 Revision 1.11  1999/11/01 20:41:51  fca
127 Added protections against using the wrong version of FRAME
128
129 Revision 1.10  1999/09/29 09:24:35  fca
130 Introduction of the Copyright and cvs Log
131
132 */
133
134 ///////////////////////////////////////////////////////////////////////////////
135 //                                                                           //
136 //  Transition Radiation Detector version 1 -- slow simulator                //
137 //                                                                           //
138 //Begin_Html
139 /*
140 <img src="picts/AliTRDfullClass.gif">
141 */
142 //End_Html
143 //                                                                           //
144 //                                                                           //
145 ///////////////////////////////////////////////////////////////////////////////
146
147 #include <stdlib.h> 
148
149 #include <TMath.h>
150 #include <TVector.h>
151 #include <TRandom.h>
152 #include <TF1.h>
153 #include <TLorentzVector.h>
154
155 #include "AliRun.h"
156 #include "AliConst.h"
157
158 #include "AliTRDv1.h"
159 #include "AliTRDhit.h"
160 #include "AliTRDmatrix.h"
161 #include "AliTRDgeometry.h"
162 #include "AliTRDsim.h"
163
164 ClassImp(AliTRDv1)
165  
166 //_____________________________________________________________________________
167 AliTRDv1::AliTRDv1():AliTRD()
168 {
169   //
170   // Default constructor
171   //
172
173   fSensSelect      =  0;
174   fSensPlane       = -1;
175   fSensChamber     = -1;
176   fSensSector      = -1;
177   fSensSectorRange =  0;
178
179   fDeltaE          = NULL;
180   fTR              = NULL;
181
182 }
183
184 //_____________________________________________________________________________
185 AliTRDv1::AliTRDv1(const char *name, const char *title) 
186          :AliTRD(name, title) 
187 {
188   //
189   // Standard constructor for Transition Radiation Detector version 1
190   //
191
192   fSensSelect      =  0;
193   fSensPlane       = -1;
194   fSensChamber     = -1;
195   fSensSector      = -1;
196   fSensSectorRange =  0;
197
198   fDeltaE          = NULL;
199   fTR              = NULL;
200
201   SetBufferSize(128000);
202
203 }
204
205 //_____________________________________________________________________________
206 AliTRDv1::AliTRDv1(const AliTRDv1 &trd)
207 {
208   //
209   // Copy constructor
210   //
211
212   ((AliTRDv1 &) trd).Copy(*this);
213
214 }
215
216 //_____________________________________________________________________________
217 AliTRDv1::~AliTRDv1()
218 {
219   //
220   // AliTRDv1 destructor
221   //
222
223   if (fDeltaE) delete fDeltaE;
224   if (fTR)     delete fTR;
225
226 }
227  
228 //_____________________________________________________________________________
229 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
230 {
231   //
232   // Assignment operator
233   //
234
235   if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
236   return *this;
237
238 }
239  
240 //_____________________________________________________________________________
241 void AliTRDv1::Copy(TObject &trd)
242 {
243   //
244   // Copy function
245   //
246
247   ((AliTRDv1 &) trd).fSensSelect      = fSensSelect;
248   ((AliTRDv1 &) trd).fSensPlane       = fSensPlane;
249   ((AliTRDv1 &) trd).fSensChamber     = fSensChamber;
250   ((AliTRDv1 &) trd).fSensSector      = fSensSector;
251   ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange;
252
253   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
254   fTR->Copy(*((AliTRDv1 &) trd).fTR);
255
256 }
257
258 //_____________________________________________________________________________
259 void AliTRDv1::CreateGeometry()
260 {
261   //
262   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
263   // This version covers the full azimuth. 
264   //
265
266   // Check that FRAME is there otherwise we have no place where to put the TRD
267   AliModule* frame = gAlice->GetModule("FRAME");
268   if (!frame) return;
269
270   // Define the chambers
271   AliTRD::CreateGeometry();
272
273 }
274
275 //_____________________________________________________________________________
276 void AliTRDv1::CreateMaterials()
277 {
278   //
279   // Create materials for the Transition Radiation Detector version 1
280   //
281
282   AliTRD::CreateMaterials();
283
284 }
285
286 //_____________________________________________________________________________
287 void AliTRDv1::CreateTRhit(Int_t det)
288 {
289   //
290   // Creates an electron cluster from a TR photon.
291   // The photon is assumed to be created a the end of the radiator. The 
292   // distance after which it deposits its energy takes into account the 
293   // absorbtion of the entrance window and of the gas mixture in drift
294   // volume.
295   //
296
297   // PDG code electron
298   const Int_t   kPdgElectron = 11;
299
300   // Ionization energy
301   const Float_t kWion        = 22.04;
302
303   // Maximum number of TR photons per track
304   const Int_t   kNTR         = 50;
305
306   TLorentzVector mom, pos;
307
308   // Create TR at the entrance of the chamber
309   if (gMC->IsTrackEntering()) {
310
311     // Create TR only for electrons 
312     Int_t iPdg = gMC->TrackPid();
313     if (TMath::Abs(iPdg) != kPdgElectron) return;
314
315     Float_t eTR[kNTR];
316     Int_t   nTR;
317
318     // Create TR photons
319     gMC->TrackMomentum(mom);
320     Float_t pTot = mom.Rho();
321     fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
322     if (nTR > kNTR) {
323       printf("AliTRDv1::CreateTRhit -- ");
324       printf("Boundary error: nTR = %d, kNTR = %d\n",nTR,kNTR);
325       exit(1);
326     }
327
328     // Loop through the TR photons
329     for (Int_t iTR = 0; iTR < nTR; iTR++) {
330
331       Float_t energyMeV = eTR[iTR] * 0.001;
332       Float_t energyeV  = eTR[iTR] * 1000.0;
333       Float_t absLength = 0;
334       Float_t sigma     = 0;
335
336       // Take the absorbtion in the entrance window into account
337       Double_t muMy = fTR->GetMuMy(energyMeV);
338       sigma = muMy * fFoilDensity;
339       absLength = gRandom->Exp(-sigma);
340       if (absLength < AliTRDgeometry::MyThick()) continue;
341
342       // The absorbtion cross sections in the drift gas
343       if (fGasMix == 1) {
344         // Gas-mixture (Xe/CO2)
345         Double_t muXe = fTR->GetMuXe(energyMeV);
346         Double_t muCO = fTR->GetMuCO(energyMeV);
347         sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity;
348       }
349       else {
350         // Gas-mixture (Xe/Isobutane) 
351         Double_t muXe = fTR->GetMuXe(energyMeV);
352         Double_t muBu = fTR->GetMuBu(energyMeV);
353         sigma = (0.97 * muXe + 0.03 * muBu) * fGasDensity;
354       }
355
356       // The distance after which the energy of the TR photon
357       // is deposited.
358       absLength = gRandom->Exp(-sigma);
359       if (absLength > AliTRDgeometry::DrThick()) continue;
360
361       // The position of the absorbtion
362       Float_t posHit[3];
363       gMC->TrackPosition(pos);
364       posHit[0] = pos[0] + mom[0] / pTot * absLength;
365       posHit[1] = pos[1] + mom[1] / pTot * absLength;
366       posHit[2] = pos[2] + mom[2] / pTot * absLength;      
367
368       // Create the charge 
369       Int_t q = ((Int_t) (energyeV / kWion));
370
371       // Add the hit to the array. TR photon hits are marked 
372       // by negative charge
373       AddHit(gAlice->CurrentTrack(),det,posHit,-q,kTRUE); 
374
375     }
376
377   }
378
379 }
380
381 //_____________________________________________________________________________
382 void AliTRDv1::Init() 
383 {
384   //
385   // Initialise Transition Radiation Detector after geometry has been built.
386   //
387
388   AliTRD::Init();
389
390   if(fDebug) printf("%s: Slow simulator\n",ClassName());
391   if (fSensSelect) {
392     if (fSensPlane   >= 0)
393       printf("          Only plane %d is sensitive\n",fSensPlane);
394     if (fSensChamber >= 0)   
395       printf("          Only chamber %d is sensitive\n",fSensChamber);
396     if (fSensSector  >= 0) {
397       Int_t sens1  = fSensSector;
398       Int_t sens2  = fSensSector + fSensSectorRange;
399             sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
400                    * AliTRDgeometry::Nsect();
401       printf("          Only sectors %d - %d are sensitive\n",sens1,sens2-1);
402     }
403   }
404   if (fTR) 
405     printf("%s: TR simulation on\n",ClassName());
406   else
407     printf("%s: TR simulation off\n",ClassName());
408   printf("\n");
409
410   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
411   const Float_t kPoti = 12.1;
412   // Maximum energy (50 keV);
413   const Float_t kEend = 50000.0;
414   // Ermilova distribution for the delta-ray spectrum
415   Float_t poti = TMath::Log(kPoti);
416   Float_t eEnd = TMath::Log(kEend);
417   fDeltaE = new TF1("deltae",Ermilova,poti,eEnd,0);
418
419   if(fDebug) {
420     printf("%s: ",ClassName());
421     for (Int_t i = 0; i < 80; i++) printf("*");
422     printf("\n");
423   }
424
425 }
426
427 //_____________________________________________________________________________
428 AliTRDsim *AliTRDv1::CreateTR()
429 {
430   //
431   // Enables the simulation of TR
432   //
433
434   fTR = new AliTRDsim();
435   return fTR;
436
437 }
438
439 //_____________________________________________________________________________
440 void AliTRDv1::SetSensPlane(Int_t iplane)
441 {
442   //
443   // Defines the hit-sensitive plane (0-5)
444   //
445
446   if ((iplane < 0) || (iplane > 5)) {
447     printf("Wrong input value: %d\n",iplane);
448     printf("Use standard setting\n");
449     fSensPlane  = -1;
450     fSensSelect =  0;
451     return;
452   }
453
454   fSensSelect = 1;
455   fSensPlane  = iplane;
456
457 }
458
459 //_____________________________________________________________________________
460 void AliTRDv1::SetSensChamber(Int_t ichamber)
461 {
462   //
463   // Defines the hit-sensitive chamber (0-4)
464   //
465
466   if ((ichamber < 0) || (ichamber > 4)) {
467     printf("Wrong input value: %d\n",ichamber);
468     printf("Use standard setting\n");
469     fSensChamber = -1;
470     fSensSelect  =  0;
471     return;
472   }
473
474   fSensSelect  = 1;
475   fSensChamber = ichamber;
476
477 }
478
479 //_____________________________________________________________________________
480 void AliTRDv1::SetSensSector(Int_t isector)
481 {
482   //
483   // Defines the hit-sensitive sector (0-17)
484   //
485
486   SetSensSector(isector,1);
487
488 }
489
490 //_____________________________________________________________________________
491 void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
492 {
493   //
494   // Defines a range of hit-sensitive sectors. The range is defined by
495   // <isector> (0-17) as the starting point and <nsector> as the number 
496   // of sectors to be included.
497   //
498
499   if ((isector < 0) || (isector > 17)) {
500     printf("Wrong input value <isector>: %d\n",isector);
501     printf("Use standard setting\n");
502     fSensSector      = -1;
503     fSensSectorRange =  0;
504     fSensSelect      =  0;
505     return;
506   }
507
508   if ((nsector < 1) || (nsector > 18)) {
509     printf("Wrong input value <nsector>: %d\n",nsector);
510     printf("Use standard setting\n");
511     fSensSector      = -1;
512     fSensSectorRange =  0;
513     fSensSelect      =  0;
514     return;
515   }
516
517   fSensSelect      = 1;
518   fSensSector      = isector;
519   fSensSectorRange = nsector;
520
521 }
522
523 //_____________________________________________________________________________
524 void AliTRDv1::StepManager()
525 {
526   //
527   // Slow simulator. Every charged track produces electron cluster as hits 
528   // along its path across the drift volume. The step size is set acording
529   // to Bethe-Bloch. The energy distribution of the delta electrons follows
530   // a spectrum taken from Ermilova et al.
531   //
532
533   Int_t    pla = 0;
534   Int_t    cha = 0;
535   Int_t    sec = 0;
536   Int_t    det = 0;
537   Int_t    iPdg;
538   Int_t    qTot;
539
540   Float_t  hits[3];
541   Double_t  random[1];
542   Float_t  charge;
543   Float_t  aMass;
544
545   Double_t pTot = 0;
546   Double_t eDelta;
547   Double_t betaGamma, pp;
548   Double_t stepSize;
549
550   Bool_t   drRegion = kFALSE;
551   Bool_t   amRegion = kFALSE;
552
553   TString  cIdCurrent;
554   TString  cIdSensDr = "J";
555   TString  cIdSensAm = "K";
556   Char_t   cIdChamber[3];
557            cIdChamber[2] = 0;
558
559   TLorentzVector pos, mom;
560
561   const Int_t    kNplan       = AliTRDgeometry::Nplan();
562   const Int_t    kNcham       = AliTRDgeometry::Ncham();
563   const Int_t    kNdetsec     = kNplan * kNcham;
564
565   const Double_t kBig         = 1.0E+12;
566
567   // Ionization energy
568   const Float_t  kWion        = 22.04;
569   // Maximum momentum for e+ e- g 
570   const Float_t  kPTotMaxEl   = 0.002;
571   // Minimum energy for the step size adjustment
572   const Float_t  kEkinMinStep = 1.0e-5;
573   // Plateau value of the energy-loss for electron in xenon
574   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
575   //const Double_t kPlateau = 1.70;
576   // the averaged value (26/3/99)
577   const Float_t  kPlateau     = 1.55;
578   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
579   const Float_t  kPrim        = 48.0;
580   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
581   const Float_t  kPoti        = 12.1;
582
583   // PDG code electron
584   const Int_t    kPdgElectron = 11;
585
586   // Set the maximum step size to a very large number for all 
587   // neutral particles and those outside the driftvolume
588   gMC->SetMaxStep(kBig); 
589
590   // Use only charged tracks 
591   if (( gMC->TrackCharge()       ) &&
592       (!gMC->IsTrackStop()       ) && 
593       (!gMC->IsTrackDisappeared())) {
594
595     // Inside a sensitive volume?
596     drRegion = kFALSE;
597     amRegion = kFALSE;
598     cIdCurrent = gMC->CurrentVolName();
599     if (cIdSensDr == cIdCurrent[1]) {
600       drRegion = kTRUE;
601     }
602     if (cIdSensAm == cIdCurrent[1]) {
603       amRegion = kTRUE;
604     }
605     if (drRegion || amRegion) {
606
607       // The hit coordinates and charge
608       gMC->TrackPosition(pos);
609       hits[0] = pos[0];
610       hits[1] = pos[1];
611       hits[2] = pos[2];
612
613       // The sector number (0 - 17)
614       // The numbering goes clockwise and starts at y = 0
615       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
616       if (phi < 90.) 
617         phi = phi + 270.;
618       else
619         phi = phi -  90.;
620       sec = ((Int_t) (phi / 20));
621
622       // The plane and chamber number
623       cIdChamber[0] = cIdCurrent[2];
624       cIdChamber[1] = cIdCurrent[3];
625       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
626       cha = ((Int_t) idChamber / kNplan);
627       pla = ((Int_t) idChamber % kNplan);
628
629       // Check on selected volumes
630       Int_t addthishit = 1;
631       if (fSensSelect) {
632         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
633         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
634         if (fSensSector  >= 0) {
635           Int_t sens1  = fSensSector;
636           Int_t sens2  = fSensSector + fSensSectorRange;
637                 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
638                        * AliTRDgeometry::Nsect();
639           if (sens1 < sens2) {
640             if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
641           }
642           else {
643             if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
644           }
645         }
646       }
647
648       // Add this hit
649       if (addthishit) {
650
651         // The detector number
652         det = fGeometry->GetDetector(pla,cha,sec);
653
654         // Special hits and TR photons only in the drift region
655         if (drRegion) {
656
657           // Create a track reference at the entrance and
658           // exit of each chamber that contain the 
659           // momentum components of the particle
660           if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
661             gMC->TrackMomentum(mom);
662             AddTrackReference(gAlice->CurrentTrack());
663           }
664
665           // Create the hits from TR photons
666           if (fTR) CreateTRhit(det);
667
668         }
669
670         // Calculate the energy of the delta-electrons
671         eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
672         eDelta = TMath::Max(eDelta,0.0);
673
674         // The number of secondary electrons created
675         qTot = ((Int_t) (eDelta / kWion) + 1);
676
677         // Create a new dEdx hit
678         if (drRegion) {
679           AddHit(gAlice->CurrentTrack(),det,hits,qTot,kTRUE);       
680         }
681         else {
682           AddHit(gAlice->CurrentTrack(),det,hits,qTot,kFALSE);      
683         }
684
685         // Calculate the maximum step size for the next tracking step
686         // Produce only one hit if Ekin is below cutoff 
687         aMass = gMC->TrackMass();
688         if ((gMC->Etot() - aMass) > kEkinMinStep) {
689
690           // The energy loss according to Bethe Bloch
691           iPdg  = TMath::Abs(gMC->TrackPid());
692           if ( (iPdg != kPdgElectron) ||
693               ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
694             gMC->TrackMomentum(mom);
695             pTot      = mom.Rho();
696             betaGamma = pTot / aMass;
697             pp        = kPrim * BetheBloch(betaGamma);
698             // Take charge > 1 into account
699             charge = gMC->TrackCharge();
700             if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
701           }
702           // Electrons above 20 Mev/c are at the plateau
703           else {
704             pp = kPrim * kPlateau;
705           }
706       
707           if (pp > 0) {
708             do 
709             gMC->GetRandom()->RndmArray(1, random);
710             while ((random[0] == 1.) || (random[0] == 0.));
711             stepSize = - TMath::Log(random[0]) / pp; 
712             gMC->SetMaxStep(stepSize);
713           }
714
715         }
716
717       }
718
719     }
720
721   }
722
723 }
724
725 //_____________________________________________________________________________
726 Double_t AliTRDv1::BetheBloch(Double_t bg) 
727 {
728   //
729   // Parametrization of the Bethe-Bloch-curve
730   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
731   //
732
733   // This parameters have been adjusted to averaged values from GEANT
734   const Double_t kP1 = 7.17960e-02;
735   const Double_t kP2 = 8.54196;
736   const Double_t kP3 = 1.38065e-06;
737   const Double_t kP4 = 5.30972;
738   const Double_t kP5 = 2.83798;
739
740   // This parameters have been adjusted to Xe-data found in:
741   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
742   //const Double_t kP1 = 0.76176E-1;
743   //const Double_t kP2 = 10.632;
744   //const Double_t kP3 = 3.17983E-6;
745   //const Double_t kP4 = 1.8631;
746   //const Double_t kP5 = 1.9479;
747
748   // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
749   const Double_t kBgMin = 0.8;
750   const Double_t kBBMax = 6.83298;
751   //const Double_t kBgMin = 0.6;
752   //const Double_t kBBMax = 17.2809;
753   //const Double_t kBgMin = 0.4;
754   //const Double_t kBBMax = 82.0;
755
756   if (bg > kBgMin) {
757     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
758     Double_t aa = TMath::Power(yy,kP4);
759     Double_t bb = TMath::Power((1./bg),kP5);
760              bb = TMath::Log(kP3 + bb);
761     return ((kP2 - aa - bb)*kP1 / aa);
762   }
763   else {
764     return kBBMax;
765   }
766
767 }
768
769 //_____________________________________________________________________________
770 Double_t Ermilova(Double_t *x, Double_t *)
771 {
772   //
773   // Calculates the delta-ray energy distribution according to Ermilova.
774   // Logarithmic scale !
775   //
776
777   Double_t energy;
778   Double_t dpos;
779   Double_t dnde;
780
781   Int_t    pos1, pos2;
782
783   const Int_t kNv = 31;
784
785   Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120  
786                      , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
787                      , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
788                      , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
789                      , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
790                      , 9.4727, 9.9035,10.3735,10.5966,10.8198
791                      ,11.5129 };
792
793   Float_t vye[kNv] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
794                      , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
795                      ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
796                      ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
797                      ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
798                      ,  0.04 ,  0.023,  0.015,  0.011,  0.01
799                      ,  0.004 };
800
801   energy = x[0];
802
803   // Find the position 
804   pos1 = pos2 = 0;
805   dpos = 0;
806   do {
807     dpos = energy - vxe[pos2++];
808   } 
809   while (dpos > 0);
810   pos2--; 
811   if (pos2 > kNv) pos2 = kNv - 1;
812   pos1 = pos2 - 1;
813
814   // Differentiate between the sampling points
815   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
816
817   return dnde;
818
819 }