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