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