X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDv1.cxx;h=3669e67034a2d8a24bf43193d73b48ce48ff5137;hb=3dfe174d961a34352b253c8be7d322d75e7bb5d5;hp=c845aa3015b06d6abd96bdde3b4627d338e6c529;hpb=bd0f8685bbc4d1379f15d2facafe9c61c20952a7;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDv1.cxx b/TRD/AliTRDv1.cxx index c845aa3015b..3669e67034a 100644 --- a/TRD/AliTRDv1.cxx +++ b/TRD/AliTRDv1.cxx @@ -15,95 +15,70 @@ /* $Id$ */ -/////////////////////////////////////////////////////////////////////////////// -// // -// Transition Radiation Detector version 1 -- slow simulator // -// // -//Begin_Html -/* - -*/ -//End_Html -// // -// // -/////////////////////////////////////////////////////////////////////////////// - -#include - -#include +//////////////////////////////////////////////////////////////////////////// +// // +// Transition Radiation Detector version 1 -- slow simulator // +// // +//////////////////////////////////////////////////////////////////////////// + #include #include #include -#include #include +#include +#include +#include -#include "AliConst.h" -#include "AliLog.h" +#include "AliTrackReference.h" #include "AliMC.h" #include "AliRun.h" +#include "AliGeomManager.h" + #include "AliTRDgeometry.h" -#include "AliTRDhit.h" -#include "AliTRDsim.h" +#include "AliTRDCommonParam.h" +#include "AliTRDsimTR.h" #include "AliTRDv1.h" ClassImp(AliTRDv1) //_____________________________________________________________________________ -AliTRDv1::AliTRDv1():AliTRD() +AliTRDv1::AliTRDv1() + :AliTRD() + ,fTRon(kTRUE) + ,fTR(NULL) + ,fStepSize(0) + ,fWion(0) { // // Default constructor // - fSensSelect = 0; - fSensPlane = -1; - fSensChamber = -1; - fSensSector = -1; - fSensSectorRange = 0; - - fDeltaE = NULL; - fDeltaG = NULL; - fTR = NULL; - fTRon = kFALSE; - - fStepSize = 0.1; - fTypeOfStepManager = 1; - } //_____________________________________________________________________________ AliTRDv1::AliTRDv1(const char *name, const char *title) - :AliTRD(name, title) + :AliTRD(name,title) + ,fTRon(kTRUE) + ,fTR(NULL) + ,fStepSize(0.1) + ,fWion(0) { // // Standard constructor for Transition Radiation Detector version 1 // - fSensSelect = 0; - fSensPlane = -1; - fSensChamber = -1; - fSensSector = -1; - fSensSectorRange = 0; - - fDeltaE = NULL; - fDeltaG = NULL; - fTR = NULL; - fTRon = kTRUE; - fStepSize = 0.1; - fTypeOfStepManager = 1; - SetBufferSize(128000); -} - -//_____________________________________________________________________________ -AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd) -{ - // - // Copy constructor - // - - ((AliTRDv1 &) trd).Copy(*this); + if (AliTRDCommonParam::Instance()->IsXenon()) { + fWion = 23.53; // Ionization energy XeCO2 (85/15) + } + else if (AliTRDCommonParam::Instance()->IsArgon()) { + fWion = 27.21; // Ionization energy ArCO2 (82/18) + } + else { + AliFatal("Wrong gas mixture"); + exit(1); + } } @@ -114,46 +89,128 @@ AliTRDv1::~AliTRDv1() // AliTRDv1 destructor // - if (fDeltaE) delete fDeltaE; - if (fDeltaG) delete fDeltaG; - if (fTR) delete fTR; + if (fTR) { + delete fTR; + fTR = 0; + } } //_____________________________________________________________________________ -AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd) +void AliTRDv1::AddAlignableVolumes() const { // - // Assignment operator + // Create entries for alignable volumes associating the symbolic volume + // name with the corresponding volume path. Needs to be syncronized with + // eventual changes in the geometry. // - if (this != &trd) ((AliTRDv1 &) trd).Copy(*this); - return *this; + TString volPath; + TString symName; + + TString vpStr = "ALIC_1/B077_1/BSEGMO"; + TString vpApp1 = "_1/BTRD"; + TString vpApp2 = "_1"; + TString vpApp3a = "/UTR1_1/UTS1_1/UTI1_1/UT"; + TString vpApp3b = "/UTR2_1/UTS2_1/UTI2_1/UT"; + TString vpApp3c = "/UTR3_1/UTS3_1/UTI3_1/UT"; + + TString snStr = "TRD/sm"; + TString snApp1 = "/st"; + TString snApp2 = "/pl"; -} - -//_____________________________________________________________________________ -void AliTRDv1::Copy(TObject &trd) const -{ - printf("void AliTRDv1::Copy(TObject &trd) const\n"); // - // Copy function + // The super modules + // The symbolic names are: TRD/sm00 + // ... + // TRD/sm17 // + for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) { - ((AliTRDv1 &) trd).fSensSelect = fSensSelect; - ((AliTRDv1 &) trd).fSensPlane = fSensPlane; - ((AliTRDv1 &) trd).fSensChamber = fSensChamber; - ((AliTRDv1 &) trd).fSensSector = fSensSector; - ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange; + volPath = vpStr; + volPath += isector; + volPath += vpApp1; + volPath += isector; + volPath += vpApp2; - ((AliTRDv1 &) trd).fTypeOfStepManager = fTypeOfStepManager; - ((AliTRDv1 &) trd).fStepSize = fStepSize; + symName = snStr; + symName += Form("%02d",isector); - ((AliTRDv1 &) trd).fTRon = fTRon; + gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data()); + + } - fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE); - fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG); - fTR->Copy(*((AliTRDv1 &) trd).fTR); + // + // The readout chambers + // The symbolic names are: TRD/sm00/st0/pl0 + // ... + // TRD/sm17/st4/pl5 + // + AliGeomManager::ELayerID idTRD1 = AliGeomManager::kTRD1; + Int_t layer, modUID; + + for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) { + + if (fGeometry->GetSMstatus(isector) == 0) continue; + + for (Int_t istack = 0; istack < AliTRDgeometry::Nstack(); istack++) { + for (Int_t ilayer = 0; ilayer < AliTRDgeometry::Nlayer(); ilayer++) { + + layer = idTRD1 + ilayer; + modUID = AliGeomManager::LayerToVolUIDSafe(layer,isector*5+istack); + + Int_t idet = AliTRDgeometry::GetDetectorSec(ilayer,istack); + + volPath = vpStr; + volPath += isector; + volPath += vpApp1; + volPath += isector; + volPath += vpApp2; + switch (isector) { + case 13: + case 14: + case 15: + if (istack == 2) { + continue; + } + volPath += vpApp3c; + break; + case 11: + case 12: + volPath += vpApp3b; + break; + default: + volPath += vpApp3a; + }; + volPath += Form("%02d",idet); + volPath += vpApp2; + + symName = snStr; + symName += Form("%02d",isector); + symName += snApp1; + symName += istack; + symName += snApp2; + symName += ilayer; + + TGeoPNEntry *alignableEntry = + gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data(),modUID); + + // Add the tracking to local matrix following the TPC example + if (alignableEntry) { + TGeoHMatrix *globMatrix = alignableEntry->GetGlobalOrig(); + Double_t sectorAngle = 20.0 * (isector % 18) + 10.0; + TGeoHMatrix *t2lMatrix = new TGeoHMatrix(); + t2lMatrix->RotateZ(sectorAngle); + t2lMatrix->MultiplyLeft(&(globMatrix->Inverse())); + alignableEntry->SetMatrix(t2lMatrix); + } + else { + AliError(Form("Alignable entry %s is not valid!",symName.Data())); + } + + } + } + } } @@ -167,7 +224,10 @@ void AliTRDv1::CreateGeometry() // Check that FRAME is there otherwise we have no place where to put the TRD AliModule* frame = gAlice->GetModule("FRAME"); - if (!frame) return; + if (!frame) { + AliError("TRD needs FRAME to be present\n"); + return; + } // Define the chambers AliTRD::CreateGeometry(); @@ -196,88 +256,89 @@ void AliTRDv1::CreateTRhit(Int_t det) // volume. // - // PDG code electron - const Int_t kPdgElectron = 11; - - // Ionization energy - const Float_t kWion = 23.53; - // Maximum number of TR photons per track const Int_t kNTR = 50; - TLorentzVector mom, pos; - - // Create TR at the entrance of the chamber - if (gMC->IsTrackEntering()) { - - // Create TR only for electrons - Int_t iPdg = gMC->TrackPid(); - if (TMath::Abs(iPdg) != kPdgElectron) return; + TLorentzVector mom; + TLorentzVector pos; - Float_t eTR[kNTR]; - Int_t nTR; + Float_t eTR[kNTR]; + Int_t nTR; - // Create TR photons - gMC->TrackMomentum(mom); - Float_t pTot = mom.Rho(); - fTR->CreatePhotons(iPdg,pTot,nTR,eTR); - if (nTR > kNTR) { - AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR)); - } + // Create TR photons + gMC->TrackMomentum(mom); + Float_t pTot = mom.Rho(); + fTR->CreatePhotons(11,pTot,nTR,eTR); + if (nTR > kNTR) { + AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR)); + } - // Loop through the TR photons - for (Int_t iTR = 0; iTR < nTR; iTR++) { + // Loop through the TR photons + for (Int_t iTR = 0; iTR < nTR; iTR++) { - Float_t energyMeV = eTR[iTR] * 0.001; - Float_t energyeV = eTR[iTR] * 1000.0; - Float_t absLength = 0; - Float_t sigma = 0; + Float_t energyMeV = eTR[iTR] * 0.001; + Float_t energyeV = eTR[iTR] * 1000.0; + Float_t absLength = 0.0; + Float_t sigma = 0.0; - // Take the absorbtion in the entrance window into account - Double_t muMy = fTR->GetMuMy(energyMeV); - sigma = muMy * fFoilDensity; - if (sigma > 0.0) { - absLength = gRandom->Exp(1.0/sigma); - if (absLength < AliTRDgeometry::MyThick()) continue; - } - else { + // Take the absorbtion in the entrance window into account + Double_t muMy = fTR->GetMuMy(energyMeV); + sigma = muMy * fFoilDensity; + if (sigma > 0.0) { + absLength = gRandom->Exp(1.0/sigma); + if (absLength < AliTRDgeometry::MyThick()) { continue; } + } + else { + continue; + } - // The absorbtion cross sections in the drift gas - // Gas-mixture (Xe/CO2) - Double_t muXe = fTR->GetMuXe(energyMeV); - Double_t muCO = fTR->GetMuCO(energyMeV); - sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp(); - - // The distance after which the energy of the TR photon - // is deposited. - if (sigma > 0.0) { - absLength = gRandom->Exp(1.0/sigma); - if (absLength > (AliTRDgeometry::DrThick() - + AliTRDgeometry::AmThick())) { - continue; - } - } - else { + // The absorbtion cross sections in the drift gas + // Gas-mixture (Xe/CO2) + Double_t muNo = 0.0; + if (AliTRDCommonParam::Instance()->IsXenon()) { + muNo = fTR->GetMuXe(energyMeV); + } + else if (AliTRDCommonParam::Instance()->IsArgon()) { + muNo = fTR->GetMuAr(energyMeV); + } + Double_t muCO = fTR->GetMuCO(energyMeV); + sigma = (fGasNobleFraction * muNo + (1.0 - fGasNobleFraction) * muCO) + * fGasDensity + * fTR->GetTemp(); + + // The distance after which the energy of the TR photon + // is deposited. + if (sigma > 0.0) { + absLength = gRandom->Exp(1.0/sigma); + if (absLength > (AliTRDgeometry::DrThick() + + AliTRDgeometry::AmThick())) { continue; } - - // The position of the absorbtion - Float_t posHit[3]; - gMC->TrackPosition(pos); - posHit[0] = pos[0] + mom[0] / pTot * absLength; - posHit[1] = pos[1] + mom[1] / pTot * absLength; - posHit[2] = pos[2] + mom[2] / pTot * absLength; - - // Create the charge - Int_t q = ((Int_t) (energyeV / kWion)); - - // Add the hit to the array. TR photon hits are marked - // by negative charge - AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE); - } + else { + continue; + } + + // The position of the absorbtion + Float_t posHit[3]; + gMC->TrackPosition(pos); + posHit[0] = pos[0] + mom[0] / pTot * absLength; + posHit[1] = pos[1] + mom[1] / pTot * absLength; + posHit[2] = pos[2] + mom[2] / pTot * absLength; + + // Create the charge + Int_t q = ((Int_t) (energyeV / fWion)); + + // Add the hit to the array. TR photon hits are marked + // by negative charge + AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber() + ,det + ,posHit + ,-q + ,gMC->TrackTime()*1.0e06 + ,kTRUE); } @@ -293,608 +354,37 @@ void AliTRDv1::Init() AliTRD::Init(); AliDebug(1,"Slow simulator\n"); - if (fSensSelect) { - if (fSensPlane >= 0) - AliInfo(Form("Only plane %d is sensitive")); - if (fSensChamber >= 0) - AliInfo(Form("Only chamber %d is sensitive",fSensChamber)); - if (fSensSector >= 0) { - Int_t sens1 = fSensSector; - Int_t sens2 = fSensSector + fSensSectorRange; - sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) - * AliTRDgeometry::Nsect(); - AliInfo(Form("Only sectors %d - %d are sensitive\n",sens1,sens2-1)); - } - } // Switch on TR simulation as default if (!fTRon) { AliInfo("TR simulation off"); } else { - fTR = new AliTRDsim(); + fTR = new AliTRDsimTR(); } - // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2) - const Float_t kPoti = 12.1; - // Maximum energy (50 keV); - const Float_t kEend = 50000.0; - // Ermilova distribution for the delta-ray spectrum - Float_t poti = TMath::Log(kPoti); - Float_t eEnd = TMath::Log(kEend); - - // Ermilova distribution for the delta-ray spectrum - fDeltaE = new TF1("deltae" ,Ermilova ,poti,eEnd,0); - - // Geant3 distribution for the delta-ray spectrum - fDeltaG = new TF1("deltag",IntSpecGeant,2.421257,28.536469,0); - AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); } -//_____________________________________________________________________________ -void AliTRDv1::SetSensPlane(Int_t iplane) -{ - // - // Defines the hit-sensitive plane (0-5) - // - - if ((iplane < 0) || (iplane > 5)) { - AliWarning(Form("Wrong input value:%d",iplane)); - AliWarning("Use standard setting"); - fSensPlane = -1; - fSensSelect = 0; - return; - } - - fSensSelect = 1; - fSensPlane = iplane; - -} - -//_____________________________________________________________________________ -void AliTRDv1::SetSensChamber(Int_t ichamber) -{ - // - // Defines the hit-sensitive chamber (0-4) - // - - if ((ichamber < 0) || (ichamber > 4)) { - AliWarning(Form("Wrong input value: %d",ichamber)); - AliWarning("Use standard setting"); - fSensChamber = -1; - fSensSelect = 0; - return; - } - - fSensSelect = 1; - fSensChamber = ichamber; - -} - -//_____________________________________________________________________________ -void AliTRDv1::SetSensSector(Int_t isector) -{ - // - // Defines the hit-sensitive sector (0-17) - // - - SetSensSector(isector,1); - -} - -//_____________________________________________________________________________ -void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector) -{ - // - // Defines a range of hit-sensitive sectors. The range is defined by - // (0-17) as the starting point and as the number - // of sectors to be included. - // - - if ((isector < 0) || (isector > 17)) { - AliWarning(Form("Wrong input value : %d",isector)); - AliWarning("Use standard setting"); - fSensSector = -1; - fSensSectorRange = 0; - fSensSelect = 0; - return; - } - - if ((nsector < 1) || (nsector > 18)) { - AliWarning(Form("Wrong input value : %d",nsector)); - AliWarning("Use standard setting"); - fSensSector = -1; - fSensSectorRange = 0; - fSensSelect = 0; - return; - } - - fSensSelect = 1; - fSensSector = isector; - fSensSectorRange = nsector; - -} - //_____________________________________________________________________________ void AliTRDv1::StepManager() -{ - // - // 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 = 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 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=((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; - } - - Int_t nsteps = 0; - do {nsteps = gRandom->Poisson(pp);} while(!nsteps); - stepSize = 1./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; - 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 = 23.53; // 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() { // // Slow simulator. Every charged track produces electron cluster as hits // along its path across the drift volume. The step size is fixed in // this version of the step manager. // + // Works for Xe/CO2 as well as Ar/CO2 + // + + // PDG code electron + const Int_t kPdgElectron = 11; - Int_t pla = 0; - Int_t cha = 0; - Int_t sec = 0; - Int_t det = 0; + Int_t layer = 0; + Int_t stack = 0; + Int_t sector = 0; + Int_t det = 0; Int_t qTot; Float_t hits[3]; @@ -903,39 +393,50 @@ void AliTRDv1::StepManagerFixedStep() 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, mom; + cIdChamber[2] = 0; - const Int_t kNplan = AliTRDgeometry::Nplan(); - const Int_t kNcham = AliTRDgeometry::Ncham(); - const Int_t kNdetsec = kNplan * kNcham; + TLorentzVector pos; + TLorentzVector mom; - const Double_t kBig = 1.0E+12; + const Int_t kNlayer = AliTRDgeometry::Nlayer(); + const Int_t kNstack = AliTRDgeometry::Nstack(); + const Int_t kNdetsec = kNlayer * kNstack; - const Float_t kWion = 23.53; // Ionization energy + const Double_t kBig = 1.0e+12; const Float_t kEkinMinStep = 1.0e-5; // Minimum energy for the step size adjustment // Set the maximum step size to a very large number for all // neutral particles and those outside the driftvolume - gMC->SetMaxStep(kBig); + if (!fPrimaryIonisation) gMC->SetMaxStep(kBig); // If not charged track or already stopped or disappeared, just return. if ((!gMC->TrackCharge()) || - gMC->IsTrackStop() || - gMC->IsTrackDisappeared()) return; + gMC->IsTrackDisappeared()) { + return; + } // Inside a sensitive volume? cIdCurrent = gMC->CurrentVolName(); - if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE; - if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE; + if (cIdSensDr == cIdCurrent[1]) { + drRegion = kTRUE; + } + if (cIdSensAm == cIdCurrent[1]) { + amRegion = kTRUE; + } - if ((!drRegion) && (!amRegion)) return; + if ((!drRegion) && + (!amRegion)) { + return; + } // The hit coordinates and charge gMC->TrackPosition(pos); @@ -943,373 +444,74 @@ void AliTRDv1::StepManagerFixedStep() 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 += 270.; - else phi -= 90.; - sec = ((Int_t) (phi / 20.)); + // The sector number (0 - 17), according to standard coordinate system + cIdPath = gGeoManager->GetPath(); + cIdSector[0] = cIdPath[21]; + cIdSector[1] = cIdPath[22]; + sector = atoi(cIdSector); // The plane and chamber number - cIdChamber[0] = cIdCurrent[2]; - cIdChamber[1] = cIdCurrent[3]; + 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; - } - } - } + stack = ((Int_t) idChamber / kNlayer); + layer = ((Int_t) idChamber % kNlayer); - if (!addthishit) return; + // The detector number + det = fGeometry->GetDetector(layer,stack,sector); - det = fGeometry->GetDetector(pla,cha,sec); // The detector number - - Int_t trkStat = 0; // 0: InFlight 1:Entering 2:Exiting + // 0: InFlight 1:Entering 2:Exiting + Int_t trkStat = 0; // 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 ((drRegion) && + (gMC->IsTrackEntering())) { - if (gMC->IsTrackEntering()) { - gMC->TrackMomentum(mom); - AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber()); - trkStat = 1; - } - if (gMC->IsTrackExiting()) { - gMC->TrackMomentum(mom); - AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber()); - trkStat = 2; + // 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(), AliTrackReference::kTRD); + trkStat = 1; + + // Create the hits from TR photons if electron/positron is + // entering the drift volume + if ((fTR) && + (fTRon) && + (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) { + CreateTRhit(det); } - // Create the hits from TR photons - if (fTR) 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(), AliTrackReference::kTRD); + trkStat = 2; } // Calculate the charge according to GEANT Edep // Create a new dEdx hit eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09; - qTot = (Int_t) (eDep / kWion); - AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber() - ,det,hits,qTot,drRegion); + qTot = (Int_t) (eDep / fWion); + if ((qTot) || + (trkStat)) { + AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber() + ,det + ,hits + ,qTot + ,gMC->TrackTime()*1.0e06 + ,drRegion); + } // Set Maximum Step Size // Produce only one hit if Ekin is below cutoff - if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return; - gMC->SetMaxStep(fStepSize); - -} - -//_____________________________________________________________________________ -Double_t AliTRDv1::BetheBloch(Double_t bg) -{ - // - // Parametrization of the Bethe-Bloch-curve - // The parametrization is the same as for the TPC and is taken from Lehrhaus. - // - - // This parameters have been adjusted to averaged values from GEANT - const Double_t kP1 = 7.17960e-02; - const Double_t kP2 = 8.54196; - const Double_t kP3 = 1.38065e-06; - const Double_t kP4 = 5.30972; - const Double_t kP5 = 2.83798; - - // This parameters have been adjusted to Xe-data found in: - // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253 - //const Double_t kP1 = 0.76176E-1; - //const Double_t kP2 = 10.632; - //const Double_t kP3 = 3.17983E-6; - //const Double_t kP4 = 1.8631; - //const Double_t kP5 = 1.9479; - - // Lower cutoff of the Bethe-Bloch-curve to limit step sizes - const Double_t kBgMin = 0.8; - const Double_t kBBMax = 6.83298; - //const Double_t kBgMin = 0.6; - //const Double_t kBBMax = 17.2809; - //const Double_t kBgMin = 0.4; - //const Double_t kBBMax = 82.0; - - if (bg > kBgMin) { - Double_t yy = bg / TMath::Sqrt(1. + bg*bg); - Double_t aa = TMath::Power(yy,kP4); - Double_t bb = TMath::Power((1./bg),kP5); - bb = TMath::Log(kP3 + bb); - return ((kP2 - aa - bb)*kP1 / aa); - } - else { - return kBBMax; + if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) { + return; } - -} - -//_____________________________________________________________________________ -Double_t AliTRDv1::BetheBlochGeant(Double_t bg) -{ - // - // Return dN/dx (number of primary collisions per centimeter) - // for given beta*gamma factor. - // - // Implemented by K.Oyama according to GEANT 3 parametrization shown in - // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html - // This must be used as a set with IntSpecGeant. - // - - Double_t arr_g[20] = { - 1.100000, 1.200000, 1.300000, 1.500000, - 1.800000, 2.000000, 2.500000, 3.000000, - 4.000000, 7.000000, 10.000000, 20.000000, - 40.000000, 70.000000, 100.000000, 300.000000, - 600.000000, 1000.000000, 3000.000000, 10000.000000 }; - - Double_t arr_nc[20] = { - 75.009056, 45.508083, 35.299252, 27.116327, - 22.734999, 21.411915, 19.934095, 19.449375, - 19.344431, 20.185553, 21.027925, 22.912676, - 24.933352, 26.504053, 27.387468, 29.566597, - 30.353779, 30.787134, 31.129285, 31.157350 }; - - // betagamma to gamma - Double_t g = TMath::Sqrt( 1. + bg*bg ); - - // Find the index just before the point we need. - int i; - for( i = 0 ; i < 18 ; i++ ) - if( arr_g[i] < g && arr_g[i+1] > g ) - break; - - // Simple interpolation. - Double_t pp = ((arr_nc[i+1] - arr_nc[i]) / - (arr_g[i+1]-arr_g[i])) * (g-arr_g[i]) + arr_nc[i]; - - return pp; //arr_nc[8]; - -} - -//_____________________________________________________________________________ -void AliTRDv1::Stepping() -{ -// Stepping info -// --- - - cout << "X(cm) " - << "Y(cm) " - << "Z(cm) " - << "KinE(MeV) " - << "dE(MeV) " - << "Step(cm) " - << "TrackL(cm) " - << "Volume " - << "Process " - << endl; - - // Position - // - Double_t x, y, z; - gMC->TrackPosition(x, y, z); - cout << setw(8) << setprecision(3) << x << " " - << setw(8) << setprecision(3) << y << " " - << setw(8) << setprecision(3) << z << " "; - - // Kinetic energy - // - Double_t px, py, pz, etot; - gMC->TrackMomentum(px, py, pz, etot); - Double_t ekin = etot - gMC->TrackMass(); - cout << setw(9) << setprecision(4) << ekin*1e03 << " "; - - // Energy deposit - // - cout << setw(9) << setprecision(4) << gMC->Edep()*1e03 << " "; - - // Step length - // - cout << setw(8) << setprecision(3) << gMC->TrackStep() << " "; - - // Track length - // - cout << setw(8) << setprecision(3) << gMC->TrackLength() << " "; - - // Volume - // - if (gMC->CurrentVolName() != 0) - cout << setw(4) << gMC->CurrentVolName() << " "; - else - cout << setw(4) << "None" << " "; - - // Process - // - TArrayI processes; - Int_t nofProcesses = gMC->StepProcesses(processes); - for(int ip=0;ip 0); - pos2--; - if (pos2 > kNv) pos2 = kNv - 1; - pos1 = pos2 - 1; - - // Differentiate between the sampling points - dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]); - - return dnde; - -} - -//_____________________________________________________________________________ -Double_t IntSpecGeant(Double_t *x, Double_t *) -{ - // - // Integrated spectrum from Geant3 - // - - const Int_t npts = 83; - Double_t arre[npts] = { - 2.421257, 2.483278, 2.534301, 2.592230, - 2.672067, 2.813299, 3.015059, 3.216819, - 3.418579, 3.620338, 3.868209, 3.920198, - 3.978284, 4.063923, 4.186264, 4.308605, - 4.430946, 4.553288, 4.724261, 4.837736, - 4.999842, 5.161949, 5.324056, 5.486163, - 5.679688, 5.752998, 5.857728, 5.962457, - 6.067185, 6.171914, 6.315653, 6.393674, - 6.471694, 6.539689, 6.597658, 6.655627, - 6.710957, 6.763648, 6.816338, 6.876198, - 6.943227, 7.010257, 7.106285, 7.252151, - 7.460531, 7.668911, 7.877290, 8.085670, - 8.302979, 8.353585, 8.413120, 8.483500, - 8.541030, 8.592857, 8.668865, 8.820485, - 9.037086, 9.253686, 9.470286, 9.686887, - 9.930838, 9.994655, 10.085822, 10.176990, - 10.268158, 10.359325, 10.503614, 10.627565, - 10.804637, 10.981709, 11.158781, 11.335854, - 11.593397, 11.781165, 12.049404, 12.317644, - 12.585884, 12.854123, 14.278421, 16.975889, - 20.829416, 24.682943, 28.536469 - }; - /* - Double_t arrdndx[npts] = { - 19.344431, 18.664679, 18.136106, 17.567745, - 16.836426, 15.677382, 14.281277, 13.140237, - 12.207677, 11.445510, 10.697049, 10.562296, - 10.414673, 10.182341, 9.775256, 9.172330, - 8.240271, 6.898587, 4.808303, 3.889751, - 3.345288, 3.093431, 2.897347, 2.692470, - 2.436222, 2.340029, 2.208579, 2.086489, - 1.975535, 1.876519, 1.759626, 1.705024, - 1.656374, 1.502638, 1.330566, 1.200697, - 1.101168, 1.019323, 0.943867, 0.851951, - 0.755229, 0.671576, 0.570675, 0.449672, - 0.326722, 0.244225, 0.188225, 0.149608, - 0.121529, 0.116289, 0.110636, 0.103490, - 0.096147, 0.089191, 0.079780, 0.063927, - 0.047642, 0.036341, 0.028250, 0.022285, - 0.017291, 0.016211, 0.014802, 0.013533, - 0.012388, 0.011352, 0.009803, 0.008537, - 0.007039, 0.005829, 0.004843, 0.004034, - 0.003101, 0.002564, 0.001956, 0.001494, - 0.001142, 0.000873, 0.000210, 0.000014, - 0.000000, 0.000000, 0.000000 - }; - */ - // Differentiate - // dnde = (arrdndx[i-1] - arrdndx[i]) / (arre[i] - arre[i-1]); - - Double_t arrdnde[npts] = { - 10.960000, 10.960000, 10.359500, 9.811340, - 9.1601500, 8.206670, 6.919630, 5.655430, - 4.6221300, 3.777610, 3.019560, 2.591950, - 2.5414600, 2.712920, 3.327460, 4.928240, - 7.6185300, 10.966700, 12.225800, 8.094750, - 3.3586900, 1.553650, 1.209600, 1.263840, - 1.3241100, 1.312140, 1.255130, 1.165770, - 1.0594500, 0.945450, 0.813231, 0.699837, - 0.6235580, 2.260990, 2.968350, 2.240320, - 1.7988300, 1.553300, 1.432070, 1.535520, - 1.4429900, 1.247990, 1.050750, 0.829549, - 0.5900280, 0.395897, 0.268741, 0.185320, - 0.1292120, 0.103545, 0.0949525, 0.101535, - 0.1276380, 0.134216, 0.123816, 0.104557, - 0.0751843, 0.0521745, 0.0373546, 0.0275391, - 0.0204713, 0.0169234, 0.0154552, 0.0139194, - 0.0125592, 0.0113638, 0.0107354, 0.0102137, - 0.00845984, 0.00683338, 0.00556836, 0.00456874, - 0.0036227, 0.00285991, 0.00226664, 0.00172234, - 0.00131226, 0.00100284, 0.000465492, 7.26607e-05, - 3.63304e-06, 0.0000000, 0.0000000 - }; - - Int_t i; - Double_t energy = x[0]; - - for( i = 0 ; i < npts ; i++ ) - if( energy < arre[i] ) break; - - if( i == 0 ) - AliErrorGeneral("AliTRDv1","Given energy value is too small or zero"); - - return arrdnde[i]; + if (!fPrimaryIonisation) gMC->SetMaxStep(fStepSize); }