[u/mrichter/AliRoot.git] / RICH / AliRICHResponseV0.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.                  *
 **************************************************************************/

/* $Id$ */

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

#include "AliRICHResponseV0.h"
#include "AliRun.h"
#include "AliSegmentation.h"

//___________________________________________
ClassImp(AliRICHResponseV0)

AliRICHResponseV0::AliRICHResponseV0()
{
   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
}//AliRICHResponseV0::ctor()
Float_t AliRICHResponseV0::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)
      {
	//printf("Voltage:%d, Yhit:%f\n",fVoltage, yhit);

	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; 
	//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;
}

Float_t AliRICHResponseV0::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;
}


// -------------------------------------------
Float_t AliRICHResponseV0::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;       
    
}

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