[u/mrichter/AliRoot.git] / RICH / AliRICHResponse.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


#include <TMath.h>
#include <TParticle.h>
#include <TRandom.h>
#include <TVirtualMC.h>

#include "AliRICHResponse.h"
#include "AliRun.h"
#include "AliSegmentation.h"
#include "AliMC.h"

ClassImp(AliRICHResponse)
//__________________________________________________________________________________________________
AliRICHResponse::AliRICHResponse()
{
   SetSigmaIntegration(5.);
   SetChargeSlope(27.);
   SetChargeSpread(0.18, 0.18);
   SetMaxAdc(4096);
   SetAlphaFeedback(0.036);
   SetEIonisation(26.e-9);
   SetSqrtKx3(0.77459667);
   SetKx2(0.962);
   SetKx4(0.379);
   SetSqrtKy3(0.77459667);
   SetKy2(0.962);
   SetKy4(0.379);
   SetPitch(0.25);
   SetWireSag(1);		      // 1->On, 0->Off
   SetVoltage(2150);		      // Should only be 2000, 2050, 2100 or 2150
}//AliRICHResponse::ctor()
//__________________________________________________________________________________________________
Float_t AliRICHResponse::IntPH(Float_t eloss, Float_t yhit)
{
    // Get number of electrons and return charge
    
    Int_t nel;
    nel= Int_t(eloss/fEIonisation);
    
    Float_t charge=0;
    Double_t gain_var=1;

    if (nel == 0) nel=1;

    if (fWireSag)
      {
	if (fVoltage==2150)
	  {
	    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;
	  }
	if (fVoltage==2100)
	    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;
	if (fVoltage==2050)
	    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;
	if (fVoltage==2000)
	    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;
		
	gain_var = gain_var/100;
	//printf("Yhit:%f, Gain variation:%f\n",yhit,gain_var);

	Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9; 
	//printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);

	for (Int_t i=1;i<=nel;i++) {
	  charge -= gain*TMath::Log(gRandom->Rndm());    
	}
      }
    else
      {
	for (Int_t i=1;i<=nel;i++) {
	  charge -= fChargeSlope*TMath::Log(gRandom->Rndm());    
	}
      }

    return charge;
}//InitPH()
//__________________________________________________________________________________________________
Float_t AliRICHResponse::IntPH(Float_t yhit)
{//  Get number of electrons and return charge, for a single photon

  Float_t charge=0;
  Double_t gain_var=1;

   if (fWireSag)
      {
	if (fVoltage==2150)
	  {
	    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;
	    //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;
	  }
	if (fVoltage==2100)
	    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;
	if (fVoltage==2050)
	    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;
	if (fVoltage==2000)
	    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;

	gain_var = gain_var/100;
	//printf(" Yhit:%f, Gain variation:%f\n",yhit, gain_var);
	
	Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9; 
	
	charge -= gain*TMath::Log(gRandom->Rndm());
	//printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
      }
   else
     {
       charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
     }
    return charge;
}//IntPH()
//__________________________________________________________________________________________________
Float_t AliRICHResponse::IntXY(AliSegmentation * segmentation)
{
    
    const Float_t kInversePitch = 1/fPitch;
    Float_t response;

    //  Integration limits defined by segmentation model
    
    Float_t xi1, xi2, yi1, yi2;
    segmentation->IntegrationLimits(xi1,xi2,yi1,yi2);

    xi1=xi1*kInversePitch;
    xi2=xi2*kInversePitch;
    yi1=yi1*kInversePitch;
    yi2=yi2*kInversePitch;

    //printf("Integration Limits: %f-%f, %f-%f\n",xi1,xi2,yi1,yi2);
    
    //printf("KInversePitch:%f\n",kInversePitch);

    //
// The Mathieson function 
    Double_t ux1=fSqrtKx3*TMath::TanH(fKx2*xi1);
    Double_t ux2=fSqrtKx3*TMath::TanH(fKx2*xi2);
    
    Double_t uy1=fSqrtKy3*TMath::TanH(fKy2*yi1);
    Double_t uy2=fSqrtKy3*TMath::TanH(fKy2*yi2);

    //printf("Integration Data: %f-%f, %f-%f\n",ux1,ux2,uy1,uy2);
    
    //printf("%f %f %f %f\n",fSqrtKx3,fKx2,fKy4,fKx4);
    
    response=4.*fKx4*(TMath::ATan(ux2)-TMath::ATan(ux1))*fKy4*(TMath::ATan(uy2)-TMath::ATan(uy1));

    //printf("Response:%f\n",response);

    return response;       
    
}//IntXY
//__________________________________________________________________________________________________
Int_t AliRICHResponse::FeedBackPhotons(Float_t *source, Float_t qtot)
{  // Generate FeedBack photons
  Int_t j, ipart, nt;
    
  Int_t sNfeed=0;
  
  
  // Local variables 
  Double_t ranf[2];
  Float_t cthf, phif, enfp = 0, sthf;
  Int_t i, ifeed;
  Float_t e1[3], e2[3], e3[3];
  Float_t vmod, uswop;
  Float_t fp;
  Double_t random;
  Float_t dir[3], phi;
  Int_t nfp;
  Float_t pol[3], mom[4];
  TLorentzVector position;
  //
  // Determine number of feedback photons

  //  Get weight of current particle
  TParticle *current = (TParticle*) 
    (*gAlice->GetMCApp()->Particles())[gAlice->GetMCApp()->GetCurrentTrackNumber()];
    
  ifeed = Int_t(current->GetWeight()/100+0.5);
  ipart = gMC->TrackPid();
  fp = fAlphaFeedback * qtot;
  nfp = gRandom->Poisson(fp);
  
  // This call to fill the time of flight
  gMC->TrackPosition(position);
  //printf("Track position: %f %f %f %15.12f\n", position[0],position[1],position[2],position[3]);
  //
  // Generate photons
  for (i = 0; i <nfp; i++) {
	
    // Direction
    gMC->GetRandom()->RndmArray(2,ranf);
    cthf = ranf[0] * 2 - 1.;
    if (cthf < 0)  continue;
    sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
    phif = ranf[1] * 2 * TMath::Pi();
    //
    //gMC->Rndm(&random, 1);
    gMC->GetRandom()->RndmArray(1, &random);
    if (random <= .57) {
      enfp = 7.5e-9;
    } else if (random <= .7) {
      enfp = 6.4e-9;
    } else {
      enfp = 7.9e-9;
    }

    dir[0] = sthf * TMath::Sin(phif);
    dir[1] = cthf;
    dir[2] = sthf * TMath::Cos(phif);
    gMC->Gdtom(dir, mom, 2);
    mom[0]*=enfp;
    mom[1]*=enfp;
    mom[2]*=enfp;
    mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
    //printf("Dir %f %f %f\n",dir[0],dir[1],dir[2]);
    //printf("Momentum %15.12f %15.12f %15.12f\n",mom[0],mom[1],mom[2]);
    //printf("Energy %e\n", mom[3]);
    
    // Polarisation
    e1[0] = 0;
    e1[1] = -dir[2];
    e1[2] = dir[1];
    
    e2[0] = -dir[1];
    e2[1] = dir[0];
    e2[2] = 0;
    
    e3[0] = dir[1];
    e3[1] = 0;
    e3[2] = -dir[0];
    
    vmod=0;
    for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
    if (!vmod) for(j=0;j<3;j++) {
      uswop=e1[j];
      e1[j]=e3[j];
      e3[j]=uswop;
    }
    vmod=0;
    for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
    if (!vmod) for(j=0;j<3;j++) {
      uswop=e2[j];
      e2[j]=e3[j];
      e3[j]=uswop;
    }
    
    vmod=0;
    for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
    vmod=TMath::Sqrt(1/vmod);
    for(j=0;j<3;j++) e1[j]*=vmod;
    
    vmod=0;
    for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
    vmod=TMath::Sqrt(1/vmod);
    for(j=0;j<3;j++) e2[j]*=vmod;
    
    //gMC->Rndm(ranf, 1);
    gMC->GetRandom()->RndmArray(1,ranf);
    phi = ranf[0] * 2 * TMath::Pi();
    for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
    gMC->Gdtom(pol, pol, 2);
    
    // Put photon on the stack and label it as feedback (51, 52) 
    ++sNfeed;

    gAlice->GetMCApp()->PushTrack(Int_t(1), gAlice->GetMCApp()->GetCurrentTrackNumber(), Int_t(50000051),
		     mom[0],mom[1],mom[2],mom[3],source[0],source[1],source[2],position[3],pol[0],pol[1],pol[2],
		     kPFeedBackPhoton, nt, 1.);
    
  }
  return(sNfeed);
}//FeedBackPhotons()
Commit	Line	Data
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
c712cb2f	16
	17	#include <TMath.h>
	18	#include <TParticle.h>
	19	#include <TRandom.h>
	20	#include <TVirtualMC.h>
