Modify services of top most layer (inside TOF acceptance)
[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 /* $Id$ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //  Transition Radiation Detector version 1 -- slow simulator             //
21 //                                                                        //
22 ////////////////////////////////////////////////////////////////////////////
23
24 #include <stdlib.h> 
25
26 #include <TF1.h>
27 #include <TLorentzVector.h>
28 #include <TMath.h>
29 #include <TRandom.h>
30 #include <TVector.h>
31 #include <TVirtualMC.h>
32 #include <TGeoManager.h>
33
34 #include "AliConst.h"
35 #include "AliLog.h"
36 #include "AliMC.h"
37 #include "AliRun.h"
38
39 #include "AliTRDgeometry.h"
40 #include "AliTRDhit.h"
41 #include "AliTRDsim.h"
42 #include "AliTRDv1.h"
43
44 ClassImp(AliTRDv1)
45  
46 //_____________________________________________________________________________
47 AliTRDv1::AliTRDv1()
48   :AliTRD()
49   ,fTRon(kFALSE)
50   ,fTR(NULL)
51   ,fTypeOfStepManager(0)
52   ,fStepSize(0)
53   ,fDeltaE(NULL)
54   ,fDeltaG(NULL)
55   ,fTrackLength0(0)
56   ,fPrimaryTrackPid(0)
57 {
58   //
59   // Default constructor
60   //
61
62 }
63
64 //_____________________________________________________________________________
65 AliTRDv1::AliTRDv1(const char *name, const char *title) 
66   :AliTRD(name,title) 
67   ,fTRon(kTRUE)
68   ,fTR(NULL)
69   ,fTypeOfStepManager(2)
70   ,fStepSize(0.1)
71   ,fDeltaE(NULL)
72   ,fDeltaG(NULL)
73   ,fTrackLength0(0)
74   ,fPrimaryTrackPid(0)
75 {
76   //
77   // Standard constructor for Transition Radiation Detector version 1
78   //
79
80   SetBufferSize(128000);
81
82 }
83
84 //_____________________________________________________________________________
85 AliTRDv1::AliTRDv1(const AliTRDv1 &trd)
86   :AliTRD(trd)
87   ,fTRon(trd.fTRon)
88   ,fTR(NULL)
89   ,fTypeOfStepManager(trd.fTypeOfStepManager)
90   ,fStepSize(trd.fStepSize)
91   ,fDeltaE(NULL)
92   ,fDeltaG(NULL)
93   ,fTrackLength0(trd.fTrackLength0)
94   ,fPrimaryTrackPid(trd.fPrimaryTrackPid)
95 {
96   //
97   // Copy constructor
98   //
99
100   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
101   fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG);
102   fTR->Copy(*((AliTRDv1 &) trd).fTR);
103
104 }
105
106 //_____________________________________________________________________________
107 AliTRDv1::~AliTRDv1()
108 {
109   //
110   // AliTRDv1 destructor
111   //
112
113   if (fDeltaE) {
114     delete fDeltaE;
115     fDeltaE = 0;
116   }
117
118   if (fDeltaG) {
119     delete fDeltaG;
120     fDeltaG = 0;
121   }
122
123   if (fTR) {
124     delete fTR;
125     fTR     = 0;
126   }
127
128 }
129  
130 //_____________________________________________________________________________
131 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
132 {
133   //
134   // Assignment operator
135   //
136
137   if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
138
139   return *this;
140
141 }
142  
143 //_____________________________________________________________________________
144 void AliTRDv1::Copy(TObject &trd) const
145 {
146   //
147   // Copy function
148   //
149
150   ((AliTRDv1 &) trd).fTypeOfStepManager = fTypeOfStepManager;
151   ((AliTRDv1 &) trd).fStepSize          = fStepSize;
152   ((AliTRDv1 &) trd).fTRon              = fTRon;
153   ((AliTRDv1 &) trd).fTrackLength0      = fTrackLength0;
154   ((AliTRDv1 &) trd).fPrimaryTrackPid   = fPrimaryTrackPid;
155
156   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
157   fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG);
158   fTR->Copy(*((AliTRDv1 &) trd).fTR);
159
160 }
161
162 //_____________________________________________________________________________
163 void AliTRDv1::AddAlignableVolumes() const
164 {
165   //
166   // Create entries for alignable volumes associating the symbolic volume
167   // name with the corresponding volume path. Needs to be syncronized with
168   // eventual changes in the geometry.
169   //
170
171   TString volPath;
172   TString symName;
173
174   TString vpStr  = "ALIC_1/B077_1/BSEGMO";
175   TString vpApp1 = "_1/BTRD";
176   TString vpApp2 = "_1";
177   TString vpApp3 = "/UTR1_1/UTS1_1/UTI1_1/UT";
178
179   TString snStr  = "TRD/sm";
180   TString snApp1 = "/st";
181   TString snApp2 = "/pl";
182
183   //
184   // The super modules
185   // The symbolic names are: TRD/sm00
186   //                           ...
187   //                         TRD/sm17
188   //
189   for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
190
191     volPath  = vpStr;
192     volPath += isect;
193     volPath += vpApp1;
194     volPath += isect;
195     volPath += vpApp2;
196
197     symName  = snStr;
198     symName += Form("%02d",isect);
199
200     gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
201
202   }
203
204   //
205   // The readout chambers
206   // The symbolic names are: TRD/sm00/st0/pl0
207   //                           ...
208   //                         TRD/sm17/st4/pl5
209   //
210   for (Int_t isect = 0; isect < AliTRDgeometry::Nsect(); isect++) {
211     for (Int_t icham = 0; icham < AliTRDgeometry::Ncham(); icham++) {
212       for (Int_t iplan = 0; iplan < AliTRDgeometry::Nplan(); iplan++) {
213
214         Int_t idet = AliTRDgeometry::GetDetectorSec(iplan,icham);
215
216         volPath  = vpStr;
217         volPath += isect;
218         volPath += vpApp1;
219         volPath += isect;
220         volPath += vpApp2;
221         volPath += vpApp3;
222         volPath += Form("%02d",idet);
223         volPath += vpApp2;
224
225         symName  = snStr;
226         symName += Form("%02d",isect);
227         symName += snApp1;
228         symName += icham;
229         symName += snApp2;
230         symName += iplan;
231
232         gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
233
234       }
235     }
236   }
237
238 }
239
240 //_____________________________________________________________________________
241 void AliTRDv1::CreateGeometry()
242 {
243   //
244   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
245   // This version covers the full azimuth. 
246   //
247
248   // Check that FRAME is there otherwise we have no place where to put the TRD
249   AliModule* frame = gAlice->GetModule("FRAME");
250   if (!frame) {
251     AliError("TRD needs FRAME to be present\n");
252     return;
253   }
254
255   // Define the chambers
256   AliTRD::CreateGeometry();
257
258 }
259
260 //_____________________________________________________________________________
261 void AliTRDv1::CreateMaterials()
262 {
263   //
264   // Create materials for the Transition Radiation Detector version 1
265   //
266
267   AliTRD::CreateMaterials();
268
269 }
270
271 //_____________________________________________________________________________
272 void AliTRDv1::CreateTRhit(Int_t det)
273 {
274   //
275   // Creates an electron cluster from a TR photon.
276   // The photon is assumed to be created a the end of the radiator. The 
277   // distance after which it deposits its energy takes into account the 
278   // absorbtion of the entrance window and of the gas mixture in drift
279   // volume.
280   //
281
282   // PDG code electron
283   const Int_t   kPdgElectron = 11;
284
285   // Ionization energy
286   const Float_t kWion        = 23.53;
287
288   // Maximum number of TR photons per track
289   const Int_t   kNTR         = 50;
290
291   TLorentzVector mom;
292   TLorentzVector pos;
293
294   // Create TR at the entrance of the chamber
295   if (gMC->IsTrackEntering()) {
296
297     // Create TR only for electrons 
298     Int_t iPdg = gMC->TrackPid();
299     if (TMath::Abs(iPdg) != kPdgElectron) {
300       return;
301     }
302
303     Float_t eTR[kNTR];
304     Int_t   nTR;
305
306     // Create TR photons
307     gMC->TrackMomentum(mom);
308     Float_t pTot = mom.Rho();
309     fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
310     if (nTR > kNTR) {
311       AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
312     }
313
314     // Loop through the TR photons
315     for (Int_t iTR = 0; iTR < nTR; iTR++) {
316
317       Float_t energyMeV = eTR[iTR] * 0.001;
318       Float_t energyeV  = eTR[iTR] * 1000.0;
319       Float_t absLength = 0.0;
320       Float_t sigma     = 0.0;
321
322       // Take the absorbtion in the entrance window into account
323       Double_t muMy = fTR->GetMuMy(energyMeV);
324       sigma         = muMy * fFoilDensity;
325       if (sigma > 0.0) {
326         absLength = gRandom->Exp(1.0/sigma);
327         if (absLength < AliTRDgeometry::MyThick()) {
328           continue;
329         }
330       }
331       else {
332         continue;
333       }
334
335       // The absorbtion cross sections in the drift gas
336       // Gas-mixture (Xe/CO2)
337       Double_t muXe = fTR->GetMuXe(energyMeV);
338       Double_t muCO = fTR->GetMuCO(energyMeV);
339       sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
340
341       // The distance after which the energy of the TR photon
342       // is deposited.
343       if (sigma > 0.0) {
344         absLength = gRandom->Exp(1.0/sigma);
345         if (absLength > (AliTRDgeometry::DrThick()
346                        + AliTRDgeometry::AmThick())) {
347           continue;
348         }
349       }
350       else {
351         continue;
352       }
353
354       // The position of the absorbtion
355       Float_t posHit[3];
356       gMC->TrackPosition(pos);
357       posHit[0] = pos[0] + mom[0] / pTot * absLength;
358       posHit[1] = pos[1] + mom[1] / pTot * absLength;
359       posHit[2] = pos[2] + mom[2] / pTot * absLength;
360
361       // Create the charge 
362       Int_t q = ((Int_t) (energyeV / kWion));
363
364       // Add the hit to the array. TR photon hits are marked 
365       // by negative charge
366       AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
367             ,det
368             ,posHit
369             ,-q
370             ,kTRUE); 
371
372     }
373
374   }
375
376 }
377
378 //_____________________________________________________________________________
379 void AliTRDv1::Init() 
380 {
381   //
382   // Initialise Transition Radiation Detector after geometry has been built.
383   //
384
385   AliTRD::Init();
386
387   AliDebug(1,"Slow simulator\n");
388
389   // Switch on TR simulation as default
390   if (!fTRon) {
391     AliInfo("TR simulation off");
392   }
393   else {
394     fTR = new AliTRDsim();
395   }
396
397   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
398   const Float_t kPoti = 12.1;
399   // Maximum energy (50 keV);
400   const Float_t kEend = 50000.0;
401   // Ermilova distribution for the delta-ray spectrum
402   Float_t poti        = TMath::Log(kPoti);
403   Float_t eEnd        = TMath::Log(kEend);
404
405   // Ermilova distribution for the delta-ray spectrum
406   fDeltaE = new TF1("deltae" ,Ermilova ,poti,eEnd,0);
407
408   // Geant3 distribution for the delta-ray spectrum
409   fDeltaG = new TF1("deltag",IntSpecGeant,2.421257,28.536469,0);
410
411   AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
412
413 }
414
415 //_____________________________________________________________________________
416 void AliTRDv1::StepManager()
417 {
418   //
419   // Slow simulator. Every charged track produces electron cluster as hits
420   // along its path across the drift volume. 
421   //
422
423   switch (fTypeOfStepManager) {
424    case 0: 
425     StepManagerErmilova();
426     break;  
427    case 1: 
428     StepManagerGeant();  
429     break;  
430    case 2: 
431     StepManagerFixedStep();
432     break;  
433    default: 
434     AliWarning("Not a valid Step Manager.");
435   }
436
437 }
438
439 //_____________________________________________________________________________
440 void AliTRDv1::SelectStepManager(Int_t t)
441 {
442   //
443   // Selects a step manager type:
444   //   0 - Ermilova
445   //   1 - Geant3
446   //   2 - Fixed step size
447   //
448
449   fTypeOfStepManager = t;
450   AliInfo(Form("Step Manager type %d was selected",fTypeOfStepManager));
451
452 }
453
454 //_____________________________________________________________________________
455 void AliTRDv1::StepManagerGeant()
456 {
457   //
458   // Slow simulator. Every charged track produces electron cluster as hits
459   // along its path across the drift volume. The step size is set acording
460   // to Bethe-Bloch. The energy distribution of the delta electrons follows
461   // a spectrum taken from Geant3.
462   //
463   // Version by A. Bercuci
464   //
465
466   Int_t    pla = 0;
467   Int_t    cha = 0;
468   Int_t    sec = 0;
469   Int_t    det = 0;
470   Int_t    iPdg;
471   Int_t    qTot;
472
473   Float_t  hits[3];
474   Float_t  charge;
475   Float_t  aMass;
476
477   Double_t pTot     = 0;
478   Double_t eDelta;
479   Double_t betaGamma;
480   Double_t pp;
481   Double_t stepSize = 0;
482
483   Bool_t   drRegion = kFALSE;
484   Bool_t   amRegion = kFALSE;
485
486   TString  cIdCurrent;
487   TString  cIdSensDr = "J";
488   TString  cIdSensAm = "K";
489   Char_t   cIdChamber[3];
490            cIdChamber[2] = 0;
491
492   TLorentzVector pos;
493   TLorentzVector mom;
494
495   TArrayI        processes;
496
497   const Int_t    kNplan       = AliTRDgeometry::Nplan();
498   const Int_t    kNcham       = AliTRDgeometry::Ncham();
499   const Int_t    kNdetsec     = kNplan * kNcham;
500
501   const Double_t kBig         = 1.0e+12; // Infinitely big
502   const Float_t  kWion        = 23.53;   // Ionization energy
503   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g
504
505   // Minimum energy for the step size adjustment
506   const Float_t  kEkinMinStep = 1.0e-5;
507   // energy threshold for production of delta electrons
508   const Float_t  kECut        = 1.0e4;
509   // Parameters entering the parametrized range for delta electrons
510   const Float_t  kRa          = 5.37e-4;
511   const Float_t  kRb          = 0.9815;
512   const Float_t  kRc          = 3.123e-3;
513   // Gas density -> To be made user adjustable !
514   // [0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
515   const Float_t  kRho         = 0.004945 ; 
516
517   // Plateau value of the energy-loss for electron in xenon
518   // The averaged value (26/3/99)
519   const Float_t  kPlateau     = 1.55;
520   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
521   const Float_t  kPrim        = 19.34;  
522   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
523   const Float_t  kPoti        = 12.1;
524   // PDG code electron
525   const Int_t    kPdgElectron = 11;  
526
527   // Set the maximum step size to a very large number for all
528   // neutral particles and those outside the driftvolume
529   gMC->SetMaxStep(kBig);
530
531   // Use only charged tracks
532   if (( gMC->TrackCharge()       ) &&
533       (!gMC->IsTrackDisappeared())) {
534
535     // Inside a sensitive volume?
536     drRegion = kFALSE;
537     amRegion = kFALSE;
538     cIdCurrent = gMC->CurrentVolName();
539     if (cIdSensDr == cIdCurrent[1]) {
540       drRegion = kTRUE;
541     }
542     if (cIdSensAm == cIdCurrent[1]) {
543       amRegion = kTRUE;
544     }
545     if (drRegion || amRegion) {
546
547       // The hit coordinates and charge
548       gMC->TrackPosition(pos);
549       hits[0] = pos[0];
550       hits[1] = pos[1];
551       hits[2] = pos[2];
552
553       // The sector number (0 - 17)
554       // The numbering goes clockwise and starts at y = 0
555       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
556       if (phi < 90.0) {
557         phi = phi + 270.0;
558       }
559       else {
560         phi = phi -  90.0;
561       }
562       sec = ((Int_t) (phi / 20.0));
563
564       // The plane and chamber number
565       cIdChamber[0]   = cIdCurrent[2];
566       cIdChamber[1]   = cIdCurrent[3];
567       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
568       cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
569       pla = ((Int_t) idChamber % kNplan);
570
571       // Check on selected volumes
572       Int_t addthishit = 1;
573
574       // Add this hit
575       if (addthishit) {
576
577         // The detector number
578         det = fGeometry->GetDetector(pla,cha,sec);
579
580         // Special hits only in the drift region
581         if (drRegion) {
582
583           // Create a track reference at the entrance and
584           // exit of each chamber that contain the
585           // momentum components of the particle
586           if (gMC->IsTrackEntering() || 
587               gMC->IsTrackExiting()) {
588             gMC->TrackMomentum(mom);
589             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
590           }
591
592           if (gMC->IsTrackEntering() && 
593               !gMC->IsNewTrack()) {
594             // determine if hit belong to primary track 
595             fPrimaryTrackPid = gAlice->GetMCApp()->GetCurrentTrackNumber();
596             // determine track length when entering the detector
597             fTrackLength0    = gMC->TrackLength();
598           }
599                                         
600           // Create the hits from TR photons
601           if (fTR) CreateTRhit(det);
602
603         }
604
605         // Calculate the energy of the delta-electrons
606         // modified by Alex Bercuci (A.Bercuci@gsi.de) on 26.01.06
607         // take into account correlation with the underlying GEANT tracking
608         // mechanism. see
609         // http://www-linux.gsi.de/~abercuci/Contributions/TRD/index.html
610         //
611         // determine the most significant process (last on the processes list)
612         // which caused this hit
613         gMC->StepProcesses(processes);
614         Int_t nofprocesses = processes.GetSize();
615         Int_t pid;
616         if (!nofprocesses) {
617           pid = 0;
618         }
619         else {
620           pid = processes[nofprocesses-1];              
621         }               
622                 
623         // generate Edep according to GEANT parametrisation
624         eDelta = TMath::Exp(fDeltaG->GetRandom()) - kPoti;
625         eDelta = TMath::Max(eDelta,0.0);
626         Float_t prRange = 0.0;
627         Float_t range   = gMC->TrackLength() - fTrackLength0;
628         // merge GEANT tracker information with locally cooked one
629         if (gAlice->GetMCApp()->GetCurrentTrackNumber() == fPrimaryTrackPid) {
630           if      (pid == 27) { 
631             if (eDelta >= kECut) {                
632               prRange = kRa * eDelta * 0.001
633                       * (1.0 - kRb / (1.0 + kRc * eDelta * 0.001)) / kRho;
634               if (prRange >= (3.7 - range)) {
635                 eDelta *= 0.1;
636               }
637             }
638           } 
639           else if (pid ==  1) { 
640             if (eDelta <  kECut) {
641               eDelta *= 0.5;
642             }
643             else {                
644               prRange = kRa * eDelta * 0.001
645                       * (1.0 - kRb / (1.0 + kRc * eDelta * 0.001)) / kRho;
646               if (prRange >= ((AliTRDgeometry::DrThick()
647                              + AliTRDgeometry::AmThick()) - range)) {
648                 eDelta *= 0.05;
649               }
650               else {
651                 eDelta *= 0.5;
652               }
653             }
654           } 
655           else {
656             eDelta = 0.0;
657           }     
658         } 
659         else {
660           eDelta = 0.0;
661         }
662
663         // Generate the electron cluster size
664         if (eDelta == 0.0) {
665           qTot = 0;
666         }
667         else {
668           qTot = ((Int_t) (eDelta / kWion) + 1);
669         }
670
671         // Create a new dEdx hit
672         AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
673               ,det
674               ,hits
675               ,qTot
676               ,drRegion);
677                                 
678         // Calculate the maximum step size for the next tracking step
679         // Produce only one hit if Ekin is below cutoff
680         aMass = gMC->TrackMass();
681         if ((gMC->Etot() - aMass) > kEkinMinStep) {
682
683           // The energy loss according to Bethe Bloch
684           iPdg = TMath::Abs(gMC->TrackPid());
685           if ((iPdg != kPdgElectron) ||
686               ((iPdg == kPdgElectron) && 
687                (pTot  < kPTotMaxEl))) {
688             gMC->TrackMomentum(mom);
689             pTot      = mom.Rho();
690             betaGamma = pTot / aMass;
691             pp        = BetheBlochGeant(betaGamma);
692             // Take charge > 1 into account
693             charge     = gMC->TrackCharge();
694             if (TMath::Abs(charge) > 1) {
695               pp = pp * charge*charge;
696             }
697           } 
698           else { 
699             // Electrons above 20 Mev/c are at the plateau
700             pp = kPrim * kPlateau;
701           }
702
703           Int_t nsteps = 0;
704           do {
705             nsteps = gRandom->Poisson(pp);
706           } while(!nsteps);
707           stepSize = 1.0 / nsteps;
708           gMC->SetMaxStep(stepSize);
709
710         }
711
712       }
713
714     }
715
716   }
717
718 }
719
720 //_____________________________________________________________________________
721 void AliTRDv1::StepManagerErmilova()
722 {
723   //
724   // Slow simulator. Every charged track produces electron cluster as hits 
725   // along its path across the drift volume. The step size is set acording
726   // to Bethe-Bloch. The energy distribution of the delta electrons follows
727   // a spectrum taken from Ermilova et al.
728   //
729
730   Int_t    pla = 0;
731   Int_t    cha = 0;
732   Int_t    sec = 0;
733   Int_t    det = 0;
734   Int_t    iPdg;
735   Int_t    qTot;
736
737   Float_t  hits[3];
738   Double_t random[1];
739   Float_t  charge;
740   Float_t  aMass;
741
742   Double_t pTot     = 0.0;
743   Double_t eDelta;
744   Double_t betaGamma;
745   Double_t pp;
746   Double_t stepSize;
747
748   Bool_t   drRegion = kFALSE;
749   Bool_t   amRegion = kFALSE;
750
751   TString  cIdCurrent;
752   TString  cIdSensDr = "J";
753   TString  cIdSensAm = "K";
754   Char_t   cIdChamber[3];
755            cIdChamber[2] = 0;
756
757   TLorentzVector pos;
758   TLorentzVector mom;
759
760   const Int_t    kNplan       = AliTRDgeometry::Nplan();
761   const Int_t    kNcham       = AliTRDgeometry::Ncham();
762   const Int_t    kNdetsec     = kNplan * kNcham;
763
764   const Double_t kBig         = 1.0e+12; // Infinitely big
765   const Float_t  kWion        = 23.53;   // Ionization energy
766   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g 
767
768   // Minimum energy for the step size adjustment
769   const Float_t  kEkinMinStep = 1.0e-5;
770
771   // Plateau value of the energy-loss for electron in xenon
772   // The averaged value (26/3/99)
773   const Float_t  kPlateau     = 1.55;
774   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
775   const Float_t  kPrim        = 48.0;  
776   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
777   const Float_t  kPoti        = 12.1;
778   // PDG code electron
779   const Int_t    kPdgElectron = 11;  
780
781   // Set the maximum step size to a very large number for all 
782   // neutral particles and those outside the driftvolume
783   gMC->SetMaxStep(kBig); 
784
785   // Use only charged tracks 
786   if (( gMC->TrackCharge()       ) &&
787       (!gMC->IsTrackDisappeared())) {
788
789     // Inside a sensitive volume?
790     drRegion = kFALSE;
791     amRegion = kFALSE;
792     cIdCurrent = gMC->CurrentVolName();
793     if (cIdSensDr == cIdCurrent[1]) {
794       drRegion = kTRUE;
795     }
796     if (cIdSensAm == cIdCurrent[1]) {
797       amRegion = kTRUE;
798     }
799     if (drRegion || amRegion) {
800
801       // The hit coordinates and charge
802       gMC->TrackPosition(pos);
803       hits[0] = pos[0];
804       hits[1] = pos[1];
805       hits[2] = pos[2];
806
807       // The sector number (0 - 17)
808       // The numbering goes clockwise and starts at y = 0
809       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
810       if (phi < 90.0) { 
811         phi = phi + 270.0;
812       }
813       else {
814         phi = phi -  90.0;
815       }
816       sec = ((Int_t) (phi / 20.0));
817
818       // The plane and chamber number
819       cIdChamber[0] = cIdCurrent[2];
820       cIdChamber[1] = cIdCurrent[3];
821       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
822       cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
823       pla = ((Int_t) idChamber % kNplan);
824
825       // Check on selected volumes
826       Int_t addthishit = 1;
827
828       // Add this hit
829       if (addthishit) {
830
831         // The detector number
832         det = fGeometry->GetDetector(pla,cha,sec);
833
834         // Special hits only in the drift region
835         if (drRegion) {
836
837           // Create a track reference at the entrance and
838           // exit of each chamber that contain the 
839           // momentum components of the particle
840           if (gMC->IsTrackEntering() || 
841               gMC->IsTrackExiting()) {
842             gMC->TrackMomentum(mom);
843             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
844           }
845           // Create the hits from TR photons
846           if (fTR) {
847             CreateTRhit(det);
848           }
849
850         }
851
852         // Calculate the energy of the delta-electrons
853         eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
854         eDelta = TMath::Max(eDelta,0.0);
855         // Generate the electron cluster size
856         if (eDelta == 0.0) {
857           qTot = 0;
858         }
859         else {
860           qTot = ((Int_t) (eDelta / kWion) + 1);
861         }
862
863         // Create a new dEdx hit
864         if (drRegion) {
865           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
866                 ,det
867                 ,hits
868                 ,qTot
869                 ,kTRUE);
870         }
871         else {
872           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
873                 ,det
874                 ,hits
875                 ,qTot
876                 ,kFALSE);
877         }
878
879         // Calculate the maximum step size for the next tracking step
880         // Produce only one hit if Ekin is below cutoff 
881         aMass = gMC->TrackMass();
882         if ((gMC->Etot() - aMass) > kEkinMinStep) {
883
884           // The energy loss according to Bethe Bloch
885           iPdg  = TMath::Abs(gMC->TrackPid());
886           if ((iPdg != kPdgElectron) ||
887               ((iPdg == kPdgElectron) && 
888                (pTot  < kPTotMaxEl))) {
889             gMC->TrackMomentum(mom);
890             pTot      = mom.Rho();
891             betaGamma = pTot / aMass;
892             pp        = kPrim * BetheBloch(betaGamma);
893             // Take charge > 1 into account
894             charge = gMC->TrackCharge();
895             if (TMath::Abs(charge) > 1) {
896               pp = pp * charge*charge;
897             }
898           } 
899           else { 
900             // Electrons above 20 Mev/c are at the plateau
901             pp = kPrim * kPlateau;
902           }
903       
904           if (pp > 0.0) {
905             do {
906               gMC->GetRandom()->RndmArray(1,random);
907             }
908             while ((random[0] == 1.0) || 
909                    (random[0] == 0.0));
910             stepSize = - TMath::Log(random[0]) / pp; 
911             gMC->SetMaxStep(stepSize);
912           }
913
914         }
915
916       }
917
918     }
919
920   }
921
922 }
923
924 //_____________________________________________________________________________
925 void AliTRDv1::StepManagerFixedStep()
926 {
927   //
928   // Slow simulator. Every charged track produces electron cluster as hits 
929   // along its path across the drift volume. The step size is fixed in
930   // this version of the step manager.
931   //
932
933   Int_t    pla = 0;
934   Int_t    cha = 0;
935   Int_t    sec = 0;
936   Int_t    det = 0;
937   Int_t    qTot;
938
939   Float_t  hits[3];
940   Double_t eDep;
941
942   Bool_t   drRegion = kFALSE;
943   Bool_t   amRegion = kFALSE;
944
945   TString  cIdCurrent;
946   TString  cIdSensDr = "J";
947   TString  cIdSensAm = "K";
948   Char_t   cIdChamber[3];
949   cIdChamber[2] = 0;
950
951   TLorentzVector pos;
952   TLorentzVector mom;
953
954   const Int_t    kNplan       = AliTRDgeometry::Nplan();
955   const Int_t    kNcham       = AliTRDgeometry::Ncham();
956   const Int_t    kNdetsec     = kNplan * kNcham;
957
958   const Double_t kBig         = 1.0e+12;
959
960   const Float_t  kWion        = 23.53;   // Ionization energy
961   const Float_t  kEkinMinStep = 1.0e-5;  // Minimum energy for the step size adjustment
962
963   // Set the maximum step size to a very large number for all 
964   // neutral particles and those outside the driftvolume
965   gMC->SetMaxStep(kBig); 
966
967   // If not charged track or already stopped or disappeared, just return.
968   if ((!gMC->TrackCharge()) || 
969         gMC->IsTrackDisappeared()) return;
970
971   // Inside a sensitive volume?
972   cIdCurrent = gMC->CurrentVolName();
973
974   if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
975   if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
976
977   if ((!drRegion) && 
978       (!amRegion)) {
979     return;
980   }
981
982   // The hit coordinates and charge
983   gMC->TrackPosition(pos);
984   hits[0] = pos[0];
985   hits[1] = pos[1];
986   hits[2] = pos[2];
987
988   // The sector number (0 - 17)
989   // The numbering goes clockwise and starts at y = 0
990   Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
991   if (phi < 90.0) {
992     phi = phi + 270.0;
993   }
994   else {          
995     phi = phi -  90.0;
996   }
997   sec = ((Int_t) (phi / 20.0));
998
999   // The plane and chamber number
1000   cIdChamber[0]   = cIdCurrent[2];
1001   cIdChamber[1]   = cIdCurrent[3];
1002   Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
1003   cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
1004   pla = ((Int_t) idChamber % kNplan);
1005
1006   // Check on selected volumes
1007   Int_t addthishit = 1;
1008
1009   if (!addthishit) {
1010     return;
1011   }
1012
1013   // The detector number
1014   det = fGeometry->GetDetector(pla,cha,sec);
1015
1016   // 0: InFlight 1:Entering 2:Exiting
1017   Int_t trkStat = 0;
1018
1019   // Special hits only in the drift region
1020   if (drRegion) {
1021
1022     // Create a track reference at the entrance and exit of each
1023     // chamber that contain the momentum components of the particle
1024
1025     if (gMC->IsTrackEntering()) {
1026       gMC->TrackMomentum(mom);
1027       AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
1028       trkStat = 1;
1029     }
1030     if (gMC->IsTrackExiting()) {
1031       gMC->TrackMomentum(mom);
1032       AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
1033       trkStat = 2;
1034     }
1035
1036     // Create the hits from TR photons
1037     if (fTR) {
1038       CreateTRhit(det);    
1039     }
1040
1041   }
1042   
1043   // Calculate the charge according to GEANT Edep
1044   // Create a new dEdx hit
1045   eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
1046   qTot = (Int_t) (eDep / kWion);
1047   AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
1048         ,det
1049         ,hits
1050         ,qTot
1051         ,drRegion);
1052
1053   // Set Maximum Step Size
1054   // Produce only one hit if Ekin is below cutoff
1055   if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) {
1056     return;
1057   }
1058   gMC->SetMaxStep(fStepSize);
1059
1060 }
1061
1062 //_____________________________________________________________________________
1063 Double_t AliTRDv1::BetheBloch(Double_t bg) 
1064 {
1065   //
1066   // Parametrization of the Bethe-Bloch-curve
1067   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
1068   //
1069
1070   // This parameters have been adjusted to averaged values from GEANT
1071   const Double_t kP1    = 7.17960e-02;
1072   const Double_t kP2    = 8.54196;
1073   const Double_t kP3    = 1.38065e-06;
1074   const Double_t kP4    = 5.30972;
1075   const Double_t kP5    = 2.83798;
1076
1077   // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
1078   const Double_t kBgMin = 0.8;
1079   const Double_t kBBMax = 6.83298;
1080
1081   if (bg > kBgMin) {
1082     Double_t yy = bg / TMath::Sqrt(1.0 + bg*bg);
1083     Double_t aa = TMath::Power(yy,kP4);
1084     Double_t bb = TMath::Power((1.0/bg),kP5);
1085              bb = TMath::Log(kP3 + bb);
1086     return ((kP2 - aa - bb) * kP1 / aa);
1087   }
1088   else {
1089     return kBBMax;
1090   }
1091
1092 }
1093
1094 //_____________________________________________________________________________
1095 Double_t AliTRDv1::BetheBlochGeant(Double_t bg)
1096 {
1097   //
1098   // Return dN/dx (number of primary collisions per centimeter)
1099   // for given beta*gamma factor.
1100   //
1101   // Implemented by K.Oyama according to GEANT 3 parametrization shown in
1102   // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
1103   // This must be used as a set with IntSpecGeant.
1104   //
1105
1106   Int_t    i = 0;
1107
1108   Double_t arrG[20]  = {     1.100000,     1.200000,     1.300000,     1.500000
1109                        ,     1.800000,     2.000000,     2.500000,     3.000000
1110                        ,     4.000000,     7.000000,    10.000000,    20.000000
1111                        ,    40.000000,    70.000000,   100.000000,   300.000000
1112                        ,   600.000000,  1000.000000,  3000.000000, 10000.000000 };
1113
1114   Double_t arrNC[20] = {    75.009056,    45.508083,    35.299252,    27.116327
1115                        ,    22.734999,    21.411915,    19.934095,    19.449375
1116                        ,    19.344431,    20.185553,    21.027925,    22.912676
1117                        ,    24.933352,    26.504053,    27.387468,    29.566597
1118                        ,    30.353779,    30.787134,    31.129285,    31.157350 };
1119
1120   // Betagamma to gamma
1121   Double_t g = TMath::Sqrt(1.0 + bg*bg);
1122
1123   // Find the index just before the point we need.
1124   for (i = 0; i < 18; i++) {
1125     if ((arrG[i]   < g) && 
1126         (arrG[i+1] > g)) {
1127       break;
1128     }
1129   }
1130
1131   // Simple interpolation.
1132   Double_t pp = ((arrNC[i+1] - arrNC[i]) / (arrG[i+1]  - arrG[i])) 
1133               * (g - arrG[i]) + arrNC[i];
1134
1135   return pp;
1136
1137 }
1138
1139 //_____________________________________________________________________________
1140 Double_t Ermilova(Double_t *x, Double_t *)
1141 {
1142   //
1143   // Calculates the delta-ray energy distribution according to Ermilova.
1144   // Logarithmic scale !
1145   //
1146
1147   Double_t energy;
1148   Double_t dpos;
1149   Double_t dnde;
1150
1151   Int_t    pos1;
1152   Int_t    pos2;
1153
1154   const Int_t kNv = 31;
1155
1156   Float_t vxe[kNv] = {  2.3026,  2.9957,  3.4012,  3.6889,  3.9120  
1157                      ,  4.0943,  4.2485,  4.3820,  4.4998,  4.6052
1158                      ,  4.7005,  5.0752,  5.2983,  5.7038,  5.9915
1159                      ,  6.2146,  6.5221,  6.9078,  7.3132,  7.6009
1160                      ,  8.0064,  8.5172,  8.6995,  8.9872,  9.2103
1161                      ,  9.4727,  9.9035, 10.3735, 10.5966, 10.8198
1162                      , 11.5129 };
1163
1164   Float_t vye[kNv] = { 80.0,    31.0,    23.3,    21.1,    21.0
1165                      , 20.9,    20.8,    20.0,    16.0,    11.0
1166                      ,  8.0,     6.0,     5.2,     4.6,     4.0
1167                      ,  3.5,     3.0,     1.4,     0.67,    0.44
1168                      ,  0.3,     0.18,    0.12,    0.08,    0.056
1169                      ,  0.04,    0.023,   0.015,   0.011,    0.01
1170                      ,  0.004  };
1171
1172   energy = x[0];
1173
1174   // Find the position 
1175   pos1 = 0;
1176   pos2 = 0;
1177   dpos = 0;
1178   do {
1179     dpos = energy - vxe[pos2++];
1180   } 
1181   while (dpos > 0);
1182   pos2--; 
1183   if (pos2 > kNv) {
1184     pos2 = kNv - 1;
1185   }
1186   pos1 = pos2 - 1;
1187
1188   // Differentiate between the sampling points
1189   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
1190
1191   return dnde;
1192
1193 }
1194
1195 //_____________________________________________________________________________
1196 Double_t IntSpecGeant(Double_t *x, Double_t *)
1197 {
1198   //
1199   // Integrated spectrum from Geant3
1200   //
1201
1202   const Int_t npts = 83;
1203   Double_t arre[npts]    = {  2.421257,     2.483278,    2.534301,     2.592230
1204                            ,  2.672067,     2.813299,    3.015059,     3.216819
1205                            ,  3.418579,     3.620338,    3.868209,     3.920198
1206                            ,  3.978284,     4.063923,    4.186264,     4.308605
1207                            ,  4.430946,     4.553288,    4.724261,     4.837736
1208                            ,  4.999842,     5.161949,    5.324056,     5.486163
1209                            ,  5.679688,     5.752998,    5.857728,     5.962457
1210                            ,  6.067185,     6.171914,    6.315653,     6.393674
1211                            ,  6.471694,     6.539689,    6.597658,     6.655627
1212                            ,  6.710957,     6.763648,    6.816338,     6.876198
1213                            ,  6.943227,     7.010257,    7.106285,     7.252151
1214                            ,  7.460531,     7.668911,    7.877290,     8.085670
1215                            ,  8.302979,     8.353585,    8.413120,     8.483500
1216                            ,  8.541030,     8.592857,    8.668865,     8.820485
1217                            ,  9.037086,     9.253686,    9.470286,     9.686887
1218                            ,  9.930838,     9.994655,   10.085822,    10.176990
1219                            , 10.268158,    10.359325,   10.503614,    10.627565
1220                            , 10.804637,    10.981709,   11.158781,    11.335854
1221                            , 11.593397,    11.781165,   12.049404,    12.317644
1222                            , 12.585884,    12.854123,   14.278421,    16.975889
1223                            , 20.829416,    24.682943,   28.536469 };
1224
1225   Double_t arrdnde[npts] = { 10.960000,    10.960000,   10.359500,     9.811340
1226                            ,  9.1601500,    8.206670,    6.919630,     5.655430
1227                            ,  4.6221300,    3.777610,    3.019560,     2.591950
1228                            ,  2.5414600,    2.712920,    3.327460,     4.928240
1229                            ,  7.6185300,   10.966700,   12.225800,     8.094750
1230                            ,  3.3586900,    1.553650,    1.209600,     1.263840
1231                            ,  1.3241100,    1.312140,    1.255130,     1.165770
1232                            ,  1.0594500,    0.945450,    0.813231,     0.699837
1233                            ,  0.6235580,    2.260990,    2.968350,     2.240320
1234                            ,  1.7988300,    1.553300,    1.432070,     1.535520
1235                            ,  1.4429900,    1.247990,    1.050750,     0.829549
1236                            ,  0.5900280,    0.395897,    0.268741,     0.185320
1237                            ,  0.1292120,    0.103545,    0.0949525,    0.101535
1238                            ,  0.1276380,    0.134216,    0.123816,     0.104557
1239                            ,  0.0751843,    0.0521745,   0.0373546,    0.0275391
1240                            ,  0.0204713,    0.0169234,   0.0154552,    0.0139194
1241                            ,  0.0125592,    0.0113638,   0.0107354,    0.0102137
1242                            ,  0.00845984,   0.00683338,  0.00556836,   0.00456874
1243                            ,  0.0036227,    0.00285991,  0.00226664,   0.00172234
1244                            ,  0.00131226,   0.00100284,  0.000465492,  7.26607e-05
1245                            ,  3.63304e-06,  0.0000000,   0.0000000   };
1246
1247   Int_t    i;
1248   Double_t energy = x[0];
1249
1250   for (i = 0; i < npts; i++) {
1251     if (energy < arre[i]) {
1252       break;
1253     }
1254   }
1255
1256   if (i == 0) {
1257     AliErrorGeneral("AliTRDv1::IntSpecGeant","Given energy value is too small or zero");
1258   }
1259
1260   return arrdnde[i];
1261
1262 }