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