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