X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDv1.cxx;h=c845aa3015b06d6abd96bdde3b4627d338e6c529;hb=4a378d7471c111b7380f83923d1b7cbbe845bc1f;hp=a09b25beb58d732627ba021c39e6e26ef564f890;hpb=b9d0a01d7a0723a09071b0b56200d72f59a9c2b6;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDv1.cxx b/TRD/AliTRDv1.cxx index a09b25beb58..c845aa3015b 100644 --- a/TRD/AliTRDv1.cxx +++ b/TRD/AliTRDv1.cxx @@ -13,105 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.33.6.2 2002/07/24 10:09:31 alibrary -Updating VirtualMC - -Revision 1.34 2002/06/13 08:11:56 cblume -Add the track references - -Revision 1.33 2002/02/20 14:01:40 hristov -Compare a TString with a string, otherwise the conversion cannot be done on Sun - -Revision 1.32 2002/02/13 16:58:37 cblume -Bug fix reported by Jiri. Make atoi input zero terminated in StepManager() - -Revision 1.31 2002/02/11 14:25:27 cblume -Geometry update, compressed hit structure - -Revision 1.30 2001/05/21 16:45:47 hristov -Last minute changes (C.Blume) - -Revision 1.29 2001/05/16 14:57:28 alibrary -New files for folders and Stack - -Revision 1.28 2001/05/07 08:03:22 cblume -Generate also hits in the amplification region - -Revision 1.27 2001/03/30 14:40:15 cblume -Update of the digitization parameter - -Revision 1.26 2000/11/30 17:38:08 cblume -Changes to get in line with new STEER and EVGEN - -Revision 1.25 2000/11/15 14:30:16 cblume -Fixed bug in calculating detector no. of extra hit - -Revision 1.24 2000/11/10 14:58:36 cblume -Introduce additional hit with amplitude 0 at the chamber borders - -Revision 1.23 2000/11/01 14:53:21 cblume -Merge with TRD-develop - -Revision 1.17.2.5 2000/10/15 23:40:01 cblume -Remove AliTRDconst - -Revision 1.17.2.4 2000/10/06 16:49:46 cblume -Made Getters const - -Revision 1.17.2.3 2000/10/04 16:34:58 cblume -Replace include files by forward declarations - -Revision 1.17.2.2 2000/09/18 13:50:17 cblume -Include TR photon generation and adapt to new AliTRDhit - -Revision 1.22 2000/06/27 13:08:50 cblume -Changed to Copy(TObject &A) to appease the HP-compiler - -Revision 1.21 2000/06/09 11:10:07 cblume -Compiler warnings and coding conventions, next round - -Revision 1.20 2000/06/08 18:32:58 cblume -Make code compliant to coding conventions - -Revision 1.19 2000/06/07 16:27:32 cblume -Try to remove compiler warnings on Sun and HP - -Revision 1.18 2000/05/08 16:17:27 cblume -Merge TRD-develop - -Revision 1.17.2.1 2000/05/08 14:59:16 cblume -Made inline function non-virtual. Bug fix in setting sensitive chamber - -Revision 1.17 2000/02/28 19:10:26 cblume -Include the new TRD classes - -Revision 1.16.4.1 2000/02/28 18:04:35 cblume -Change to new hit version, introduce geometry class, and move digitization and clustering to AliTRDdigitizer/AliTRDclusterizerV1 - -Revision 1.16 1999/11/05 22:50:28 fca -Do not use Atan, removed from ROOT too - -Revision 1.15 1999/11/02 17:20:19 fca -initialise nbytes before using it - -Revision 1.14 1999/11/02 17:15:54 fca -Correct ansi scoping not accepted by HP compilers - -Revision 1.13 1999/11/02 17:14:51 fca -Correct ansi scoping not accepted by HP compilers - -Revision 1.12 1999/11/02 16:35:56 fca -New version of TRD introduced - -Revision 1.11 1999/11/01 20:41:51 fca -Added protections against using the wrong version of FRAME - -Revision 1.10 1999/09/29 09:24:35 fca -Introduction of the Copyright and cvs Log - -*/ +/* $Id$ */ /////////////////////////////////////////////////////////////////////////////// // // @@ -128,21 +30,21 @@ Introduction of the Copyright and cvs Log #include -#include -#include -#include #include #include +#include +#include +#include +#include -#include "AliRun.h" -#include "AliMC.h" #include "AliConst.h" - -#include "AliTRDv1.h" -#include "AliTRDhit.h" -#include "AliTRDmatrix.h" +#include "AliLog.h" +#include "AliMC.h" +#include "AliRun.h" #include "AliTRDgeometry.h" +#include "AliTRDhit.h" #include "AliTRDsim.h" +#include "AliTRDv1.h" ClassImp(AliTRDv1) @@ -153,14 +55,19 @@ AliTRDv1::AliTRDv1():AliTRD() // Default constructor // - fSensSelect = 0; - fSensPlane = -1; - fSensChamber = -1; - fSensSector = -1; - fSensSectorRange = 0; + fSensSelect = 0; + fSensPlane = -1; + fSensChamber = -1; + fSensSector = -1; + fSensSectorRange = 0; + + fDeltaE = NULL; + fDeltaG = NULL; + fTR = NULL; + fTRon = kFALSE; - fDeltaE = NULL; - fTR = NULL; + fStepSize = 0.1; + fTypeOfStepManager = 1; } @@ -172,21 +79,25 @@ AliTRDv1::AliTRDv1(const char *name, const char *title) // Standard constructor for Transition Radiation Detector version 1 // - fSensSelect = 0; - fSensPlane = -1; - fSensChamber = -1; - fSensSector = -1; - fSensSectorRange = 0; + fSensSelect = 0; + fSensPlane = -1; + fSensChamber = -1; + fSensSector = -1; + fSensSectorRange = 0; - fDeltaE = NULL; - fTR = NULL; + fDeltaE = NULL; + fDeltaG = NULL; + fTR = NULL; + fTRon = kTRUE; + fStepSize = 0.1; + fTypeOfStepManager = 1; SetBufferSize(128000); } //_____________________________________________________________________________ -AliTRDv1::AliTRDv1(const AliTRDv1 &trd) +AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd) { // // Copy constructor @@ -204,6 +115,7 @@ AliTRDv1::~AliTRDv1() // if (fDeltaE) delete fDeltaE; + if (fDeltaG) delete fDeltaG; if (fTR) delete fTR; } @@ -221,19 +133,26 @@ AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd) } //_____________________________________________________________________________ -void AliTRDv1::Copy(TObject &trd) +void AliTRDv1::Copy(TObject &trd) const { + printf("void AliTRDv1::Copy(TObject &trd) const\n"); // // Copy function // - ((AliTRDv1 &) trd).fSensSelect = fSensSelect; - ((AliTRDv1 &) trd).fSensPlane = fSensPlane; - ((AliTRDv1 &) trd).fSensChamber = fSensChamber; - ((AliTRDv1 &) trd).fSensSector = fSensSector; - ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange; + ((AliTRDv1 &) trd).fSensSelect = fSensSelect; + ((AliTRDv1 &) trd).fSensPlane = fSensPlane; + ((AliTRDv1 &) trd).fSensChamber = fSensChamber; + ((AliTRDv1 &) trd).fSensSector = fSensSector; + ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange; + + ((AliTRDv1 &) trd).fTypeOfStepManager = fTypeOfStepManager; + ((AliTRDv1 &) trd).fStepSize = fStepSize; + + ((AliTRDv1 &) trd).fTRon = fTRon; fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE); + fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG); fTR->Copy(*((AliTRDv1 &) trd).fTR); } @@ -281,7 +200,7 @@ void AliTRDv1::CreateTRhit(Int_t det) const Int_t kPdgElectron = 11; // Ionization energy - const Float_t kWion = 22.04; + const Float_t kWion = 23.53; // Maximum number of TR photons per track const Int_t kNTR = 50; @@ -303,9 +222,7 @@ void AliTRDv1::CreateTRhit(Int_t det) Float_t pTot = mom.Rho(); fTR->CreatePhotons(iPdg,pTot,nTR,eTR); if (nTR > kNTR) { - printf("AliTRDv1::CreateTRhit -- "); - printf("Boundary error: nTR = %d, kNTR = %d\n",nTR,kNTR); - exit(1); + AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR)); } // Loop through the TR photons @@ -319,41 +236,46 @@ void AliTRDv1::CreateTRhit(Int_t det) // Take the absorbtion in the entrance window into account Double_t muMy = fTR->GetMuMy(energyMeV); sigma = muMy * fFoilDensity; - absLength = gRandom->Exp(sigma); - if (absLength < AliTRDgeometry::MyThick()) continue; - - // The absorbtion cross sections in the drift gas - if (fGasMix == 1) { - // Gas-mixture (Xe/CO2) - Double_t muXe = fTR->GetMuXe(energyMeV); - Double_t muCO = fTR->GetMuCO(energyMeV); - sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity; + if (sigma > 0.0) { + absLength = gRandom->Exp(1.0/sigma); + if (absLength < AliTRDgeometry::MyThick()) continue; } else { - // Gas-mixture (Xe/Isobutane) - Double_t muXe = fTR->GetMuXe(energyMeV); - Double_t muBu = fTR->GetMuBu(energyMeV); - sigma = (0.97 * muXe + 0.03 * muBu) * fGasDensity; + 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. - absLength = gRandom->Exp(sigma); - if (absLength > AliTRDgeometry::DrThick()) continue; + if (sigma > 0.0) { + absLength = gRandom->Exp(1.0/sigma); + if (absLength > (AliTRDgeometry::DrThick() + + AliTRDgeometry::AmThick())) { + continue; + } + } + 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; + 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->CurrentTrack(),det,posHit,-q,kTRUE); + AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE); } @@ -370,25 +292,28 @@ void AliTRDv1::Init() AliTRD::Init(); - if(fDebug) printf("%s: Slow simulator\n",ClassName()); + AliDebug(1,"Slow simulator\n"); if (fSensSelect) { if (fSensPlane >= 0) - printf(" Only plane %d is sensitive\n",fSensPlane); + AliInfo(Form("Only plane %d is sensitive")); if (fSensChamber >= 0) - printf(" Only chamber %d is sensitive\n",fSensChamber); + 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(); - printf(" Only sectors %d - %d are sensitive\n",sens1,sens2-1); + AliInfo(Form("Only sectors %d - %d are sensitive\n",sens1,sens2-1)); } } - if (fTR) - printf("%s: TR simulation on\n",ClassName()); - else - printf("%s: TR simulation off\n",ClassName()); - printf("\n"); + + // Switch on TR simulation as default + if (!fTRon) { + AliInfo("TR simulation off"); + } + else { + fTR = new AliTRDsim(); + } // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2) const Float_t kPoti = 12.1; @@ -397,25 +322,14 @@ void AliTRDv1::Init() // Ermilova distribution for the delta-ray spectrum Float_t poti = TMath::Log(kPoti); Float_t eEnd = TMath::Log(kEend); - fDeltaE = new TF1("deltae",Ermilova,poti,eEnd,0); - if(fDebug) { - printf("%s: ",ClassName()); - for (Int_t i = 0; i < 80; i++) printf("*"); - printf("\n"); - } + // 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); -//_____________________________________________________________________________ -AliTRDsim *AliTRDv1::CreateTR() -{ - // - // Enables the simulation of TR - // - - fTR = new AliTRDsim(); - return fTR; + AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"); } @@ -427,8 +341,8 @@ void AliTRDv1::SetSensPlane(Int_t iplane) // if ((iplane < 0) || (iplane > 5)) { - printf("Wrong input value: %d\n",iplane); - printf("Use standard setting\n"); + AliWarning(Form("Wrong input value:%d",iplane)); + AliWarning("Use standard setting"); fSensPlane = -1; fSensSelect = 0; return; @@ -447,8 +361,8 @@ void AliTRDv1::SetSensChamber(Int_t ichamber) // if ((ichamber < 0) || (ichamber > 4)) { - printf("Wrong input value: %d\n",ichamber); - printf("Use standard setting\n"); + AliWarning(Form("Wrong input value: %d",ichamber)); + AliWarning("Use standard setting"); fSensChamber = -1; fSensSelect = 0; return; @@ -480,8 +394,8 @@ void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector) // if ((isector < 0) || (isector > 17)) { - printf("Wrong input value : %d\n",isector); - printf("Use standard setting\n"); + AliWarning(Form("Wrong input value : %d",isector)); + AliWarning("Use standard setting"); fSensSector = -1; fSensSectorRange = 0; fSensSelect = 0; @@ -489,8 +403,8 @@ void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector) } if ((nsector < 1) || (nsector > 18)) { - printf("Wrong input value : %d\n",nsector); - printf("Use standard setting\n"); + AliWarning(Form("Wrong input value : %d",nsector)); + AliWarning("Use standard setting"); fSensSector = -1; fSensSectorRange = 0; fSensSelect = 0; @@ -505,6 +419,270 @@ void AliTRDv1::SetSensSector(Int_t isector, Int_t 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 @@ -521,7 +699,7 @@ void AliTRDv1::StepManager() Int_t qTot; Float_t hits[3]; - Double_t random[1]; + Double_t random[1]; Float_t charge; Float_t aMass; @@ -542,26 +720,34 @@ void AliTRDv1::StepManager() TLorentzVector pos, mom; const Int_t kNplan = AliTRDgeometry::Nplan(); - const Double_t kBig = 1.0E+12; + 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)] - // Ionization energy - const Float_t kWion = 22.04; - // Maximum momentum for e+ e- g - const Float_t kPTotMaxEl = 0.002; // 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; - // dN1/dx|min for the gas mixture (90% Xe + 10% CO2) - const Float_t kPrim = 48.0; + + 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; - // PDG code electron - const Int_t kPdgElectron = 11; + 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 @@ -602,8 +788,8 @@ void AliTRDv1::StepManager() // The plane and chamber number cIdChamber[0] = cIdCurrent[2]; cIdChamber[1] = cIdCurrent[3]; - Int_t idChamber = atoi(cIdChamber); - cha = ((Int_t) idChamber / kNplan); + Int_t idChamber = (atoi(cIdChamber) % kNdetsec); + cha = kNcham - ((Int_t) idChamber / kNplan) - 1; pla = ((Int_t) idChamber % kNplan); // Check on selected volumes @@ -618,11 +804,11 @@ void AliTRDv1::StepManager() * AliTRDgeometry::Nsect(); if (sens1 < sens2) { if ((sec < sens1) || (sec >= sens2)) addthishit = 0; - } + } else { if ((sec < sens1) && (sec >= sens2)) addthishit = 0; - } - } + } + } } // Add this hit @@ -631,7 +817,7 @@ void AliTRDv1::StepManager() // The detector number det = fGeometry->GetDetector(pla,cha,sec); - // Special hits and TR photons only in the drift region + // Special hits only in the drift region if (drRegion) { // Create a track reference at the entrance and @@ -639,28 +825,29 @@ void AliTRDv1::StepManager() // momentum components of the particle if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) { gMC->TrackMomentum(mom); - AddTrackReference(gAlice->CurrentTrack(),mom,pos); + 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); - // The number of secondary electrons created - qTot = ((Int_t) (eDelta / kWion) + 1); - - // Create a new dEdx hit + // Create a new dEdx hit if (drRegion) { - AddHit(gAlice->CurrentTrack(),det,hits,qTot,kTRUE); - } + AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber() + ,det,hits,qTot, kTRUE); + } else { - AddHit(gAlice->CurrentTrack(),det,hits,qTot,kFALSE); - } + 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 @@ -670,7 +857,7 @@ void AliTRDv1::StepManager() // The energy loss according to Bethe Bloch iPdg = TMath::Abs(gMC->TrackPid()); if ( (iPdg != kPdgElectron) || - ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) { + ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) { gMC->TrackMomentum(mom); pTot = mom.Rho(); betaGamma = pTot / aMass; @@ -678,10 +865,8 @@ void AliTRDv1::StepManager() // Take charge > 1 into account charge = gMC->TrackCharge(); if (TMath::Abs(charge) > 1) pp = pp * charge*charge; - } - // Electrons above 20 Mev/c are at the plateau - else { - pp = kPrim * kPlateau; + } else { // Electrons above 20 Mev/c are at the plateau + pp = kPrim * kPlateau; } if (pp > 0) { @@ -690,15 +875,145 @@ void AliTRDv1::StepManager() 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. + // + + Int_t pla = 0; + Int_t cha = 0; + Int_t sec = 0; + Int_t det = 0; + Int_t qTot; + + Float_t hits[3]; + Double_t eDep; + + 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; + + const Float_t kWion = 23.53; // Ionization energy + 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 not charged track or already stopped or disappeared, just return. + if ((!gMC->TrackCharge()) || + gMC->IsTrackStop() || + gMC->IsTrackDisappeared()) return; + + // Inside a sensitive volume? + cIdCurrent = gMC->CurrentVolName(); + + if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE; + if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE; + + if ((!drRegion) && (!amRegion)) return; + + // 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 += 270.; + else 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; } + } + } + + if (!addthishit) return; + + det = fGeometry->GetDetector(pla,cha,sec); // The detector number + + Int_t trkStat = 0; // 0: InFlight 1:Entering 2:Exiting + + // 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->TrackMomentum(mom); + AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber()); + trkStat = 1; + } + if (gMC->IsTrackExiting()) { + gMC->TrackMomentum(mom); + AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber()); + trkStat = 2; } + // Create the hits from TR photons + if (fTR) CreateTRhit(det); + } + + // 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); + + // Set Maximum Step Size + // Produce only one hit if Ekin is below cutoff + if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return; + gMC->SetMaxStep(fStepSize); } @@ -746,6 +1061,111 @@ Double_t AliTRDv1::BetheBloch(Double_t bg) } +//_____________________________________________________________________________ +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