1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // Transition Radiation Detector version 1 -- slow simulator //
24 <img src="picts/AliTRDfullClass.gif">
29 ///////////////////////////////////////////////////////////////////////////////
34 #include <TLorentzVector.h>
38 #include <TVirtualMC.h>
44 #include "AliTRDgeometry.h"
45 #include "AliTRDhit.h"
46 #include "AliTRDsim.h"
51 //_____________________________________________________________________________
52 AliTRDv1::AliTRDv1():AliTRD()
55 // Default constructor
70 fTypeOfStepManager = 1;
74 //_____________________________________________________________________________
75 AliTRDv1::AliTRDv1(const char *name, const char *title)
79 // Standard constructor for Transition Radiation Detector version 1
93 fTypeOfStepManager = 1;
95 SetBufferSize(128000);
99 //_____________________________________________________________________________
100 AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd)
106 ((AliTRDv1 &) trd).Copy(*this);
110 //_____________________________________________________________________________
111 AliTRDv1::~AliTRDv1()
114 // AliTRDv1 destructor
117 if (fDeltaE) delete fDeltaE;
118 if (fDeltaG) delete fDeltaG;
123 //_____________________________________________________________________________
124 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
127 // Assignment operator
130 if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
135 //_____________________________________________________________________________
136 void AliTRDv1::Copy(TObject &trd) const
138 printf("void AliTRDv1::Copy(TObject &trd) const\n");
143 ((AliTRDv1 &) trd).fSensSelect = fSensSelect;
144 ((AliTRDv1 &) trd).fSensPlane = fSensPlane;
145 ((AliTRDv1 &) trd).fSensChamber = fSensChamber;
146 ((AliTRDv1 &) trd).fSensSector = fSensSector;
147 ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange;
149 ((AliTRDv1 &) trd).fTypeOfStepManager = fTypeOfStepManager;
150 ((AliTRDv1 &) trd).fStepSize = fStepSize;
152 ((AliTRDv1 &) trd).fTRon = fTRon;
154 fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
155 fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG);
156 fTR->Copy(*((AliTRDv1 &) trd).fTR);
160 //_____________________________________________________________________________
161 void AliTRDv1::CreateGeometry()
164 // Create the GEANT geometry for the Transition Radiation Detector - Version 1
165 // This version covers the full azimuth.
168 // Check that FRAME is there otherwise we have no place where to put the TRD
169 AliModule* frame = gAlice->GetModule("FRAME");
172 // Define the chambers
173 AliTRD::CreateGeometry();
177 //_____________________________________________________________________________
178 void AliTRDv1::CreateMaterials()
181 // Create materials for the Transition Radiation Detector version 1
184 AliTRD::CreateMaterials();
188 //_____________________________________________________________________________
189 void AliTRDv1::CreateTRhit(Int_t det)
192 // Creates an electron cluster from a TR photon.
193 // The photon is assumed to be created a the end of the radiator. The
194 // distance after which it deposits its energy takes into account the
195 // absorbtion of the entrance window and of the gas mixture in drift
200 const Int_t kPdgElectron = 11;
203 const Float_t kWion = 23.53;
205 // Maximum number of TR photons per track
206 const Int_t kNTR = 50;
208 TLorentzVector mom, pos;
210 // Create TR at the entrance of the chamber
211 if (gMC->IsTrackEntering()) {
213 // Create TR only for electrons
214 Int_t iPdg = gMC->TrackPid();
215 if (TMath::Abs(iPdg) != kPdgElectron) return;
221 gMC->TrackMomentum(mom);
222 Float_t pTot = mom.Rho();
223 fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
225 AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
228 // Loop through the TR photons
229 for (Int_t iTR = 0; iTR < nTR; iTR++) {
231 Float_t energyMeV = eTR[iTR] * 0.001;
232 Float_t energyeV = eTR[iTR] * 1000.0;
233 Float_t absLength = 0;
236 // Take the absorbtion in the entrance window into account
237 Double_t muMy = fTR->GetMuMy(energyMeV);
238 sigma = muMy * fFoilDensity;
240 absLength = gRandom->Exp(1.0/sigma);
241 if (absLength < AliTRDgeometry::MyThick()) continue;
247 // The absorbtion cross sections in the drift gas
248 // Gas-mixture (Xe/CO2)
249 Double_t muXe = fTR->GetMuXe(energyMeV);
250 Double_t muCO = fTR->GetMuCO(energyMeV);
251 sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
253 // The distance after which the energy of the TR photon
256 absLength = gRandom->Exp(1.0/sigma);
257 if (absLength > (AliTRDgeometry::DrThick()
258 + AliTRDgeometry::AmThick())) {
266 // The position of the absorbtion
268 gMC->TrackPosition(pos);
269 posHit[0] = pos[0] + mom[0] / pTot * absLength;
270 posHit[1] = pos[1] + mom[1] / pTot * absLength;
271 posHit[2] = pos[2] + mom[2] / pTot * absLength;
274 Int_t q = ((Int_t) (energyeV / kWion));
276 // Add the hit to the array. TR photon hits are marked
277 // by negative charge
278 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE);
286 //_____________________________________________________________________________
287 void AliTRDv1::Init()
290 // Initialise Transition Radiation Detector after geometry has been built.
295 AliDebug(1,"Slow simulator\n");
298 AliInfo(Form("Only plane %d is sensitive"));
299 if (fSensChamber >= 0)
300 AliInfo(Form("Only chamber %d is sensitive",fSensChamber));
301 if (fSensSector >= 0) {
302 Int_t sens1 = fSensSector;
303 Int_t sens2 = fSensSector + fSensSectorRange;
304 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
305 * AliTRDgeometry::Nsect();
306 AliInfo(Form("Only sectors %d - %d are sensitive\n",sens1,sens2-1));
310 // Switch on TR simulation as default
312 AliInfo("TR simulation off");
315 fTR = new AliTRDsim();
318 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
319 const Float_t kPoti = 12.1;
320 // Maximum energy (50 keV);
321 const Float_t kEend = 50000.0;
322 // Ermilova distribution for the delta-ray spectrum
323 Float_t poti = TMath::Log(kPoti);
324 Float_t eEnd = TMath::Log(kEend);
326 // Ermilova distribution for the delta-ray spectrum
327 fDeltaE = new TF1("deltae" ,Ermilova ,poti,eEnd,0);
329 // Geant3 distribution for the delta-ray spectrum
330 fDeltaG = new TF1("deltag",IntSpecGeant,2.421257,28.536469,0);
332 AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
336 //_____________________________________________________________________________
337 void AliTRDv1::SetSensPlane(Int_t iplane)
340 // Defines the hit-sensitive plane (0-5)
343 if ((iplane < 0) || (iplane > 5)) {
344 AliWarning(Form("Wrong input value:%d",iplane));
345 AliWarning("Use standard setting");
356 //_____________________________________________________________________________
357 void AliTRDv1::SetSensChamber(Int_t ichamber)
360 // Defines the hit-sensitive chamber (0-4)
363 if ((ichamber < 0) || (ichamber > 4)) {
364 AliWarning(Form("Wrong input value: %d",ichamber));
365 AliWarning("Use standard setting");
372 fSensChamber = ichamber;
376 //_____________________________________________________________________________
377 void AliTRDv1::SetSensSector(Int_t isector)
380 // Defines the hit-sensitive sector (0-17)
383 SetSensSector(isector,1);
387 //_____________________________________________________________________________
388 void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
391 // Defines a range of hit-sensitive sectors. The range is defined by
392 // <isector> (0-17) as the starting point and <nsector> as the number
393 // of sectors to be included.
396 if ((isector < 0) || (isector > 17)) {
397 AliWarning(Form("Wrong input value <isector>: %d",isector));
398 AliWarning("Use standard setting");
400 fSensSectorRange = 0;
405 if ((nsector < 1) || (nsector > 18)) {
406 AliWarning(Form("Wrong input value <nsector>: %d",nsector));
407 AliWarning("Use standard setting");
409 fSensSectorRange = 0;
415 fSensSector = isector;
416 fSensSectorRange = nsector;
420 //_____________________________________________________________________________
421 void AliTRDv1::StepManager()
424 // Slow simulator. Every charged track produces electron cluster as hits
425 // along its path across the drift volume.
428 switch (fTypeOfStepManager) {
429 case 0 : StepManagerErmilova(); break; // 0 is Ermilova
430 case 1 : StepManagerGeant(); break; // 1 is Geant
431 case 2 : StepManagerFixedStep(); break; // 2 is fixed step
432 default : AliWarning("Not a valid Step Manager.");
437 //_____________________________________________________________________________
438 void AliTRDv1::SelectStepManager(Int_t t)
441 // Selects a step manager type:
444 // 2 - Fixed step size
447 fTypeOfStepManager = t;
448 AliInfo(Form("Step Manager type %d was selected",fTypeOfStepManager));
452 //_____________________________________________________________________________
453 void AliTRDv1::StepManagerGeant()
456 // Slow simulator. Every charged track produces electron cluster as hits
457 // along its path across the drift volume. The step size is set acording
458 // to Bethe-Bloch. The energy distribution of the delta electrons follows
459 // a spectrum taken from Geant3.
461 // Version by A. Bercuci
477 Double_t betaGamma, pp;
478 Double_t stepSize = 0;
480 Bool_t drRegion = kFALSE;
481 Bool_t amRegion = kFALSE;
484 TString cIdSensDr = "J";
485 TString cIdSensAm = "K";
486 Char_t cIdChamber[3];
489 TLorentzVector pos, mom;
493 const Int_t kNplan = AliTRDgeometry::Nplan();
494 const Int_t kNcham = AliTRDgeometry::Ncham();
495 const Int_t kNdetsec = kNplan * kNcham;
497 const Double_t kBig = 1.0E+12; // Infinitely big
498 const Float_t kWion = 23.53; // Ionization energy
499 const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
501 // Minimum energy for the step size adjustment
502 const Float_t kEkinMinStep = 1.0e-5;
503 // energy threshold for production of delta electrons
504 const Float_t kECut = 1.0e4;
505 // Parameters entering the parametrized range for delta electrons
506 const Float_t kRa = 5.37E-4;
507 const Float_t kRb = 0.9815;
508 const Float_t kRc = 3.123E-3;
509 // Gas density -> To be made user adjustable !
510 const Float_t kRho = 0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
512 // Plateau value of the energy-loss for electron in xenon
513 // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
514 //const Double_t kPlateau = 1.70;
515 // the averaged value (26/3/99)
516 const Float_t kPlateau = 1.55;
518 const Float_t kPrim = 19.34; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
519 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
520 const Float_t kPoti = 12.1;
522 const Int_t kPdgElectron = 11; // PDG code electron
524 // Set the maximum step size to a very large number for all
525 // neutral particles and those outside the driftvolume
526 gMC->SetMaxStep(kBig);
528 // Use only charged tracks
529 if (( gMC->TrackCharge() ) &&
530 (!gMC->IsTrackDisappeared())) {
532 // Inside a sensitive volume?
535 cIdCurrent = gMC->CurrentVolName();
536 if (cIdSensDr == cIdCurrent[1]) {
539 if (cIdSensAm == cIdCurrent[1]) {
542 if (drRegion || amRegion) {
544 // The hit coordinates and charge
545 gMC->TrackPosition(pos);
550 // The sector number (0 - 17)
551 // The numbering goes clockwise and starts at y = 0
552 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
557 sec = ((Int_t) (phi / 20));
559 // The plane and chamber number
560 cIdChamber[0] = cIdCurrent[2];
561 cIdChamber[1] = cIdCurrent[3];
562 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
563 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
564 pla = ((Int_t) idChamber % kNplan);
566 // Check on selected volumes
567 Int_t addthishit = 1;
569 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
570 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
571 if (fSensSector >= 0) {
572 Int_t sens1 = fSensSector;
573 Int_t sens2 = fSensSector + fSensSectorRange;
574 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
575 * AliTRDgeometry::Nsect();
577 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
580 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
588 // The detector number
589 det = fGeometry->GetDetector(pla,cha,sec);
591 // Special hits only in the drift region
594 // Create a track reference at the entrance and
595 // exit of each chamber that contain the
596 // momentum components of the particle
597 if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
598 gMC->TrackMomentum(mom);
599 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
602 if (gMC->IsTrackEntering() && !gMC->IsNewTrack()) {
603 // determine if hit belong to primary track
604 fPrimaryTrackPid = gAlice->GetMCApp()->GetCurrentTrackNumber();
605 // determine track length when entering the detector
606 fTrackLength0 = gMC->TrackLength();
609 // Create the hits from TR photons
610 if (fTR) CreateTRhit(det);
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
618 // http://www-linux.gsi.de/~abercuci/Contributions/TRD/index.html
620 // determine the most significant process (last on the processes list)
621 // which caused this hit
622 gMC->StepProcesses(processes);
623 Int_t nofprocesses = processes.GetSize();
629 pid = processes[nofprocesses-1];
632 // generate Edep according to GEANT parametrisation
633 eDelta = TMath::Exp(fDeltaG->GetRandom()) - kPoti;
634 eDelta = TMath::Max(eDelta,0.0);
635 Float_t pr_range = 0.0;
636 Float_t range = gMC->TrackLength() - fTrackLength0;
637 // merge GEANT tracker information with locally cooked one
638 if (gAlice->GetMCApp()->GetCurrentTrackNumber() == fPrimaryTrackPid) {
639 // printf("primary pid=%d eDelta=%f\n",pid,eDelta);
641 if (eDelta >= kECut) {
642 pr_range = kRa*eDelta*0.001*(1.0-kRb/(1.0+kRc*eDelta*0.001))/kRho;
643 if (pr_range >= (3.7-range)) eDelta *= 0.1;
647 if (eDelta < kECut) {
651 pr_range = kRa*eDelta*0.001*(1.0-kRb/(1.0+kRc*eDelta*0.001))/kRho;
652 if (pr_range >= ((AliTRDgeometry::DrThick()
653 + AliTRDgeometry::AmThick()) - range)) {
669 // Generate the electron cluster size
674 qTot = ((Int_t) (eDelta / kWion) + 1);
677 // Create a new dEdx hit
678 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot, drRegion);
680 // Calculate the maximum step size for the next tracking step
681 // Produce only one hit if Ekin is below cutoff
682 aMass = gMC->TrackMass();
683 if ((gMC->Etot() - aMass) > kEkinMinStep) {
685 // The energy loss according to Bethe Bloch
686 iPdg = TMath::Abs(gMC->TrackPid());
687 if ((iPdg != kPdgElectron) ||
688 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
689 gMC->TrackMomentum(mom);
691 betaGamma = pTot / aMass;
692 pp = BetheBlochGeant(betaGamma);
693 // Take charge > 1 into account
694 charge = gMC->TrackCharge();
695 if (TMath::Abs(charge) > 1) {
696 pp = pp * charge*charge;
700 // Electrons above 20 Mev/c are at the plateau
701 pp = kPrim * kPlateau;
705 do {nsteps = gRandom->Poisson(pp);} while(!nsteps);
706 stepSize = 1./nsteps;
707 gMC->SetMaxStep(stepSize);
719 //_____________________________________________________________________________
720 void AliTRDv1::StepManagerErmilova()
723 // Slow simulator. Every charged track produces electron cluster as hits
724 // along its path across the drift volume. The step size is set acording
725 // to Bethe-Bloch. The energy distribution of the delta electrons follows
726 // a spectrum taken from Ermilova et al.
743 Double_t betaGamma, pp;
746 Bool_t drRegion = kFALSE;
747 Bool_t amRegion = kFALSE;
750 TString cIdSensDr = "J";
751 TString cIdSensAm = "K";
752 Char_t cIdChamber[3];
755 TLorentzVector pos, mom;
757 const Int_t kNplan = AliTRDgeometry::Nplan();
758 const Int_t kNcham = AliTRDgeometry::Ncham();
759 const Int_t kNdetsec = kNplan * kNcham;
761 const Double_t kBig = 1.0E+12; // Infinitely big
762 const Float_t kWion = 23.53; // Ionization energy
763 const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
765 // energy threshold for production of delta electrons
766 //const Float_t kECut = 1.0e4;
767 // Parameters entering the parametrized range for delta electrons
768 //const float ra=5.37E-4, rb=0.9815, rc=3.123E-3;
769 // Gas density -> To be made user adjustable !
770 //const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
772 // Minimum energy for the step size adjustment
773 const Float_t kEkinMinStep = 1.0e-5;
775 // Plateau value of the energy-loss for electron in xenon
776 // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
777 //const Double_t kPlateau = 1.70;
778 // the averaged value (26/3/99)
779 const Float_t kPlateau = 1.55;
781 const Float_t kPrim = 48.0; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
782 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
783 const Float_t kPoti = 12.1;
785 const Int_t kPdgElectron = 11; // PDG code electron
787 // Set the maximum step size to a very large number for all
788 // neutral particles and those outside the driftvolume
789 gMC->SetMaxStep(kBig);
791 // Use only charged tracks
792 if (( gMC->TrackCharge() ) &&
793 (!gMC->IsTrackDisappeared())) {
795 // Inside a sensitive volume?
798 cIdCurrent = gMC->CurrentVolName();
799 if (cIdSensDr == cIdCurrent[1]) {
802 if (cIdSensAm == cIdCurrent[1]) {
805 if (drRegion || amRegion) {
807 // The hit coordinates and charge
808 gMC->TrackPosition(pos);
813 // The sector number (0 - 17)
814 // The numbering goes clockwise and starts at y = 0
815 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
820 sec = ((Int_t) (phi / 20));
822 // The plane and chamber number
823 cIdChamber[0] = cIdCurrent[2];
824 cIdChamber[1] = cIdCurrent[3];
825 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
826 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
827 pla = ((Int_t) idChamber % kNplan);
829 // Check on selected volumes
830 Int_t addthishit = 1;
832 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
833 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
834 if (fSensSector >= 0) {
835 Int_t sens1 = fSensSector;
836 Int_t sens2 = fSensSector + fSensSectorRange;
837 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
838 * AliTRDgeometry::Nsect();
840 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
843 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
851 // The detector number
852 det = fGeometry->GetDetector(pla,cha,sec);
854 // Special hits only in the drift region
857 // Create a track reference at the entrance and
858 // exit of each chamber that contain the
859 // momentum components of the particle
860 if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
861 gMC->TrackMomentum(mom);
862 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
864 // Create the hits from TR photons
865 if (fTR) CreateTRhit(det);
869 // Calculate the energy of the delta-electrons
870 eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
871 eDelta = TMath::Max(eDelta,0.0);
872 // Generate the electron cluster size
873 if(eDelta==0.) qTot=0;
874 else qTot = ((Int_t) (eDelta / kWion) + 1);
876 // Create a new dEdx hit
878 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
879 ,det,hits,qTot, kTRUE);
882 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
883 ,det,hits,qTot,kFALSE);
886 // Calculate the maximum step size for the next tracking step
887 // Produce only one hit if Ekin is below cutoff
888 aMass = gMC->TrackMass();
889 if ((gMC->Etot() - aMass) > kEkinMinStep) {
891 // The energy loss according to Bethe Bloch
892 iPdg = TMath::Abs(gMC->TrackPid());
893 if ( (iPdg != kPdgElectron) ||
894 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
895 gMC->TrackMomentum(mom);
897 betaGamma = pTot / aMass;
898 pp = kPrim * BetheBloch(betaGamma);
899 // Take charge > 1 into account
900 charge = gMC->TrackCharge();
901 if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
902 } else { // Electrons above 20 Mev/c are at the plateau
903 pp = kPrim * kPlateau;
908 gMC->GetRandom()->RndmArray(1, random);
909 while ((random[0] == 1.) || (random[0] == 0.));
910 stepSize = - TMath::Log(random[0]) / pp;
911 gMC->SetMaxStep(stepSize);
919 //_____________________________________________________________________________
920 void AliTRDv1::StepManagerFixedStep()
923 // Slow simulator. Every charged track produces electron cluster as hits
924 // along its path across the drift volume. The step size is fixed in
925 // this version of the step manager.
937 Bool_t drRegion = kFALSE;
938 Bool_t amRegion = kFALSE;
941 TString cIdSensDr = "J";
942 TString cIdSensAm = "K";
943 Char_t cIdChamber[3];
946 TLorentzVector pos, mom;
948 const Int_t kNplan = AliTRDgeometry::Nplan();
949 const Int_t kNcham = AliTRDgeometry::Ncham();
950 const Int_t kNdetsec = kNplan * kNcham;
952 const Double_t kBig = 1.0E+12;
954 const Float_t kWion = 23.53; // Ionization energy
955 const Float_t kEkinMinStep = 1.0e-5; // Minimum energy for the step size adjustment
957 // Set the maximum step size to a very large number for all
958 // neutral particles and those outside the driftvolume
959 gMC->SetMaxStep(kBig);
961 // If not charged track or already stopped or disappeared, just return.
962 if ((!gMC->TrackCharge()) ||
963 gMC->IsTrackDisappeared()) return;
965 // Inside a sensitive volume?
966 cIdCurrent = gMC->CurrentVolName();
968 if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
969 if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
971 if ((!drRegion) && (!amRegion)) return;
973 // The hit coordinates and charge
974 gMC->TrackPosition(pos);
979 // The sector number (0 - 17)
980 // The numbering goes clockwise and starts at y = 0
981 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
982 if (phi < 90.) phi += 270.;
984 sec = ((Int_t) (phi / 20.));
986 // The plane and chamber number
987 cIdChamber[0] = cIdCurrent[2];
988 cIdChamber[1] = cIdCurrent[3];
989 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
990 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
991 pla = ((Int_t) idChamber % kNplan);
993 // Check on selected volumes
994 Int_t addthishit = 1;
996 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
997 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
998 if (fSensSector >= 0) {
999 Int_t sens1 = fSensSector;
1000 Int_t sens2 = fSensSector + fSensSectorRange;
1001 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) * AliTRDgeometry::Nsect();
1002 if (sens1 < sens2) {
1003 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
1006 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
1011 if (!addthishit) return;
1013 det = fGeometry->GetDetector(pla,cha,sec); // The detector number
1015 Int_t trkStat = 0; // 0: InFlight 1:Entering 2:Exiting
1017 // Special hits only in the drift region
1020 // Create a track reference at the entrance and exit of each
1021 // chamber that contain the momentum components of the particle
1023 if (gMC->IsTrackEntering()) {
1024 gMC->TrackMomentum(mom);
1025 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
1028 if (gMC->IsTrackExiting()) {
1029 gMC->TrackMomentum(mom);
1030 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
1034 // Create the hits from TR photons
1035 if (fTR) CreateTRhit(det);
1039 // Calculate the charge according to GEANT Edep
1040 // Create a new dEdx hit
1041 eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
1042 qTot = (Int_t) (eDep / kWion);
1043 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
1044 ,det,hits,qTot,drRegion);
1046 // Set Maximum Step Size
1047 // Produce only one hit if Ekin is below cutoff
1048 if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return;
1049 gMC->SetMaxStep(fStepSize);
1053 //_____________________________________________________________________________
1054 Double_t AliTRDv1::BetheBloch(Double_t bg)
1057 // Parametrization of the Bethe-Bloch-curve
1058 // The parametrization is the same as for the TPC and is taken from Lehrhaus.
1061 // This parameters have been adjusted to averaged values from GEANT
1062 const Double_t kP1 = 7.17960e-02;
1063 const Double_t kP2 = 8.54196;
1064 const Double_t kP3 = 1.38065e-06;
1065 const Double_t kP4 = 5.30972;
1066 const Double_t kP5 = 2.83798;
1068 // This parameters have been adjusted to Xe-data found in:
1069 // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
1070 //const Double_t kP1 = 0.76176E-1;
1071 //const Double_t kP2 = 10.632;
1072 //const Double_t kP3 = 3.17983E-6;
1073 //const Double_t kP4 = 1.8631;
1074 //const Double_t kP5 = 1.9479;
1076 // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
1077 const Double_t kBgMin = 0.8;
1078 const Double_t kBBMax = 6.83298;
1079 //const Double_t kBgMin = 0.6;
1080 //const Double_t kBBMax = 17.2809;
1081 //const Double_t kBgMin = 0.4;
1082 //const Double_t kBBMax = 82.0;
1085 Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
1086 Double_t aa = TMath::Power(yy,kP4);
1087 Double_t bb = TMath::Power((1./bg),kP5);
1088 bb = TMath::Log(kP3 + bb);
1089 return ((kP2 - aa - bb)*kP1 / aa);
1097 //_____________________________________________________________________________
1098 Double_t AliTRDv1::BetheBlochGeant(Double_t bg)
1101 // Return dN/dx (number of primary collisions per centimeter)
1102 // for given beta*gamma factor.
1104 // Implemented by K.Oyama according to GEANT 3 parametrization shown in
1105 // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
1106 // This must be used as a set with IntSpecGeant.
1109 Double_t arr_g[20] = {
1110 1.100000, 1.200000, 1.300000, 1.500000,
1111 1.800000, 2.000000, 2.500000, 3.000000,
1112 4.000000, 7.000000, 10.000000, 20.000000,
1113 40.000000, 70.000000, 100.000000, 300.000000,
1114 600.000000, 1000.000000, 3000.000000, 10000.000000 };
1116 Double_t arr_nc[20] = {
1117 75.009056, 45.508083, 35.299252, 27.116327,
1118 22.734999, 21.411915, 19.934095, 19.449375,
1119 19.344431, 20.185553, 21.027925, 22.912676,
1120 24.933352, 26.504053, 27.387468, 29.566597,
1121 30.353779, 30.787134, 31.129285, 31.157350 };
1123 // betagamma to gamma
1124 Double_t g = TMath::Sqrt( 1. + bg*bg );
1126 // Find the index just before the point we need.
1128 for( i = 0 ; i < 18 ; i++ )
1129 if( arr_g[i] < g && arr_g[i+1] > g )
1132 // Simple interpolation.
1133 Double_t pp = ((arr_nc[i+1] - arr_nc[i]) /
1134 (arr_g[i+1]-arr_g[i])) * (g-arr_g[i]) + arr_nc[i];
1136 return pp; //arr_nc[8];
1140 //_____________________________________________________________________________
1141 void AliTRDv1::Stepping()
1160 gMC->TrackPosition(x, y, z);
1161 cout << setw(8) << setprecision(3) << x << " "
1162 << setw(8) << setprecision(3) << y << " "
1163 << setw(8) << setprecision(3) << z << " ";
1167 Double_t px, py, pz, etot;
1168 gMC->TrackMomentum(px, py, pz, etot);
1169 Double_t ekin = etot - gMC->TrackMass();
1170 cout << setw(9) << setprecision(4) << ekin*1e03 << " ";
1174 cout << setw(9) << setprecision(4) << gMC->Edep()*1e03 << " ";
1178 cout << setw(8) << setprecision(3) << gMC->TrackStep() << " ";
1182 cout << setw(8) << setprecision(3) << gMC->TrackLength() << " ";
1186 if (gMC->CurrentVolName() != 0)
1187 cout << setw(4) << gMC->CurrentVolName() << " ";
1189 cout << setw(4) << "None" << " ";
1194 Int_t nofProcesses = gMC->StepProcesses(processes);
1195 for(int ip=0;ip<nofProcesses; ip++)
1196 cout << TMCProcessName[processes[ip]]<<" / ";
1202 //_____________________________________________________________________________
1203 Double_t Ermilova(Double_t *x, Double_t *)
1206 // Calculates the delta-ray energy distribution according to Ermilova.
1207 // Logarithmic scale !
1216 const Int_t kNv = 31;
1218 Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
1219 , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
1220 , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
1221 , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
1222 , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
1223 , 9.4727, 9.9035,10.3735,10.5966,10.8198
1226 Float_t vye[kNv] = { 80.0 , 31.0 , 23.3 , 21.1 , 21.0
1227 , 20.9 , 20.8 , 20.0 , 16.0 , 11.0
1228 , 8.0 , 6.0 , 5.2 , 4.6 , 4.0
1229 , 3.5 , 3.0 , 1.4 , 0.67 , 0.44
1230 , 0.3 , 0.18 , 0.12 , 0.08 , 0.056
1231 , 0.04 , 0.023, 0.015, 0.011, 0.01
1236 // Find the position
1240 dpos = energy - vxe[pos2++];
1244 if (pos2 > kNv) pos2 = kNv - 1;
1247 // Differentiate between the sampling points
1248 dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
1254 //_____________________________________________________________________________
1255 Double_t IntSpecGeant(Double_t *x, Double_t *)
1258 // Integrated spectrum from Geant3
1261 const Int_t npts = 83;
1262 Double_t arre[npts] = {
1263 2.421257, 2.483278, 2.534301, 2.592230,
1264 2.672067, 2.813299, 3.015059, 3.216819,
1265 3.418579, 3.620338, 3.868209, 3.920198,
1266 3.978284, 4.063923, 4.186264, 4.308605,
1267 4.430946, 4.553288, 4.724261, 4.837736,
1268 4.999842, 5.161949, 5.324056, 5.486163,
1269 5.679688, 5.752998, 5.857728, 5.962457,
1270 6.067185, 6.171914, 6.315653, 6.393674,
1271 6.471694, 6.539689, 6.597658, 6.655627,
1272 6.710957, 6.763648, 6.816338, 6.876198,
1273 6.943227, 7.010257, 7.106285, 7.252151,
1274 7.460531, 7.668911, 7.877290, 8.085670,
1275 8.302979, 8.353585, 8.413120, 8.483500,
1276 8.541030, 8.592857, 8.668865, 8.820485,
1277 9.037086, 9.253686, 9.470286, 9.686887,
1278 9.930838, 9.994655, 10.085822, 10.176990,
1279 10.268158, 10.359325, 10.503614, 10.627565,
1280 10.804637, 10.981709, 11.158781, 11.335854,
1281 11.593397, 11.781165, 12.049404, 12.317644,
1282 12.585884, 12.854123, 14.278421, 16.975889,
1283 20.829416, 24.682943, 28.536469
1286 Double_t arrdndx[npts] = {
1287 19.344431, 18.664679, 18.136106, 17.567745,
1288 16.836426, 15.677382, 14.281277, 13.140237,
1289 12.207677, 11.445510, 10.697049, 10.562296,
1290 10.414673, 10.182341, 9.775256, 9.172330,
1291 8.240271, 6.898587, 4.808303, 3.889751,
1292 3.345288, 3.093431, 2.897347, 2.692470,
1293 2.436222, 2.340029, 2.208579, 2.086489,
1294 1.975535, 1.876519, 1.759626, 1.705024,
1295 1.656374, 1.502638, 1.330566, 1.200697,
1296 1.101168, 1.019323, 0.943867, 0.851951,
1297 0.755229, 0.671576, 0.570675, 0.449672,
1298 0.326722, 0.244225, 0.188225, 0.149608,
1299 0.121529, 0.116289, 0.110636, 0.103490,
1300 0.096147, 0.089191, 0.079780, 0.063927,
1301 0.047642, 0.036341, 0.028250, 0.022285,
1302 0.017291, 0.016211, 0.014802, 0.013533,
1303 0.012388, 0.011352, 0.009803, 0.008537,
1304 0.007039, 0.005829, 0.004843, 0.004034,
1305 0.003101, 0.002564, 0.001956, 0.001494,
1306 0.001142, 0.000873, 0.000210, 0.000014,
1307 0.000000, 0.000000, 0.000000
1311 // dnde = (arrdndx[i-1] - arrdndx[i]) / (arre[i] - arre[i-1]);
1313 Double_t arrdnde[npts] = {
1314 10.960000, 10.960000, 10.359500, 9.811340,
1315 9.1601500, 8.206670, 6.919630, 5.655430,
1316 4.6221300, 3.777610, 3.019560, 2.591950,
1317 2.5414600, 2.712920, 3.327460, 4.928240,
1318 7.6185300, 10.966700, 12.225800, 8.094750,
1319 3.3586900, 1.553650, 1.209600, 1.263840,
1320 1.3241100, 1.312140, 1.255130, 1.165770,
1321 1.0594500, 0.945450, 0.813231, 0.699837,
1322 0.6235580, 2.260990, 2.968350, 2.240320,
1323 1.7988300, 1.553300, 1.432070, 1.535520,
1324 1.4429900, 1.247990, 1.050750, 0.829549,
1325 0.5900280, 0.395897, 0.268741, 0.185320,
1326 0.1292120, 0.103545, 0.0949525, 0.101535,
1327 0.1276380, 0.134216, 0.123816, 0.104557,
1328 0.0751843, 0.0521745, 0.0373546, 0.0275391,
1329 0.0204713, 0.0169234, 0.0154552, 0.0139194,
1330 0.0125592, 0.0113638, 0.0107354, 0.0102137,
1331 0.00845984, 0.00683338, 0.00556836, 0.00456874,
1332 0.0036227, 0.00285991, 0.00226664, 0.00172234,
1333 0.00131226, 0.00100284, 0.000465492, 7.26607e-05,
1334 3.63304e-06, 0.0000000, 0.0000000
1338 Double_t energy = x[0];
1340 for( i = 0 ; i < npts ; i++ )
1341 if( energy < arre[i] ) break;
1344 AliErrorGeneral("AliTRDv1","Given energy value is too small or zero");