MollerPrimaryGenAction_phasespace_motts.cc

Go to the documentation of this file.
00001 #include "G4Event.hh"
00002 #include "G4ParticleGun.hh"
00003 //#include "G4ParticleTable.hh"
00004 //#include "G4ParticleDefinition.hh"
00005 #include "G4Electron.hh"
00006 #include "globals.hh"
00007 //#include "Randomize.hh"
00008 
00009 #include "MollerPrimaryGenAction.hh"
00010 #include "RootAnalysis.hh"
00011 
00012 
00013 void MollerPrimaryGenAction::GeneratePrimaries_phasespace_motts(G4Event* anEvent)
00014 {
00015   //This generator throws events in phase space
00016 
00017   G4int numsectors = 7; 
00018   G4double sectorhalfopeningangle = pi/numsectors/2.;
00019   G4double etheta=0., ephi=0., zpos=0.;
00020 
00021   G4int numthetas=17, numphis=13, numzpos=2, numphicentre=7;
00022   G4double thetamin=0.004, thetamax=0.020;
00023   G4double zposmin=-0.75, zposmax=0.75;
00024   G4double phicentre=0;
00025   G4double phimin=0, phimax=2*pi;
00026   G4double phicentremin=0, phicentremax=12./7*pi;
00027 
00028   int thetanum, ephinum, zposnum, phicentrenum;
00029   double thetafrac, ephifrac, zposfrac, phicentrefrac;
00030   // step through theta then phi then target position z then sector   
00031   thetanum = eventnumber%(numthetas);  // % gives remainder after division
00032   ephinum = (eventnumber/(numthetas))%(numphis);
00033   zposnum = (eventnumber/(numthetas*numphis))%(numzpos);
00034   phicentrenum = (eventnumber/(numthetas*numphis*numzpos))%(numphicentre);
00035 
00036   thetafrac = ((double)thetanum)/(numthetas-1);
00037   ephifrac = ((double)ephinum)/(numphis-1);
00038   zposfrac = ((double)zposnum)/(numzpos-1);
00039   phicentrefrac = ((double)phicentrenum)/(numphicentre-1);
00040 
00041   etheta = thetamin + thetafrac*(thetamax-thetamin);
00042   phicentre = phicentremin + phicentrefrac*(phicentremax-phicentremin);
00043   phimin = phicentre - sectorhalfopeningangle*0.999;
00044   phimax = phicentre + sectorhalfopeningangle*0.999;
00045   ephi = phimin + ephifrac*(phimax-phimin);
00046   zpos = zposmin + zposfrac*(zposmax-zposmin);
00047 
00048   G4double eenergy;
00049     eenergy = ebeam - ebeam/(1 + mp/(2*ebeam*sin(etheta/2.)*sin(etheta/2.))); // ep elastic electrons
00050 
00051   G4double xdir = sin(etheta)*cos(ephi);
00052   G4double ydir = sin(etheta)*sin(ephi);
00053   G4double zdir = cos(etheta);
00054   printf("Primaries:%5i phi:%2i %.2f zpos:%2i %.2f; phicentre:%2i %.2f; theta:%2i %.2f =%.3f deg =%4.1f mrad  eenergy=%.3f GeV  phi=%6.1f deg, xdir=%.3f, ydir=%.3f, zdir=%.3f\n", 
00055        eventnumber, ephinum, ephifrac, zposnum, zposfrac, phicentrenum, phicentrefrac, 
00056        thetanum, thetafrac, etheta/degree, etheta/rad*1000, eenergy, ephi/degree, xdir, ydir, zdir);
00057   
00058   particleGun->SetParticleEnergy(eenergy*GeV);
00059   particleGun->SetParticleMomentumDirection(G4ThreeVector(xdir,ydir,zdir));
00060   particleGun->SetParticlePosition(G4ThreeVector(0.*cm, 0.*cm, zpos*m));
00061   particleGun->GeneratePrimaryVertex(anEvent);
00062   eventnumber++;
00063 
00064   gRootAnalysis->SetMomentum0(11000,0,0,11000);
00065   gRootAnalysis->SetMomentum2(eenergy*MeV,eenergy*MeV*xdir,eenergy*MeV*ydir,eenergy*MeV*zdir);
00066   gRootAnalysis->SetMomentum1(eenergy*MeV,eenergy*MeV*xdir,eenergy*MeV*ydir,eenergy*MeV*zdir);
00067   gRootAnalysis->SetProcess(1);
00068 
00069   if (eventnumber==(numphis*numthetas*numzpos*7)) {
00070     printf("\nFully completed one set of optics rays.\n\n");
00071   }
00072 }
00073 
00074 
00075 
00076 /* emacs
00077  * Local Variables:
00078  * mode:C++
00079  * mode:font-lock
00080  * c-file-style: "stroustrup"
00081  * tab-width: 4
00082  * End:
00083  */