Difference between revisions of "LUND Spread LH2 0Phi.C"
Jump to navigation
Jump to search
(Created page with "#include <math.h> #include <TRandom3.h> #include <stdio.h> #include <stdlib.h> double XSect(double Mol_CM_E,double Mol_CM_Theta) { double angle,dSigma_dOmega, w; double alph...") |
|||
| Line 1: | Line 1: | ||
| + | <pre> | ||
#include <math.h> | #include <math.h> | ||
#include <TRandom3.h> | #include <TRandom3.h> | ||
| Line 296: | Line 297: | ||
}//End main | }//End main | ||
| + | </pre> | ||
Latest revision as of 22:27, 2 November 2017
#include <math.h>
#include <TRandom3.h>
#include <stdio.h>
#include <stdlib.h>
double XSect(double Mol_CM_E,double Mol_CM_Theta)
{
double angle,dSigma_dOmega, w;
double alpha=7.2973525664e-3;
//Convert MeV to eV
Mol_CM_E=Mol_CM_E*1e6;
//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
dSigma_dOmega=dSigma_dOmega/2; //Adjust for 2 electrons
//printf("Angle, XSect");
//printf(" %f\t%f\n", angle*180/3.14, dSigma_dOmega);
w=dSigma_dOmega;
//printf("weight=");
//printf(" %f\n",w);
return w;
}
void LUND_Spread_LH2_0Phi()
{
TLorentzVector Fnl_Lab_4Mom,Mol_Lab_4Mom;
TLorentzVector Fnl_CM_4Mom,Mol_CM_4Mom;
TLorentzVector Init_e_Lab_4Mom, Init_Mol_Lab_4Mom;
TLorentzVector Init_e_CM_4Mom, Init_Mol_CM_4Mom;
TLorentzVector Lab2CMS,CMS2Lab;
double theta_wire[113];
double Fnl_CM[3],Fnl_CM_P,Fnl_CM_E;
double Fnl_Mol_CM[3],Fnl_Mol_CM_P,Fnl_Mol_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_Lab_P;
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,s;
double Theta;
double random_phi;
double Energy;
double weight;
double wire_number;
double t,a1,a2,B1Right,B2Right,B3Right,B4Right,B5Right,B6Right,B1Left,B2Left,B3Left,B4Left,B5Left,B6Left;
double wireNumRight,wireNumLeft;
int trigger=0;
double low_limit=5.000;
double high_limit=5.001;
TFile *f = new TFile("LUND_Spread_LH2_0Phi.root","RECREATE");
TTree *tree = new TTree("Moller","Moller data from ascii file");
FILE *f0;
f0=fopen("LUND_Spread_LH2_0Phi.LUND","w");
TH1F* FullMolThetaCM=new TH1F("FullMolThetaCM","Moller Theta CM Frame",900001,90,180);
TH1F* FullMolThetaLab=new TH1F("FullMolThetaLab","Moller Theta Lab Frame",901,0,90);
TH1F* FullMolThetaCMweighted=new TH1F("FullMolThetaCMweighted","Moller Theta CM Frame Weighted",900001,90,180);
TH1F* FullMolThetaLabweighted=new TH1F("FullMolThetaLabweighted","Moller Theta Lab Frame Weighted",901,0,90);
TH1D* PhiDistribution=new TH1D("PhiDistribution","Random Phi Distribution",601,-30,30);
TH1F* MolThetaCM=new TH1F("MolThetaCM","Moller Theta CM Frame",900001,90,180);
TH1F* MolThetaCMweighted=new TH1F("MolThetaCMweighted","Moller Theta CM Frame Weighted",900001,90,180);
TH1F* MolThetaCMweightedReBin=new TH1F("MolThetaCMweightedReBin","Moller Theta CM Frame Weighted Re-binned",901,90,180);
TH1F* MolThetaLab=new TH1F("MolThetaLab","Moller Theta Lab Frame",501,0,50);
TH1F* MolThetaLabweighted=new TH1F("MolThetaLabweighted","Moller Theta Lab Frame Weighted",901,0,90);
//Loop over Moller Theta in CM (90-180)
Mol_Theta_CM=89.99999;
for(int theta=0;theta<8999999;theta++) //Theta 90-179.9999 degrees, by .0001 degree increment
{
random_phi=0.0;
Mol_Theta_CM=Mol_Theta_CM+.00001;
random_phi=random_phi*3.14159265359;
random_phi=random_phi/180;
//
//Define 4 Momentum vector componets for Initial Lab frame electrons(11GeV electron, stationary electron)
//
Init_e_Lab_P=sqrt(11000*11000-0.510998910*0.510998910);
Init_e_Lab_4Mom.SetPxPyPzE(0,0,Init_e_Lab_P,11000);
Init_Mol_Lab_4Mom.SetPxPyPzE(0,0,0,0.510998910);
Lab2CMS=Init_e_Lab_4Mom+Init_Mol_Lab_4Mom;
Total_Lab_E=Lab2CMS.E();
Total_Lab_P=Lab2CMS.P();
//Define the Mandelstram variable s
s=sqrt(Total_Lab_E*Total_Lab_E-Total_Lab_P*Total_Lab_P);
//Boost to CM before collision
Init_e_CM_4Mom=Init_e_Lab_4Mom;
Init_Mol_CM_4Mom=Init_Mol_Lab_4Mom;
Init_Mol_CM_4Mom=Init_Mol_Lab_4Mom;
Init_e_CM_4Mom.Boost(-Lab2CMS.BoostVector());
Init_Mol_CM_4Mom.Boost(-Lab2CMS.BoostVector());
//
//Cycle through Theta and E for Final CM Frame of Moller Electron
//
//Convert Moller CM Theta to radians
Mol_Theta_CM=Mol_Theta_CM*3.14159265359;
Mol_Theta_CM=Mol_Theta_CM/180;
// Moller CM Energy
Mol_CM_E=s/2;
//Moller CM Total Momentum
Mol_CM_P=sqrt(Mol_CM_E*Mol_CM_E-0.510998910*0.510998910);
//Moller CM Pz
Mol_CM[2]=cos(Mol_Theta_CM);
Mol_CM[2]=Mol_CM_P*Mol_CM[2];
//Moller CM Px
Mol_CM[0]=cos(random_phi);
Mol_CM[0]=Mol_CM[0]*sqrt((Mol_CM_P*Mol_CM_P)-(Mol_CM[2]*Mol_CM[2]));
//Adjust sign for Px based on quadrant position
if(random_phi>=1.57079632679 && random_phi<=3.14159265359 && Mol_CM[0]>0)
Mol_CM[0]=-Mol_CM[0];
if(random_phi>=-1.57079632679 && random_phi<=-3.14159265359 && Mol_CM[0]>0)
Mol_CM[0]=-Mol_CM[0];
//Moller CM Py
Mol_CM[1]=sqrt(Mol_CM_P*Mol_CM_P-Mol_CM[2]*Mol_CM[2]-Mol_CM[0]*Mol_CM[0]);
//Adjust sign for Py based on quadrant position
if((random_phi<=0) && (random_phi>=-3.14159265359) && Mol_CM[1]>0)
Mol_CM[1]=-Mol_CM[1];
//Define Moller Lab Final 4 vector
Mol_CM_4Mom.SetPxPyPzE(Mol_CM[0],Mol_CM[1],Mol_CM[2],Mol_CM_E);
Fnl_CM_4Mom.SetPxPyPzE(-Mol_CM[0],-Mol_CM[1],-Mol_CM_4Mom.Pz(),Mol_CM_4Mom.E());
//Define CM total frame
CMS2Lab=Mol_CM_4Mom+Fnl_CM_4Mom;
//Boost to Lab Frame
Mol_Lab_4Mom=Mol_CM_4Mom;
Fnl_Lab_4Mom=Fnl_CM_4Mom;
Mol_Lab_4Mom.Boost(Lab2CMS.BoostVector());
Fnl_Lab_4Mom.Boost(Lab2CMS.BoostVector());
//
//Define the LUND components
//
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!
//
//
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_4Mom.Px()/1000;
Py=Fnl_Lab_4Mom.Py()/1000;
Pz=Fnl_Lab_4Mom.Pz()/1000;
px=Mol_Lab_4Mom.Px()/1000;
py=Mol_Lab_4Mom.Py()/1000;
pz=Mol_Lab_4Mom.Pz()/1000;
KE=Fnl_Lab_4Mom.E()/1000;
ke=Mol_Lab_4Mom.E()/1000;
Theta=Mol_Lab_4Mom.Theta()/3.14159265359*180;
Energy=Mol_Lab_4Mom.E()/1000; //GeV
FullMolThetaCM->Fill(Mol_CM_4Mom.Theta()*180/3.14159265359); //account for 2 electrons of LH2
FullMolThetaLab->Fill(Mol_Lab_4Mom.Theta()*180/3.14159265359);
FullMolThetaCMweighted->Fill(Mol_CM_4Mom.Theta()*180/3.14159265359,x_bj); //account for 2 electrons of LH2
FullMolThetaLabweighted->Fill(Mol_Lab_4Mom.Theta()*180/3.14159265359,x_bj/1.19e-6); //account for LH2 density,4e7 incident e
if((Theta>low_limit) && (Theta<high_limit) && (high_limit<40.001))
{
low_limit=low_limit+.001;
high_limit=high_limit+.001;
//printf("ThetaLab=%f\n",Theta);
// for( int phi_step=0;phi_step<301;phi_step++) //Phi -30.0-30.0 degrees, by .2 degree increment
// {
// random_phi=-30.0+(phi_step*.2);
random_phi=random_phi*3.14159265359;
random_phi=random_phi/180;
//Moller CM Px
Mol_CM[0]=cos(random_phi);
Mol_CM[0]=Mol_CM[0]*sqrt((Mol_CM_P*Mol_CM_P)-(Mol_CM[2]*Mol_CM[2]));
//Adjust sign for Px based on quadrant position
if(random_phi>=1.57079632679 && random_phi<=3.14159265359 && Mol_CM[0]>0)
Mol_CM[0]=-Mol_CM[0];
if(random_phi>=-1.57079632679 && random_phi<=-3.14159265359 && Mol_CM[0]>0)
Mol_CM[0]=-Mol_CM[0];
//Moller CM Py
if(random_phi==0)
Mol_CM[1]=0;
if(random_phi!=0)
Mol_CM[1]=sqrt(Mol_CM_P*Mol_CM_P-Mol_CM[2]*Mol_CM[2]-Mol_CM[0]*Mol_CM[0]);
//Adjust sign for Py based on quadrant position
if((random_phi<=0) && (random_phi>=-3.14159265359) && Mol_CM[1]>0)
Mol_CM[1]=-Mol_CM[1];
Mol_CM[1]=Mol_CM[1]/1000;
Mol_CM[0]=Mol_CM[0]/1000;
//printf("Phi=%f",random_phi*180/3.14159265359);
//printf(" Low Limit=%f",low_limit);
//printf(" MomentumTheta=%f\n",atan(Mol_CM[0]/pz)*180/3.14159265359);
//printf(" Theta=%f\n",Theta);
//printf(" PX=%.15f",Mol_CM[0]);
//printf(" PY=%.15f\n",Mol_CM[1]);
fprintf(f0,"%d\t%g\t%g\t%g\t%g\t%g\t%.12f\t%.15f\t%.15f\t%.15f\t%.15f\n",2,1.,1.,1.,1.,1.,x_bj,y,W,Qsqrd,Nu);
fprintf(f0, "%d\t%g\t%g\t%g\t%g\t%g\t%.15f\t%.15f\t%.15f\t%.15f\t%g\t%g\t%g\t%g\n",
1,-1.,1.,11.,0.,0.,-Mol_CM[0],-Mol_CM[1],Pz,KE,0.000511, 0.,0.,0.);
fprintf(f0, "%d\t%g\t%g\t%g\t%g\t%g\t%.15f\t%.15f\t%.15f\t%.15f\t%g\t%g\t%g\t%g\n",
2,-1.,1.,11.,0.,0.,Mol_CM[0],Mol_CM[1],pz,ke,0.000511, 0.,0.,0.);
PhiDistribution->Fill(random_phi*180/3.14159265359);
MolThetaCM->Fill(Mol_CM_4Mom.Theta()*180/3.14159265359); //account for 2 e$
MolThetaLab->Fill(Mol_Lab_4Mom.Theta()*180/3.14159265359);
MolThetaCMweighted->Fill(Mol_CM_4Mom.Theta()*180/3.14159265359,x_bj); //account for 2 electrons of LH2
MolThetaCMweightedReBin->Fill(Mol_CM_4Mom.Theta()*180/3.14159265359,x_bj); //account for 2 electrons of LH2
MolThetaLabweighted->Fill(Mol_Lab_4Mom.Theta()*180/3.14159265359,x_bj/1.19e-6); //account for LH2 density, 4e7 $
// }//End Phi loop
}//End bin condition
//Convert Moller Lab Theta to back to degrees for loop
Mol_Theta_CM=Mol_Theta_CM/3.14159265359;
Mol_Theta_CM=Mol_Theta_CM*180;
};//End Theta loop
tree->Print();
tree->Write();
f->Write();
delete tree;
delete f;
}//End main