New gain variation diagnostics.
[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.4  2000/12/01 17:37:44  morsch
19   Replace in argument of SetTrack(..) string constant by  kPFeedBackPhoton.
20
21   Revision 1.3  2000/10/03 21:44:09  morsch
22   Use AliSegmentation and AliHit abstract base classes.
23
24   Revision 1.2  2000/10/02 21:28:12  fca
25   Removal of useless dependecies via forward declarations
26
27   Revision 1.1  2000/06/12 15:29:37  jbarbosa
28   Cleaned up version.
29
30 */
31
32 #include "AliRICHResponseV0.h"
33 #include "AliSegmentation.h"
34 #include "AliRun.h"
35 #include "AliMC.h"
36
37 #include <TMath.h>
38 #include <TRandom.h>
39 #include <TParticle.h>
40 //___________________________________________
41 ClassImp(AliRICHResponseV0)
42
43 Float_t AliRICHResponseV0::IntPH(Float_t eloss, Float_t yhit)
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;
51     Double_t gain_var=1;
52
53     if (nel == 0) nel=1;
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
88     return charge;
89 }
90
91 Float_t AliRICHResponseV0::IntPH(Float_t yhit)
92 {
93
94 //  Get number of electrons and return charge, for a single photon
95
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      }
125     return charge;
126 }
127
128
129
130 // -------------------------------------------
131 Float_t AliRICHResponseV0::IntXY(AliSegmentation * segmentation)
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
172 Int_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],
282                      kPFeedBackPhoton, nt, 1.);
283   }
284   return(sNfeed);
285 }
286
287
288
289