Old style include needed on HP
[u/mrichter/AliRoot.git] / RICH / AliRICHResponseV0.cxx
CommitLineData
237c933d 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17 $Log$
d9042b8a 18 Revision 1.7 2001/05/10 12:32:27 jbarbosa
19 Changed call to SetTrack.
20
3c5ffeb4 21 Revision 1.6 2001/02/23 17:39:02 jbarbosa
22 Removed verbose output.
23
d65c4a20 24 Revision 1.5 2001/02/23 17:25:08 jbarbosa
25 Re-definition of IntPH() to accomodate for wire sag effect.
26
733b4fa4 27 Revision 1.4 2000/12/01 17:37:44 morsch
28 Replace in argument of SetTrack(..) string constant by kPFeedBackPhoton.
29
9e4670f0 30 Revision 1.3 2000/10/03 21:44:09 morsch
31 Use AliSegmentation and AliHit abstract base classes.
32
a2f7eaf6 33 Revision 1.2 2000/10/02 21:28:12 fca
34 Removal of useless dependecies via forward declarations
35
94de3818 36 Revision 1.1 2000/06/12 15:29:37 jbarbosa
37 Cleaned up version.
38
237c933d 39*/
40
41#include "AliRICHResponseV0.h"
a2f7eaf6 42#include "AliSegmentation.h"
237c933d 43#include "AliRun.h"
94de3818 44#include "AliMC.h"
237c933d 45
46#include <TMath.h>
47#include <TRandom.h>
48#include <TParticle.h>
49//___________________________________________
50ClassImp(AliRICHResponseV0)
51
733b4fa4 52Float_t AliRICHResponseV0::IntPH(Float_t eloss, Float_t yhit)
237c933d 53{
54 // Get number of electrons and return charge
55
56 Int_t nel;
57 nel= Int_t(eloss/fEIonisation);
58
59 Float_t charge=0;
733b4fa4 60 Double_t gain_var=1;
61
237c933d 62 if (nel == 0) nel=1;
733b4fa4 63
64 if (fWireSag)
65 {
66 //printf("Voltage:%d, Yhit:%f\n",fVoltage, yhit);
67
68 if (fVoltage==2150)
69 {
70 gain_var = 9e-6*TMath::Power(yhit,4) + 2e-7*TMath::Power(yhit,3) - 0.0316*TMath::Power(yhit,2) - 3e-4*yhit + 25.367;
71 //gain_var = 9e-5*TMath::Power(yhit,4) + 2e-6*TMath::Power(yhit,3) - 0.316*TMath::Power(yhit,2) - 3e-3*yhit + 253.67;
72 }
73 if (fVoltage==2100)
74 gain_var = 8e-6*TMath::Power(yhit,4) + 2e-7*TMath::Power(yhit,3) - 0.0283*TMath::Power(yhit,2) - 2e-4*yhit + 23.015;
75 if (fVoltage==2050)
76 gain_var = 7e-6*TMath::Power(yhit,4) + 1e-7*TMath::Power(yhit,3) - 0.0254*TMath::Power(yhit,2) - 2e-4*yhit + 20.888;
77 if (fVoltage==2000)
78 gain_var = 6e-6*TMath::Power(yhit,4) + 8e-8*TMath::Power(yhit,3) - 0.0227*TMath::Power(yhit,2) - 1e-4*yhit + 18.961;
79
80 gain_var = gain_var/100;
81 //printf("Yhit:%f, Gain variation:%f\n",yhit,gain_var);
82
83 Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9;
d65c4a20 84 //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
733b4fa4 85
86 for (Int_t i=1;i<=nel;i++) {
87 charge -= gain*TMath::Log(gRandom->Rndm());
88 }
89 }
90 else
91 {
92 for (Int_t i=1;i<=nel;i++) {
93 charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
94 }
95 }
96
237c933d 97 return charge;
98}
99
733b4fa4 100Float_t AliRICHResponseV0::IntPH(Float_t yhit)
237c933d 101{
102
103// Get number of electrons and return charge, for a single photon
104
733b4fa4 105 Float_t charge=0;
106 Double_t gain_var=1;
107
108 if (fWireSag)
109 {
110 if (fVoltage==2150)
111 {
112 gain_var = 9e-6*TMath::Power(yhit,4) + 2e-7*TMath::Power(yhit,3) - 0.0316*TMath::Power(yhit,2) - 3e-4*yhit + 25.367;
113 //gain_var = 9e-5*TMath::Power(yhit,4) + 2e-6*TMath::Power(yhit,3) - 0.316*TMath::Power(yhit,2) - 3e-3*yhit + 253.67;
114 }
115 if (fVoltage==2100)
116 gain_var = 8e-6*TMath::Power(yhit,4) + 2e-7*TMath::Power(yhit,3) - 0.0283*TMath::Power(yhit,2) - 2e-4*yhit + 23.015;
117 if (fVoltage==2050)
118 gain_var = 7e-6*TMath::Power(yhit,4) + 1e-7*TMath::Power(yhit,3) - 0.0254*TMath::Power(yhit,2) - 2e-4*yhit + 20.888;
119 if (fVoltage==2000)
120 gain_var = 6e-6*TMath::Power(yhit,4) + 8e-8*TMath::Power(yhit,3) - 0.0227*TMath::Power(yhit,2) - 1e-4*yhit + 18.961;
121
122 gain_var = gain_var/100;
123 //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain_var);
124
125 Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9;
126
127 charge -= gain*TMath::Log(gRandom->Rndm());
d65c4a20 128 //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
733b4fa4 129 }
130 else
131 {
132 charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
133 }
237c933d 134 return charge;
135}
136
137
138
139// -------------------------------------------
a2f7eaf6 140Float_t AliRICHResponseV0::IntXY(AliSegmentation * segmentation)
237c933d 141{
142
143 const Float_t kInversePitch = 1/fPitch;
144 Float_t response;
145//
146// Integration limits defined by segmentation model
147//
148
149 Float_t xi1, xi2, yi1, yi2;
150 segmentation->IntegrationLimits(xi1,xi2,yi1,yi2);
151
152 xi1=xi1*kInversePitch;
153 xi2=xi2*kInversePitch;
154 yi1=yi1*kInversePitch;
155 yi2=yi2*kInversePitch;
156
157 //printf("Integration Limits: %f-%f, %f-%f\n",xi1,xi2,yi1,yi2);
158
159 //printf("KInversePitch:%f\n",kInversePitch);
160
161 //
162// The Mathieson function
163 Double_t ux1=fSqrtKx3*TMath::TanH(fKx2*xi1);
164 Double_t ux2=fSqrtKx3*TMath::TanH(fKx2*xi2);
165
166 Double_t uy1=fSqrtKy3*TMath::TanH(fKy2*yi1);
167 Double_t uy2=fSqrtKy3*TMath::TanH(fKy2*yi2);
168
169 //printf("Integration Data: %f-%f, %f-%f\n",ux1,ux2,uy1,uy2);
170
171 //printf("%f %f %f %f\n",fSqrtKx3,fKx2,fKy4,fKx4);
172
173 response=4.*fKx4*(TMath::ATan(ux2)-TMath::ATan(ux1))*fKy4*(TMath::ATan(uy2)-TMath::ATan(uy1));
174
175 //printf("Response:%f\n",response);
176
177 return response;
178
179}
180
181Int_t AliRICHResponseV0::FeedBackPhotons(Float_t *source, Float_t qtot)
182{
183 //
184 // Generate FeedBack photons
185 //
186 Int_t j, ipart, nt;
187
188 Int_t sNfeed=0;
189
190
191 // Local variables
192 Float_t cthf, ranf[2], phif, enfp = 0, sthf;
193 Int_t i, ifeed;
194 Float_t e1[3], e2[3], e3[3];
195 Float_t vmod, uswop;
196 Float_t fp, random;
197 Float_t dir[3], phi;
198 Int_t nfp;
d9042b8a 199 Float_t pol[3], mom[4];
237c933d 200 TLorentzVector position;
201 //
202 // Determine number of feedback photons
203
204 // Get weight of current particle
205 TParticle *current = (TParticle*)
206 (*gAlice->Particles())[gAlice->CurrentTrack()];
207
208 ifeed = Int_t(current->GetWeight()/100+0.5);
209 ipart = gMC->TrackPid();
210 fp = fAlphaFeedback * qtot;
211 nfp = gRandom->Poisson(fp);
212
213 // This call to fill the time of flight
214 gMC->TrackPosition(position);
3c5ffeb4 215 //printf("Track position: %f %f %f %15.12f\n", position[0],position[1],position[2],position[3]);
237c933d 216 //
217 // Generate photons
218 for (i = 0; i <nfp; i++) {
219
220 // Direction
221 gMC->Rndm(ranf, 2);
222 cthf = ranf[0] * 2 - 1.;
223 if (cthf < 0) continue;
224 sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
225 phif = ranf[1] * 2 * TMath::Pi();
226 //
227 gMC->Rndm(&random, 1);
228 if (random <= .57) {
229 enfp = 7.5e-9;
230 } else if (random <= .7) {
231 enfp = 6.4e-9;
232 } else {
233 enfp = 7.9e-9;
234 }
235
236 dir[0] = sthf * TMath::Sin(phif);
237 dir[1] = cthf;
238 dir[2] = sthf * TMath::Cos(phif);
239 gMC->Gdtom(dir, mom, 2);
240 mom[0]*=enfp;
241 mom[1]*=enfp;
242 mom[2]*=enfp;
d9042b8a 243 mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
3c5ffeb4 244 //printf("Dir %f %f %f\n",dir[0],dir[1],dir[2]);
245 //printf("Momentum %15.12f %15.12f %15.12f\n",mom[0],mom[1],mom[2]);
d9042b8a 246 //printf("Energy %e\n", mom[3]);
237c933d 247
248 // Polarisation
249 e1[0] = 0;
250 e1[1] = -dir[2];
251 e1[2] = dir[1];
252
253 e2[0] = -dir[1];
254 e2[1] = dir[0];
255 e2[2] = 0;
256
257 e3[0] = dir[1];
258 e3[1] = 0;
259 e3[2] = -dir[0];
260
261 vmod=0;
262 for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
263 if (!vmod) for(j=0;j<3;j++) {
264 uswop=e1[j];
265 e1[j]=e3[j];
266 e3[j]=uswop;
267 }
268 vmod=0;
269 for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
270 if (!vmod) for(j=0;j<3;j++) {
271 uswop=e2[j];
272 e2[j]=e3[j];
273 e3[j]=uswop;
274 }
275
276 vmod=0;
277 for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
278 vmod=TMath::Sqrt(1/vmod);
279 for(j=0;j<3;j++) e1[j]*=vmod;
280
281 vmod=0;
282 for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
283 vmod=TMath::Sqrt(1/vmod);
284 for(j=0;j<3;j++) e2[j]*=vmod;
285
286 gMC->Rndm(ranf, 1);
287 phi = ranf[0] * 2 * TMath::Pi();
288 for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
289 gMC->Gdtom(pol, pol, 2);
290
291 // Put photon on the stack and label it as feedback (51, 52)
292 ++sNfeed;
293
294 gAlice->SetTrack(Int_t(1), gAlice->CurrentTrack(), Int_t(50000051),
3c5ffeb4 295 mom[0],mom[1],mom[2],mom[3],source[0],source[1],source[2],position[3],pol[0],pol[1],pol[2],
9e4670f0 296 kPFeedBackPhoton, nt, 1.);
3c5ffeb4 297
298 //printf("Adding feedback with tof %f and going to %f %f %f\n",position[3],mom[0],mom[1],mom[2]);
237c933d 299 }
3c5ffeb4 300 //if(sNfeed)
301 //printf("feedbacks produced:%d\n",sNfeed);
237c933d 302 return(sNfeed);
303}
304
305
306
307