Difference between revisions of "LUND Spread LH2 0Phi.C"

From New IAC Wiki
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