// $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03 #include #include #include #include #include #include #include #include #include "PamVMCPrimaryGenerator.h" using namespace TMath; 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.), frandom(0) { // Standard constructor ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.); ftheta->SetNpx(1000); fprimColl = new TClonesArray("PamVMCPrimary"); 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), frandom(0) { // Default constructor //Default primary proton ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.); ftheta->SetNpx(1000); 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 ftheta; delete fprimColl; } // private methods 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) //printf("generateprimary check fprimP0 = %f\n",fprim.fP0); Double_t kinEnergy = MomentumToKinE(fprim.fP0); Double_t e = fmass + kinEnergy; // Particle momentum Double_t px, py, pz; 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=frandom->Uniform(KinEToMomentum(PEmin),KinEToMomentum(PEmax)); } else{ fprim.fP0=frandom->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)+frandom->Uniform(0.,1.)*(Log(PEmax)-Log(PEmin))); } else { if(PEmin==0.) PEmin=1.E-10; fprim.fP0=Power((frandom->Uniform(0.,1.)*(Power(PEmax,alpha)-Power(PEmin,alpha))+Power(PEmin,alpha)),1./alpha); } cout<<"GenSpe fprim.fP0= "<Uniform(PEmin,PEmax); wurfy = frandom->Uniform(funct_min,funct_max); if( wurfy<(function3par(wurfP,a,b,c) )) { // this is ok! fprim.fP0=wurfP; found=1; } } //printf("exit+++++++++++++++++++ %f %f \n",wurfP,fprim.fP0); } // cecilia pizzolotto void PamVMCPrimaryGenerator::GenSpe_Flat(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy) { // Generates a flat spectrum from PEmin to PElim. Then a power law Double_t PElim = 1.; //Double_t alpha = 1.+gamma; //integral spectral index Bool_t okflag=0.; Double_t throw_x =0.; Double_t throw_y =0.; while(okflag==0) { throw_x=frandom->Uniform(PEmin,PEmax); // cout<<" x "<Uniform(0.,1.); if( throw_y<(1*pow(throw_x,gamma))) { okflag=1.; } } } fprim.fP0=throw_x; //h->Fill(fprimf.P0); okflag=0.; // reset if(isEnergy) fprim.fP0=KinEToMomentum(fprim.fP0); } // Spherical distribution -- Test by Cecilia P july 2009 ---- // flusso isotropo su 2pi void PamVMCPrimaryGenerator::GenSphericalPhiThe() { // Generate phi theta Double_t theta=0.; Double_t phi=0.; Double_t xcos = sqrt( frandom->Uniform(0.,1.) ); theta = acos(xcos); //RAD phi = frandom->Uniform(0.,2.*Pi()); SetDirection(theta, phi); return; } void PamVMCPrimaryGenerator::GenSphPhiThe(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Double_t zmin, Double_t zmax) { Bool_t trkGood = kFALSE; Double_t theta = 999.; Double_t phi = 0.; Double_t x2,y2,x3,y3; Double_t x0,y0,z0; //static const Double_t rad2deg = 57.2958; // S21 and S31 position/size taken as reference (z on top of det) // constraint: must pass throuth these planes static const Double_t s2_xmax=9.05, s2_ymax=7.55, s2_z=73.439; // z on top of det static const Double_t s3_xmax=9.05, s3_ymax=7.55, s3_z=26.093; // z on top of det //Double_t thetamax=3.14; //thetamax = atan((xmax+s3_xmax)/(zmax-s3_z)); //cout<<" Quanto รจ il theta max? "<Uniform(xmin,xmax); y0= frandom->Uniform(ymin,ymax); z0= frandom->Uniform(zmin,zmax); // Generate phi theta theta=999.; // init while (theta>=0.65) // take only theta smaller than 37deg=0.65rad { Double_t xcos = sqrt( frandom->Uniform(0.,1.) ); theta = acos(xcos); //RAD } phi = frandom->Uniform(0.,2.*Pi()); // Calculate xy at the constraint Double_t fact2 = (s2_z-z0)/cos(theta); x2 = x0 + fabs(fact2) * sin(theta) * cos(phi); y2 = y0 + fabs(fact2) * sin(theta) * sin(phi); Double_t fact3 = (s3_z-z0)/cos(theta); x3 = x0 + fabs(fact3) * sin(theta) * cos(phi); y3 = y0 + fabs(fact3) * sin(theta) * sin(phi); //cout<<" x/y0= "<