+
/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-/*
-$Log$
-Revision 1.9 2001/05/21 16:45:47 hristov
-Last minute changes (C.Blume)
-
-Revision 1.8 2001/01/26 19:56:57 hristov
-Major upgrade of AliRoot code
-
-Revision 1.7 2000/12/20 13:00:45 cblume
-Modifications for the HP-compiler
-
-Revision 1.6 2000/12/12 10:20:10 cblume
-Initialize fSepctrum = 0 in ctors
-
-Revision 1.5 2000/10/15 23:40:01 cblume
-Remove AliTRDconst
-
-Revision 1.4 2000/10/06 16:49:46 cblume
-Made Getters const
-
-Revision 1.3.2.1 2000/09/18 13:45:30 cblume
-New class AliTRDsim that simulates TR photons
-
-Revision 1.2 1999/09/29 09:24:35 fca
-Introduction of the Copyright and cvs Log
-
-*/
+/* $Id$ */
///////////////////////////////////////////////////////////////////////////////
// //
// 17.07.1998 - A.Andronic //
// 08.12.1998 - simplified version //
// 11.07.2000 - Adapted code to aliroot environment (C.Blume) //
+// 04.06.2004 - Momentum dependent parameters implemented (CBL) //
// //
///////////////////////////////////////////////////////////////////////////////
#include <TRandom.h>
#include <TMath.h>
#include <TParticle.h>
+#include <TVirtualMC.h>
+#include <TVirtualMCStack.h>
#include "AliModule.h"
+#include "AliLog.h"
+#include "AliMC.h"
#include "AliTRDsim.h"
ClassImp(AliTRDsim)
//_____________________________________________________________________________
-AliTRDsim::AliTRDsim():TObject()
+AliTRDsim::AliTRDsim()
+ :TObject()
+ ,fNFoilsDim(0)
+ ,fNFoils(0)
+ ,fNFoilsUp(0)
+ ,fFoilThick(0)
+ ,fGapThick(0)
+ ,fFoilDens(0)
+ ,fGapDens(0)
+ ,fFoilOmega(0)
+ ,fGapOmega()
+ ,fFoilZ(0)
+ ,fGapZ(0)
+ ,fFoilA(0)
+ ,fGapA(0)
+ ,fTemp(0)
+ ,fSpNBins(0)
+ ,fSpRange(0)
+ ,fSpBinWidth(0)
+ ,fSpLower(0)
+ ,fSpUpper(0)
+ ,fSigma(0)
+ ,fSpectrum(0)
{
//
// AliTRDsim default constructor
//
- fSpectrum = 0;
- fSigma = 0;
-
Init();
}
//_____________________________________________________________________________
AliTRDsim::AliTRDsim(AliModule *mod, Int_t foil, Int_t gap)
+ :TObject()
+ ,fNFoilsDim(0)
+ ,fNFoils(0)
+ ,fNFoilsUp(0)
+ ,fFoilThick(0)
+ ,fGapThick(0)
+ ,fFoilDens(0)
+ ,fGapDens(0)
+ ,fFoilOmega(0)
+ ,fGapOmega()
+ ,fFoilZ(0)
+ ,fGapZ(0)
+ ,fFoilA(0)
+ ,fGapA(0)
+ ,fTemp(0)
+ ,fSpNBins(0)
+ ,fSpRange(0)
+ ,fSpBinWidth(0)
+ ,fSpLower(0)
+ ,fSpUpper(0)
+ ,fSigma(0)
+ ,fSpectrum(0)
{
//
// AliTRDsim constructor. Takes the material properties of the radiator
// thickness of the gaps is 500 mu.
//
- Float_t aFoil, zFoil, rhoFoil;
- Float_t aGap, zGap, rhoGap;
- Float_t rad, abs;
- Char_t name[21];
+ Float_t aFoil;
+ Float_t zFoil;
+ Float_t rhoFoil;
- fSpectrum = 0;
- fSigma = 0;
+ Float_t aGap;
+ Float_t zGap;
+ Float_t rhoGap;
+
+ Float_t rad;
+ Float_t abs;
+
+ Char_t name[21];
Init();
//_____________________________________________________________________________
AliTRDsim::AliTRDsim(const AliTRDsim &s)
+ :TObject(s)
+ ,fNFoilsDim(s.fNFoilsDim)
+ ,fNFoils(0)
+ ,fNFoilsUp(0)
+ ,fFoilThick(s.fFoilThick)
+ ,fGapThick(s.fGapThick)
+ ,fFoilDens(s.fFoilDens)
+ ,fGapDens(s.fGapDens)
+ ,fFoilOmega(s.fFoilOmega)
+ ,fGapOmega(s.fGapOmega)
+ ,fFoilZ(s.fFoilZ)
+ ,fGapZ(s.fGapZ)
+ ,fFoilA(s.fFoilA)
+ ,fGapA(s.fGapA)
+ ,fTemp(s.fTemp)
+ ,fSpNBins(s.fSpNBins)
+ ,fSpRange(s.fSpRange)
+ ,fSpBinWidth(s.fSpBinWidth)
+ ,fSpLower(s.fSpLower)
+ ,fSpUpper(s.fSpUpper)
+ ,fSigma(0)
+ ,fSpectrum(0)
{
//
// AliTRDsim copy constructor
//
- ((AliTRDsim &) s).Copy(*this);
+ if (((AliTRDsim &) s).fNFoils) {
+ delete [] ((AliTRDsim &) s).fNFoils;
+ }
+ ((AliTRDsim &) s).fNFoils = new Int_t[fNFoilsDim];
+ for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
+ ((AliTRDsim &) s).fNFoils[iFoil] = fNFoils[iFoil];
+ }
+
+ if (((AliTRDsim &) s).fNFoilsUp) {
+ delete [] ((AliTRDsim &) s).fNFoilsUp;
+ }
+ ((AliTRDsim &) s).fNFoilsUp = new Double_t[fNFoilsDim];
+ for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
+ ((AliTRDsim &) s).fNFoilsUp[iFoil] = fNFoilsUp[iFoil];
+ }
+
+ if (((AliTRDsim &) s).fSigma) {
+ delete [] ((AliTRDsim &) s).fSigma;
+ }
+ ((AliTRDsim &) s).fSigma = new Double_t[fSpNBins];
+ for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
+ ((AliTRDsim &) s).fSigma[iBin] = fSigma[iBin];
+ }
+
+ fSpectrum->Copy(*((AliTRDsim &) s).fSpectrum);
}
// AliTRDsim destructor
//
- // if (fSpectrum) delete fSpectrum;
- if (fSigma) delete [] fSigma;
+ if (fSigma) {
+ delete [] fSigma;
+ fSigma = 0;
+ }
+
+ if (fNFoils) {
+ delete [] fNFoils;
+ fNFoils = 0;
+ }
+
+ if (fNFoilsUp) {
+ delete [] fNFoilsUp;
+ fNFoilsUp = 0;
+ }
}
//
if (this != &s) ((AliTRDsim &) s).Copy(*this);
+
return *this;
}
//_____________________________________________________________________________
-void AliTRDsim::Copy(TObject &s)
+void AliTRDsim::Copy(TObject &s) const
{
//
// Copy function
//
- ((AliTRDsim &) s).fNFoils = fNFoils;
((AliTRDsim &) s).fFoilThick = fFoilThick;
((AliTRDsim &) s).fFoilDens = fFoilDens;
((AliTRDsim &) s).fFoilOmega = fFoilOmega;
((AliTRDsim &) s).fSpLower = fSpLower;
((AliTRDsim &) s).fSpUpper = fSpUpper;
- if (((AliTRDsim &) s).fSigma) delete [] ((AliTRDsim &) s).fSigma;
- ((AliTRDsim &) s).fSigma = new Double_t[fSpNBins];
+ if (((AliTRDsim &) s).fNFoils) {
+ delete [] ((AliTRDsim &) s).fNFoils;
+ }
+ ((AliTRDsim &) s).fNFoils = new Int_t[fNFoilsDim];
+ for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
+ ((AliTRDsim &) s).fNFoils[iFoil] = fNFoils[iFoil];
+ }
+
+ if (((AliTRDsim &) s).fNFoilsUp) {
+ delete [] ((AliTRDsim &) s).fNFoilsUp;
+ }
+ ((AliTRDsim &) s).fNFoilsUp = new Double_t[fNFoilsDim];
+ for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
+ ((AliTRDsim &) s).fNFoilsUp[iFoil] = fNFoilsUp[iFoil];
+ }
+
+ if (((AliTRDsim &) s).fSigma) {
+ delete [] ((AliTRDsim &) s).fSigma;
+ }
+ ((AliTRDsim &) s).fSigma = new Double_t[fSpNBins];
for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
- ((AliTRDsim &) s).fSigma[iBin] = fSigma[iBin];
+ ((AliTRDsim &) s).fSigma[iBin] = fSigma[iBin];
}
fSpectrum->Copy(*((AliTRDsim &) s).fSpectrum);
{
//
// Initialization
- // The default radiator are 100 prolypropilene foils of 13 mu thickness
- // with gaps of 60 mu filled with CO2.
+ // The default radiator are prolypropilene foils of 10 mu thickness
+ // with gaps of 80 mu filled with N2.
//
- fNFoils = 100;
+ fNFoilsDim = 7;
+
+ if (fNFoils) {
+ delete [] fNFoils;
+ }
+ fNFoils = new Int_t[fNFoilsDim];
+ fNFoils[0] = 170;
+ fNFoils[1] = 225;
+ fNFoils[2] = 275;
+ fNFoils[3] = 305;
+ fNFoils[4] = 325;
+ fNFoils[5] = 340;
+ fNFoils[6] = 350;
+
+ if (fNFoilsUp) {
+ delete [] fNFoilsUp;
+ }
+ fNFoilsUp = new Double_t[fNFoilsDim];
+ fNFoilsUp[0] = 1.25;
+ fNFoilsUp[1] = 1.75;
+ fNFoilsUp[2] = 2.50;
+ fNFoilsUp[3] = 3.50;
+ fNFoilsUp[4] = 4.50;
+ fNFoilsUp[5] = 5.50;
+ fNFoilsUp[6] = 10000.0;
fFoilThick = 0.0013;
fFoilDens = 0.92;
fFoilOmega = Omega(fFoilDens,fFoilZ,fFoilA);
fGapThick = 0.0060;
- fGapDens = 0.001977;
- fGapZ = 7.45455;
- fGapA = 14.9091;
+ fGapDens = 0.00125;
+ fGapZ = 7.0;
+ fGapA = 14.00674;
fGapOmega = Omega(fGapDens ,fGapZ ,fGapA );
fTemp = 293.16;
if (fSpectrum) delete fSpectrum;
fSpectrum = new TH1D("TRspectrum","TR spectrum",fSpNBins,fSpLower,fSpUpper);
+ fSpectrum->SetDirectory(0);
// Set the sigma values
SetSigma();
break;
};
- // Calculate gamma
- Double_t gamma = TMath::Sqrt(p*p + mass*mass) / mass;
-
// Calculate the TR photons
- return TrPhotons(gamma, nPhoton, ePhoton);
+ return TrPhotons(p, mass, nPhoton, ePhoton);
}
//_____________________________________________________________________________
-Int_t AliTRDsim::TrPhotons(Double_t gamma, Int_t &nPhoton, Float_t *ePhoton)
+Int_t AliTRDsim::TrPhotons(Float_t p, Float_t mass
+ , Int_t &nPhoton, Float_t *ePhoton)
{
//
- // Produces TR photons.
+ // Produces TR photons using a parametric model for regular radiator. Photons
+ // with energy larger than 15 keV are included in the MC stack and tracked by VMC
+ // machinary.
+ //
+ // Input parameters:
+ // p - parent momentum [GeV/c]
+ // mass - parent mass
//
+ // Output :
+ // nPhoton - number of photons which have to be processed by custom code
+ // ePhoton - energy of this photons in keV.
+ //
const Double_t kAlpha = 0.0072973;
- const Int_t kSumMax = 10;
+ const Int_t kSumMax = 30;
+
+ Double_t tau = fGapThick / fFoilThick;
+
+ // Calculate gamma
+ Double_t gamma = TMath::Sqrt(p*p + mass*mass) / mass;
- Double_t kappa = fGapThick / fFoilThick;
+ // Select the number of foils corresponding to momentum
+ Int_t foils = SelectNFoils(p);
fSpectrum->Reset();
// The TR spectrum
- Double_t stemp = 0;
- for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
-
- // keV -> eV
- Double_t energyeV = (fSpBinWidth * iBin + 1.0) * 1e3;
-
- Double_t csFoil = fFoilOmega / energyeV;
- Double_t csGap = fGapOmega / energyeV;
-
- Double_t rho1 = energyeV * fFoilThick * 1e4 * 2.5
- * (1.0 / (gamma*gamma) + csFoil*csFoil);
- Double_t rho2 = energyeV * fFoilThick * 1e4 * 2.5
- * (1.0 / (gamma*gamma) + csGap *csGap);
+ Double_t csi1;
+ Double_t csi2;
+ Double_t rho1;
+ Double_t rho2;
+ Double_t fSigma;
+ Double_t sum;
+ Double_t nEqu;
+ Double_t thetaN;
+ Double_t aux;
+ Double_t energyeV;
+ Double_t energykeV;
+ for (Int_t iBin = 1; iBin <= fSpNBins; iBin++) {
+
+ energykeV = fSpectrum->GetBinCenter(iBin);
+ energyeV = energykeV * 1.0e3;
+
+ fSigma = Sigma(energykeV);
+
+ csi1 = fFoilOmega / energyeV;
+ csi2 = fGapOmega / energyeV;
+
+ rho1 = 2.5 * energyeV * fFoilThick * 1.0e4
+ * (1.0 / (gamma*gamma) + csi1*csi1);
+ rho2 = 2.5 * energyeV * fFoilThick * 1.0e4
+ * (1.0 / (gamma*gamma) + csi2 *csi2);
// Calculate the sum
- Double_t sum = 0;
- for (Int_t iSum = 0; iSum < kSumMax; iSum++) {
- Double_t tetan = (TMath::Pi() * 2.0 * (iSum+1) - (rho1 + kappa * rho2))
- / (kappa + 1.0);
- if (tetan < 0.0) tetan = 0.0;
- Double_t aux = 1.0 / (rho1 + tetan) - 1.0 / (rho2 + tetan);
- sum += tetan * (aux*aux) * (1.0 - TMath::Cos(rho1 + tetan));
+ sum = 0.0;
+ for (Int_t n = 1; n <= kSumMax; n++) {
+ thetaN = (TMath::Pi() * 2.0 * n - (rho1 + tau * rho2)) / (1.0 + tau);
+ if (thetaN < 0.0) {
+ thetaN = 0.0;
+ }
+ aux = 1.0 / (rho1 + thetaN) - 1.0 / (rho2 + thetaN);
+ sum += thetaN * (aux*aux) * (1.0 - TMath::Cos(rho1 + thetaN));
}
- // Absorbtion
- Double_t conv = 1.0 - TMath::Exp(-fNFoils * fSigma[iBin]);
-
- // eV -> keV
- Float_t energykeV = energyeV * 0.001;
+ // Equivalent number of foils
+ nEqu = (1.0 - TMath::Exp(-foils * fSigma)) / (1.0 - TMath::Exp(-fSigma));
// dN / domega
- Double_t wn = kAlpha * 4.0 / (fSigma[iBin] * (kappa + 1.0))
- * conv * sum / energykeV;
- fSpectrum->SetBinContent(iBin,wn);
-
- stemp += wn;
+ fSpectrum->SetBinContent(iBin,4.0 * kAlpha * nEqu * sum / (energykeV * (1.0 + tau)));
}
// <nTR> (binsize corr.)
- Float_t ntr = stemp * fSpBinWidth;
- // Number of TR photons from Poisson distribution with mean <ntr>
- nPhoton = gRandom->Poisson(ntr);
- // Energy of the TR photons
- for (Int_t iPhoton = 0; iPhoton < nPhoton; iPhoton++) {
- ePhoton[iPhoton] = fSpectrum->GetRandom();
+ Float_t nTr = fSpBinWidth * fSpectrum->Integral();
+ // Number of TR photons from Poisson distribution with mean <nTr>
+ Int_t nPhCand = gRandom->Poisson(nTr);
+
+ // Link the MC stack and get info about parent electron
+ TVirtualMCStack *stack = gMC->GetStack();
+ TParticle *trGenerator = stack->GetCurrentTrack();
+ Int_t track = stack->GetCurrentTrackNumber();
+ Double_t px = trGenerator->Px() / trGenerator->P();
+ Double_t py = trGenerator->Py() / trGenerator->P();
+ Double_t pz = trGenerator->Pz() / trGenerator->P();
+
+ // Current position of electron
+ Double_t x;
+ Double_t y;
+ Double_t z;
+ gMC->TrackPosition(x,y,z);
+
+ // Counter for TR analysed in custom code (e < 15keV)
+ nPhoton = 0;
+
+ for (Int_t iPhoton = 0; iPhoton < nPhCand; iPhoton++) {
+
+ // Energy of the TR photon
+ Double_t e = fSpectrum->GetRandom();
+
+ // Put TR photon on particle stack
+ if (e > 15.0 ) {
+
+ e *= 1.0e-6; // Convert it to GeV
+
+ Int_t phtrack;
+ stack-> PushTrack(1 // Must be 1
+ ,track // Identifier of the parent track, -1 for a primary
+ ,22 // Particle code.
+ ,px*e // 4 momentum (The photon is generated on the same
+ ,py*e // direction as the parent. For irregular radiator one
+ ,pz*e // can calculate also the angle but this is a secondary
+ ,e // order effect)
+ ,x,y,z,0.0 // 4 vertex
+ ,0.0,0.0,0.0 // Polarisation
+ ,kPFeedBackPhoton // Production mechanism (there is no TR in G3 so one
+ // has to make some convention)
+ ,phtrack // On output the number of the track stored
+ ,1.0
+ ,1);
+
+ }
+ // Custom treatment of TR photons
+ else {
+
+ ePhoton[nPhoton++] = e;
+
+ }
+
}
return 1;
// Sets the absorbtion crosssection for the energies of the TR spectrum
//
- if (fSigma) delete [] fSigma;
+ if (fSigma) {
+ delete [] fSigma;
+ }
fSigma = new Double_t[fSpNBins];
+
for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
Double_t energykeV = iBin * fSpBinWidth + 1.0;
fSigma[iBin] = Sigma(energykeV);
- //printf("SetSigma(): iBin = %d fSigma %g\n",iBin,fSigma[iBin]);
}
}
// Calculates the absorbtion crosssection for a one-foil-one-gap-radiator
//
- // Gas at 0 C
- const Double_t kTemp0 = 273.16;
-
// keV -> MeV
Double_t energyMeV = energykeV * 0.001;
if (energyMeV >= 0.001) {
- return(GetMuPo(energyMeV) * fFoilDens * fFoilThick +
- GetMuCO(energyMeV) * fGapDens * fGapThick * fTemp/kTemp0);
+ return(GetMuPo(energyMeV) * fFoilDens * fFoilThick +
+ GetMuAi(energyMeV) * fGapDens * fGapThick * GetTemp());
}
else {
- return 1e6;
+ return 1.0e6;
}
}
}
+//_____________________________________________________________________________
+Double_t AliTRDsim::GetMuAi(Double_t energyMeV)
+{
+ //
+ // Returns the photon absorbtion cross section for air
+ // Implemented by Oliver Busch
+ //
+
+ const Int_t kN = 38;
+
+ Double_t mu[kN] = { 0.35854E+04, 0.11841E+04, 0.52458E+03,
+ 0.16143E+03, 0.14250E+03, 0.15722E+03,
+ 0.77538E+02, 0.40099E+02, 0.23313E+02,
+ 0.98816E+01, 0.51000E+01, 0.16079E+01,
+ 0.77536E+00, 0.35282E+00, 0.24790E+00,
+ 0.20750E+00, 0.18703E+00, 0.16589E+00,
+ 0.15375E+00, 0.13530E+00, 0.12311E+00,
+ 0.10654E+00, 0.95297E-01, 0.86939E-01,
+ 0.80390E-01, 0.70596E-01, 0.63452E-01,
+ 0.56754E-01, 0.51644E-01, 0.44382E-01,
+ 0.35733E-01, 0.30721E-01, 0.27450E-01,
+ 0.25171E-01, 0.22205E-01, 0.20399E-01,
+ 0.18053E-01, 0.18057E-01 };
+
+
+
+ Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02,
+ 0.30000E-02, 0.32029E-02, 0.32029E-02,
+ 0.40000E-02, 0.50000E-02, 0.60000E-02,
+ 0.80000E-02, 0.10000E-01, 0.15000E-01,
+ 0.20000E-01, 0.30000E-01, 0.40000E-01,
+ 0.50000E-01, 0.60000E-01, 0.80000E-01,
+ 0.10000E+00, 0.15000E+00, 0.20000E+00,
+ 0.30000E+00, 0.40000E+00, 0.50000E+00,
+ 0.60000E+00, 0.80000E+00, 0.10000E+01,
+ 0.12500E+01, 0.15000E+01, 0.20000E+01,
+ 0.30000E+01, 0.40000E+01, 0.50000E+01,
+ 0.60000E+01, 0.80000E+01, 0.10000E+02,
+ 0.15000E+02, 0.20000E+02 };
+
+ return Interpolate(energyMeV,en,mu,kN);
+
+}
+
//_____________________________________________________________________________
Double_t AliTRDsim::Interpolate(Double_t energyMeV
, Double_t *en, Double_t *mu, Int_t n)
// Locates a point (xval) in a 1-dim grid (xv(n))
//
- if (xval >= xv[n-1]) return 1;
- if (xval < xv[0]) return -1;
+ if (xval >= xv[n-1]) {
+ return 1;
+ }
+ if (xval < xv[0]) {
+ return -1;
+ }
Int_t km;
Int_t kh = n - 1;
kl = 0;
while (kh - kl > 1) {
- if (xval < xv[km = (kl+kh)/2]) kh = km;
- else kl = km;
+ if (xval < xv[km = (kl+kh)/2]) {
+ kh = km;
+ }
+ else {
+ kl = km;
+ }
}
- if (xval < xv[kl] || xval > xv[kl+1] || kl >= n-1) {
- printf("Locate failed xv[%d] %f xval %f xv[%d] %f!!!\n"
- ,kl,xv[kl],xval,kl+1,xv[kl+1]);
+ if ((xval < xv[kl]) ||
+ (xval > xv[kl+1]) ||
+ (kl >= n-1)) {
+ AliFatal(Form("Locate failed xv[%d] %f xval %f xv[%d] %f!!!\n"
+ ,kl,xv[kl],xval,kl+1,xv[kl+1]));
exit(1);
}
return 0;
}
+
+//_____________________________________________________________________________
+Int_t AliTRDsim::SelectNFoils(Float_t p)
+{
+ //
+ // Selects the number of foils corresponding to the momentum
+ //
+
+ Int_t foils = fNFoils[fNFoilsDim-1];
+
+ for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
+ if (p < fNFoilsUp[iFoil]) {
+ foils = fNFoils[iFoil];
+ break;
+ }
+ }
+
+ return foils;
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff