New gain variation diagnostics.
[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$
733b4fa4 18 Revision 1.4 2000/12/01 17:37:44 morsch
19 Replace in argument of SetTrack(..) string constant by kPFeedBackPhoton.
20
9e4670f0 21 Revision 1.3 2000/10/03 21:44:09 morsch
22 Use AliSegmentation and AliHit abstract base classes.
23
a2f7eaf6 24 Revision 1.2 2000/10/02 21:28:12 fca
25 Removal of useless dependecies via forward declarations
26
94de3818 27 Revision 1.1 2000/06/12 15:29:37 jbarbosa
28 Cleaned up version.
29
237c933d 30*/
31
32#include "AliRICHResponseV0.h"
a2f7eaf6 33#include "AliSegmentation.h"
237c933d 34#include "AliRun.h"
94de3818 35#include "AliMC.h"
237c933d 36
37#include <TMath.h>
38#include <TRandom.h>
39#include <TParticle.h>
40//___________________________________________
41ClassImp(AliRICHResponseV0)
42
733b4fa4 43Float_t AliRICHResponseV0::IntPH(Float_t eloss, Float_t yhit)
237c933d 44{
45 // Get number of electrons and return charge
46
47 Int_t nel;
48 nel= Int_t(eloss/fEIonisation);
49
50 Float_t charge=0;
733b4fa4 51 Double_t gain_var=1;
52
237c933d 53 if (nel == 0) nel=1;
733b4fa4 54
55 if (fWireSag)
56 {
57 //printf("Voltage:%d, Yhit:%f\n",fVoltage, yhit);
58
59 if (fVoltage==2150)
60 {
61 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;
62 //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;
63 }
64 if (fVoltage==2100)
65 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;
66 if (fVoltage==2050)
67 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;
68 if (fVoltage==2000)
69 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;
70
71 gain_var = gain_var/100;
72 //printf("Yhit:%f, Gain variation:%f\n",yhit,gain_var);
73
74 Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9;
75 printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
76
77 for (Int_t i=1;i<=nel;i++) {
78 charge -= gain*TMath::Log(gRandom->Rndm());
79 }
80 }
81 else
82 {
83 for (Int_t i=1;i<=nel;i++) {
84 charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
85 }
86 }
87
237c933d 88 return charge;
89}
90
733b4fa4 91Float_t AliRICHResponseV0::IntPH(Float_t yhit)
237c933d 92{
93
94// Get number of electrons and return charge, for a single photon
95
733b4fa4 96 Float_t charge=0;
97 Double_t gain_var=1;
98
99 if (fWireSag)
100 {
101 if (fVoltage==2150)
102 {
103 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;
104 //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;
105 }
106 if (fVoltage==2100)
107 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;
108 if (fVoltage==2050)
109 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;
110 if (fVoltage==2000)
111 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;
112
113 gain_var = gain_var/100;
114 //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain_var);
115
116 Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9;
117
118 charge -= gain*TMath::Log(gRandom->Rndm());
119 printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
120 }
121 else
122 {
123 charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
124 }
237c933d 125 return charge;
126}
127
128
129
130// -------------------------------------------
a2f7eaf6 131Float_t AliRICHResponseV0::IntXY(AliSegmentation * segmentation)
237c933d 132{
133
134 const Float_t kInversePitch = 1/fPitch;
135 Float_t response;
136//
137// Integration limits defined by segmentation model
138//
139
140 Float_t xi1, xi2, yi1, yi2;
141 segmentation->IntegrationLimits(xi1,xi2,yi1,yi2);
142
143 xi1=xi1*kInversePitch;
144 xi2=xi2*kInversePitch;
145 yi1=yi1*kInversePitch;
146 yi2=yi2*kInversePitch;
147
148 //printf("Integration Limits: %f-%f, %f-%f\n",xi1,xi2,yi1,yi2);
149
150 //printf("KInversePitch:%f\n",kInversePitch);
151
152 //
153// The Mathieson function
154 Double_t ux1=fSqrtKx3*TMath::TanH(fKx2*xi1);
155 Double_t ux2=fSqrtKx3*TMath::TanH(fKx2*xi2);
156
157 Double_t uy1=fSqrtKy3*TMath::TanH(fKy2*yi1);
158 Double_t uy2=fSqrtKy3*TMath::TanH(fKy2*yi2);
159
160 //printf("Integration Data: %f-%f, %f-%f\n",ux1,ux2,uy1,uy2);
161
162 //printf("%f %f %f %f\n",fSqrtKx3,fKx2,fKy4,fKx4);
163
164 response=4.*fKx4*(TMath::ATan(ux2)-TMath::ATan(ux1))*fKy4*(TMath::ATan(uy2)-TMath::ATan(uy1));
165
166 //printf("Response:%f\n",response);
167
168 return response;
169
170}
171
172Int_t AliRICHResponseV0::FeedBackPhotons(Float_t *source, Float_t qtot)
173{
174 //
175 // Generate FeedBack photons
176 //
177 Int_t j, ipart, nt;
178
179 Int_t sNfeed=0;
180
181
182 // Local variables
183 Float_t cthf, ranf[2], phif, enfp = 0, sthf;
184 Int_t i, ifeed;
185 Float_t e1[3], e2[3], e3[3];
186 Float_t vmod, uswop;
187 Float_t fp, random;
188 Float_t dir[3], phi;
189 Int_t nfp;
190 Float_t pol[3], mom[3];
191 TLorentzVector position;
192 //
193 // Determine number of feedback photons
194
195 // Get weight of current particle
196 TParticle *current = (TParticle*)
197 (*gAlice->Particles())[gAlice->CurrentTrack()];
198
199 ifeed = Int_t(current->GetWeight()/100+0.5);
200 ipart = gMC->TrackPid();
201 fp = fAlphaFeedback * qtot;
202 nfp = gRandom->Poisson(fp);
203
204 // This call to fill the time of flight
205 gMC->TrackPosition(position);
206 //
207 // Generate photons
208 for (i = 0; i <nfp; i++) {
209
210 // Direction
211 gMC->Rndm(ranf, 2);
212 cthf = ranf[0] * 2 - 1.;
213 if (cthf < 0) continue;
214 sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
215 phif = ranf[1] * 2 * TMath::Pi();
216 //
217 gMC->Rndm(&random, 1);
218 if (random <= .57) {
219 enfp = 7.5e-9;
220 } else if (random <= .7) {
221 enfp = 6.4e-9;
222 } else {
223 enfp = 7.9e-9;
224 }
225
226 dir[0] = sthf * TMath::Sin(phif);
227 dir[1] = cthf;
228 dir[2] = sthf * TMath::Cos(phif);
229 gMC->Gdtom(dir, mom, 2);
230 mom[0]*=enfp;
231 mom[1]*=enfp;
232 mom[2]*=enfp;
233
234 // Polarisation
235 e1[0] = 0;
236 e1[1] = -dir[2];
237 e1[2] = dir[1];
238
239 e2[0] = -dir[1];
240 e2[1] = dir[0];
241 e2[2] = 0;
242
243 e3[0] = dir[1];
244 e3[1] = 0;
245 e3[2] = -dir[0];
246
247 vmod=0;
248 for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
249 if (!vmod) for(j=0;j<3;j++) {
250 uswop=e1[j];
251 e1[j]=e3[j];
252 e3[j]=uswop;
253 }
254 vmod=0;
255 for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
256 if (!vmod) for(j=0;j<3;j++) {
257 uswop=e2[j];
258 e2[j]=e3[j];
259 e3[j]=uswop;
260 }
261
262 vmod=0;
263 for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
264 vmod=TMath::Sqrt(1/vmod);
265 for(j=0;j<3;j++) e1[j]*=vmod;
266
267 vmod=0;
268 for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
269 vmod=TMath::Sqrt(1/vmod);
270 for(j=0;j<3;j++) e2[j]*=vmod;
271
272 gMC->Rndm(ranf, 1);
273 phi = ranf[0] * 2 * TMath::Pi();
274 for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
275 gMC->Gdtom(pol, pol, 2);
276
277 // Put photon on the stack and label it as feedback (51, 52)
278 ++sNfeed;
279
280 gAlice->SetTrack(Int_t(1), gAlice->CurrentTrack(), Int_t(50000051),
281 mom,source,pol,position[3],
9e4670f0 282 kPFeedBackPhoton, nt, 1.);
237c933d 283 }
284 return(sNfeed);
285}
286
287
288
289