803d1ab0	21
237c933d	22	#include "AliRICHResponse.h"
c712cb2f	23	#include "AliRun.h"
	24	#include "AliSegmentation.h"
	25	#include "AliMC.h"
237c933d	26
237c933d	27	ClassImp(AliRICHResponse)
c712cb2f	28	//__________________________________________________________________________________________________
	29	AliRICHResponse::AliRICHResponse()
	30	{
	31	SetSigmaIntegration(5.);
	32	SetChargeSlope(27.);
	33	SetChargeSpread(0.18, 0.18);
	34	SetMaxAdc(4096);
	35	SetAlphaFeedback(0.036);
	36	SetEIonisation(26.e-9);
	37	SetSqrtKx3(0.77459667);
	38	SetKx2(0.962);
	39	SetKx4(0.379);
	40	SetSqrtKy3(0.77459667);
	41	SetKy2(0.962);
	42	SetKy4(0.379);
	43	SetPitch(0.25);
	44	SetWireSag(1); // 1->On, 0->Off
	45	SetVoltage(2150); // Should only be 2000, 2050, 2100 or 2150
	46	}//AliRICHResponse::ctor()
	47	//__________________________________________________________________________________________________
	48	Float_t AliRICHResponse::IntPH(Float_t eloss, Float_t yhit)
	49	{
	50	// Get number of electrons and return charge
	51
	52	Int_t nel;
	53	nel= Int_t(eloss/fEIonisation);
	54
	55	Float_t charge=0;
	56	Double_t gain_var=1;
	57
	58	if (nel == 0) nel=1;
	59
	60	if (fWireSag)
	61	{
	62	if (fVoltage==2150)
	63	{
	64	gain_var = 9e-6TMath::Power(yhit,4) + 2e-7TMath::Power(yhit,3) - 0.0316TMath::Power(yhit,2) - 3e-4yhit + 25.367;
	65	}
	66	if (fVoltage==2100)
	67	gain_var = 8e-6TMath::Power(yhit,4) + 2e-7TMath::Power(yhit,3) - 0.0283TMath::Power(yhit,2) - 2e-4yhit + 23.015;
	68	if (fVoltage==2050)
	69	gain_var = 7e-6TMath::Power(yhit,4) + 1e-7TMath::Power(yhit,3) - 0.0254TMath::Power(yhit,2) - 2e-4yhit + 20.888;
	70	if (fVoltage==2000)
	71	gain_var = 6e-6TMath::Power(yhit,4) + 8e-8TMath::Power(yhit,3) - 0.0227TMath::Power(yhit,2) - 1e-4yhit + 18.961;
	72
	73	gain_var = gain_var/100;
	74	//printf("Yhit:%f, Gain variation:%f\n",yhit,gain_var);
	75
	76	Float_t gain = (fChargeSlope + fChargeSlopegain_var).9;
	77	//printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
	78
	79	for (Int_t i=1;i<=nel;i++) {
	80	charge -= gain*TMath::Log(gRandom->Rndm());
	81	}
	82	}
	83	else
	84	{
	85	for (Int_t i=1;i<=nel;i++) {
	86	charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
	87	}
	88	}
	89
	90	return charge;
	91	}//InitPH()
