Update of TR simulation by Alexandru
[u/mrichter/AliRoot.git] / TRD / AliTRDv1.cxx
CommitLineData
4c039060 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
88cb7938 16/* $Id$ */
4c039060 17
fe4da5cc 18///////////////////////////////////////////////////////////////////////////////
19// //
769257f4 20// Transition Radiation Detector version 1 -- slow simulator //
fe4da5cc 21// //
22//Begin_Html
23/*
5c7f4665 24<img src="picts/AliTRDfullClass.gif">
fe4da5cc 25*/
26//End_Html
27// //
28// //
29///////////////////////////////////////////////////////////////////////////////
30
769257f4 31#include <stdlib.h>
32
793ff80c 33#include <TF1.h>
1819f4bb 34#include <TLorentzVector.h>
88cb7938 35#include <TMath.h>
36#include <TRandom.h>
37#include <TVector.h>
38#include <TVirtualMC.h>
fe4da5cc 39
d3f347ff 40#include "AliConst.h"
45160b1f 41#include "AliLog.h"
42#include "AliMC.h"
88cb7938 43#include "AliRun.h"
44#include "AliTRDgeometry.h"
793ff80c 45#include "AliTRDhit.h"
793ff80c 46#include "AliTRDsim.h"
88cb7938 47#include "AliTRDv1.h"
851d3db9 48
fe4da5cc 49ClassImp(AliTRDv1)
8230f242 50
51//_____________________________________________________________________________
52AliTRDv1::AliTRDv1():AliTRD()
53{
54 //
55 // Default constructor
56 //
57
a328fff9 58 fSensSelect = 0;
59 fSensPlane = -1;
60 fSensChamber = -1;
61 fSensSector = -1;
62 fSensSectorRange = 0;
8230f242 63
a328fff9 64 fDeltaE = NULL;
65 fDeltaG = NULL;
66 fTR = NULL;
67
68 fStepSize = 0.1;
69 fTypeOfStepManager = 2;
8230f242 70
71}
72
fe4da5cc 73//_____________________________________________________________________________
74AliTRDv1::AliTRDv1(const char *name, const char *title)
75 :AliTRD(name, title)
76{
77 //
851d3db9 78 // Standard constructor for Transition Radiation Detector version 1
fe4da5cc 79 //
82bbf98a 80
a328fff9 81 fSensSelect = 0;
82 fSensPlane = -1;
83 fSensChamber = -1;
84 fSensSector = -1;
85 fSensSectorRange = 0;
5c7f4665 86
a328fff9 87 fDeltaE = NULL;
88 fDeltaG = NULL;
89 fTR = NULL;
90 fStepSize = 0.1;
91 fTypeOfStepManager = 2;
5c7f4665 92
93 SetBufferSize(128000);
94
95}
96
97//_____________________________________________________________________________
73ae7b59 98AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd)
8230f242 99{
100 //
101 // Copy constructor
102 //
103
dd9a6ee3 104 ((AliTRDv1 &) trd).Copy(*this);
8230f242 105
106}
107
108//_____________________________________________________________________________
5c7f4665 109AliTRDv1::~AliTRDv1()
110{
dd9a6ee3 111 //
112 // AliTRDv1 destructor
113 //
82bbf98a 114
5c7f4665 115 if (fDeltaE) delete fDeltaE;
a328fff9 116 if (fDeltaG) delete fDeltaG;
793ff80c 117 if (fTR) delete fTR;
82bbf98a 118
fe4da5cc 119}
120
dd9a6ee3 121//_____________________________________________________________________________
122AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
123{
124 //
125 // Assignment operator
126 //
127
128 if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
129 return *this;
130
131}
8230f242 132
133//_____________________________________________________________________________
e0d47c25 134void AliTRDv1::Copy(TObject &trd) const
8230f242 135{
136 //
137 // Copy function
138 //
139
a328fff9 140 ((AliTRDv1 &) trd).fSensSelect = fSensSelect;
141 ((AliTRDv1 &) trd).fSensPlane = fSensPlane;
142 ((AliTRDv1 &) trd).fSensChamber = fSensChamber;
143 ((AliTRDv1 &) trd).fSensSector = fSensSector;
144 ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange;
145
146 ((AliTRDv1 &) trd).fTypeOfStepManager = fTypeOfStepManager;
147 ((AliTRDv1 &) trd).fStepSize = fStepSize;
8230f242 148
793ff80c 149 fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
a328fff9 150 fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG);
793ff80c 151 fTR->Copy(*((AliTRDv1 &) trd).fTR);
8230f242 152
153}
154
fe4da5cc 155//_____________________________________________________________________________
156void AliTRDv1::CreateGeometry()
157{
158 //
851d3db9 159 // Create the GEANT geometry for the Transition Radiation Detector - Version 1
5c7f4665 160 // This version covers the full azimuth.
d3f347ff 161 //
162
82bbf98a 163 // Check that FRAME is there otherwise we have no place where to put the TRD
8230f242 164 AliModule* frame = gAlice->GetModule("FRAME");
165 if (!frame) return;
d3f347ff 166
82bbf98a 167 // Define the chambers
168 AliTRD::CreateGeometry();
d3f347ff 169
fe4da5cc 170}
171
172//_____________________________________________________________________________
173void AliTRDv1::CreateMaterials()
174{
175 //
851d3db9 176 // Create materials for the Transition Radiation Detector version 1
fe4da5cc 177 //
82bbf98a 178
d3f347ff 179 AliTRD::CreateMaterials();
82bbf98a 180
fe4da5cc 181}
182
183//_____________________________________________________________________________
793ff80c 184void AliTRDv1::CreateTRhit(Int_t det)
185{
186 //
187 // Creates an electron cluster from a TR photon.
188 // The photon is assumed to be created a the end of the radiator. The
189 // distance after which it deposits its energy takes into account the
190 // absorbtion of the entrance window and of the gas mixture in drift
191 // volume.
192 //
193
194 // PDG code electron
195 const Int_t kPdgElectron = 11;
196
197 // Ionization energy
198 const Float_t kWion = 22.04;
199
200 // Maximum number of TR photons per track
201 const Int_t kNTR = 50;
202
203 TLorentzVector mom, pos;
793ff80c 204
793ff80c 205 // Create TR at the entrance of the chamber
206 if (gMC->IsTrackEntering()) {
207
f73816f5 208 // Create TR only for electrons
209 Int_t iPdg = gMC->TrackPid();
210 if (TMath::Abs(iPdg) != kPdgElectron) return;
211
793ff80c 212 Float_t eTR[kNTR];
213 Int_t nTR;
214
215 // Create TR photons
216 gMC->TrackMomentum(mom);
217 Float_t pTot = mom.Rho();
218 fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
219 if (nTR > kNTR) {
45160b1f 220 AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
793ff80c 221 }
222
223 // Loop through the TR photons
224 for (Int_t iTR = 0; iTR < nTR; iTR++) {
225
226 Float_t energyMeV = eTR[iTR] * 0.001;
227 Float_t energyeV = eTR[iTR] * 1000.0;
228 Float_t absLength = 0;
229 Float_t sigma = 0;
230
231 // Take the absorbtion in the entrance window into account
232 Double_t muMy = fTR->GetMuMy(energyMeV);
233 sigma = muMy * fFoilDensity;
842287f2 234 if (sigma > 0.0) {
235 absLength = gRandom->Exp(1.0/sigma);
236 if (absLength < AliTRDgeometry::MyThick()) continue;
237 }
238 else {
239 continue;
240 }
793ff80c 241
242 // The absorbtion cross sections in the drift gas
3dac2b2d 243 // Gas-mixture (Xe/CO2)
244 Double_t muXe = fTR->GetMuXe(energyMeV);
245 Double_t muCO = fTR->GetMuCO(energyMeV);
246 sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
793ff80c 247
248 // The distance after which the energy of the TR photon
249 // is deposited.
842287f2 250 if (sigma > 0.0) {
251 absLength = gRandom->Exp(1.0/sigma);
a328fff9 252 if (absLength > (AliTRDgeometry::DrThick()
253 + AliTRDgeometry::AmThick())) {
254 continue;
255 }
842287f2 256 }
257 else {
258 continue;
259 }
793ff80c 260
261 // The position of the absorbtion
262 Float_t posHit[3];
263 gMC->TrackPosition(pos);
264 posHit[0] = pos[0] + mom[0] / pTot * absLength;
265 posHit[1] = pos[1] + mom[1] / pTot * absLength;
266 posHit[2] = pos[2] + mom[2] / pTot * absLength;
267
268 // Create the charge
269 Int_t q = ((Int_t) (energyeV / kWion));
270
271 // Add the hit to the array. TR photon hits are marked
272 // by negative charge
5d12ce38 273 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE);
793ff80c 274
275 }
276
277 }
278
279}
280
281//_____________________________________________________________________________
5c7f4665 282void AliTRDv1::Init()
283{
284 //
285 // Initialise Transition Radiation Detector after geometry has been built.
5c7f4665 286 //
287
288 AliTRD::Init();
289
45160b1f 290 AliDebug(1,"Slow simulator\n");
851d3db9 291 if (fSensSelect) {
292 if (fSensPlane >= 0)
45160b1f 293 AliInfo(Form("Only plane %d is sensitive"));
851d3db9 294 if (fSensChamber >= 0)
45160b1f 295 AliInfo(Form("Only chamber %d is sensitive",fSensChamber));
9d0b222b 296 if (fSensSector >= 0) {
297 Int_t sens1 = fSensSector;
298 Int_t sens2 = fSensSector + fSensSectorRange;
793ff80c 299 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
300 * AliTRDgeometry::Nsect();
45160b1f 301 AliInfo(Form("Only sectors %d - %d are sensitive\n",sens1,sens2-1));
9d0b222b 302 }
851d3db9 303 }
793ff80c 304 if (fTR)
45160b1f 305 AliInfo("TR simulation on")
793ff80c 306 else
45160b1f 307 AliInfo("TR simulation off");
5c7f4665 308
309 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
310 const Float_t kPoti = 12.1;
311 // Maximum energy (50 keV);
312 const Float_t kEend = 50000.0;
313 // Ermilova distribution for the delta-ray spectrum
8230f242 314 Float_t poti = TMath::Log(kPoti);
315 Float_t eEnd = TMath::Log(kEend);
a328fff9 316
317 // Ermilova distribution for the delta-ray spectrum
318 fDeltaE = new TF1("deltae" ,Ermilova ,poti,eEnd,0);
319
320 // Geant3 distribution for the delta-ray spectrum
321 fDeltaG = new TF1("deltaeg",IntSpecGeant,poti,eEnd,0);
5c7f4665 322
45160b1f 323 AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
5c7f4665 324
fe4da5cc 325}
326
327//_____________________________________________________________________________
793ff80c 328AliTRDsim *AliTRDv1::CreateTR()
329{
330 //
331 // Enables the simulation of TR
332 //
333
334 fTR = new AliTRDsim();
335 return fTR;
336
337}
338
339//_____________________________________________________________________________
5c7f4665 340void AliTRDv1::SetSensPlane(Int_t iplane)
341{
342 //
851d3db9 343 // Defines the hit-sensitive plane (0-5)
5c7f4665 344 //
82bbf98a 345
851d3db9 346 if ((iplane < 0) || (iplane > 5)) {
45160b1f 347 AliWarning(Form("Wrong input value:%d",iplane));
348 AliWarning("Use standard setting");
851d3db9 349 fSensPlane = -1;
350 fSensSelect = 0;
5c7f4665 351 return;
352 }
82bbf98a 353
5c7f4665 354 fSensSelect = 1;
355 fSensPlane = iplane;
82bbf98a 356
5c7f4665 357}
358
359//_____________________________________________________________________________
360void AliTRDv1::SetSensChamber(Int_t ichamber)
361{
362 //
851d3db9 363 // Defines the hit-sensitive chamber (0-4)
5c7f4665 364 //
365
851d3db9 366 if ((ichamber < 0) || (ichamber > 4)) {
45160b1f 367 AliWarning(Form("Wrong input value: %d",ichamber));
368 AliWarning("Use standard setting");
851d3db9 369 fSensChamber = -1;
370 fSensSelect = 0;
5c7f4665 371 return;
372 }
373
374 fSensSelect = 1;
375 fSensChamber = ichamber;
376
377}
378
379//_____________________________________________________________________________
380void AliTRDv1::SetSensSector(Int_t isector)
381{
382 //
851d3db9 383 // Defines the hit-sensitive sector (0-17)
5c7f4665 384 //
385
9d0b222b 386 SetSensSector(isector,1);
387
388}
389
390//_____________________________________________________________________________
391void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
392{
393 //
394 // Defines a range of hit-sensitive sectors. The range is defined by
395 // <isector> (0-17) as the starting point and <nsector> as the number
396 // of sectors to be included.
397 //
398
851d3db9 399 if ((isector < 0) || (isector > 17)) {
45160b1f 400 AliWarning(Form("Wrong input value <isector>: %d",isector));
401 AliWarning("Use standard setting");
9d0b222b 402 fSensSector = -1;
403 fSensSectorRange = 0;
404 fSensSelect = 0;
5c7f4665 405 return;
406 }
407
9d0b222b 408 if ((nsector < 1) || (nsector > 18)) {
45160b1f 409 AliWarning(Form("Wrong input value <nsector>: %d",nsector));
410 AliWarning("Use standard setting");
9d0b222b 411 fSensSector = -1;
412 fSensSectorRange = 0;
413 fSensSelect = 0;
414 return;
415 }
416
417 fSensSelect = 1;
418 fSensSector = isector;
419 fSensSectorRange = nsector;
5c7f4665 420
421}
422
423//_____________________________________________________________________________
424void AliTRDv1::StepManager()
425{
426 //
5c7f4665 427 // Slow simulator. Every charged track produces electron cluster as hits
a328fff9 428 // along its path across the drift volume.
429 //
430
431 switch (fTypeOfStepManager) {
432 case 0 : StepManagerErmilova(); break; // 0 is Ermilova
433 case 1 : StepManagerGeant(); break; // 1 is Geant
434 case 2 : StepManagerFixedStep(); break; // 2 is fixed step
45160b1f 435 default : AliWarning("Not a valid Step Manager.");
a328fff9 436 }
437
438}
439
440//_____________________________________________________________________________
441void AliTRDv1::SelectStepManager(Int_t t)
442{
443 //
444 // Selects a step manager type:
445 // 0 - Ermilova
446 // 1 - Geant3
447 // 2 - Fixed step size
448 //
449
450 if (t == 1) {
45160b1f 451 AliWarning("Sorry, Geant parametrization step manager is not implemented yet. Please ask K.Oyama for detail.");
a328fff9 452 }
453
454 fTypeOfStepManager = t;
45160b1f 455 AliInfo(Form("Step Manager type %d was selected",fTypeOfStepManager));
a328fff9 456
457}
458
459//_____________________________________________________________________________
460void AliTRDv1::StepManagerGeant()
461{
462 //
463 // Slow simulator. Every charged track produces electron cluster as hits
464 // along its path across the drift volume. The step size is set acording
465 // to Bethe-Bloch. The energy distribution of the delta electrons follows
466 // a spectrum taken from Geant3.
467 //
468
45160b1f 469 AliWarning("Not implemented yet.");
a328fff9 470
471}
472
473//_____________________________________________________________________________
474void AliTRDv1::StepManagerErmilova()
475{
476 //
477 // Slow simulator. Every charged track produces electron cluster as hits
5c7f4665 478 // along its path across the drift volume. The step size is set acording
479 // to Bethe-Bloch. The energy distribution of the delta electrons follows
480 // a spectrum taken from Ermilova et al.
481 //
482
851d3db9 483 Int_t pla = 0;
484 Int_t cha = 0;
485 Int_t sec = 0;
793ff80c 486 Int_t det = 0;
851d3db9 487 Int_t iPdg;
793ff80c 488 Int_t qTot;
5c7f4665 489
793ff80c 490 Float_t hits[3];
a5cadd36 491 Double_t random[1];
5c7f4665 492 Float_t charge;
493 Float_t aMass;
494
f73816f5 495 Double_t pTot = 0;
5c7f4665 496 Double_t eDelta;
497 Double_t betaGamma, pp;
f73816f5 498 Double_t stepSize;
5c7f4665 499
332e9569 500 Bool_t drRegion = kFALSE;
501 Bool_t amRegion = kFALSE;
502
503 TString cIdCurrent;
504 TString cIdSensDr = "J";
505 TString cIdSensAm = "K";
593a9fc3 506 Char_t cIdChamber[3];
507 cIdChamber[2] = 0;
332e9569 508
5c7f4665 509 TLorentzVector pos, mom;
82bbf98a 510
332e9569 511 const Int_t kNplan = AliTRDgeometry::Nplan();
e644678a 512 const Int_t kNcham = AliTRDgeometry::Ncham();
513 const Int_t kNdetsec = kNplan * kNcham;
514
a328fff9 515 const Double_t kBig = 1.0E+12; // Infinitely big
516 const Float_t kWion = 22.04; // Ionization energy
517 const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
5c7f4665 518
f73816f5 519 // Minimum energy for the step size adjustment
520 const Float_t kEkinMinStep = 1.0e-5;
a328fff9 521
5c7f4665 522 // Plateau value of the energy-loss for electron in xenon
523 // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
524 //const Double_t kPlateau = 1.70;
525 // the averaged value (26/3/99)
a3c76cdc 526 const Float_t kPlateau = 1.55;
a328fff9 527
528 const Float_t kPrim = 48.0; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
5c7f4665 529 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
a3c76cdc 530 const Float_t kPoti = 12.1;
851d3db9 531
a328fff9 532 const Int_t kPdgElectron = 11; // PDG code electron
5c7f4665 533
534 // Set the maximum step size to a very large number for all
535 // neutral particles and those outside the driftvolume
536 gMC->SetMaxStep(kBig);
537
538 // Use only charged tracks
539 if (( gMC->TrackCharge() ) &&
540 (!gMC->IsTrackStop() ) &&
541 (!gMC->IsTrackDisappeared())) {
fe4da5cc 542
5c7f4665 543 // Inside a sensitive volume?
332e9569 544 drRegion = kFALSE;
545 amRegion = kFALSE;
546 cIdCurrent = gMC->CurrentVolName();
e6674585 547 if (cIdSensDr == cIdCurrent[1]) {
332e9569 548 drRegion = kTRUE;
549 }
e6674585 550 if (cIdSensAm == cIdCurrent[1]) {
332e9569 551 amRegion = kTRUE;
552 }
553 if (drRegion || amRegion) {
fe4da5cc 554
5c7f4665 555 // The hit coordinates and charge
556 gMC->TrackPosition(pos);
557 hits[0] = pos[0];
558 hits[1] = pos[1];
559 hits[2] = pos[2];
5c7f4665 560
851d3db9 561 // The sector number (0 - 17)
562 // The numbering goes clockwise and starts at y = 0
e15eb584 563 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
851d3db9 564 if (phi < 90.)
565 phi = phi + 270.;
566 else
567 phi = phi - 90.;
568 sec = ((Int_t) (phi / 20));
82bbf98a 569
332e9569 570 // The plane and chamber number
571 cIdChamber[0] = cIdCurrent[2];
572 cIdChamber[1] = cIdCurrent[3];
e644678a 573 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
a5cadd36 574 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
332e9569 575 pla = ((Int_t) idChamber % kNplan);
82bbf98a 576
5c7f4665 577 // Check on selected volumes
578 Int_t addthishit = 1;
579 if (fSensSelect) {
6f1e466d 580 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
581 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
9d0b222b 582 if (fSensSector >= 0) {
583 Int_t sens1 = fSensSector;
584 Int_t sens2 = fSensSector + fSensSectorRange;
793ff80c 585 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
586 * AliTRDgeometry::Nsect();
9d0b222b 587 if (sens1 < sens2) {
588 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
589 }
590 else {
591 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
592 }
593 }
5c7f4665 594 }
595
596 // Add this hit
597 if (addthishit) {
598
f73816f5 599 // The detector number
793ff80c 600 det = fGeometry->GetDetector(pla,cha,sec);
601
a328fff9 602 // Special hits only in the drift region
332e9569 603 if (drRegion) {
f73816f5 604
c61f1a66 605 // Create a track reference at the entrance and
606 // exit of each chamber that contain the
607 // momentum components of the particle
f73816f5 608 if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
609 gMC->TrackMomentum(mom);
5d12ce38 610 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
f73816f5 611 }
612
613 // Create the hits from TR photons
614 if (fTR) CreateTRhit(det);
615
616 }
617
618 // Calculate the energy of the delta-electrons
619 eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
620 eDelta = TMath::Max(eDelta,0.0);
621
622 // The number of secondary electrons created
623 qTot = ((Int_t) (eDelta / kWion) + 1);
624
625 // Create a new dEdx hit
332e9569 626 if (drRegion) {
a328fff9 627 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
628 ,det,hits,qTot,kTRUE);
f73816f5 629 }
5c7f4665 630 else {
a328fff9 631 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
632 ,det,hits,qTot,kFALSE);
f73816f5 633 }
634
5c7f4665 635 // Calculate the maximum step size for the next tracking step
f73816f5 636 // Produce only one hit if Ekin is below cutoff
637 aMass = gMC->TrackMass();
638 if ((gMC->Etot() - aMass) > kEkinMinStep) {
639
640 // The energy loss according to Bethe Bloch
641 iPdg = TMath::Abs(gMC->TrackPid());
642 if ( (iPdg != kPdgElectron) ||
643 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
644 gMC->TrackMomentum(mom);
645 pTot = mom.Rho();
646 betaGamma = pTot / aMass;
647 pp = kPrim * BetheBloch(betaGamma);
648 // Take charge > 1 into account
649 charge = gMC->TrackCharge();
650 if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
651 }
652 // Electrons above 20 Mev/c are at the plateau
653 else {
654 pp = kPrim * kPlateau;
655 }
656
657 if (pp > 0) {
658 do
b9d0a01d 659 gMC->GetRandom()->RndmArray(1, random);
f73816f5 660 while ((random[0] == 1.) || (random[0] == 0.));
661 stepSize = - TMath::Log(random[0]) / pp;
662 gMC->SetMaxStep(stepSize);
663 }
664
5c7f4665 665 }
666
667 }
d3f347ff 668
669 }
670
5c7f4665 671 }
672
673}
674
675//_____________________________________________________________________________
a328fff9 676void AliTRDv1::StepManagerFixedStep()
677{
678 //
679 // Slow simulator. Every charged track produces electron cluster as hits
680 // along its path across the drift volume. The step size is fixed in
681 // this version of the step manager.
682 //
683
684 Int_t pla = 0;
685 Int_t cha = 0;
686 Int_t sec = 0;
687 Int_t det = 0;
688 Int_t qTot;
689
690 Float_t hits[3];
691 Double_t eDep;
692
693 Bool_t drRegion = kFALSE;
694 Bool_t amRegion = kFALSE;
695
696 TString cIdCurrent;
697 TString cIdSensDr = "J";
698 TString cIdSensAm = "K";
699 Char_t cIdChamber[3];
700 cIdChamber[2] = 0;
701
702 TLorentzVector pos, mom;
703
704 const Int_t kNplan = AliTRDgeometry::Nplan();
705 const Int_t kNcham = AliTRDgeometry::Ncham();
706 const Int_t kNdetsec = kNplan * kNcham;
707
708 const Double_t kBig = 1.0E+12;
709
710 const Float_t kWion = 22.04; // Ionization energy
711 const Float_t kEkinMinStep = 1.0e-5; // Minimum energy for the step size adjustment
712
713 // Set the maximum step size to a very large number for all
714 // neutral particles and those outside the driftvolume
715 gMC->SetMaxStep(kBig);
716
717 // If not charged track or already stopped or disappeared, just return.
718 if ((!gMC->TrackCharge()) ||
719 gMC->IsTrackStop() ||
720 gMC->IsTrackDisappeared()) return;
721
722 // Inside a sensitive volume?
723 cIdCurrent = gMC->CurrentVolName();
724
725 if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
726 if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
727
728 if ((!drRegion) && (!amRegion)) return;
729
730 // The hit coordinates and charge
731 gMC->TrackPosition(pos);
732 hits[0] = pos[0];
733 hits[1] = pos[1];
734 hits[2] = pos[2];
735
736 // The sector number (0 - 17)
737 // The numbering goes clockwise and starts at y = 0
738 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
739 if (phi < 90.) phi += 270.;
740 else phi -= 90.;
741 sec = ((Int_t) (phi / 20.));
742
743 // The plane and chamber number
744 cIdChamber[0] = cIdCurrent[2];
745 cIdChamber[1] = cIdCurrent[3];
746 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
a5cadd36 747 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
a328fff9 748 pla = ((Int_t) idChamber % kNplan);
e0d47c25 749
a328fff9 750 // Check on selected volumes
751 Int_t addthishit = 1;
752 if(fSensSelect) {
753 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
754 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
755 if (fSensSector >= 0) {
756 Int_t sens1 = fSensSector;
757 Int_t sens2 = fSensSector + fSensSectorRange;
758 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) * AliTRDgeometry::Nsect();
759 if (sens1 < sens2) {
760 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
761 }
762 else {
763 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
764 }
765 }
766 }
767
768 if (!addthishit) return;
769
770 det = fGeometry->GetDetector(pla,cha,sec); // The detector number
771
772 Int_t trkStat = 0; // 0: InFlight 1:Entering 2:Exiting
773
774 // Special hits only in the drift region
775 if (drRegion) {
776
777 // Create a track reference at the entrance and exit of each
778 // chamber that contain the momentum components of the particle
779
780 if (gMC->IsTrackEntering()) {
781 gMC->TrackMomentum(mom);
782 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
783 trkStat = 1;
784 }
785 if (gMC->IsTrackExiting()) {
786 gMC->TrackMomentum(mom);
787 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
788 trkStat = 2;
789 }
790
791 // Create the hits from TR photons
792 if (fTR) CreateTRhit(det);
793
794 }
795
796 // Calculate the charge according to GEANT Edep
797 // Create a new dEdx hit
798 eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
799 qTot = (Int_t) (eDep / kWion);
800 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot,drRegion);
801
802 // Set Maximum Step Size
803 // Produce only one hit if Ekin is below cutoff
804 if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return;
805 gMC->SetMaxStep(fStepSize);
806
807}
808
809//_____________________________________________________________________________
5c7f4665 810Double_t AliTRDv1::BetheBloch(Double_t bg)
811{
812 //
813 // Parametrization of the Bethe-Bloch-curve
814 // The parametrization is the same as for the TPC and is taken from Lehrhaus.
815 //
816
817 // This parameters have been adjusted to averaged values from GEANT
818 const Double_t kP1 = 7.17960e-02;
819 const Double_t kP2 = 8.54196;
820 const Double_t kP3 = 1.38065e-06;
821 const Double_t kP4 = 5.30972;
822 const Double_t kP5 = 2.83798;
823
824 // This parameters have been adjusted to Xe-data found in:
825 // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
826 //const Double_t kP1 = 0.76176E-1;
827 //const Double_t kP2 = 10.632;
828 //const Double_t kP3 = 3.17983E-6;
829 //const Double_t kP4 = 1.8631;
830 //const Double_t kP5 = 1.9479;
831
f73816f5 832 // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
833 const Double_t kBgMin = 0.8;
834 const Double_t kBBMax = 6.83298;
835 //const Double_t kBgMin = 0.6;
836 //const Double_t kBBMax = 17.2809;
837 //const Double_t kBgMin = 0.4;
838 //const Double_t kBBMax = 82.0;
839
840 if (bg > kBgMin) {
5c7f4665 841 Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
842 Double_t aa = TMath::Power(yy,kP4);
843 Double_t bb = TMath::Power((1./bg),kP5);
844 bb = TMath::Log(kP3 + bb);
845 return ((kP2 - aa - bb)*kP1 / aa);
846 }
f73816f5 847 else {
848 return kBBMax;
849 }
d3f347ff 850
fe4da5cc 851}
5c7f4665 852
853//_____________________________________________________________________________
a328fff9 854Double_t BetheBlochGeant(Double_t bg)
855{
856 //
857 // Return dN/dx (number of primary collisions per centimeter)
858 // for given beta*gamma factor.
859 //
860 // Implemented by K.Oyama according to GEANT 3 parametrization shown in
861 // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
862 // This must be used as a set with IntSpecGeant.
863 //
864
865 Double_t arr_g[20] = {
866 1.100000, 1.200000, 1.300000, 1.500000,
867 1.800000, 2.000000, 2.500000, 3.000000,
868 4.000000, 7.000000, 10.000000, 20.000000,
869 40.000000, 70.000000, 100.000000, 300.000000,
870 600.000000, 1000.000000, 3000.000000, 10000.000000 };
871
872 Double_t arr_nc[20] = {
873 75.009056, 45.508083, 35.299252, 27.116327,
874 22.734999, 21.411915, 19.934095, 19.449375,
875 19.344431, 20.185553, 21.027925, 22.912676,
876 24.933352, 26.504053, 27.387468, 29.566597,
877 30.353779, 30.787134, 31.129285, 31.157350 };
878
879 // betagamma to gamma
880 Double_t g = TMath::Sqrt( 1. + bg*bg );
881
882 // Find the index just before the point we need.
883 int i;
884 for( i = 0 ; i < 18 ; i++ )
885 if( arr_g[i] < g && arr_g[i+1] > g )
886 break;
887
888 // Simple interpolation.
889 Double_t pp = ((arr_nc[i+1] - arr_nc[i]) /
890 (arr_g[i+1]-arr_g[i])) * (g-arr_g[i]) + arr_nc[i];
891
892 return pp;
893
894}
895
896//_____________________________________________________________________________
5c7f4665 897Double_t Ermilova(Double_t *x, Double_t *)
898{
899 //
900 // Calculates the delta-ray energy distribution according to Ermilova.
901 // Logarithmic scale !
902 //
903
904 Double_t energy;
905 Double_t dpos;
906 Double_t dnde;
907
908 Int_t pos1, pos2;
909
8230f242 910 const Int_t kNv = 31;
5c7f4665 911
8230f242 912 Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
913 , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
914 , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
915 , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
916 , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
917 , 9.4727, 9.9035,10.3735,10.5966,10.8198
918 ,11.5129 };
5c7f4665 919
8230f242 920 Float_t vye[kNv] = { 80.0 , 31.0 , 23.3 , 21.1 , 21.0
921 , 20.9 , 20.8 , 20.0 , 16.0 , 11.0
922 , 8.0 , 6.0 , 5.2 , 4.6 , 4.0
923 , 3.5 , 3.0 , 1.4 , 0.67 , 0.44
924 , 0.3 , 0.18 , 0.12 , 0.08 , 0.056
925 , 0.04 , 0.023, 0.015, 0.011, 0.01
926 , 0.004 };
5c7f4665 927
928 energy = x[0];
929
930 // Find the position
931 pos1 = pos2 = 0;
932 dpos = 0;
933 do {
934 dpos = energy - vxe[pos2++];
935 }
936 while (dpos > 0);
937 pos2--;
f73816f5 938 if (pos2 > kNv) pos2 = kNv - 1;
5c7f4665 939 pos1 = pos2 - 1;
940
941 // Differentiate between the sampling points
942 dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
943
944 return dnde;
945
946}
a328fff9 947
948//_____________________________________________________________________________
949Double_t IntSpecGeant(Double_t *x, Double_t *)
950{
951 //
952 // Integrated spectrum from Geant3
953 //
954
955 const Int_t n_pts = 83;
956 Double_t arr_e[n_pts] = {
957 2.421257, 2.483278, 2.534301, 2.592230,
958 2.672067, 2.813299, 3.015059, 3.216819,
959 3.418579, 3.620338, 3.868209, 3.920198,
960 3.978284, 4.063923, 4.186264, 4.308605,
961 4.430946, 4.553288, 4.724261, 4.837736,
962 4.999842, 5.161949, 5.324056, 5.486163,
963 5.679688, 5.752998, 5.857728, 5.962457,
964 6.067185, 6.171914, 6.315653, 6.393674,
965 6.471694, 6.539689, 6.597658, 6.655627,
966 6.710957, 6.763648, 6.816338, 6.876198,
967 6.943227, 7.010257, 7.106285, 7.252151,
968 7.460531, 7.668911, 7.877290, 8.085670,
969 8.302979, 8.353585, 8.413120, 8.483500,
970 8.541030, 8.592857, 8.668865, 8.820485,
971 9.037086, 9.253686, 9.470286, 9.686887,
972 9.930838, 9.994655, 10.085822, 10.176990,
973 10.268158, 10.359325, 10.503614, 10.627565,
974 10.804637, 10.981709, 11.158781, 11.335854,
975 11.593397, 11.781165, 12.049404, 12.317644,
976 12.585884, 12.854123, 14.278421, 16.975889,
977 20.829416, 24.682943, 28.536469
978 };
979 Double_t arr_dndx[n_pts] = {
980 19.344431, 18.664679, 18.136106, 17.567745,
981 16.836426, 15.677382, 14.281277, 13.140237,
982 12.207677, 11.445510, 10.697049, 10.562296,
983 10.414673, 10.182341, 9.775256, 9.172330,
984 8.240271, 6.898587, 4.808303, 3.889751,
985 3.345288, 3.093431, 2.897347, 2.692470,
986 2.436222, 2.340029, 2.208579, 2.086489,
987 1.975535, 1.876519, 1.759626, 1.705024,
988 1.656374, 1.502638, 1.330566, 1.200697,
989 1.101168, 1.019323, 0.943867, 0.851951,
990 0.755229, 0.671576, 0.570675, 0.449672,
991 0.326722, 0.244225, 0.188225, 0.149608,
992 0.121529, 0.116289, 0.110636, 0.103490,
993 0.096147, 0.089191, 0.079780, 0.063927,
994 0.047642, 0.036341, 0.028250, 0.022285,
995 0.017291, 0.016211, 0.014802, 0.013533,
996 0.012388, 0.011352, 0.009803, 0.008537,
997 0.007039, 0.005829, 0.004843, 0.004034,
998 0.003101, 0.002564, 0.001956, 0.001494,
999 0.001142, 0.000873, 0.000210, 0.000014,
1000 0.000000, 0.000000, 0.000000
1001 };
1002
1003 Int_t i;
1004 Double_t energy = x[0];
1005 Double_t dnde;
1006
1007 for( i = 0 ; i < n_pts ; i++ )
1008 if( energy < arr_e[i] ) break;
1009
1010 if( i == 0 )
45160b1f 1011 AliErrorGeneral("AliTRDv1","Given energy value is too small or zero");
a328fff9 1012
1013 // Interpolate
1014 dnde = (arr_dndx[i-1] - arr_dndx[i]) / (arr_e[i] - arr_e[i-1]);
1015
1016 return dnde;
1017
1018}