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