--- PamVMC/src/PamVMCPrimaryGenerator.cxx 2007/06/28 07:16:57 1.1 +++ PamVMC/src/PamVMCPrimaryGenerator.cxx 2010/09/15 07:01:57 1.6 @@ -1,5 +1,5 @@ // $Id: PamVMCPrimaryGenerator.cxx,v 1.0 2006/06/03 -// + #include #include @@ -8,48 +8,87 @@ #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), - fPdg(kProton), - fKinEnergy(100.0),//100 GeV - fDirX(0.), - fDirY(0.), - fDirZ(-1.), - fPolAngle(0.), - fNofPrimaries(1) + 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), - fPdg(0), - fKinEnergy(0.), - fDirX(0.), - fDirY(0.), - fDirZ(0.), - fPolAngle(0.), - fNofPrimaries(0) + fevno(0), + fmass(0.), + fcharge(0.), + fprimColl(0), + frandom(0) { -// Default constructor + // 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 +// Destructor + delete ftheta; + delete fprimColl; } // private methods -#include void PamVMCPrimaryGenerator::GeneratePrimary() { @@ -62,71 +101,231 @@ Int_t toBeDone = 1; // Particle type - Int_t pdg = fPdg; - TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(fPdg); - + Int_t pdg = fprim.fPDG; + + Double_t fvx, fvy, fvz; + fvx=fprim.fX0; + fvy=fprim.fY0; + fvz=fprim.fZ0; + // Position - Double_t vx = 0.; - Double_t vy = 0.; - Double_t vz = 110.; + Double_t tof = 0.; // Energy (in GeV) - Double_t kinEnergy = fKinEnergy; - Double_t mass = particlePDG->Mass(); - Double_t e = mass + kinEnergy; + //printf("generateprimary check fprimP0 = %f\n",fprim.fP0); + Double_t kinEnergy = MomentumToKinE(fprim.fP0); + Double_t e = fmass + kinEnergy; // Particle momentum - Double_t p0, px, py, pz; - p0 = sqrt(e*e - mass*mass); - px = p0 * fDirX; - py = p0 * fDirY; - pz = p0 * fDirZ; + 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; -// if ( fPdg == 50000050 ) { -// TVector3 normal (1., 0., 0.); -// TVector3 kphoton = TVector3(fDirX, fDirY, fDirZ); -// TVector3 product = normal.Cross(kphoton); -// Double_t modul2 = product*product; - -// TVector3 e_perpend (0., 0., 1.); -// if (modul2 > 0.) e_perpend = (1./sqrt(modul2))*product; -// TVector3 e_paralle = e_perpend.Cross(kphoton); - -/// polar = TMath::Cos(fPolAngle*TMath::DegToRad())*e_paralle - // + TMath::Sin(fPolAngle*TMath::DegToRad())*e_perpend; - // } - //else - // Warning("GeneratePrimary", - // "The primary particle is not an opticalphoton"); // Add particle to stack - fStack->PushTrack(toBeDone, -1, pdg, px, py, pz, e, vx, vy, vz, tof, + 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; } -// public methods -void PamVMCPrimaryGenerator::GeneratePrimaries() -{ -// Fill the user stack (derived from TVirtualMCStack) with primary particle +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); + } - for (Int_t i=0; iUniform(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()); - fDirX = dirX/norm; - fDirY = dirY/norm; - fDirZ = dirZ/norm; -} + 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= "<