[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.                  *
 **************************************************************************/

/*
  $Log$
  Revision 1.10  2002/10/14 14:57:32  hristov
  Merging the VirtualMC branch to the main development branch (HEAD)

  Revision 1.9.6.1  2002/06/10 15:12:46  hristov
  Merged with v3-08-02

  Revision 1.9  2001/11/08 20:10:45  dibari
  ctor with initialising of parmeters added

  Revision 1.8  2001/09/05 09:09:38  hristov
  The energy of feedback photons calculated correctly

  Revision 1.7  2001/05/10 12:32:27  jbarbosa
  Changed call to SetTrack.

  Revision 1.6  2001/02/23 17:39:02  jbarbosa
  Removed verbose output.

  Revision 1.5  2001/02/23 17:25:08  jbarbosa
  Re-definition of IntPH() to accomodate for wire sag effect.

  Revision 1.4  2000/12/01 17:37:44  morsch
  Replace in argument of SetTrack(..) string constant by  kPFeedBackPhoton.

  Revision 1.3  2000/10/03 21:44:09  morsch
  Use AliSegmentation and AliHit abstract base classes.

  Revision 1.2  2000/10/02 21:28:12  fca
  Removal of useless dependecies via forward declarations

  Revision 1.1  2000/06/12 15:29:37  jbarbosa
  Cleaned up version.

*/

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

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