Creating uniform LUND files
Jump to navigation
Jump to search
#include <math.h>
#include <TRandom3.h>
#include <stdio.h>
#include <stdlib.h>
double XSect(double Mol_CM_E,double Mol_CM_Theta)
{
float angle,dSigma_dOmega, w;
double alpha=7.2973525664e-3;
//Convert MeV to eV
Mol_CM_E=53e6;
//printf("Mol_CM_E=");
//printf(" %f\n",Mol_CM_E);
//printf("Mol_CM_Theta=");
//printf(" %f\n",Mol_CM_Theta);
angle=Mol_CM_Theta;
/* using equation from Landau-Lifshitz & Azfar */
/* units in eV */
dSigma_dOmega=(alpha*(3+cos(angle)*cos(angle)));
dSigma_dOmega=dSigma_dOmega*dSigma_dOmega;
dSigma_dOmega=dSigma_dOmega/(4*Mol_CM_E*Mol_CM_E*sin(angle)*sin(angle)*sin(angle)*sin(angle));
dSigma_dOmega=dSigma_dOmega*1e18;
dSigma_dOmega=dSigma_dOmega*.3892e-3;//barns
//printf("Angle, XSect");
//printf(" %f\t%f\n", angle, dSigma_dOmega);
w=dSigma_dOmega;
//printf("weight=");
//printf(" %f\n",w);
return w;
}
void LUND_Spread()
{
TLorentzVector Fnl_Lab_4Mom,Fnl_Lab_4Mom2;
TLorentzVector Mol_Lab_4Mom,Mol_Lab_4Mom2;
TLorentzVector Fnl_CM_4Mom;
TLorentzVector Mol_CM_4Mom;
TLorentzVector CMS;
double Fnl_CM[3],Fnl_CM_P,Fnl_CM_E;
double Mol_CM[3],Mol_CM_P,Mol_CM_E;
double Fnl_Lab[3],Fnl_Lab_P,Fnl_Lab_E;
double Mol_Lab[3],Mol_Lab_P,Mol_Lab_E;
double Total_Lab_E, Total_CM_E, Total_Lab_P, Total_CM_P, Total_Lab_Mass;
double Total_Lab[3];
double Fnl_Theta_CM, Fnl_Theta_Lab, Fnl_Phi_CM, Fnl_Phi_Lab;
double Mol_Theta_CM, Mol_Theta_Lab, Mol_Phi_CM, Mol_Phi_Lab;
double Init_e_Lab_E, Init_e_Lab_P;
double Nu,Qsqrd,x_bj,y,W;
double weight;
TFile *f = new TFile("LUND_spread.root","RECREATE");
TTree *tree = new TTree("Moller","Moller data from ascii file");
FILE *f0;
f0=fopen("LUND_spread.LUND","w");
TH1F* FnlMomLab2=new TH1F("FnlMomLab2","Final Scattered Electron Momentum Lab Frame 2",3000,0,15000);
TH1F* MolMomLab2=new TH1F("MolMomLab2","Moller Momentum Lab Frame 2", 30000,0,60);
TH1F* FnlMomCM2=new TH1F("FnlMomCM2","Final Scattered Electron Momentum CM Frame 2",1000,40,60);
TH1F* MolMomCM2=new TH1F("MolMomCM2","Moller Momentum CM Frame 2",1000,40,60);
TH1F* FnlThetaLab2=new TH1F("FinalThetaLab2","Final Scattered Electron Theta Lab Frame 2",200,0,1);
TH1F* MolThetaLab2=new TH1F("MolThetaLab2","Moller Theta Lab Frame 2",2000,0,60);
TH1F* FnlThetaCM2=new TH1F("FnlThetaCM2","Final Scattered Electron Theta CM Frame 2",360,0,90);
TH1F* MolThetaCM2=new TH1F("MolThetaCM2","Moller Theta CM Frame 2",360,90,180);
TH1F* FnlMomRangeLab2=new TH1F("FnlMomRangeLab2","Final Scattered Electron Momentum Range Lab Frame 2",3000,0,15000);
TH1F* MolMomRangeLab2=new TH1F("MolMomRangeLab2","Moller Momentum Lab Range Frame 2", 30000,0,60);
TH1F* FnlMomRangeCM2=new TH1F("FnlMomRangeCM2","Final Scattered Electron Momentum Range CM Frame 2",1000,40,60);
TH1F* MolMomRangeCM2=new TH1F("MolMomRangeCM2","Moller Momentum Range CM Frame 2",1000,40,60);
TH1F* FnlThetaRangeLab2=new TH1F("FnlThetaRangeLab2","Final Scattered Electron Theta Range Lab Frame 2",200,0,1);
TH1F* MolThetaRangeLa2b=new TH1F("MolThetaRangeLab2","Moller Theta Range Lab Frame 2",2000,0,60);
TH1F* FnlThetaRangeCM2=new TH1F("FnlThetaRangeCM2","Final Scattered Electron Theta Range CM Frame 2",360,0,90);
TH1F* MolThetaRangeCM2=new TH1F("MolThetaRangeCM2","Moller Theta Range CM Frame 2",360,90,180);
// tree->Branch("evt",&evt.event,"event/I:FnlPx/F:FnlPy:FnlPz:MolPx:MolPy:MolPz");
//Loop over Moller Energy 2MeV-5500MeV
Mol_Lab_E=2; //MeV
for(int e=1;e<550;e++) //Energy 2-5500 MeV by 10 MeV increments
{
Mol_Lab_E=Mol_Lab_E+10;
Mol_Theta_Lab=4.5; //5 degrees=.0872664626 radians
for(int t=0;t<70;t++) //Theta 5-40 degrees, by .5 degree increments
{
//printf("Mol_Lab_E=");
//printf(" %f\n",Mol_Lab_E);
Mol_Theta_Lab=Mol_Theta_Lab+0.5;
//printf("Mol_Theta_Lab=");
//printf(" %f\n",Mol_Theta_Lab);
//Convert Moller Lab Theta to radians
Mol_Theta_Lab=Mol_Theta_Lab*3.14;
Mol_Theta_Lab=Mol_Theta_Lab/180;
//Moller Lab Total Momentum
Mol_Lab_P=sqrt(Mol_Lab_E*Mol_Lab_E-.511*.511);
//printf("Mol_Lab_P=");
//printf(" %f\n",Mol_Lab_P);
//Moller Lab Pz
Mol_Lab[2]=cos(Mol_Theta_Lab);
Mol_Lab[2]=Mol_Lab_P*Mol_Lab[2];
//printf("Mol_Lab[2]=");
//printf(" %f\n",Mol_Lab[2]);
//Assign a Phi Angle (10 degrees)
float random_phi=0.174532925199;
//printf("Random Phi=");
//printf(" %f\n",random_phi);
//Moller Lab Px
Mol_Lab[0]=cos(random_phi);
Mol_Lab[0]=Mol_Lab[0]*sqrt((Mol_Lab_P*Mol_Lab_P)-(Mol_Lab[2]*Mol_Lab[2]));
if(random_phi>=1.57079632679 && random_phi<=3.14159265359)
Mol_Lab[0]=-Mol_Lab[0];
if(random_phi>=-1.57079632679 && random_phi<=-3.14159265359)
Mol_Lab[0]=-Mol_Lab[0];
//printf("Mol_Lab[0]=");
//printf( " %f\n",Mol_Lab[0]);
//Moller Lab Py
Mol_Lab[1]=sqrt(Mol_Lab_P*Mol_Lab_P-Mol_Lab[2]*Mol_Lab[2]-Mol_Lab[0]*Mol_Lab[0]);
if((random_phi<=0) && (random_phi>=-1.57079632679))
Mol_Lab[1]=-Mol_Lab[1];
if(random_phi>=-1.57079632679 && random_phi<=-3.14159265359)
Mol_Lab[1]=-Mol_Lab[1];
//printf("Mol_Lab[1]=");
//printf( " %f\n",Mol_Lab[1]);
//printf( "\n");
//Define Moller Lab 4 vector
Mol_Lab_4Mom.SetPxPyPzE(Mol_Lab[0],Mol_Lab[1],Mol_Lab[2],Mol_Lab_E);
//printf("Mol_Lab_4Mom.E=");
//printf(" %f\n",Mol_Lab_4Mom.E());
//printf("Mol_Lab_4Mom.P=");
//printf(" %f\n",Mol_Lab_4Mom.P());
//printf("Mol_Lab_4Mom.Px=");
//printf(" %f\n",Mol_Lab_4Mom.Px());
//printf("Mol_Lab_4Mom.Py=");
//printf(" %f\n",Mol_Lab_4Mom.Py());
//printf("Mol_Lab_4Mom.Pz=");
//printf(" %f\n",Mol_Lab_4Mom.Pz());
//printf("Mol_Lab_4Mom.Theta=");
//printf(" %f\n",Mol_Lab_4Mom.Theta());
//printf(" \n");
//Find 4 Momentum vector componets for Lab frame electron
//Fnl_Lab_E
Fnl_Lab_E=11000-Mol_Lab_4Mom.E();
//printf("Fnl_Lab_E=");
//printf( " %f\n",Fnl_Lab_E);
//P_Lab(total)
//Fnl_Lab_P=sqrt((Fnl_Lab_E*Fnl_Lab_E)-(.511*.511));
//printf("Fnl_Lab_P=");
//printf( " %f\n",Fnl_Lab_P);
//Px
//Fnl_Lab[0]=-Mol_Lab[0];
//printf("Fnl_Lab[0]=");
//printf(" %f\n",Fnl_Lab[0]);
//Py
//Fnl_Lab[1]=-Mol_Lab_4Mom.Py();
//printf("Fnl_Lab[1]=");
//printf(" %f\n",Fnl_Lab[1]);
//Pz
Fnl_Lab[2]=((Fnl_Lab_E*Fnl_Lab_E)-(Mol_Lab_4Mom.Px()*Mol_Lab_4Mom.Px())-(Mol_Lab_4Mom.Py()*Mol_Lab_4Mom.Py()));
Fnl_Lab[2]=sqrt(Fnl_Lab[2]);
//printf("Fnl_Lab[2]=");
//printf(" %f\n",Fnl_Lab[2]);
//Define 4vector
Fnl_Lab_4Mom.SetPxPyPzE(-Mol_Lab[0],-Mol_Lab[1],Fnl_Lab[2],Fnl_Lab_E);
//printf("Fnl_Lab_4Mom.Px=");
//printf( " %f\n",Fnl_Lab_4Mom.Px());
//printf("Fnl_Lab_4Mom.Py=");
//printf( " %f\n",Fnl_Lab_4Mom.Py());
//printf("Fnl_Lab_4Mom.E=");
//printf( " %f\n",Fnl_Lab_4Mom.E());
//printf("Fnl_Lab_4Mom.P=");
//printf( " %f\n",Fnl_Lab_4Mom.P());
//printf("Fnl_Lab_4Mom.Pz=");
//printf( " %f\n",Fnl_Lab_4Mom.Pz());
//printf("Fnl_Lab_4Mom.Theta=");
//printf( " %f\n",Fnl_Lab_4Mom.Theta());
//printf( "\n");
//Boost vectors
CMS=Fnl_Lab_4Mom+Mol_Lab_4Mom;
//printf("CMS_P=");
//printf(" %f\n",CMS.P());
//printf("CMS_E=");
//printf(" %f\n",CMS.E());
//printf("CMS.Px=");
//printf( " %f\n",CMS.Px());
//printf("CMS.Py=");
//printf( " %f\n",CMS.Py());
//printf("CMS.Pz=");
//printf( " %f\n",CMS.Pz());
//printf("\n");
//Boost to CM Frame
Fnl_Lab_4Mom2=Fnl_Lab_4Mom;
Mol_Lab_4Mom2=Mol_Lab_4Mom;
Fnl_Lab_4Mom.Boost(-CMS.BoostVector());
Fnl_CM_4Mom=Fnl_Lab_4Mom;
Mol_Lab_4Mom.Boost(-CMS.BoostVector());
Mol_CM_4Mom=Mol_Lab_4Mom;
//printf("Mol_CM_P=");
//printf( "%f\n",Mol_CM_4Mom.P());
//printf("Mol_CM_E=");
//printf( "%f\n",Mol_CM_4Mom.E());
//printf("\n");
//Find 4 Momentum vector components for CM Moller electron
//Mol_CM_E
Mol_CM_E=Mol_Lab_4Mom.E();
//printf("Mol_CM_E=");
//printf(" %f\n",Mol_CM_4Mom.E());
//Mol_CM_P
Mol_CM_P=sqrt((Mol_CM_E*Mol_CM_E)-(.511*.511));
//printf("Mol_CM_P=");
//printf( " %f\n",Mol_CM_P);
//Px
Mol_CM[0]=Mol_Lab[0];
//printf("Mol_CM[0]=");
//printf( " %f\n",Mol_CM[0]);
//Py
Mol_CM[1]=Mol_Lab[1];
//printf("Mol_CM[1]=");
//printf( " %f\n",Mol_CM[1]);
//Pz
Mol_CM[2]=sqrt(Mol_CM_P*Mol_CM_P-Mol_CM[0]*Mol_CM[0]-Mol_CM[1]*Mol_CM[1]);
//Mol_CM[2]=-Mol_CM[2];
//printf("Mol_CM[2]=");
//printf( " %f\n",Mol_CM[2]);
//Theta
Mol_Theta_CM=acos(Mol_CM[2]/Mol_CM_P);
//printf("Mol_Theta_CM=");
//printf( " %f\n",Mol_Theta_CM);
//printf( "\n");
//Define 4vector
//Mol_CM_4Mom.SetPxPyPzE(Mol_CM[0],Mol_CM[1],Mol_CM[2],Mol_CM_E);
//printf("Mol_CM_4Mom.E=");
//printf(" %f\n",Mol_CM_4Mom.E());
//printf("Mol_CM_4Mom.P=");
//printf(" %f\n",Mol_CM_4Mom.P());
//printf("Mol_CM_4Mom.Px=");
//printf(" %f\n",Mol_CM_4Mom.Px());
//printf("Mol_CM_4Mom.Py=");
//printf(" %f\n",Mol_CM_4Mom.Py());
//printf("Mol_CM_4Mom.Pz=");
//printf(" %f\n",Mol_CM_4Mom.Pz());
//printf("Mol_CM_4Mom.Theta=");
//printf(" %f\n",Mol_CM_4Mom.Theta());
//printf(" \n");
//printf("nlines=");
//printf(" %i\n",nlines);
//Find 4 Momentum vector components for CM electron
//Fnl_CM_E
Fnl_CM_E=Mol_CM_E; //Mev
//printf("Fnl_CM_E=");
//printf( " %f\n",Fnl_CM_E);
//Fnl_CM_P
Fnl_CM_P=Mol_CM_P;
//printf("Fnl_CM_P=");
//printf( " %f\n",Fnl_CM_P);
//Px
Fnl_CM[0]=-Mol_CM[0];
//printf("Fnl_CM[0]=");
//printf( " %f\n",Fnl_CM[0]);
//Py
Fnl_CM[1]=-Mol_CM[1];
//printf("Fnl_CM[1]=");
//printf( " %f\n",Fnl_CM[1]);
//Pz
Fnl_CM[2]=Mol_CM[2];
//printf("Fnl_CM[2]=");
//printf( " %f\n",Fnl_CM[2]);
//Fnl_Theta_CM
Fnl_Theta_CM=-Mol_Theta_CM+3.14159265359;
//printf("Fnl_Theta_CM=");
//printf( " %f\n",Fnl_Theta_CM);
//printf( "\n");
//Define 4vector
//Fnl_CM_4Mom.SetPxPyPzE(Fnl_CM[0],Fnl_CM[1],Fnl_CM[2],Fnl_CM_E);
//printf("Fnl_CM_4Mom.E=");
//printf( " %f\n",Fnl_CM_4Mom.E());
//printf("Fnl_CM_4Mom.P=");
//printf( " %f\n",Fnl_CM_4Mom.P());
//printf("Fnl_CM_4Mom.Px=");
//printf( " %f\n",Fnl_CM_4Mom.Px());
//printf("Fnl_CM_4Mom.Py=");
//printf( " %f\n",Fnl_CM_4Mom.Py());
//printf("Fnl_CM_4Mom.Pz=");
//printf( " %f\n",Fnl_CM_4Mom.Pz());
//printf("Fnl_CM_4Mom.Theta=");
//printf( " %f\n",Fnl_CM_4Mom.Theta());
//printf( "\n");
//Create Lab Frame Histograms
//FnlMomLab2->Fill(Fnl_Lab_4Mom.P());
//MolMomLab2->Fill(Mol_Lab_4Mom.P());
//FnlThetaLab2->Fill(Fnl_Lab_4Mom.Theta()*180/3.14159265359);
//MolThetaLab2->Fill(Mol_Lab_4Mom.Theta()*180/3.14159265359);
//Create CM Histograms
//FnlMomCM2->Fill(Fnl_CM_4Mom.P());
//MolMomCM2->Fill(Mol_CM_4Mom.P());
//FnlThetaCM2->Fill(Fnl_CM_4Mom.Theta()*180/3.14159265359);
//MolThetaCM2->Fill(Mol_CM_4Mom.Theta()*180/3.14159265359);
Nu=Mol_Lab_4Mom.E()-Fnl_Lab_4Mom.E();
Qsqrd=4*Mol_Lab_4Mom.E()*Fnl_Lab_4Mom.E()*(1-Fnl_Lab_4Mom.Pz()/Fnl_Lab_4Mom.E())/2; /* should be final momentum and not final energy*/
x_bj=Qsqrd/(2*0.938*Nu);
y=Nu/Mol_Lab_4Mom.E();
W=0.938*0.938+2*0.938*Nu-Qsqrd;
if(W>0)
W=sqrt(W);
//
// MeV MUST be converted to GeV for LUND format! Use your MeV values, but limit them to fake 1-10GeV range (.1-10MeV ) in name only.
double px,py,pz;
double Px,Py,Pz;
double KE,ke;
weight=XSect(Mol_CM_4Mom.E(),Mol_CM_4Mom.Theta());
x_bj=weight;
Px=Fnl_Lab_4Mom2.Px()/1000;
Py=Fnl_Lab_4Mom2.Py()/1000;
Pz=Fnl_Lab_4Mom2.Pz()/1000;
px=Mol_Lab_4Mom2.Px()/1000;
py=Mol_Lab_4Mom2.Py()/1000;
pz=Mol_Lab_4Mom2.Pz()/1000;
KE=Fnl_Lab_4Mom2.E()/1000;
ke=Mol_Lab_4Mom2.E()/1000;
fprintf(f0,"%d\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n",2,1.,1.,1.,1.,x_bj,y,W,Qsqrd,Nu);
fprintf(f0, "%d\t%g\t%g\t%g\t%g\t%g\t%f\t%f\t%f\t%f\t%g\t%g\t%g\t%g\n",
1,-1.,1.,11.,0.,0.,Px,Py,Pz,KE,0.000511, 0.,0.,0.);
fprintf(f0, "%d\t%g\t%g\t%g\t%g\t%g\t%f\t%f\t%f\t%f\t%g\t%g\t%g\t%g\n",
2,-1.,1.,11.,0.,0.,px,py,pz,ke,0.000511, 0.,0.,0.);
//Convert Moller Lab Theta to degrees
Mol_Theta_Lab=Mol_Theta_Lab/3.14;
Mol_Theta_Lab=Mol_Theta_Lab*180;
}//End Theta loop
}//End Energy loop
tree->Print();
tree->Write();
f->Write();
delete tree;
delete f;
}//End main