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