-{
- //
- // Slow simulator. Every charged track produces electron cluster as hits
- // along its path across the drift volume.
- //
-
- switch (fTypeOfStepManager) {
- case 0 : StepManagerErmilova(); break; // 0 is Ermilova
- case 1 : StepManagerGeant(); break; // 1 is Geant
- case 2 : StepManagerFixedStep(); break; // 2 is fixed step
- 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
- //
-
-/* if (t == 1) {
- AliWarning("Sorry, Geant parametrization step manager is not implemented yet. Please ask K.Oyama for detail.");
- }
-*/
- 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.
- //
- 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, pp;
- Double_t stepSize=0;
-
- Bool_t drRegion = kFALSE;
- Bool_t amRegion = kFALSE;
-
- TString cIdCurrent;
- TString cIdSensDr = "J";
- TString cIdSensAm = "K";
- Char_t cIdChamber[3];
- cIdChamber[2] = 0;
-
- TLorentzVector pos, 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 = 22.57; // 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 ra=5.37E-4, rb=0.9815, rc=3.123E-3;
- // Gas density -> To be made user adjustable !
- const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
-
- // Plateau value of the energy-loss for electron in xenon
- // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
- //const Double_t kPlateau = 1.70;
- // the averaged value (26/3/99)
- const Float_t kPlateau = 1.55;
-
- const Float_t kPrim = 19.34; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
- // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
- const Float_t kPoti = 12.1;
-
- const Int_t kPdgElectron = 11; // PDG code electron
-
- // 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->IsTrackStop() ) &&
- (!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)
- // The numbering goes clockwise and starts at y = 0
- Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
- if (phi < 90.)
- phi = phi + 270.;
- else
- phi = phi - 90.;
- sec = ((Int_t) (phi / 20));
-
- // 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);
-
- // Check on selected volumes
- Int_t addthishit = 1;
- if (fSensSelect) {
- if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
- if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
- if (fSensSector >= 0) {
- Int_t sens1 = fSensSector;
- Int_t sens2 = fSensSector + fSensSectorRange;
- sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
- * AliTRDgeometry::Nsect();
- if (sens1 < sens2) {
- if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
- }
- else {
- if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
- }
- }
- }
-
- // Add this hit
- if (addthishit) {
-
- // The detector number
- det = fGeometry->GetDetector(pla,cha,sec);
-
- // Special hits only in the drift region
- if (drRegion) {
- // Create a track reference at the entrance and
- // exit of each chamber that contain the
- // momentum components of the particle
- if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
- gMC->TrackMomentum(mom);
- AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
- }
-
- if (gMC->IsTrackEntering() && !gMC->IsNewTrack()) {
- // determine if hit belong to primary track
- fPrimaryTrackPid=gAlice->GetMCApp()->GetCurrentTrackNumber();
- //determine track length when entering the detector
- fTrackLength0=gMC->TrackLength();
- }
-
- // Create the hits from TR photons
- if (fTR) CreateTRhit(det);
- }
-
- // 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
- TArrayI processes;
- gMC->StepProcesses(processes);
- int nofprocesses=processes.GetSize(), 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 pr_range=0.;
- float range=gMC->TrackLength()-fTrackLength0;
- // merge GEANT tracker information with localy cooked one
- if(gAlice->GetMCApp()->GetCurrentTrackNumber()==fPrimaryTrackPid) {
-// printf("primary pid=%d eDelta=%f\n",pid,eDelta);
- if(pid==27){
- if(eDelta>=kECut){
- pr_range=ra*eDelta*.001*(1.-rb/(1.+rc*eDelta*0.001))/rho;
- if(pr_range>=(3.7-range)) eDelta*=.1;
- }
- } else if(pid==1){
- if(eDelta<kECut) eDelta*=.5;
- else {
- pr_range=ra*eDelta*.001*(1.-rb/(1.+rc*eDelta*0.001))/rho;
- if(pr_range>=((AliTRDgeometry::DrThick()
- + AliTRDgeometry::AmThick())-range)) eDelta*=.05;
- else eDelta*=.5;
- }
- } else eDelta=0.;
- } else eDelta=0.;
-
- // Generate the electron cluster size
- if(eDelta==0.) qTot=0;
- else qTot = ((Int_t) (eDelta / kWion) + 1);
- // Create a new dEdx hit
- AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot, 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;
- }
-
- stepSize = 1./gRandom->Poisson(pp);
- 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;
- Double_t eDelta;
- Double_t betaGamma, pp;
- Double_t stepSize;
-
- Bool_t drRegion = kFALSE;
- Bool_t amRegion = kFALSE;
-
- TString cIdCurrent;
- TString cIdSensDr = "J";
- TString cIdSensAm = "K";
- Char_t cIdChamber[3];
- cIdChamber[2] = 0;
-
- TLorentzVector pos, 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 = 22.57; // Ionization energy
- const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
-
- // energy threshold for production of delta electrons
- //const Float_t kECut = 1.0e4;
- // Parameters entering the parametrized range for delta electrons
- //const float ra=5.37E-4, rb=0.9815, rc=3.123E-3;
- // Gas density -> To be made user adjustable !
- //const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
-
- // Minimum energy for the step size adjustment
- const Float_t kEkinMinStep = 1.0e-5;
-
- // Plateau value of the energy-loss for electron in xenon
- // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
- //const Double_t kPlateau = 1.70;
- // the averaged value (26/3/99)
- const Float_t kPlateau = 1.55;
-
- const Float_t kPrim = 48.0; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
- // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
- const Float_t kPoti = 12.1;
-
- const Int_t kPdgElectron = 11; // PDG code electron
-
- // 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->IsTrackStop() ) &&
- (!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)
- // The numbering goes clockwise and starts at y = 0
- Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
- if (phi < 90.)
- phi = phi + 270.;
- else
- phi = phi - 90.;
- sec = ((Int_t) (phi / 20));
-
- // 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);
-
- // Check on selected volumes
- Int_t addthishit = 1;
- if (fSensSelect) {
- if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
- if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
- if (fSensSector >= 0) {
- Int_t sens1 = fSensSector;
- Int_t sens2 = fSensSector + fSensSectorRange;
- sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
- * AliTRDgeometry::Nsect();
- if (sens1 < sens2) {
- if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
- }
- else {
- if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
- }
- }
- }
-
- // Add this hit
- if (addthishit) {
-
- // The detector number
- det = fGeometry->GetDetector(pla,cha,sec);
-
- // Special hits only in the drift region
- if (drRegion) {
-
- // Create a track reference at the entrance and
- // exit of each chamber that contain the
- // momentum components of the particle
- if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
- gMC->TrackMomentum(mom);
- AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
- }
- // Create the hits from TR photons
- if (fTR) CreateTRhit(det);
- }
-
-
- // 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.) qTot=0;
- else qTot = ((Int_t) (eDelta / kWion) + 1);
-
- // Create a new dEdx hit
- if (drRegion) {
- AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
- ,det,hits,qTot, kTRUE);
- }
- else {
- AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
- ,det,hits,qTot,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) {
- do
- gMC->GetRandom()->RndmArray(1, random);
- while ((random[0] == 1.) || (random[0] == 0.));
- stepSize = - TMath::Log(random[0]) / pp;
- gMC->SetMaxStep(stepSize);
- }
- }
- }
- }
- }
-}
-
-//_____________________________________________________________________________
-void AliTRDv1::StepManagerFixedStep()