1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 #include <TParticle.h>
21 #include <TVirtualMC.h>
23 #include "AliRICHResponseV0.h"
25 #include "AliSegmentation.h"
27 //___________________________________________
28 ClassImp(AliRICHResponseV0)
30 AliRICHResponseV0::AliRICHResponseV0()
32 SetSigmaIntegration(5.);
34 SetChargeSpread(0.18, 0.18);
36 SetAlphaFeedback(0.036);
37 SetEIonisation(26.e-9);
38 SetSqrtKx3(0.77459667);
41 SetSqrtKy3(0.77459667);
45 SetWireSag(1); // 1->On, 0->Off
46 SetVoltage(2150); // Should only be 2000, 2050, 2100 or 2150
47 }//AliRICHResponseV0::ctor()
48 Float_t AliRICHResponseV0::IntPH(Float_t eloss, Float_t yhit)
50 // Get number of electrons and return charge
53 nel= Int_t(eloss/fEIonisation);
62 //printf("Voltage:%d, Yhit:%f\n",fVoltage, yhit);
66 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;
67 //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;
70 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;
72 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;
74 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;
76 gain_var = gain_var/100;
77 //printf("Yhit:%f, Gain variation:%f\n",yhit,gain_var);
79 Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9;
80 //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
82 for (Int_t i=1;i<=nel;i++) {
83 charge -= gain*TMath::Log(gRandom->Rndm());
88 for (Int_t i=1;i<=nel;i++) {
89 charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
96 Float_t AliRICHResponseV0::IntPH(Float_t yhit)
99 // Get number of electrons and return charge, for a single photon
108 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;
109 //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;
112 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;
114 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;
116 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;
118 gain_var = gain_var/100;
119 //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain_var);
121 Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9;
123 charge -= gain*TMath::Log(gRandom->Rndm());
124 //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
128 charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
135 // -------------------------------------------
136 Float_t AliRICHResponseV0::IntXY(AliSegmentation * segmentation)
139 const Float_t kInversePitch = 1/fPitch;
142 // Integration limits defined by segmentation model
145 Float_t xi1, xi2, yi1, yi2;
146 segmentation->IntegrationLimits(xi1,xi2,yi1,yi2);
148 xi1=xi1*kInversePitch;
149 xi2=xi2*kInversePitch;
150 yi1=yi1*kInversePitch;
151 yi2=yi2*kInversePitch;
153 //printf("Integration Limits: %f-%f, %f-%f\n",xi1,xi2,yi1,yi2);
155 //printf("KInversePitch:%f\n",kInversePitch);
158 // The Mathieson function
159 Double_t ux1=fSqrtKx3*TMath::TanH(fKx2*xi1);
160 Double_t ux2=fSqrtKx3*TMath::TanH(fKx2*xi2);
162 Double_t uy1=fSqrtKy3*TMath::TanH(fKy2*yi1);
163 Double_t uy2=fSqrtKy3*TMath::TanH(fKy2*yi2);
165 //printf("Integration Data: %f-%f, %f-%f\n",ux1,ux2,uy1,uy2);
167 //printf("%f %f %f %f\n",fSqrtKx3,fKx2,fKy4,fKx4);
169 response=4.*fKx4*(TMath::ATan(ux2)-TMath::ATan(ux1))*fKy4*(TMath::ATan(uy2)-TMath::ATan(uy1));
171 //printf("Response:%f\n",response);
177 Int_t AliRICHResponseV0::FeedBackPhotons(Float_t *source, Float_t qtot)
180 // Generate FeedBack photons
189 Float_t cthf, phif, enfp = 0, sthf;
191 Float_t e1[3], e2[3], e3[3];
197 Float_t pol[3], mom[4];
198 TLorentzVector position;
200 // Determine number of feedback photons
202 // Get weight of current particle
203 TParticle *current = (TParticle*)
204 (*gAlice->Particles())[gAlice->GetCurrentTrackNumber()];
206 ifeed = Int_t(current->GetWeight()/100+0.5);
207 ipart = gMC->TrackPid();
208 fp = fAlphaFeedback * qtot;
209 nfp = gRandom->Poisson(fp);
211 // This call to fill the time of flight
212 gMC->TrackPosition(position);
213 //printf("Track position: %f %f %f %15.12f\n", position[0],position[1],position[2],position[3]);
216 for (i = 0; i <nfp; i++) {
219 gMC->GetRandom()->RndmArray(2,ranf);
220 cthf = ranf[0] * 2 - 1.;
221 if (cthf < 0) continue;
222 sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
223 phif = ranf[1] * 2 * TMath::Pi();
225 //gMC->Rndm(&random, 1);
226 gMC->GetRandom()->RndmArray(1, &random);
229 } else if (random <= .7) {
235 dir[0] = sthf * TMath::Sin(phif);
237 dir[2] = sthf * TMath::Cos(phif);
238 gMC->Gdtom(dir, mom, 2);
242 mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
243 //printf("Dir %f %f %f\n",dir[0],dir[1],dir[2]);
244 //printf("Momentum %15.12f %15.12f %15.12f\n",mom[0],mom[1],mom[2]);
245 //printf("Energy %e\n", mom[3]);
261 for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
262 if (!vmod) for(j=0;j<3;j++) {
268 for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
269 if (!vmod) for(j=0;j<3;j++) {
276 for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
277 vmod=TMath::Sqrt(1/vmod);
278 for(j=0;j<3;j++) e1[j]*=vmod;
281 for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
282 vmod=TMath::Sqrt(1/vmod);
283 for(j=0;j<3;j++) e2[j]*=vmod;
285 //gMC->Rndm(ranf, 1);
286 gMC->GetRandom()->RndmArray(1,ranf);
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);
291 // Put photon on the stack and label it as feedback (51, 52)
294 gAlice->PushTrack(Int_t(1), gAlice->GetCurrentTrackNumber(), Int_t(50000051),
295 mom[0],mom[1],mom[2],mom[3],source[0],source[1],source[2],position[3],pol[0],pol[1],pol[2],
296 kPFeedBackPhoton, nt, 1.);
298 //printf("Adding feedback with tof %f and going to %f %f %f\n",position[3],mom[0],mom[1],mom[2]);
301 //printf("feedbacks produced:%d\n",sNfeed);