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