-{
- //
- // Slow simulator. Every charged track produces electron cluster as hits
- // along its path across the drift volume.
- //
-
- switch (fTypeOfStepManager) {
- case 0:
- StepManagerErmilova();
- break;
- case 1:
- StepManagerGeant();
- break;
- case 2:
- StepManagerFixedStep();
- break;
- default:
- AliWarning("Not a valid Step Manager.");
- }
-
-}
-
-//_____________________________________________________________________________
-void AliTRDv1::SelectStepManager(Int_t t)
-{
- //
- // Selects a step manager type:
- // 0 - Ermilova
- // 1 - Geant3
- // 2 - Fixed step size
- //
-
- fTypeOfStepManager = t;
- AliInfo(Form("Step Manager type %d was selected",fTypeOfStepManager));
-
-}
-
-//_____________________________________________________________________________
-void AliTRDv1::StepManagerGeant()
-{
- //
- // Slow simulator. Every charged track produces electron cluster as hits
- // along its path across the drift volume. The step size is set acording
- // to Bethe-Bloch. The energy distribution of the delta electrons follows
- // a spectrum taken from Geant3.
- //
- // Version by A. Bercuci
- //
-
- Int_t pla = 0;
- Int_t cha = 0;
- Int_t sec = 0;
- Int_t det = 0;
- Int_t iPdg;
- Int_t qTot;
-
- Float_t hits[3];
- Float_t charge;
- Float_t aMass;
-
- Double_t pTot = 0;
- Double_t eDelta;
- Double_t betaGamma;
- Double_t pp;
- Double_t stepSize = 0;
-
- Bool_t drRegion = kFALSE;
- Bool_t amRegion = kFALSE;
-
- TString cIdPath;
- Char_t cIdSector[3];
- cIdSector[2] = 0;
-
- TString cIdCurrent;
- TString cIdSensDr = "J";
- TString cIdSensAm = "K";
- Char_t cIdChamber[3];
- cIdChamber[2] = 0;
-
- TLorentzVector pos;
- TLorentzVector mom;
-
- TArrayI processes;
-
- const Int_t kNplan = AliTRDgeometry::Nplan();
- const Int_t kNcham = AliTRDgeometry::Ncham();
- const Int_t kNdetsec = kNplan * kNcham;
-
- const Double_t kBig = 1.0e+12; // Infinitely big
- const Float_t kWion = 23.53; // Ionization energy
- const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
-
- // Minimum energy for the step size adjustment
- const Float_t kEkinMinStep = 1.0e-5;
- // energy threshold for production of delta electrons
- const Float_t kECut = 1.0e4;
- // Parameters entering the parametrized range for delta electrons
- const Float_t kRa = 5.37e-4;
- const Float_t kRb = 0.9815;
- const Float_t kRc = 3.123e-3;
- // Gas density -> To be made user adjustable !
- // [0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
- const Float_t kRho = 0.004945 ;
-
- // Plateau value of the energy-loss for electron in xenon
- // The averaged value (26/3/99)
- const Float_t kPlateau = 1.55;
- // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
- const Float_t kPrim = 19.34;
- // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
- const Float_t kPoti = 12.1;
- // PDG code electron
- const Int_t kPdgElectron = 11;
-
- // Set the maximum step size to a very large number for all
- // neutral particles and those outside the driftvolume
- gMC->SetMaxStep(kBig);
-
- // Use only charged tracks
- if (( gMC->TrackCharge() ) &&
- (!gMC->IsTrackDisappeared())) {
-
- // Inside a sensitive volume?
- drRegion = kFALSE;
- amRegion = kFALSE;
- cIdCurrent = gMC->CurrentVolName();
- if (cIdSensDr == cIdCurrent[1]) {
- drRegion = kTRUE;
- }
- if (cIdSensAm == cIdCurrent[1]) {
- amRegion = kTRUE;
- }
- if (drRegion || amRegion) {
-
- // The hit coordinates and charge
- gMC->TrackPosition(pos);
- hits[0] = pos[0];
- hits[1] = pos[1];
- hits[2] = pos[2];
-
- // The sector number (0 - 17), according to standard coordinate system
- cIdPath = gGeoManager->GetPath();
- cIdSector[0] = cIdPath[21];
- cIdSector[1] = cIdPath[22];
- sec = atoi(cIdSector);
-
- // The plane and chamber number
- cIdChamber[0] = cIdCurrent[2];
- cIdChamber[1] = cIdCurrent[3];
- Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
- cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
- pla = ((Int_t) idChamber % kNplan);
-
- // The detector number
- det = fGeometry->GetDetector(pla,cha,sec);
-
- // Special hits only in the drift region
- if ((drRegion) &&
- (gMC->IsTrackEntering())) {
-
- // Create a track reference at the entrance of each
- // chamber that contains the momentum components of the particle
- gMC->TrackMomentum(mom);
- AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
-
- // Create the hits from TR photons if electron/positron is
- // entering the drift volume
- if ((fTR) &&
- (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
- CreateTRhit(det);
- }
-
- }
- else if ((amRegion) &&
- (gMC->IsTrackExiting())) {
-
- // Create a track reference at the exit of each
- // chamber that contains the momentum components of the particle
- gMC->TrackMomentum(mom);
- AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
-
- }
-
- // Calculate the energy of the delta-electrons
- // modified by Alex Bercuci (A.Bercuci@gsi.de) on 26.01.06
- // take into account correlation with the underlying GEANT tracking
- // mechanism. see
- // http://www-linux.gsi.de/~abercuci/Contributions/TRD/index.html
- //
- // determine the most significant process (last on the processes list)
- // which caused this hit
- gMC->StepProcesses(processes);
- Int_t nofprocesses = processes.GetSize();
- Int_t pid;
- if (!nofprocesses) {
- pid = 0;
- }
- else {
- pid = processes[nofprocesses-1];
- }
-
- // Generate Edep according to GEANT parametrisation
- eDelta = TMath::Exp(fDeltaG->GetRandom()) - kPoti;
- eDelta = TMath::Max(eDelta,0.0);
- Float_t prRange = 0.0;
- Float_t range = gMC->TrackLength() - fTrackLength0;
- // merge GEANT tracker information with locally cooked one
- if (gAlice->GetMCApp()->GetCurrentTrackNumber() == fPrimaryTrackPid) {
- if (pid == 27) {
- if (eDelta >= kECut) {
- prRange = kRa * eDelta * 0.001
- * (1.0 - kRb / (1.0 + kRc * eDelta * 0.001)) / kRho;
- if (prRange >= (3.7 - range)) {
- eDelta *= 0.1;
- }
- }
- }
- else if (pid == 1) {
- if (eDelta < kECut) {
- eDelta *= 0.5;
- }
- else {
- prRange = kRa * eDelta * 0.001
- * (1.0 - kRb / (1.0 + kRc * eDelta * 0.001)) / kRho;
- if (prRange >= ((AliTRDgeometry::DrThick()
- + AliTRDgeometry::AmThick()) - range)) {
- eDelta *= 0.05;
- }
- else {
- eDelta *= 0.5;
- }
- }
- }
- else {
- eDelta = 0.0;
- }
- }
- else {
- eDelta = 0.0;
- }
-
- // Generate the electron cluster size
- if (eDelta > 0.0) {
-
- qTot = ((Int_t) (eDelta / kWion) + 1);
-
- // Create a new dEdx hit
- AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
- ,det
- ,hits
- ,qTot
- ,gMC->TrackTime()*1.0e06
- ,drRegion);
-
- }
-
- // Calculate the maximum step size for the next tracking step
- // Produce only one hit if Ekin is below cutoff
- aMass = gMC->TrackMass();
- if ((gMC->Etot() - aMass) > kEkinMinStep) {
-
- // The energy loss according to Bethe Bloch
- iPdg = TMath::Abs(gMC->TrackPid());
- if ((iPdg != kPdgElectron) ||
- ((iPdg == kPdgElectron) &&
- (pTot < kPTotMaxEl))) {
- gMC->TrackMomentum(mom);
- pTot = mom.Rho();
- betaGamma = pTot / aMass;
- pp = BetheBlochGeant(betaGamma);
- // Take charge > 1 into account
- charge = gMC->TrackCharge();
- if (TMath::Abs(charge) > 1) {
- pp = pp * charge*charge;
- }
- }
- else {
- // Electrons above 20 Mev/c are at the plateau
- pp = kPrim * kPlateau;
- }
-
- Int_t nsteps = 0;
- do {
- nsteps = gRandom->Poisson(pp);
- } while(!nsteps);
- stepSize = 1.0 / nsteps;
- gMC->SetMaxStep(stepSize);
-
- }
-
- }
-
- }
-
-}
-
-//_____________________________________________________________________________
-void AliTRDv1::StepManagerErmilova()
-{
- //
- // Slow simulator. Every charged track produces electron cluster as hits
- // along its path across the drift volume. The step size is set acording
- // to Bethe-Bloch. The energy distribution of the delta electrons follows
- // a spectrum taken from Ermilova et al.
- //
-
- Int_t pla = 0;
- Int_t cha = 0;
- Int_t sec = 0;
- Int_t det = 0;
- Int_t iPdg;
- Int_t qTot;
-
- Float_t hits[3];
- Double_t random[1];
- Float_t charge;
- Float_t aMass;
-
- Double_t pTot = 0.0;
- Double_t eDelta;
- Double_t betaGamma;
- Double_t pp;
- Double_t stepSize;
-
- Bool_t drRegion = kFALSE;
- Bool_t amRegion = kFALSE;
-
- TString cIdPath;
- Char_t cIdSector[3];
- cIdSector[2] = 0;
-
- TString cIdCurrent;
- TString cIdSensDr = "J";
- TString cIdSensAm = "K";
- Char_t cIdChamber[3];
- cIdChamber[2] = 0;
-
- TLorentzVector pos;
- TLorentzVector mom;
-
- const Int_t kNplan = AliTRDgeometry::Nplan();
- const Int_t kNcham = AliTRDgeometry::Ncham();
- const Int_t kNdetsec = kNplan * kNcham;
-
- const Double_t kBig = 1.0e+12; // Infinitely big
- const Float_t kWion = 23.53; // Ionization energy
- const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
-
- // Minimum energy for the step size adjustment
- const Float_t kEkinMinStep = 1.0e-5;
-
- // Plateau value of the energy-loss for electron in xenon
- // The averaged value (26/3/99)
- const Float_t kPlateau = 1.55;
- // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
- const Float_t kPrim = 48.0;
- // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
- const Float_t kPoti = 12.1;
- // PDG code electron
- const Int_t kPdgElectron = 11;
-
- // Set the maximum step size to a very large number for all
- // neutral particles and those outside the driftvolume
- gMC->SetMaxStep(kBig);
-
- // Use only charged tracks
- if (( gMC->TrackCharge() ) &&
- (!gMC->IsTrackDisappeared())) {
-
- // Inside a sensitive volume?
- drRegion = kFALSE;
- amRegion = kFALSE;
- cIdCurrent = gMC->CurrentVolName();
- if (cIdSensDr == cIdCurrent[1]) {
- drRegion = kTRUE;
- }
- if (cIdSensAm == cIdCurrent[1]) {
- amRegion = kTRUE;
- }
- if (drRegion || amRegion) {
-
- // The hit coordinates and charge
- gMC->TrackPosition(pos);
- hits[0] = pos[0];
- hits[1] = pos[1];
- hits[2] = pos[2];
-
- // The sector number (0 - 17), according to standard coordinate system
- cIdPath = gGeoManager->GetPath();
- cIdSector[0] = cIdPath[21];
- cIdSector[1] = cIdPath[22];
- sec = atoi(cIdSector);
-
- // The plane and chamber number
- cIdChamber[0] = cIdCurrent[2];
- cIdChamber[1] = cIdCurrent[3];
- Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
- cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
- pla = ((Int_t) idChamber % kNplan);
-
- // The detector number
- det = fGeometry->GetDetector(pla,cha,sec);
-
- // Special hits only in the drift region
- if ((drRegion) &&
- (gMC->IsTrackEntering())) {
-
- // Create a track reference at the entrance of each
- // chamber that contains the momentum components of the particle
- gMC->TrackMomentum(mom);
- AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
-
- // Create the hits from TR photons if electron/positron is
- // entering the drift volume
- if ((fTR) &&
- (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
- CreateTRhit(det);
- }
-
- }
- else if ((amRegion) &&
- (gMC->IsTrackExiting())) {
-
- // Create a track reference at the exit of each
- // chamber that contains the momentum components of the particle
- gMC->TrackMomentum(mom);
- AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
-
- }
-
- // Calculate the energy of the delta-electrons
- eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
- eDelta = TMath::Max(eDelta,0.0);
-
- // Generate the electron cluster size
- if (eDelta > 0.0) {
-
- qTot = ((Int_t) (eDelta / kWion) + 1);
-
- // Create a new dEdx hit
- if (drRegion) {
- AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
- ,det
- ,hits
- ,qTot
- ,gMC->TrackTime()*1.0e06
- ,kTRUE);
- }
- else {
- AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
- ,det
- ,hits
- ,qTot
- ,gMC->TrackTime()*1.0e06
- ,kFALSE);
- }
-
- }
-
- // Calculate the maximum step size for the next tracking step
- // Produce only one hit if Ekin is below cutoff
- aMass = gMC->TrackMass();
- if ((gMC->Etot() - aMass) > kEkinMinStep) {
-
- // The energy loss according to Bethe Bloch
- iPdg = TMath::Abs(gMC->TrackPid());
- if ((iPdg != kPdgElectron) ||
- ((iPdg == kPdgElectron) &&
- (pTot < kPTotMaxEl))) {
- gMC->TrackMomentum(mom);
- pTot = mom.Rho();
- betaGamma = pTot / aMass;
- pp = kPrim * BetheBloch(betaGamma);
- // Take charge > 1 into account
- charge = gMC->TrackCharge();
- if (TMath::Abs(charge) > 1) {
- pp = pp * charge*charge;
- }
- }
- else {
- // Electrons above 20 Mev/c are at the plateau
- pp = kPrim * kPlateau;
- }
-
- if (pp > 0.0) {
- do {
- gMC->GetRandom()->RndmArray(1,random);
- }
- while ((random[0] == 1.0) ||
- (random[0] == 0.0));
- stepSize = - TMath::Log(random[0]) / pp;
- gMC->SetMaxStep(stepSize);
- }
-
- }
-
- }
-
- }
-
-}
-
-//_____________________________________________________________________________
-void AliTRDv1::StepManagerFixedStep()