| 1 |
// $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03 |
| 2 |
|
| 3 |
|
| 4 |
#include <TVirtualMC.h> |
| 5 |
#include <TVirtualMCStack.h> |
| 6 |
#include <TPDGCode.h> |
| 7 |
#include <TDatabasePDG.h> |
| 8 |
#include <TParticlePDG.h> |
| 9 |
#include <TVector3.h> |
| 10 |
#include <TMath.h> |
| 11 |
|
| 12 |
#include "PamVMCPrimaryGenerator.h" |
| 13 |
|
| 14 |
using TMath::Sqrt; |
| 15 |
using TMath::Sin; |
| 16 |
using TMath::Cos; |
| 17 |
using TMath::ACos; |
| 18 |
using TMath::Tan; |
| 19 |
using TMath::ATan; |
| 20 |
using TMath::ATan2; |
| 21 |
using TMath::Log; |
| 22 |
using TMath::Power; |
| 23 |
using TMath::Exp; |
| 24 |
using TMath::Pi; |
| 25 |
using TMath::Abs; |
| 26 |
|
| 27 |
ClassImp(PamVMCPrimary) |
| 28 |
|
| 29 |
PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b) |
| 30 |
{ |
| 31 |
a.fPDG=b.fPDG; |
| 32 |
a.fX0=b.fX0; |
| 33 |
a.fY0=b.fY0; |
| 34 |
a.fZ0=b.fZ0; |
| 35 |
a.fTHETA=b.fTHETA; |
| 36 |
a.fPHI=b.fPHI; |
| 37 |
a.fP0=b.fP0; |
| 38 |
a.fGOOD=b.fGOOD; |
| 39 |
|
| 40 |
return a; |
| 41 |
} |
| 42 |
|
| 43 |
|
| 44 |
ClassImp(PamVMCPrimaryGenerator) |
| 45 |
|
| 46 |
PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack, UInt_t seed) |
| 47 |
: TObject(), |
| 48 |
fStack(stack), |
| 49 |
fevno(0), |
| 50 |
fmass(0.), |
| 51 |
fcharge(0.), |
| 52 |
frnd(0) |
| 53 |
{ |
| 54 |
// Standard constructor |
| 55 |
fprimColl = new TClonesArray("PamVMCPrimary"); |
| 56 |
frnd = new TRandom3(seed); |
| 57 |
|
| 58 |
fprim.fPDG=kProton; |
| 59 |
fprim.fX0=1.; |
| 60 |
fprim.fY0=1.; |
| 61 |
fprim.fZ0=130.; |
| 62 |
fprim.fTHETA=0.; |
| 63 |
fprim.fPHI=0.; |
| 64 |
fprim.fP0=1.; //1GV |
| 65 |
|
| 66 |
} |
| 67 |
|
| 68 |
// PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(UInt_t seed) |
| 69 |
// : TObject(), |
| 70 |
// fStack(0), |
| 71 |
// fevno(0), |
| 72 |
// fmass(0.), |
| 73 |
// fcharge(0.), |
| 74 |
// fprimColl(0), |
| 75 |
// frnd(0) |
| 76 |
// { |
| 77 |
// frnd = new TRandom3(seed); |
| 78 |
// // Default constructor |
| 79 |
// //Default primary proton |
| 80 |
// fprim.fPDG=kProton; |
| 81 |
// fprim.fX0=1.; |
| 82 |
// fprim.fY0=1.; |
| 83 |
// fprim.fZ0=130.; |
| 84 |
// fprim.fTHETA=0.; |
| 85 |
// fprim.fPHI=0.; |
| 86 |
// fprim.fP0=1.; //1GV |
| 87 |
// } |
| 88 |
|
| 89 |
PamVMCPrimaryGenerator::PamVMCPrimaryGenerator() |
| 90 |
: TObject(), |
| 91 |
fStack(0), |
| 92 |
fevno(0), |
| 93 |
fmass(0.), |
| 94 |
fcharge(0.), |
| 95 |
fprimColl(0), |
| 96 |
frnd(0) |
| 97 |
{ |
| 98 |
frnd = new TRandom3(0); |
| 99 |
// Default constructor |
| 100 |
//Default primary proton |
| 101 |
fprim.fPDG=kProton; |
| 102 |
fprim.fX0=1.; |
| 103 |
fprim.fY0=1.; |
| 104 |
fprim.fZ0=130.; |
| 105 |
fprim.fTHETA=0.; |
| 106 |
fprim.fPHI=0.; |
| 107 |
fprim.fP0=1.; //1GV |
| 108 |
} |
| 109 |
|
| 110 |
PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator() |
| 111 |
{ |
| 112 |
// Destructor |
| 113 |
delete frnd; |
| 114 |
delete fprimColl; |
| 115 |
} |
| 116 |
|
| 117 |
// private methods |
| 118 |
|
| 119 |
#include <Riostream.h> |
| 120 |
|
| 121 |
void PamVMCPrimaryGenerator::GeneratePrimary() |
| 122 |
{ |
| 123 |
// Add one primary particle to the user stack (derived from TVirtualMCStack). |
| 124 |
|
| 125 |
// Track ID (filled by stack) |
| 126 |
Int_t ntr; |
| 127 |
|
| 128 |
// Option: to be tracked |
| 129 |
Int_t toBeDone = 1; |
| 130 |
|
| 131 |
// Particle type |
| 132 |
Int_t pdg = fprim.fPDG; |
| 133 |
|
| 134 |
Double_t fvx, fvy, fvz; |
| 135 |
fvx=fprim.fX0; |
| 136 |
fvy=fprim.fY0; |
| 137 |
fvz=fprim.fZ0; |
| 138 |
// Position |
| 139 |
|
| 140 |
Double_t tof = 0.; |
| 141 |
|
| 142 |
// Energy (in GeV) |
| 143 |
Double_t kinEnergy = MomentumToKinE(fprim.fP0); |
| 144 |
Double_t e = fmass + kinEnergy; |
| 145 |
|
| 146 |
// Particle momentum |
| 147 |
Double_t px, py, pz; |
| 148 |
|
| 149 |
//px = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Cos((Pi()/180.)*(fprim.fPHI)); // ritabrata |
| 150 |
//py = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Sin((Pi()/180.)*(fprim.fPHI)); |
| 151 |
//pz = -fprim.fP0*Cos((Pi()/180.)*(fprim.fTHETA)); // converting in radian |
| 152 |
px = fprim.fP0*Sin(fprim.fTHETA)*Cos(fprim.fPHI); // RADIANTS |
| 153 |
py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI); |
| 154 |
pz = -fprim.fP0*Cos(fprim.fTHETA); |
| 155 |
|
| 156 |
// Polarization |
| 157 |
TVector3 polar; |
| 158 |
|
| 159 |
// Add particle to stack |
| 160 |
fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, fvx, fvy, fvz, tof, |
| 161 |
polar.X(), polar.Y(), polar.Z(), |
| 162 |
kPPrimary, ntr, 1., 0); |
| 163 |
|
| 164 |
PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++); |
| 165 |
|
| 166 |
*pc = fprim; |
| 167 |
} |
| 168 |
|
| 169 |
|
| 170 |
void PamVMCPrimaryGenerator::SetParticle(Int_t pdg){ |
| 171 |
fprim.fPDG=pdg; |
| 172 |
//TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fprim.fPDG); |
| 173 |
fmass = (TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Mass(); |
| 174 |
fcharge = ((TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Charge())/3.; |
| 175 |
} |
| 176 |
|
| 177 |
void PamVMCPrimaryGenerator::SetMomentum( |
| 178 |
Double_t px, Double_t py, Double_t pz) |
| 179 |
{ |
| 180 |
fprim.fP0= Sqrt(px*px+py*py+pz*pz); |
| 181 |
fprim.fTHETA=ATan(Sqrt(px*px+py*py)/pz); |
| 182 |
fprim.fPHI=ATan(py/px); |
| 183 |
} |
| 184 |
|
| 185 |
void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy) |
| 186 |
{ |
| 187 |
if(isEnergy) { |
| 188 |
fprim.fP0=frnd->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax)); |
| 189 |
} else{ |
| 190 |
fprim.fP0=frnd->Uniform(PEmin,PEmax); |
| 191 |
} |
| 192 |
|
| 193 |
} |
| 194 |
|
| 195 |
|
| 196 |
|
| 197 |
|
| 198 |
void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy) |
| 199 |
{ |
| 200 |
Double_t alpha = 1.+gamma; //integral spectral index |
| 201 |
if(alpha==0.){ |
| 202 |
fprim.fP0=Exp(Log(PEmin)+frnd->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin))); |
| 203 |
} else { |
| 204 |
if(PEmin==0.) PEmin=1.E-10; |
| 205 |
fprim.fP0=Power((frnd->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha); |
| 206 |
} |
| 207 |
|
| 208 |
if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0); |
| 209 |
|
| 210 |
} |
| 211 |
|
| 212 |
|
| 213 |
// Spherical distribution -- Test by Cecilia P march 2009 -----------------------------/// |
| 214 |
void PamVMCPrimaryGenerator::GenSphPhiThe(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, |
| 215 |
Double_t zmin, Double_t zmax) |
| 216 |
{ |
| 217 |
Bool_t trkGood = kFALSE; |
| 218 |
Double_t theta = 999.; |
| 219 |
Double_t phi = 0.; |
| 220 |
Double_t x2,y2,x3,y3; |
| 221 |
|
| 222 |
//static const Double_t rad2deg = 57.2958; |
| 223 |
|
| 224 |
//static const Double_t s1_xmax=20.4, s1_ymax=16.5, s1_pz=102.8866; |
| 225 |
// calo cavity 8.07x6.57 with z2=71.6 z3=27.4 circa |
| 226 |
static const Double_t s2_xmax=7.8, s2_ymax=6.0, s2_z=73.489; |
| 227 |
static const Double_t s3_xmax=8.0, s3_ymax=6.0, s3_z=25.3159; |
| 228 |
|
| 229 |
// Generate random position unif. distr. |
| 230 |
// GenPosition(-25.,25., 25.,25., 108.0,108.0); |
| 231 |
// Double_t x0= fprim.fX0; |
| 232 |
// Double_t y0= fprim.fY0; |
| 233 |
// Double_t z0= fprim.fZ0; |
| 234 |
Double_t x0= frnd->Uniform(xmin,xmax); |
| 235 |
Double_t y0= frnd->Uniform(ymin,ymax); |
| 236 |
Double_t z0= frnd->Uniform(zmin,zmax); |
| 237 |
Int_t posGood = 0; |
| 238 |
|
| 239 |
|
| 240 |
while (trkGood!=kTRUE) |
| 241 |
{ |
| 242 |
// Generate phi theta |
| 243 |
theta=999.; // init |
| 244 |
while (theta>=0.82) //take only theta smaller than 30deg=0.52rad |
| 245 |
{ |
| 246 |
Double_t xcos = sqrt( frnd->Uniform(0.,1.) ); |
| 247 |
theta = acos(xcos); //RAD |
| 248 |
} |
| 249 |
phi = frnd->Uniform(0.,2.*Pi()); |
| 250 |
|
| 251 |
// Calculate xy at beginning/end of magnet cavity |
| 252 |
Double_t fact2 = (s2_z-z0)/cos(theta); |
| 253 |
x2 = x0 + fabs(fact2) * sin(theta) * cos(phi); |
| 254 |
y2 = y0 + fabs(fact2) * sin(theta) * sin(phi); |
| 255 |
Double_t fact3 = (s3_z-z0)/cos(theta); |
| 256 |
x3 = x0 + fabs(fact3) * sin(theta) * cos(phi); |
| 257 |
y3 = y0 + fabs(fact3) * sin(theta) * sin(phi); |
| 258 |
|
| 259 |
// Test condition on the direction |
| 260 |
if ( Abs(x2) <= Abs(s2_xmax) && Abs(y2) <= Abs(s2_ymax) && |
| 261 |
Abs(x3) <= Abs(s3_xmax) && Abs(y3) <= Abs(s3_ymax) ) { |
| 262 |
trkGood = kTRUE; |
| 263 |
//cout<<" x/y0= "<<x0<<" "<<y0<<" x/y2= "<<fact2*sin(theta)*cos(phi)<<" "<<x2<<" xy3= "<< |
| 264 |
// fact3*sin(theta)*cos(phi)<<" "<<x3<<endl; |
| 265 |
} |
| 266 |
|
| 267 |
//if from current x0y0z0 the condition are not satysfied for any angle => def new start position |
| 268 |
posGood++; |
| 269 |
if ( posGood == 100 ) |
| 270 |
{ |
| 271 |
x0= frnd->Uniform(xmin,xmax); |
| 272 |
y0= frnd->Uniform(ymin,ymax); |
| 273 |
z0= frnd->Uniform(zmin,zmax); |
| 274 |
posGood = 0; |
| 275 |
} |
| 276 |
|
| 277 |
} |
| 278 |
|
| 279 |
// Set direction and position: |
| 280 |
SetDirection(theta, phi); |
| 281 |
SetPosition(x0, y0, z0); |
| 282 |
|
| 283 |
return; |
| 284 |
} |
| 285 |
//--- end Test by Cecilia ------------------------------------/// |
| 286 |
|
| 287 |
|
| 288 |
|
| 289 |
|
| 290 |
void PamVMCPrimaryGenerator::GenSphDist(Double_t r, Double_t Thmin, Double_t Thmax, Double_t Phmin, Double_t Phmax) |
| 291 |
{ |
| 292 |
// all angles in RADIANTS |
| 293 |
Double_t theta, phi, y, f; |
| 294 |
phi = frnd->Uniform(Phmin,Phmax); |
| 295 |
|
| 296 |
do |
| 297 |
{ y = frnd->Uniform(0.,1.); |
| 298 |
theta = frnd->Uniform(Thmin,Thmax); |
| 299 |
f = Sin(theta); |
| 300 |
} while (y>f); |
| 301 |
|
| 302 |
//SetPosition((r*Sin(theta)*Cos(phi)), (r*Sin(theta)*Sin(phi)), (r*Cos(theta))); // ritabrata |
| 303 |
SetPosition(1.,1.,130.); // cecilia |
| 304 |
|
| 305 |
|
| 306 |
//random distribution of theta phi in the angle at the vertex at (0,0,r) |
| 307 |
//by the S3 max distant corners. |
| 308 |
static const Double_t s3_x=9.0, s3_y=7.5, s3_pz=25.3159; |
| 309 |
Double_t ang = ATan((Sqrt(s3_x*s3_x+s3_y*s3_y))/(r-s3_pz)); |
| 310 |
|
| 311 |
//SetDirection((frnd->Uniform((theta-ang),(theta+ang)))/(Pi()/180.), (frnd->Uniform((phi-ang),(phi+ang)))/(Pi()/180.)); |
| 312 |
SetDirection( frnd->Uniform((theta-ang),(theta+ang)) , frnd->Uniform(0.,2*Pi()) ); |
| 313 |
//SetDirection(0., 0.); |
| 314 |
} |
| 315 |
|
| 316 |
void PamVMCPrimaryGenerator::GenSphDist(Double_t r) |
| 317 |
{ |
| 318 |
//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; |
| 319 |
static const Double_t s1_x=20.4, s1_y=16.5, s1_pz=102.8866; |
| 320 |
static const Double_t s2_x=9.0, s2_y=7.5, s2_pz=73.489; |
| 321 |
static const Double_t s3_x=9.0, s3_y=7.5, s3_pz=25.3159; |
| 322 |
|
| 323 |
//calculate max theta and phi angles to be allowed (calculating wrt one |
| 324 |
//corner of s3 to opposite corner of s1) |
| 325 |
Double_t rx = s1_x+s3_x; |
| 326 |
Double_t ry = s1_y+s3_y; |
| 327 |
Double_t rz = s1_pz-s3_pz; |
| 328 |
|
| 329 |
Double_t thmax = (180./Pi())*(ACos(rz/Sqrt(rx*rx+ry*ry+rz*rz))); |
| 330 |
//Double_t phmax = (180./Pi())*(ATan2(ry,rx)); |
| 331 |
|
| 332 |
//cout << "~~~~~~Theta max Phi max : " << thmax <<", "<< phmax << endl; |
| 333 |
|
| 334 |
//generate a track check and let it go only if it is passing through all |
| 335 |
//the 3 TOF |
| 336 |
Bool_t trkGood = kFALSE; |
| 337 |
|
| 338 |
do |
| 339 |
{ |
| 340 |
//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; |
| 341 |
|
| 342 |
GenSphDist(r, 0., thmax, 0., 2*Pi() ); |
| 343 |
|
| 344 |
Double_t th = fprim.fTHETA; |
| 345 |
Double_t ph = fprim.fPHI; |
| 346 |
|
| 347 |
//cout << "~~~~~~Theta Phi : " << fprim.fTHETA <<", "<< fprim.fPHI << endl; |
| 348 |
|
| 349 |
Double_t x1, y1, x2, y2, x3, y3; |
| 350 |
|
| 351 |
x1 = s1_pz*Tan(th)*Cos(ph) - fprim.fX0; |
| 352 |
y1 = s1_pz*Tan(th)*Sin(ph) - fprim.fY0; |
| 353 |
x2 = s2_pz*Tan(th)*Cos(ph) - fprim.fX0; |
| 354 |
y2 = s2_pz*Tan(th)*Sin(ph) - fprim.fY0; |
| 355 |
x3 = s3_pz*Tan(th)*Cos(ph) - fprim.fX0; |
| 356 |
y3 = s3_pz*Tan(th)*Sin(ph) - fprim.fY0; |
| 357 |
|
| 358 |
if ( Abs(x1) <= Abs(s1_x) && Abs(y1) <= Abs(s1_y) && |
| 359 |
Abs(x2) <= Abs(s2_x) && Abs(y2) <= Abs(s2_y) && |
| 360 |
Abs(x3) <= Abs(s3_x) && Abs(y3) <= Abs(s3_y) ) trkGood = kTRUE; |
| 361 |
//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; |
| 362 |
}while (!trkGood); |
| 363 |
//cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; |
| 364 |
|
| 365 |
} |
Parent Directory
| 
Revision Log