Difference between revisions of "Creating uniform LUND files"

From New IAC Wiki
Jump to navigation Jump to search
 
Line 25: Line 25:
 
                         dSigma_dOmega=dSigma_dOmega*1e18;
 
                         dSigma_dOmega=dSigma_dOmega*1e18;
 
                         dSigma_dOmega=dSigma_dOmega*.3892e-3;//barns
 
                         dSigma_dOmega=dSigma_dOmega*.3892e-3;//barns
                        dSigma_dOmega=dSigma_dOmega/1e-6;//micro barns
 
                        dSigma_dOmega=dSigma_dOmega/10; //Adjust for 10 electrons
 
 
                         //printf("Angle, XSect");
 
                         //printf("Angle, XSect");
 
                         //printf(" %f\t%f\n", angle, dSigma_dOmega);
 
                         //printf(" %f\t%f\n", angle, dSigma_dOmega);

Latest revision as of 17:16, 26 April 2016

#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