DV Plotting XSect

From New IAC Wiki
Revision as of 20:48, 15 April 2016 by Vanwdani (talk | contribs)
Jump to navigation Jump to search

Creating a program called MollerDiffXSect.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
main()
{
float angle,dSigma_dOmega;
double alpha=7.2973525664e-3;
double E_CM=53e6;

                /* using equation from Landau-Lifshitz & Azfar */
                /* units in MeV */
                for(angle=90;angle<=180; angle=angle+.25)
                        {
                        dSigma_dOmega=(alpha*(3+cos(angle*3.14/180)*cos(angle*3.14/180)));
                        dSigma_dOmega=dSigma_dOmega*dSigma_dOmega;

dSigma_dOmega=dSigma_dOmega/(4*E_CM*E_CM*sin(angle*3.14/180)*sin(angle*3.14/180)*sin(angle*3.14/180)*sin(angle*3.14/180));
                        dSigma_dOmega=dSigma_dOmega*1e18;
                        dSigma_dOmega=dSigma_dOmega*.3892e-3;//barns
                        dSigma_dOmega=dSigma_dOmega/1e-6;//micro barns
                        fprintf(stdout," %f\t%f\n", angle, dSigma_dOmega);
                        }

}

Creating a Makefile with contents:

CC = gcc
DiffXSect: MollerDiffXSect.o
        $(CC) $(CFLAGS) MollerDiffXSect.o  -lm -o DiffXSect
MollerDiffXSect.o: MollerDiffXSect.c
        $(CC)  $(CFLAGS) -c MollerDiffXSect.c


The program can be run with the commands:

make DiffXSect
./DiffXSect >DiffXSect.dat

This creates a data file with the scattering angle theta and the corresponding theoretical Moller differential cross-section for said angle. A c++ macro can be written to input this into root

void MollerDiffXSect2Root()
{
    struct evt_t
        {
                Int_t event;
                Float_t Angle, dSigma_dOmega;
        };
  ifstream in;
  in.open("DiffXSect.dat");
  evt_t evt;
  Int_t nlines=0;

  int count=0;
  TFile *f = new TFile("DiffXSect.root","RECREATE");

  TTree *tree = new TTree("DiffXSect","Moller Diff Xsect");

  TH1F *Th=new TH1F("Theory","Theoretical Differential Cross Section",360,90,180);

  tree->Branch("evt",&evt.event,"event/F:Angle/F:dSigma_dOmega");
  while(in.good())
        {
                evt.event=nlines;
                in >> evt.Angle >> evt.dSigma_dOmega;
                        if(nlines==count)
                        {
                                printf( " %d  %f %f\n",evt.event, evt.Angle, evt.dSigma_dOmega);
                                count=count+4;
                        }
                 nlines++;
                tree->Fill();
                Th->Fill(evt.Angle,evt.dSigma_dOmega);
}
 tree->Print();
 tree->Write();
 f->Write();
 in.close();
 delete tree;
 delete f;
}