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