| 1 |
From Cecilia.Pizzolotto@ts.infn.it Fri Apr 10 09:30:53 2009 |
| 2 |
Date: Fri, 10 Apr 2009 09:30:35 +0200 |
| 3 |
From: Cecilia Pizzolotto <Cecilia.Pizzolotto@ts.infn.it> |
| 4 |
To: Emiliano Mocchiutti <Emiliano.Mocchiutti@ts.infn.it> |
| 5 |
Subject: spherical distribution |
| 6 |
|
| 7 |
// Spherical distribution -- Test by Cecilia P march 2009 |
| 8 |
-----------------------------/// |
| 9 |
-25 25 -25 25 108 108 |
| 10 |
void PamVMCPrimaryGenerator::GenSphPhiThe(Double_t xmin, Double_t xmax, |
| 11 |
Double_t ymin, Double_t ymax, |
| 12 |
Double_t zmin, Double_t zmax) |
| 13 |
{ |
| 14 |
Bool_t trkGood = kFALSE; |
| 15 |
Double_t theta = 999.; |
| 16 |
Double_t phi = 0.; |
| 17 |
Double_t x2,y2,x3,y3; |
| 18 |
|
| 19 |
//static const Double_t rad2deg = 57.2958; |
| 20 |
|
| 21 |
//static const Double_t s1_xmax=20.4, s1_ymax=16.5, s1_pz=102.8866; |
| 22 |
// calo cavity 8.07x6.57 with z2=71.6 z3=27.4 circa |
| 23 |
static const Double_t s2_xmax=7.8, s2_ymax=6.0, s2_z=73.489; |
| 24 |
static const Double_t s3_xmax=8.0, s3_ymax=6.0, s3_z=25.3159; |
| 25 |
|
| 26 |
// Generate random position unif. distr. |
| 27 |
// GenPosition(-25.,25., 25.,25., 108.0,108.0); |
| 28 |
// Double_t x0= fprim.fX0; |
| 29 |
// Double_t y0= fprim.fY0; |
| 30 |
// Double_t z0= fprim.fZ0; |
| 31 |
Double_t x0= frnd->Uniform(xmin,xmax); |
| 32 |
Double_t y0= frnd->Uniform(ymin,ymax); |
| 33 |
Double_t z0= frnd->Uniform(zmin,zmax); |
| 34 |
Int_t posGood = 0; |
| 35 |
|
| 36 |
|
| 37 |
while (trkGood!=kTRUE) |
| 38 |
{ |
| 39 |
// Generate phi theta |
| 40 |
theta=999.; // init |
| 41 |
while (theta>=0.82) //take only theta smaller than 30deg=0.52rad |
| 42 |
{ |
| 43 |
Double_t xcos = sqrt( frnd->Uniform(0.,1.) ); |
| 44 |
theta = acos(xcos); //RAD |
| 45 |
} |
| 46 |
phi = frnd->Uniform(0.,2.*Pi()); |
| 47 |
|
| 48 |
// Calculate xy at beginning/end of magnet cavity |
| 49 |
Double_t fact2 = (s2_z-z0)/cos(theta); |
| 50 |
x2 = x0 + fabs(fact2) * sin(theta) * cos(phi); |
| 51 |
y2 = y0 + fabs(fact2) * sin(theta) * sin(phi); |
| 52 |
Double_t fact3 = (s3_z-z0)/cos(theta); |
| 53 |
x3 = x0 + fabs(fact3) * sin(theta) * cos(phi); |
| 54 |
y3 = y0 + fabs(fact3) * sin(theta) * sin(phi); |
| 55 |
|
| 56 |
// Test condition on the direction |
| 57 |
if ( Abs(x2) <= Abs(s2_xmax) && Abs(y2) <= Abs(s2_ymax) && |
| 58 |
Abs(x3) <= Abs(s3_xmax) && Abs(y3) <= Abs(s3_ymax) ) { |
| 59 |
trkGood = kTRUE; |
| 60 |
//cout<<" x/y0= "<<x0<<" "<<y0<<" x/y2= "<<fact2*sin(theta)*cos(phi)<<" "<<x2<<" xy3= "<< |
| 61 |
// fact3*sin(theta)*cos(phi)<<" "<<x3<<endl; |
| 62 |
} |
| 63 |
|
| 64 |
//if from current x0y0z0 the condition are not satysfied for any angle => def new start position |
| 65 |
posGood++; |
| 66 |
if ( posGood == 100 ) |
| 67 |
{ |
| 68 |
x0= frnd->Uniform(xmin,xmax); |
| 69 |
y0= frnd->Uniform(ymin,ymax); |
| 70 |
z0= frnd->Uniform(zmin,zmax); |
| 71 |
posGood = 0; |
| 72 |
} |
| 73 |
|
| 74 |
} |
| 75 |
|
| 76 |
// Set direction and position: |
| 77 |
SetDirection(theta, phi); |
| 78 |
SetPosition(x0, y0, z0); |
| 79 |
|
| 80 |
return; |
| 81 |
} |
| 82 |
//--- end Test by Cecilia ------------------------------------/// |
| 83 |
|
Parent Directory
| 
Revision Log