]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | ////////////////////////////////////////////////////////////// |
2 | ////////////////////////////////////////////////////////////// | |
3 | // | |
4 | // Utility class to evaluate the material budget from | |
5 | // a given radius to the surface of an arbitrary cylinder | |
6 | // along radial directions from the centre: | |
7 | // | |
8 | // - radiation length | |
9 | // - Interaction length | |
10 | // - g/cm2 | |
11 | // | |
12 | // Geantinos are shot in the bins in the fNtheta bins in theta | |
13 | // and fNphi bins in phi with specified rectangular limits. | |
14 | // The statistics are accumulated per | |
15 | // fRadMin < r < fRadMax and <0 < z < fZMax | |
16 | // | |
17 | // To activate this option, you can do: | |
18 | // Root > gAlice.RunLego(); | |
19 | // or Root > .x menu.C then select menu item "RunLego" | |
20 | // Note that when calling gAlice->RunLego, an optional list | |
21 | // of arguments may be specified. | |
22 | // | |
23 | //Begin_Html | |
24 | /* | |
1439f98e | 25 | <img src="picts/alilego.gif"> |
fe4da5cc | 26 | */ |
27 | //End_Html | |
28 | // | |
29 | ////////////////////////////////////////////////////////////// | |
30 | ||
31 | #include "TMath.h" | |
1578254f | 32 | #include "AliLego.h" |
fe4da5cc | 33 | #include "AliRun.h" |
34 | #include "AliConst.h" | |
35 | ||
36 | ClassImp(AliLego) | |
37 | ||
38 | ||
39 | //___________________________________________ | |
40 | AliLego::AliLego() | |
41 | { | |
42 | fHistRadl = 0; | |
43 | fHistAbso = 0; | |
44 | fHistGcm2 = 0; | |
45 | fHistReta = 0; | |
46 | } | |
47 | ||
48 | //___________________________________________ | |
49 | AliLego::AliLego(const char *name, const char *title) | |
50 | : TNamed(name,title) | |
51 | { | |
52 | fHistRadl = 0; | |
53 | fHistAbso = 0; | |
54 | fHistGcm2 = 0; | |
55 | fHistReta = 0; | |
56 | } | |
57 | ||
58 | //___________________________________________ | |
59 | AliLego::~AliLego() | |
60 | { | |
61 | delete fHistRadl; | |
62 | delete fHistAbso; | |
63 | delete fHistGcm2; | |
64 | delete fHistReta; | |
65 | } | |
66 | ||
67 | //___________________________________________ | |
68 | void AliLego::GenerateKinematics() | |
69 | { | |
70 | // Create a geantino with kinematics corresponding to the current | |
71 | // bins in theta and phi. | |
72 | ||
1578254f | 73 | // |
74 | // Rootinos are 0 | |
75 | const Int_t mpart = 0; | |
fe4da5cc | 76 | Float_t orig[3], pmom[3]; |
77 | Float_t t, cost, sint, cosp, sinp; | |
78 | ||
79 | // --- Set to 0 radiation length, absorption length and g/cm2 --- | |
80 | fTotRadl = 0; | |
81 | fTotAbso = 0; | |
82 | fTotGcm2 = 0; | |
83 | ||
84 | fCurTheta = (fThetaMin+(fThetaBin-0.5)*(fThetaMax-fThetaMin)/fNtheta)*kDegrad; | |
85 | fCurPhi = (fPhiMin+(fPhiBin-0.5)*(fPhiMax-fPhiMin)/fNphi)*kDegrad; | |
86 | cost = TMath::Cos(fCurTheta); | |
87 | sint = TMath::Sin(fCurTheta); | |
88 | cosp = TMath::Cos(fCurPhi); | |
89 | sinp = TMath::Sin(fCurPhi); | |
90 | ||
91 | pmom[0] = cosp*sint; | |
92 | pmom[1] = sinp*sint; | |
93 | pmom[2] = cost; | |
94 | ||
95 | // --- Where to start | |
96 | orig[0] = orig[1] = orig[2] = 0; | |
97 | Float_t dalicz = 3000; | |
98 | if (fRadMin > 0) { | |
99 | t = PropagateCylinder(orig,pmom,fRadMin,dalicz); | |
100 | orig[0] = pmom[0]*t; | |
101 | orig[1] = pmom[1]*t; | |
102 | orig[2] = pmom[2]*t; | |
103 | if (TMath::Abs(orig[2]) > fZMax) return; | |
104 | } | |
105 | ||
106 | // --- We do start here | |
107 | Float_t polar[3]={0.,0.,0.}; | |
108 | Int_t ntr; | |
109 | gAlice->SetTrack(1, 0, mpart, pmom, orig, polar, 0, "LEGO ray", ntr); | |
110 | } | |
111 | ||
112 | //___________________________________________ | |
113 | void AliLego::Init(Int_t ntheta,Float_t themin,Float_t themax, | |
114 | Int_t nphi,Float_t phimin,Float_t phimax,Float_t rmin,Float_t rmax, | |
115 | Float_t zmax) | |
116 | { | |
117 | // specify the angular limits and the size of the rectangular box | |
118 | fNtheta = ntheta; | |
119 | fThetaMin = themin; | |
120 | fThetaMax = themax; | |
121 | fNphi = nphi; | |
122 | fPhiMin = phimin; | |
123 | fPhiMax = phimax; | |
124 | fRadMin = rmin; | |
125 | fRadMax = rmax; | |
126 | fZMax = zmax; | |
127 | Float_t etamin = -TMath::Log(TMath::Tan(TMath::Min((Double_t)fThetaMax*kDegrad/2,TMath::Pi()/2-1.e-10))); | |
128 | Float_t etamax = -TMath::Log(TMath::Tan(TMath::Max((Double_t)fThetaMin*kDegrad/2, 1.e-10))); | |
129 | ||
130 | fHistRadl = new TH2F("hradl","Radiation length map", nphi,phimin,phimax,ntheta,themin,themax); | |
131 | fHistAbso = new TH2F("habso","Interaction length map", nphi,phimin,phimax,ntheta,themin,themax); | |
132 | fHistGcm2 = new TH2F("hgcm2","g/cm2 length map", nphi,phimin,phimax,ntheta,themin,themax); | |
133 | fHistReta = new TH2F("hetar","Radiation length vs. eta",nphi,phimin,phimax,ntheta,etamin,etamax); | |
134 | ||
135 | } | |
136 | ||
137 | //___________________________________________ | |
138 | Float_t AliLego::PropagateCylinder(Float_t *x, Float_t *v, Float_t r, Float_t z) | |
139 | { | |
140 | // Propagate to cylinder from inside | |
141 | ||
142 | Double_t hnorm, sz, t, t1, t2, t3, sr; | |
143 | Double_t d[3]; | |
144 | const Float_t kSmall = 1e-8; | |
145 | const Float_t kSmall2 = kSmall*kSmall; | |
146 | ||
147 | // ---> Find intesection with Z planes | |
148 | d[0] = v[0]; | |
149 | d[1] = v[1]; | |
150 | d[2] = v[2]; | |
151 | hnorm = TMath::Sqrt(1/(d[0]*d[0]+d[1]*d[1]+d[2]*d[2])); | |
152 | d[0] *= hnorm; | |
153 | d[1] *= hnorm; | |
154 | d[2] *= hnorm; | |
155 | if (d[2] > kSmall) sz = (z-x[2])/d[2]; | |
156 | else if (d[2] < -kSmall) sz = -(z+x[2])/d[2]; | |
157 | else sz = 1.e10; // ---> Direction parallel to X-Y, no intersection | |
158 | ||
159 | // ---> Intersection with cylinders | |
160 | // Intersection point (x,y,z) | |
161 | // (x,y,z) is on track : x=X(1)+t*D(1) | |
162 | // y=X(2)+t*D(2) | |
163 | // z=X(3)+t*D(3) | |
164 | // (x,y,z) is on cylinder : x**2 + y**2 = R**2 | |
165 | // | |
166 | // (D(1)**2+D(2)**2)*t**2 | |
167 | // +2.*(X(1)*D(1)+X(2)*D(2))*t | |
168 | // +X(1)**2+X(2)**2-R**2=0 | |
169 | // ---> Solve second degree equation | |
170 | t1 = d[0]*d[0] + d[1]*d[1]; | |
171 | if (t1 <= kSmall2) { | |
172 | t = sz; // ---> Track parallel to the z-axis, take distance to planes | |
173 | } else { | |
174 | t2 = x[0]*d[0] + x[1]*d[1]; | |
175 | t3 = x[0]*x[0] + x[1]*x[1]; | |
176 | // ---> It should be positive, but there may be numerical problems | |
177 | sr = (t2 +TMath::Sqrt(TMath::Max(t2*t2-(t3-r*r)*t1,0.)))/t1; | |
178 | // ---> Find minimum distance between planes and cylinder | |
179 | t = TMath::Min(sz,sr); | |
180 | } | |
181 | return t; | |
182 | } | |
183 | ||
184 | //___________________________________________ | |
185 | void AliLego::Run() | |
186 | { | |
187 | // loop on phi,theta bins | |
188 | AliMC* pMC=AliMC::GetMC(); | |
189 | pMC->InitLego(); | |
190 | Float_t thed, phid, eta; | |
191 | for (fPhiBin=1; fPhiBin<=fNphi; fPhiBin++) { | |
192 | printf("AliLego Generating rays in phi bin:%d\n",fPhiBin); | |
193 | for (fThetaBin=1; fThetaBin<=fNtheta; fThetaBin++) { | |
194 | pMC->Gtrigi(); | |
195 | pMC->Gtrigc(); | |
196 | GenerateKinematics(); | |
1578254f | 197 | pMC->Gtreve_root(); |
fe4da5cc | 198 | |
199 | thed = fCurTheta*kRaddeg; | |
200 | phid = fCurPhi*kRaddeg; | |
201 | eta = -TMath::Log(TMath::Tan(TMath::Max( | |
202 | TMath::Min((Double_t)fCurTheta/2,TMath::Pi()/2-1.e-10),1.e-10))); | |
203 | fHistRadl->Fill(phid,thed,fTotRadl); | |
204 | fHistAbso->Fill(phid,thed,fTotAbso); | |
205 | fHistGcm2->Fill(phid,thed,fTotGcm2); | |
206 | fHistReta->Fill(phid,eta,fTotRadl); | |
207 | gAlice->FinishEvent(); | |
208 | } | |
209 | } | |
210 | // store histograms in current Root file | |
211 | fHistRadl->Write(); | |
212 | fHistAbso->Write(); | |
213 | fHistGcm2->Write(); | |
214 | fHistReta->Write(); | |
215 | } | |
216 | ||
217 | //___________________________________________ | |
218 | void AliLego::StepManager() | |
219 | { | |
220 | // called from AliRun::Stepmanager from gustep. | |
221 | // Accumulate the 3 parameters step by step | |
222 | AliMC* pMC = AliMC::GetMC(); | |
223 | ||
224 | Float_t t, tt; | |
225 | Float_t a,z,dens,radl,absl; | |
226 | ||
227 | Float_t step = pMC->TrackStep(); | |
228 | ||
229 | Float_t vect[3], pmom[4]; | |
230 | pMC->TrackPosition(vect); | |
231 | pMC->TrackMomentum(pmom); | |
232 | pMC->CurrentMaterial(a,z,dens,radl,absl); | |
233 | ||
234 | if (z < 1) return; | |
235 | ||
236 | // --- See if we have to stop now | |
237 | if (TMath::Abs(vect[2]) > fZMax || | |
238 | vect[0]*vect[0] +vect[1]*vect[1] > fRadMax*fRadMax) { | |
239 | pMC->StopEvent(); | |
240 | } else { | |
241 | ||
242 | // --- See how long we have to go | |
243 | t = PropagateCylinder(vect,pmom,fRadMax,fZMax); | |
244 | tt = TMath::Min(step,t); | |
245 | ||
246 | fTotAbso += tt/absl; | |
247 | fTotRadl += tt/radl; | |
248 | fTotGcm2 += tt*dens; | |
249 | } | |
250 | } | |
251 |