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