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