// $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03 #include #include #include #include #include #include #include #include "PamVMCPrimaryGenerator.h" using TMath::Sqrt; using TMath::Sin; using TMath::Cos; using TMath::ACos; using TMath::Tan; using TMath::ATan; using TMath::ATan2; using TMath::Log; using TMath::Power; using TMath::Exp; using TMath::Pi; using TMath::Abs; ClassImp(PamVMCPrimary) PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b) { a.fPDG=b.fPDG; a.fX0=b.fX0; a.fY0=b.fY0; a.fZ0=b.fZ0; a.fTHETA=b.fTHETA; a.fPHI=b.fPHI; a.fP0=b.fP0; a.fGOOD=b.fGOOD; return a; } ClassImp(PamVMCPrimaryGenerator) PamVMCPrimaryGenerator::PamVMCPrimaryGenerator(TVirtualMCStack* stack) : TObject(), fStack(stack), fevno(0), fmass(0.), fcharge(0.), frnd(0) { // Standard constructor fprimColl = new TClonesArray("PamVMCPrimary"); frnd = new TRandom3(0); fprim.fPDG=kProton; fprim.fX0=1.; fprim.fY0=1.; fprim.fZ0=130.; fprim.fTHETA=0.; fprim.fPHI=0.; fprim.fP0=1.; //1GV } PamVMCPrimaryGenerator::PamVMCPrimaryGenerator() : TObject(), fStack(0), fevno(0), fmass(0.), fcharge(0.), fprimColl(0), frnd(0) { frnd = new TRandom3(0); // Default constructor //Default primary proton fprim.fPDG=kProton; fprim.fX0=1.; fprim.fY0=1.; fprim.fZ0=130.; fprim.fTHETA=0.; fprim.fPHI=0.; fprim.fP0=1.; //1GV } PamVMCPrimaryGenerator::~PamVMCPrimaryGenerator() { // Destructor delete frnd; delete fprimColl; } // private methods #include void PamVMCPrimaryGenerator::GeneratePrimary() { // Add one primary particle to the user stack (derived from TVirtualMCStack). // Track ID (filled by stack) Int_t ntr; // Option: to be tracked Int_t toBeDone = 1; // Particle type Int_t pdg = fprim.fPDG; Double_t fvx, fvy, fvz; fvx=fprim.fX0; fvy=fprim.fY0; fvz=fprim.fZ0; // Position Double_t tof = 0.; // Energy (in GeV) Double_t kinEnergy = MomentumToKinE(fprim.fP0); Double_t e = fmass + kinEnergy; // Particle momentum Double_t px, py, pz; px = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Cos((Pi()/180.)*(fprim.fPHI)); py = -fprim.fP0*Sin((Pi()/180.)*(fprim.fTHETA))*Sin((Pi()/180.)*(fprim.fPHI)); pz = -fprim.fP0*Cos((Pi()/180.)*(fprim.fTHETA)); // converting in radian //px = fprim.fP0*Sin(fprim.fTHETA)*Cos(fprim.fPHI); //py = fprim.fP0*Sin(fprim.fTHETA)*Sin(fprim.fPHI); //pz = -fprim.fP0*Cos(fprim.fTHETA); // Polarization TVector3 polar; // Add particle to stack fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, fvx, fvy, fvz, tof, polar.X(), polar.Y(), polar.Z(), kPPrimary, ntr, 1., 0); PamVMCPrimary * pc = (PamVMCPrimary *)fprimColl->New(fevno++); *pc = fprim; } void PamVMCPrimaryGenerator::SetParticle(Int_t pdg){ fprim.fPDG=pdg; //TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fprim.fPDG); fmass = (TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Mass(); fcharge = ((TDatabasePDG::Instance()->GetParticle(fprim.fPDG))->Charge())/3.; } void PamVMCPrimaryGenerator::SetMomentum( Double_t px, Double_t py, Double_t pz) { fprim.fP0= Sqrt(px*px+py*py+pz*pz); fprim.fTHETA=ATan(Sqrt(px*px+py*py)/pz); fprim.fPHI=ATan(py/px); } void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy) { if(isEnergy) { fprim.fP0=frnd->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax)); } else{ fprim.fP0=frnd->Uniform(PEmin,PEmax); } } void PamVMCPrimaryGenerator::GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy) { Double_t alpha = 1.+gamma; //integral spectral index if(alpha==0.){ fprim.fP0=Exp(Log(PEmin)+frnd->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin))); } else { if(PEmin==0.) PEmin=1.E-10; fprim.fP0=Power((frnd->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha); } if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0); } void PamVMCPrimaryGenerator::GenSphDist(Double_t r, Double_t Thmin, Double_t Thmax, Double_t Phmin, Double_t Phmax) { Double_t theta, phi, y, f; phi = (Pi()/180.)*frnd->Uniform(Phmin,Phmax); do { y = frnd->Uniform(0.,1.); theta = (Pi()/180.)*frnd->Uniform(Thmin,Thmax); f = Sin(theta); } while (y>f); //theta = phi = (Pi()/180.)*45.; SetPosition((r*Sin(theta)*Cos(phi)), (r*Sin(theta)*Sin(phi)), (r*Cos(theta))); //random distribution of theta phi in the angle at the vertex at (0,0,r) //by the S3 max distant corners. static const Double_t s3_x=9.0, s3_y=7.5, s3_pz=25.3159; Double_t ang = ATan((Sqrt(s3_x*s3_x+s3_y*s3_y))/(r-s3_pz)); //SetDirection((frnd->Uniform((theta-ang),(theta+ang)))/(Pi()/180.), (frnd->Uniform((phi-ang),(phi+ang)))/(Pi()/180.)); SetDirection((frnd->Uniform((theta-ang),(theta+ang)))/(Pi()/180.), (frnd->Uniform(0.,2*Pi()))/(Pi()/180.)); //SetDirection(0., 0.); } void PamVMCPrimaryGenerator::GenSphDist(Double_t r) { //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; static const Double_t s1_x=20.4, s1_y=16.5, s1_pz=102.8866; static const Double_t s2_x=9.0, s2_y=7.5, s2_pz=73.489; static const Double_t s3_x=9.0, s3_y=7.5, s3_pz=25.3159; //calculate max theta and phi angles to be allowed (calculating wrt one //corner of s3 to opposite corner of s1) Double_t rx = s1_x+s3_x; Double_t ry = s1_y+s3_y; Double_t rz = s1_pz-s3_pz; Double_t thmax = (180./Pi())*(ACos(rz/Sqrt(rx*rx+ry*ry+rz*rz))); //Double_t phmax = (180./Pi())*(ATan2(ry,rx)); //cout << "~~~~~~Theta max Phi max : " << thmax <<", "<< phmax << endl; //generate a track check and let it go only if it is passing through all //the 3 TOF Bool_t trkGood = kFALSE; do { //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; GenSphDist(r, 0., thmax, 0., 360.); Double_t th = (Pi()/180.)*fprim.fTHETA; Double_t ph = (Pi()/180.)*fprim.fPHI; //cout << "~~~~~~Theta Phi : " << fprim.fTHETA <<", "<< fprim.fPHI << endl; Double_t x1, y1, x2, y2, x3, y3; x1 = s1_pz*Tan(th)*Cos(ph) - fprim.fX0; y1 = s1_pz*Tan(th)*Sin(ph) - fprim.fY0; x2 = s2_pz*Tan(th)*Cos(ph) - fprim.fX0; y2 = s2_pz*Tan(th)*Sin(ph) - fprim.fY0; x3 = s3_pz*Tan(th)*Cos(ph) - fprim.fX0; y3 = s3_pz*Tan(th)*Sin(ph) - fprim.fY0; if ( Abs(x1) <= Abs(s1_x) && Abs(y1) <= Abs(s1_y) && Abs(x2) <= Abs(s2_x) && Abs(y2) <= Abs(s2_y) && Abs(x3) <= Abs(s3_x) && Abs(y3) <= Abs(s3_y) ) trkGood = kTRUE; //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; }while (!trkGood); //cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+" << endl; }