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