| 1 |
#include "PamVMC2PiAux.h" |
| 2 |
|
| 3 |
void PamVMC2PiAux::gsurf(){ |
| 4 |
|
| 5 |
Double_t th1 = DegToRad()*17.; |
| 6 |
Double_t th2 = DegToRad()*20.; |
| 7 |
Double_t z1 = 87.4389954 + 0.150935784 + 0.4*Sin(th1) + 18.15*Cos(th1)/2.; |
| 8 |
Double_t z2 = 87.4389954 + 0.00949157123 + 0.4*Sin(th2) + 18.15*Cos(th2)/2.; |
| 9 |
Double_t xcard = GPAXCARD - 0.4*Cos(th2) + 18.15*Sin(th2)/2.; |
| 10 |
Double_t ycard = GPAYCARD - 0.4*Cos(th1) + ( z2-(87.4389954 + 0.150935784 + 0.4*Sin(th1)) )*Sin(th1); |
| 11 |
Double_t th = ATan(Sqrt(Power((xcard-XS2),2.) + Power((ycard-YS2),2.))/(z2-ZS2MAX)); |
| 12 |
Double_t ph = ATan((ycard-YS2)/(xcard-XS2)); |
| 13 |
|
| 14 |
if (TWOPIDEBUG){ |
| 15 |
cout<<"th1="<<th1<<" th2="<<th2<<" z1="<<z1<<" z2="<<z2<<endl; |
| 16 |
cout<<"xcard="<<xcard<<" ycard="<<ycard<<" th="<<th<<" ph="<<ph<<endl; |
| 17 |
} |
| 18 |
|
| 19 |
Double_t b = 2.*(XS2*Sin(th)*Cos(ph)+YS2*Sin(ph)*Sin(th) + ZS2MAX*Cos(th)+SHELL3RZ*Cos(th)); |
| 20 |
Double_t c = XS2*XS2 + YS2*YS2 + (ZS2MAX+SHELL3RZ)*(ZS2MAX+SHELL3RZ) - SHELL3RR*SHELL3RR; |
| 21 |
Double_t l = (-b+Sqrt(b*b-4.*c))/2.; |
| 22 |
|
| 23 |
Double_t ZMINS = Sqrt(Power(SHELL3RR,2.) - Power(RCY,2.)) - SHELL3RZ; |
| 24 |
fzmins = ZMINS; |
| 25 |
|
| 26 |
|
| 27 |
// if (TWOPIDEBUG) cout<<"ZMINS="<<ZMINS<<" RCY="<<RCY<<endl; |
| 28 |
|
| 29 |
Double_t h1 = (SHELL3RR-SHELL3RZ-ZMINS); |
| 30 |
Double_t SPHAR = TwoPi()*h1*SHELL3RR; |
| 31 |
|
| 32 |
|
| 33 |
Double_t ta = SPHAR; |
| 34 |
Double_t THMAX = ATan(RCY/(ZMINS+SHELL3RZ)); |
| 35 |
fthmax = THMAX; |
| 36 |
if (TWOPIDEBUG){ |
| 37 |
cout<<"Spherical surface:"<<endl; |
| 38 |
cout<<" Zcenter: "<<-SHELL3RZ<<" cm, Radius: "<<SHELL3RR |
| 39 |
<<" cm, MaxAng: "<<THMAX*RadToDeg()<<endl; |
| 40 |
cout<<" Spherical surface area: "<<SPHAR<<" cm^2"<<endl; |
| 41 |
} |
| 42 |
|
| 43 |
Double_t tgy = (ZS1-z2)/(YS1+ycard); |
| 44 |
Double_t ry = tgy/(RCY+YS1); |
| 45 |
Double_t h2 = (ZS2MAX+l*Cos(th)-ZS1); |
| 46 |
h2 = ry + h2; |
| 47 |
Double_t ZMINC = ZMINS - h2; |
| 48 |
fzminc = ZMINC; |
| 49 |
|
| 50 |
Double_t CYLAR = TwoPi()*h2*RCY; |
| 51 |
|
| 52 |
if (TWOPIDEBUG){ |
| 53 |
cout<<"Cylindrical surface:"<<endl; |
| 54 |
cout<<" Height of cylinder: "<<h2<<" cm, Radius: "<<RCY<<" cm"<<endl; |
| 55 |
cout<<" Zmin: "<<fzminc<<" cm, Zmax: "<<fzmins<<" cm"<<endl; |
| 56 |
cout<<" Cylindrical surface area: "<<CYLAR<<" cm^2"<<endl; |
| 57 |
} |
| 58 |
|
| 59 |
ta += CYLAR; |
| 60 |
|
| 61 |
if(TWOPIDEBUG) cout<<"Total generation surface area = "<<ta<<" cm^2"<<endl; |
| 62 |
|
| 63 |
SPHAR/=ta; |
| 64 |
CYLAR/=ta; |
| 65 |
|
| 66 |
if(TWOPIDEBUG) cout<<"SPHAR,CYLAR="<<SPHAR<<","<<CYLAR<<endl; |
| 67 |
|
| 68 |
fsphar = SPHAR; |
| 69 |
farea = ta; |
| 70 |
} |
| 71 |
|
| 72 |
void PamVMC2PiAux::GenPosAng(TVector3& pos, TVector2& ang){ |
| 73 |
|
| 74 |
Double_t phi = TwoPi()*frandom->Uniform(); |
| 75 |
Double_t the; |
| 76 |
Double_t ct, st, cf, sf, xbeam, ybeam, zbeam; |
| 77 |
|
| 78 |
|
| 79 |
//POSITIONS |
| 80 |
if( (frandom->Uniform()) > fsphar ){ |
| 81 |
//generation on the lateral cylindrical surface |
| 82 |
the = Pi()/2.; |
| 83 |
ct = Cos(the); |
| 84 |
st = Sin(the); |
| 85 |
cf = Cos(phi); |
| 86 |
sf = Sin(phi); |
| 87 |
xbeam = RCY * cf; |
| 88 |
ybeam = RCY * sf; |
| 89 |
zbeam = -(fzminc+frandom->Uniform()*(fzmins-fzminc)); |
| 90 |
} else { |
| 91 |
//generation from top spherical surface |
| 92 |
the = ASin( Sqrt(frandom->Uniform())*Sin(fthmax) ); //CORRECT! |
| 93 |
ct = Cos(the); |
| 94 |
st = Sin(the); |
| 95 |
cf = Cos(phi); |
| 96 |
sf = Sin(phi); |
| 97 |
xbeam = SHELL3RR * st * cf; |
| 98 |
ybeam = SHELL3RR * st * sf; |
| 99 |
zbeam = -SHELL3RR * ct + SHELL3RZ; |
| 100 |
} |
| 101 |
|
| 102 |
pos.SetXYZ(xbeam,ybeam,-zbeam); |
| 103 |
|
| 104 |
//DIRECTION COSINES |
| 105 |
// local angles. COS^2 generation (theta_min=0, theta_max=90) |
| 106 |
Double_t lthe = ASin(Sqrt(frandom->Uniform())); |
| 107 |
Double_t lphi = TwoPi() * frandom->Uniform(); |
| 108 |
|
| 109 |
// local director cosines |
| 110 |
|
| 111 |
Double_t aaa = Sin(lthe) * Cos(lphi); |
| 112 |
Double_t bbb = Sin(lthe) * Sin(lphi); |
| 113 |
Double_t ccc = Cos(lthe); |
| 114 |
|
| 115 |
// global director cosines |
| 116 |
//MY CORRECT COSINES |
| 117 |
|
| 118 |
|
| 119 |
Double_t ubeam = aaa*cf*ct - bbb*sf - ccc*cf*st; //cos(px) |
| 120 |
Double_t vbeam = aaa*sf*ct + bbb*cf - ccc*sf*st; //cos(py) |
| 121 |
Double_t wbeam = aaa*st + ccc*ct; //cos(pz) |
| 122 |
|
| 123 |
Double_t theta_dir = ACos(wbeam); |
| 124 |
Double_t phi_dir = ATan2(vbeam,ubeam); |
| 125 |
if(phi_dir<0.) phi_dir+=TwoPi(); |
| 126 |
|
| 127 |
ang.Set(theta_dir,phi_dir); |
| 128 |
} |
Parent Directory
| 
Revision Log