92	//__________________________________________________________________________________________________
93	Float_t AliRICHResponse::IntPH(Float_t yhit)
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-6TMath::Power(yhit,4) + 2e-7TMath::Power(yhit,3) - 0.0316TMath::Power(yhit,2) - 3e-4yhit + 25.367;
104	//gain_var = 9e-5TMath::Power(yhit,4) + 2e-6TMath::Power(yhit,3) - 0.316TMath::Power(yhit,2) - 3e-3yhit + 253.67;
105	}
106	if (fVoltage==2100)
107	gain_var = 8e-6TMath::Power(yhit,4) + 2e-7TMath::Power(yhit,3) - 0.0283TMath::Power(yhit,2) - 2e-4yhit + 23.015;
108	if (fVoltage==2050)
109	gain_var = 7e-6TMath::Power(yhit,4) + 1e-7TMath::Power(yhit,3) - 0.0254TMath::Power(yhit,2) - 2e-4yhit + 20.888;
110	if (fVoltage==2000)
111	gain_var = 6e-6TMath::Power(yhit,4) + 8e-8TMath::Power(yhit,3) - 0.0227TMath::Power(yhit,2) - 1e-4yhit + 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 + fChargeSlopegain_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	}//IntPH()
127	//__________________________________________________________________________________________________
128	Float_t AliRICHResponse::IntXY(AliSegmentation * segmentation)
129	{
130
131	const Float_t kInversePitch = 1/fPitch;
132	Float_t response;
133
134	// Integration limits defined by segmentation model
135
136	Float_t xi1, xi2, yi1, yi2;
137	segmentation->IntegrationLimits(xi1,xi2,yi1,yi2);
138
139	xi1=xi1*kInversePitch;
140	xi2=xi2*kInversePitch;
141	yi1=yi1*kInversePitch;
142	yi2=yi2*kInversePitch;
143
144	//printf("Integration Limits: %f-%f, %f-%f\n",xi1,xi2,yi1,yi2);
145
146	//printf("KInversePitch:%f\n",kInversePitch);
147
148	//
149	// The Mathieson function
150	Double_t ux1=fSqrtKx3TMath::TanH(fKx2xi1);
151	Double_t ux2=fSqrtKx3TMath::TanH(fKx2xi2);
152
153	Double_t uy1=fSqrtKy3TMath::TanH(fKy2yi1);
154	Double_t uy2=fSqrtKy3TMath::TanH(fKy2yi2);
155
156	//printf("Integration Data: %f-%f, %f-%f\n",ux1,ux2,uy1,uy2);
157
158	//printf("%f %f %f %f\n",fSqrtKx3,fKx2,fKy4,fKx4);
159
160	response=4.fKx4(TMath::ATan(ux2)-TMath::ATan(ux1))fKy4(TMath::ATan(uy2)-TMath::ATan(uy1));
161
162	//printf("Response:%f\n",response);
163
164	return response;
165
166	}//IntXY
167	//__________________________________________________________________________________________________
168	Int_t AliRICHResponse::FeedBackPhotons(Float_t *source, Float_t qtot)
169	{ // Generate FeedBack photons
170	Int_t j, ipart, nt;
171
172	Int_t sNfeed=0;
173
174
175	// Local variables
176	Double_t ranf[2];
177	Float_t cthf, phif, enfp = 0, sthf;
178	Int_t i, ifeed;
179	Float_t e1[3], e2[3], e3[3];
180	Float_t vmod, uswop;
181	Float_t fp;
182	Double_t random;
183	Float_t dir[3], phi;
184	Int_t nfp;
185	Float_t pol[3], mom[4];
186	TLorentzVector position;
187	//
188	// Determine number of feedback photons
189
190	// Get weight of current particle
191	TParticle current = (TParticle)
192	(*gAlice->GetMCApp()->Particles())[gAlice->GetMCApp()->GetCurrentTrackNumber()];
193
194	ifeed = Int_t(current->GetWeight()/100+0.5);
195	ipart = gMC->TrackPid();
196	fp = fAlphaFeedback * qtot;
197	nfp = gRandom->Poisson(fp);
198
199	// This call to fill the time of flight
200	gMC->TrackPosition(position);
201	//printf("Track position: %f %f %f %15.12f\n", position[0],position[1],position[2],position[3]);
202	//
203	// Generate photons
204	for (i = 0; i <nfp; i++) {
205
206	// Direction
207	gMC->GetRandom()->RndmArray(2,ranf);
208	cthf = ranf[0] * 2 - 1.;
209	if (cthf < 0) continue;
210	sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
211	phif = ranf[1] * 2 * TMath::Pi();
212	//
213	//gMC->Rndm(&random, 1);
214	gMC->GetRandom()->RndmArray(1, &random);
215	if (random <= .57) {
216	enfp = 7.5e-9;
217	} else if (random <= .7) {
218	enfp = 6.4e-9;
219	} else {
220	enfp = 7.9e-9;
221	}
222
223	dir[0] = sthf * TMath::Sin(phif);
224	dir[1] = cthf;
225	dir[2] = sthf * TMath::Cos(phif);
226	gMC->Gdtom(dir, mom, 2);
227	mom[0]*=enfp;
228	mom[1]*=enfp;
229	mom[2]*=enfp;
230	mom[3] = TMath::Sqrt(mom[0]mom[0]+mom[1]mom[1]+mom[2]*mom[2]);
231	//printf("Dir %f %f %f\n",dir[0],dir[1],dir[2]);
232	//printf("Momentum %15.12f %15.12f %15.12f\n",mom[0],mom[1],mom[2]);
233	//printf("Energy %e\n", mom[3]);
234
235	// Polarisation
236	e1[0] = 0;
237	e1[1] = -dir[2];
238	e1[2] = dir[1];
239
240	e2[0] = -dir[1];
241	e2[1] = dir[0];
242	e2[2] = 0;
243
244	e3[0] = dir[1];
245	e3[1] = 0;
246	e3[2] = -dir[0];
247
248	vmod=0;
249	for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
250	if (!vmod) for(j=0;j<3;j++) {
251	uswop=e1[j];
252	e1[j]=e3[j];
253	e3[j]=uswop;
254	}
255	vmod=0;
256	for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
257	if (!vmod) for(j=0;j<3;j++) {
258	uswop=e2[j];
259	e2[j]=e3[j];
260	e3[j]=uswop;
261	}
262
263	vmod=0;
264	for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
265	vmod=TMath::Sqrt(1/vmod);
266	for(j=0;j<3;j++) e1[j]*=vmod;
267
268	vmod=0;
269	for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
270	vmod=TMath::Sqrt(1/vmod);
271	for(j=0;j<3;j++) e2[j]*=vmod;
272
273	//gMC->Rndm(ranf, 1);
274	gMC->GetRandom()->RndmArray(1,ranf);
275	phi = ranf[0] * 2 * TMath::Pi();
276	for(j=0;j<3;j++) pol[j]=e1[j]TMath::Sin(phi)+e2[j]TMath::Cos(phi);
277	gMC->Gdtom(pol, pol, 2);
278
279	// Put photon on the stack and label it as feedback (51, 52)
280	++sNfeed;
281
282	gAlice->GetMCApp()->PushTrack(Int_t(1), gAlice->GetMCApp()->GetCurrentTrackNumber(), Int_t(50000051),
283	mom[0],mom[1],mom[2],mom[3],source[0],source[1],source[2],position[3],pol[0],pol[1],pol[2],
284	kPFeedBackPhoton, nt, 1.);
285
286	}
287	return(sNfeed);
288	}//FeedBackPhotons()