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
448 AliWarning("Sorry, Geant parametrization step manager is not implemented yet. Please ask K.Oyama for detail.");
451 fTypeOfStepManager = t;
452 AliInfo(Form("Step Manager type %d was selected",fTypeOfStepManager));
456 //_____________________________________________________________________________
457 void AliTRDv1::StepManagerGeant()
460 // Slow simulator. Every charged track produces electron cluster as hits
461 // along its path across the drift volume. The step size is set acording
462 // to Bethe-Bloch. The energy distribution of the delta electrons follows
463 // a spectrum taken from Geant3.
478 Double_t betaGamma, pp;
481 Bool_t drRegion = kFALSE;
482 Bool_t amRegion = kFALSE;
485 TString cIdSensDr = "J";
486 TString cIdSensAm = "K";
487 Char_t cIdChamber[3];
490 TLorentzVector pos, mom;
492 const Int_t kNplan = AliTRDgeometry::Nplan();
493 const Int_t kNcham = AliTRDgeometry::Ncham();
494 const Int_t kNdetsec = kNplan * kNcham;
496 const Double_t kBig = 1.0E+12; // Infinitely big
497 const Float_t kWion = 23.53; // Ionization energy
498 const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
500 // Minimum energy for the step size adjustment
501 const Float_t kEkinMinStep = 1.0e-5;
502 // energy threshold for production of delta electrons
503 const Float_t kECut = 1.0e4;
504 // Parameters entering the parametrized range for delta electrons
505 const float ra=5.37E-4, rb=0.9815, rc=3.123E-3;
506 // Gas density -> To be made user adjustable !
507 const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
509 // Plateau value of the energy-loss for electron in xenon
510 // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
511 //const Double_t kPlateau = 1.70;
512 // the averaged value (26/3/99)
513 const Float_t kPlateau = 1.55;
515 const Float_t kPrim = 19.34; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
516 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
517 const Float_t kPoti = 12.1;
519 const Int_t kPdgElectron = 11; // PDG code electron
521 // Set the maximum step size to a very large number for all
522 // neutral particles and those outside the driftvolume
523 gMC->SetMaxStep(kBig);
525 // Use only charged tracks
526 if (( gMC->TrackCharge() ) &&
527 (!gMC->IsTrackStop() ) &&
528 (!gMC->IsTrackDisappeared())) {
530 // Inside a sensitive volume?
533 cIdCurrent = gMC->CurrentVolName();
534 if (cIdSensDr == cIdCurrent[1]) {
537 if (cIdSensAm == cIdCurrent[1]) {
540 if (drRegion || amRegion) {
542 // The hit coordinates and charge
543 gMC->TrackPosition(pos);
548 // The sector number (0 - 17)
549 // The numbering goes clockwise and starts at y = 0
550 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
555 sec = ((Int_t) (phi / 20));
557 // The plane and chamber number
558 cIdChamber[0] = cIdCurrent[2];
559 cIdChamber[1] = cIdCurrent[3];
560 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
561 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
562 pla = ((Int_t) idChamber % kNplan);
564 // Check on selected volumes
565 Int_t addthishit = 1;
567 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
568 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
569 if (fSensSector >= 0) {
570 Int_t sens1 = fSensSector;
571 Int_t sens2 = fSensSector + fSensSectorRange;
572 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
573 * AliTRDgeometry::Nsect();
575 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
578 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
586 // The detector number
587 det = fGeometry->GetDetector(pla,cha,sec);
589 // Special hits only in the drift region
591 // Create a track reference at the entrance and
592 // exit of each chamber that contain the
593 // momentum components of the particle
594 if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
595 gMC->TrackMomentum(mom);
596 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
599 if (gMC->IsTrackEntering() && !gMC->IsNewTrack()) {
600 // determine if hit belong to primary track
601 fPrimaryTrackPid=gAlice->GetMCApp()->GetCurrentTrackNumber();
602 //determine track length when entering the detector
603 fTrackLength0=gMC->TrackLength();
606 // Create the hits from TR photons
607 if (fTR) CreateTRhit(det);
610 // Calculate the energy of the delta-electrons
611 // modified by Alex Bercuci (A.Bercuci@gsi.de) on 26.01.06
612 // take into account correlation with the underlying GEANT tracking
614 // http://www-linux.gsi.de/~abercuci/Contributions/TRD/index.html
616 // determine the most significant process (last on the processes list)
617 // which caused this hit
619 gMC->StepProcesses(processes);
620 int nofprocesses=processes.GetSize(), pid;
621 if(!nofprocesses) pid=0;
622 else pid= processes[nofprocesses-1];
624 // generate Edep according to GEANT parametrisation
625 eDelta =TMath::Exp(fDeltaG->GetRandom()) - kPoti;
626 eDelta=TMath::Max(eDelta,0.0);
628 float range=gMC->TrackLength()-fTrackLength0;
629 // merge GEANT tracker information with localy cooked one
630 if(gAlice->GetMCApp()->GetCurrentTrackNumber()==fPrimaryTrackPid) {
631 // printf("primary pid=%d eDelta=%f\n",pid,eDelta);
634 pr_range=ra*eDelta*.001*(1.-rb/(1.+rc*eDelta*0.001))/rho;
635 if(pr_range>=(3.7-range)) eDelta*=.1;
638 if(eDelta<kECut) eDelta*=.5;
640 pr_range=ra*eDelta*.001*(1.-rb/(1.+rc*eDelta*0.001))/rho;
641 if(pr_range>=((AliTRDgeometry::DrThick()
642 + AliTRDgeometry::AmThick())-range)) eDelta*=.05;
648 // Generate the electron cluster size
649 if(eDelta==0.) qTot=0;
650 else qTot = ((Int_t) (eDelta / kWion) + 1);
651 // Create a new dEdx hit
652 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot, drRegion);
654 // Calculate the maximum step size for the next tracking step
655 // Produce only one hit if Ekin is below cutoff
656 aMass = gMC->TrackMass();
657 if ((gMC->Etot() - aMass) > kEkinMinStep) {
659 // The energy loss according to Bethe Bloch
660 iPdg = TMath::Abs(gMC->TrackPid());
661 if ( (iPdg != kPdgElectron) ||
662 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
663 gMC->TrackMomentum(mom);
665 betaGamma = pTot / aMass;
666 pp = BetheBlochGeant(betaGamma);
667 // Take charge > 1 into account
668 charge = gMC->TrackCharge();
669 if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
670 } else { // Electrons above 20 Mev/c are at the plateau
671 pp = kPrim * kPlateau;
675 do {nsteps = gRandom->Poisson(pp);} while(!nsteps);
676 stepSize = 1./nsteps;
677 gMC->SetMaxStep(stepSize);
684 //_____________________________________________________________________________
685 void AliTRDv1::StepManagerErmilova()
688 // Slow simulator. Every charged track produces electron cluster as hits
689 // along its path across the drift volume. The step size is set acording
690 // to Bethe-Bloch. The energy distribution of the delta electrons follows
691 // a spectrum taken from Ermilova et al.
708 Double_t betaGamma, pp;
711 Bool_t drRegion = kFALSE;
712 Bool_t amRegion = kFALSE;
715 TString cIdSensDr = "J";
716 TString cIdSensAm = "K";
717 Char_t cIdChamber[3];
720 TLorentzVector pos, mom;
722 const Int_t kNplan = AliTRDgeometry::Nplan();
723 const Int_t kNcham = AliTRDgeometry::Ncham();
724 const Int_t kNdetsec = kNplan * kNcham;
726 const Double_t kBig = 1.0E+12; // Infinitely big
727 const Float_t kWion = 23.53; // Ionization energy
728 const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
730 // energy threshold for production of delta electrons
731 //const Float_t kECut = 1.0e4;
732 // Parameters entering the parametrized range for delta electrons
733 //const float ra=5.37E-4, rb=0.9815, rc=3.123E-3;
734 // Gas density -> To be made user adjustable !
735 //const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
737 // Minimum energy for the step size adjustment
738 const Float_t kEkinMinStep = 1.0e-5;
740 // Plateau value of the energy-loss for electron in xenon
741 // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
742 //const Double_t kPlateau = 1.70;
743 // the averaged value (26/3/99)
744 const Float_t kPlateau = 1.55;
746 const Float_t kPrim = 48.0; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
747 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
748 const Float_t kPoti = 12.1;
750 const Int_t kPdgElectron = 11; // PDG code electron
752 // Set the maximum step size to a very large number for all
753 // neutral particles and those outside the driftvolume
754 gMC->SetMaxStep(kBig);
756 // Use only charged tracks
757 if (( gMC->TrackCharge() ) &&
758 (!gMC->IsTrackStop() ) &&
759 (!gMC->IsTrackDisappeared())) {
761 // Inside a sensitive volume?
764 cIdCurrent = gMC->CurrentVolName();
765 if (cIdSensDr == cIdCurrent[1]) {
768 if (cIdSensAm == cIdCurrent[1]) {
771 if (drRegion || amRegion) {
773 // The hit coordinates and charge
774 gMC->TrackPosition(pos);
779 // The sector number (0 - 17)
780 // The numbering goes clockwise and starts at y = 0
781 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
786 sec = ((Int_t) (phi / 20));
788 // The plane and chamber number
789 cIdChamber[0] = cIdCurrent[2];
790 cIdChamber[1] = cIdCurrent[3];
791 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
792 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
793 pla = ((Int_t) idChamber % kNplan);
795 // Check on selected volumes
796 Int_t addthishit = 1;
798 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
799 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
800 if (fSensSector >= 0) {
801 Int_t sens1 = fSensSector;
802 Int_t sens2 = fSensSector + fSensSectorRange;
803 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
804 * AliTRDgeometry::Nsect();
806 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
809 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
817 // The detector number
818 det = fGeometry->GetDetector(pla,cha,sec);
820 // Special hits only in the drift region
823 // Create a track reference at the entrance and
824 // exit of each chamber that contain the
825 // momentum components of the particle
826 if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
827 gMC->TrackMomentum(mom);
828 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
830 // Create the hits from TR photons
831 if (fTR) CreateTRhit(det);
835 // Calculate the energy of the delta-electrons
836 eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
837 eDelta = TMath::Max(eDelta,0.0);
838 // Generate the electron cluster size
839 if(eDelta==0.) qTot=0;
840 else qTot = ((Int_t) (eDelta / kWion) + 1);
842 // Create a new dEdx hit
844 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
845 ,det,hits,qTot, kTRUE);
848 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
849 ,det,hits,qTot,kFALSE);
852 // Calculate the maximum step size for the next tracking step
853 // Produce only one hit if Ekin is below cutoff
854 aMass = gMC->TrackMass();
855 if ((gMC->Etot() - aMass) > kEkinMinStep) {
857 // The energy loss according to Bethe Bloch
858 iPdg = TMath::Abs(gMC->TrackPid());
859 if ( (iPdg != kPdgElectron) ||
860 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
861 gMC->TrackMomentum(mom);
863 betaGamma = pTot / aMass;
864 pp = kPrim * BetheBloch(betaGamma);
865 // Take charge > 1 into account
866 charge = gMC->TrackCharge();
867 if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
868 } else { // Electrons above 20 Mev/c are at the plateau
869 pp = kPrim * kPlateau;
874 gMC->GetRandom()->RndmArray(1, random);
875 while ((random[0] == 1.) || (random[0] == 0.));
876 stepSize = - TMath::Log(random[0]) / pp;
877 gMC->SetMaxStep(stepSize);
885 //_____________________________________________________________________________
886 void AliTRDv1::StepManagerFixedStep()
889 // Slow simulator. Every charged track produces electron cluster as hits
890 // along its path across the drift volume. The step size is fixed in
891 // this version of the step manager.
903 Bool_t drRegion = kFALSE;
904 Bool_t amRegion = kFALSE;
907 TString cIdSensDr = "J";
908 TString cIdSensAm = "K";
909 Char_t cIdChamber[3];
912 TLorentzVector pos, mom;
914 const Int_t kNplan = AliTRDgeometry::Nplan();
915 const Int_t kNcham = AliTRDgeometry::Ncham();
916 const Int_t kNdetsec = kNplan * kNcham;
918 const Double_t kBig = 1.0E+12;
920 const Float_t kWion = 23.53; // Ionization energy
921 const Float_t kEkinMinStep = 1.0e-5; // Minimum energy for the step size adjustment
923 // Set the maximum step size to a very large number for all
924 // neutral particles and those outside the driftvolume
925 gMC->SetMaxStep(kBig);
927 // If not charged track or already stopped or disappeared, just return.
928 if ((!gMC->TrackCharge()) ||
929 gMC->IsTrackStop() ||
930 gMC->IsTrackDisappeared()) return;
932 // Inside a sensitive volume?
933 cIdCurrent = gMC->CurrentVolName();
935 if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
936 if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
938 if ((!drRegion) && (!amRegion)) return;
940 // The hit coordinates and charge
941 gMC->TrackPosition(pos);
946 // The sector number (0 - 17)
947 // The numbering goes clockwise and starts at y = 0
948 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
949 if (phi < 90.) phi += 270.;
951 sec = ((Int_t) (phi / 20.));
953 // The plane and chamber number
954 cIdChamber[0] = cIdCurrent[2];
955 cIdChamber[1] = cIdCurrent[3];
956 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
957 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
958 pla = ((Int_t) idChamber % kNplan);
960 // Check on selected volumes
961 Int_t addthishit = 1;
963 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
964 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
965 if (fSensSector >= 0) {
966 Int_t sens1 = fSensSector;
967 Int_t sens2 = fSensSector + fSensSectorRange;
968 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) * AliTRDgeometry::Nsect();
970 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
973 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
978 if (!addthishit) return;
980 det = fGeometry->GetDetector(pla,cha,sec); // The detector number
982 Int_t trkStat = 0; // 0: InFlight 1:Entering 2:Exiting
984 // Special hits only in the drift region
987 // Create a track reference at the entrance and exit of each
988 // chamber that contain the momentum components of the particle
990 if (gMC->IsTrackEntering()) {
991 gMC->TrackMomentum(mom);
992 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
995 if (gMC->IsTrackExiting()) {
996 gMC->TrackMomentum(mom);
997 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
1001 // Create the hits from TR photons
1002 if (fTR) CreateTRhit(det);
1006 // Calculate the charge according to GEANT Edep
1007 // Create a new dEdx hit
1008 eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
1009 qTot = (Int_t) (eDep / kWion);
1010 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
1011 ,det,hits,qTot,drRegion);
1013 // Set Maximum Step Size
1014 // Produce only one hit if Ekin is below cutoff
1015 if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return;
1016 gMC->SetMaxStep(fStepSize);
1020 //_____________________________________________________________________________
1021 Double_t AliTRDv1::BetheBloch(Double_t bg)
1024 // Parametrization of the Bethe-Bloch-curve
1025 // The parametrization is the same as for the TPC and is taken from Lehrhaus.
1028 // This parameters have been adjusted to averaged values from GEANT
1029 const Double_t kP1 = 7.17960e-02;
1030 const Double_t kP2 = 8.54196;
1031 const Double_t kP3 = 1.38065e-06;
1032 const Double_t kP4 = 5.30972;
1033 const Double_t kP5 = 2.83798;
1035 // This parameters have been adjusted to Xe-data found in:
1036 // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
1037 //const Double_t kP1 = 0.76176E-1;
1038 //const Double_t kP2 = 10.632;
1039 //const Double_t kP3 = 3.17983E-6;
1040 //const Double_t kP4 = 1.8631;
1041 //const Double_t kP5 = 1.9479;
1043 // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
1044 const Double_t kBgMin = 0.8;
1045 const Double_t kBBMax = 6.83298;
1046 //const Double_t kBgMin = 0.6;
1047 //const Double_t kBBMax = 17.2809;
1048 //const Double_t kBgMin = 0.4;
1049 //const Double_t kBBMax = 82.0;
1052 Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
1053 Double_t aa = TMath::Power(yy,kP4);
1054 Double_t bb = TMath::Power((1./bg),kP5);
1055 bb = TMath::Log(kP3 + bb);
1056 return ((kP2 - aa - bb)*kP1 / aa);
1064 //_____________________________________________________________________________
1065 Double_t AliTRDv1::BetheBlochGeant(Double_t bg)
1068 // Return dN/dx (number of primary collisions per centimeter)
1069 // for given beta*gamma factor.
1071 // Implemented by K.Oyama according to GEANT 3 parametrization shown in
1072 // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
1073 // This must be used as a set with IntSpecGeant.
1076 Double_t arr_g[20] = {
1077 1.100000, 1.200000, 1.300000, 1.500000,
1078 1.800000, 2.000000, 2.500000, 3.000000,
1079 4.000000, 7.000000, 10.000000, 20.000000,
1080 40.000000, 70.000000, 100.000000, 300.000000,
1081 600.000000, 1000.000000, 3000.000000, 10000.000000 };
1083 Double_t arr_nc[20] = {
1084 75.009056, 45.508083, 35.299252, 27.116327,
1085 22.734999, 21.411915, 19.934095, 19.449375,
1086 19.344431, 20.185553, 21.027925, 22.912676,
1087 24.933352, 26.504053, 27.387468, 29.566597,
1088 30.353779, 30.787134, 31.129285, 31.157350 };
1090 // betagamma to gamma
1091 Double_t g = TMath::Sqrt( 1. + bg*bg );
1093 // Find the index just before the point we need.
1095 for( i = 0 ; i < 18 ; i++ )
1096 if( arr_g[i] < g && arr_g[i+1] > g )
1099 // Simple interpolation.
1100 Double_t pp = ((arr_nc[i+1] - arr_nc[i]) /
1101 (arr_g[i+1]-arr_g[i])) * (g-arr_g[i]) + arr_nc[i];
1103 return pp; //arr_nc[8];
1107 //_____________________________________________________________________________
1108 void AliTRDv1::Stepping()
1127 gMC->TrackPosition(x, y, z);
1128 cout << setw(8) << setprecision(3) << x << " "
1129 << setw(8) << setprecision(3) << y << " "
1130 << setw(8) << setprecision(3) << z << " ";
1134 Double_t px, py, pz, etot;
1135 gMC->TrackMomentum(px, py, pz, etot);
1136 Double_t ekin = etot - gMC->TrackMass();
1137 cout << setw(9) << setprecision(4) << ekin*1e03 << " ";
1141 cout << setw(9) << setprecision(4) << gMC->Edep()*1e03 << " ";
1145 cout << setw(8) << setprecision(3) << gMC->TrackStep() << " ";
1149 cout << setw(8) << setprecision(3) << gMC->TrackLength() << " ";
1153 if (gMC->CurrentVolName() != 0)
1154 cout << setw(4) << gMC->CurrentVolName() << " ";
1156 cout << setw(4) << "None" << " ";
1161 Int_t nofProcesses = gMC->StepProcesses(processes);
1162 for(int ip=0;ip<nofProcesses; ip++)
1163 cout << TMCProcessName[processes[ip]]<<" / ";
1169 //_____________________________________________________________________________
1170 Double_t Ermilova(Double_t *x, Double_t *)
1173 // Calculates the delta-ray energy distribution according to Ermilova.
1174 // Logarithmic scale !
1183 const Int_t kNv = 31;
1185 Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
1186 , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
1187 , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
1188 , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
1189 , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
1190 , 9.4727, 9.9035,10.3735,10.5966,10.8198
1193 Float_t vye[kNv] = { 80.0 , 31.0 , 23.3 , 21.1 , 21.0
1194 , 20.9 , 20.8 , 20.0 , 16.0 , 11.0
1195 , 8.0 , 6.0 , 5.2 , 4.6 , 4.0
1196 , 3.5 , 3.0 , 1.4 , 0.67 , 0.44
1197 , 0.3 , 0.18 , 0.12 , 0.08 , 0.056
1198 , 0.04 , 0.023, 0.015, 0.011, 0.01
1203 // Find the position
1207 dpos = energy - vxe[pos2++];
1211 if (pos2 > kNv) pos2 = kNv - 1;
1214 // Differentiate between the sampling points
1215 dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
1221 //_____________________________________________________________________________
1222 Double_t IntSpecGeant(Double_t *x, Double_t *)
1225 // Integrated spectrum from Geant3
1228 const Int_t npts = 83;
1229 Double_t arre[npts] = {
1230 2.421257, 2.483278, 2.534301, 2.592230,
1231 2.672067, 2.813299, 3.015059, 3.216819,
1232 3.418579, 3.620338, 3.868209, 3.920198,
1233 3.978284, 4.063923, 4.186264, 4.308605,
1234 4.430946, 4.553288, 4.724261, 4.837736,
1235 4.999842, 5.161949, 5.324056, 5.486163,
1236 5.679688, 5.752998, 5.857728, 5.962457,
1237 6.067185, 6.171914, 6.315653, 6.393674,
1238 6.471694, 6.539689, 6.597658, 6.655627,
1239 6.710957, 6.763648, 6.816338, 6.876198,
1240 6.943227, 7.010257, 7.106285, 7.252151,
1241 7.460531, 7.668911, 7.877290, 8.085670,
1242 8.302979, 8.353585, 8.413120, 8.483500,
1243 8.541030, 8.592857, 8.668865, 8.820485,
1244 9.037086, 9.253686, 9.470286, 9.686887,
1245 9.930838, 9.994655, 10.085822, 10.176990,
1246 10.268158, 10.359325, 10.503614, 10.627565,
1247 10.804637, 10.981709, 11.158781, 11.335854,
1248 11.593397, 11.781165, 12.049404, 12.317644,
1249 12.585884, 12.854123, 14.278421, 16.975889,
1250 20.829416, 24.682943, 28.536469
1253 Double_t arrdndx[npts] = {
1254 19.344431, 18.664679, 18.136106, 17.567745,
1255 16.836426, 15.677382, 14.281277, 13.140237,
1256 12.207677, 11.445510, 10.697049, 10.562296,
1257 10.414673, 10.182341, 9.775256, 9.172330,
1258 8.240271, 6.898587, 4.808303, 3.889751,
1259 3.345288, 3.093431, 2.897347, 2.692470,
1260 2.436222, 2.340029, 2.208579, 2.086489,
1261 1.975535, 1.876519, 1.759626, 1.705024,
1262 1.656374, 1.502638, 1.330566, 1.200697,
1263 1.101168, 1.019323, 0.943867, 0.851951,
1264 0.755229, 0.671576, 0.570675, 0.449672,
1265 0.326722, 0.244225, 0.188225, 0.149608,
1266 0.121529, 0.116289, 0.110636, 0.103490,
1267 0.096147, 0.089191, 0.079780, 0.063927,
1268 0.047642, 0.036341, 0.028250, 0.022285,
1269 0.017291, 0.016211, 0.014802, 0.013533,
1270 0.012388, 0.011352, 0.009803, 0.008537,
1271 0.007039, 0.005829, 0.004843, 0.004034,
1272 0.003101, 0.002564, 0.001956, 0.001494,
1273 0.001142, 0.000873, 0.000210, 0.000014,
1274 0.000000, 0.000000, 0.000000
1278 // dnde = (arrdndx[i-1] - arrdndx[i]) / (arre[i] - arre[i-1]);
1280 Double_t arrdnde[npts] = {
1281 10.960000, 10.960000, 10.359500, 9.811340,
1282 9.1601500, 8.206670, 6.919630, 5.655430,
1283 4.6221300, 3.777610, 3.019560, 2.591950,
1284 2.5414600, 2.712920, 3.327460, 4.928240,
1285 7.6185300, 10.966700, 12.225800, 8.094750,
1286 3.3586900, 1.553650, 1.209600, 1.263840,
1287 1.3241100, 1.312140, 1.255130, 1.165770,
1288 1.0594500, 0.945450, 0.813231, 0.699837,
1289 0.6235580, 2.260990, 2.968350, 2.240320,
1290 1.7988300, 1.553300, 1.432070, 1.535520,
1291 1.4429900, 1.247990, 1.050750, 0.829549,
1292 0.5900280, 0.395897, 0.268741, 0.185320,
1293 0.1292120, 0.103545, 0.0949525, 0.101535,
1294 0.1276380, 0.134216, 0.123816, 0.104557,
1295 0.0751843, 0.0521745, 0.0373546, 0.0275391,
1296 0.0204713, 0.0169234, 0.0154552, 0.0139194,
1297 0.0125592, 0.0113638, 0.0107354, 0.0102137,
1298 0.00845984, 0.00683338, 0.00556836, 0.00456874,
1299 0.0036227, 0.00285991, 0.00226664, 0.00172234,
1300 0.00131226, 0.00100284, 0.000465492, 7.26607e-05,
1301 3.63304e-06, 0.0000000, 0.0000000
1305 Double_t energy = x[0];
1307 for( i = 0 ; i < npts ; i++ )
1308 if( energy < arre[i] ) break;
1311 AliErrorGeneral("AliTRDv1","Given energy value is too small or zero");