#include "TRandom.h"
#include <TPolyLine.h>
#include <time.h>
#include "TTree.h"
#include <TROOT.h>
#include <TApplication.h>
#include <TRint.h>
#include <TSystem.h>
#include <TH1.h>
#include <TH2.h>
#include <TAxis.h>
#include <TGaxis.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TGraphAsymmErrors.h>
#include <TMultiGraph.h>
#include <TStyle.h>
#include <TKey.h>
#include <TLegend.h>
#include <TColor.h>
#include <TPad.h>
#include <TText.h>
#include <TPaveText.h>
#include <TBox.h>
#include <TLine.h>
#include <TMarker.h>
#include <TLatex.h>
#include <TMath.h>
#include <TF1.h>
#include <TFile.h>
#include <TClass.h>
#include "TBrowser.h"
#include "TFrame.h"
#include "TLegendEntry.h"
#include "TPave.h"
#include <fstream>
#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"
#include <TVirtualFitter.h>
using namespace std;
float massAl=23.;
float massn=1.008665;
float massfac=massAl/(massAl+massn);
//constants
float kT=20000;
const double hbar=1.05457182*pow(10,-34);
const double kb=1.380649*pow(10,-23);
const double eVtoK=11606; // 1 eV is 11606 K
const double eVtojoule=1.60218*pow(10,-19);
const double pi=3.14159265359;
const double amu=1.66054*pow(10,-27);
double mu=massfac*amu;
const double Temp=kT*eVtoK;
const double vT=sqrt(2*kb*Temp/mu);
Double_t getmacs(string filename);
void read(){
TH1F *hmacs=new TH1F("","",5000,0,4);
TH1F *hmacs0=new TH1F("","",5000,0,4);
TH1F *hmacs1=new TH1F("","",5000,0,4);
ifstream list1;
list1.open("list.dat");
string fname1;
string fname2;
string fname3;
string fname4;
string fname5;
string fname6;
if(list1.good()){
while(1){
if(!list1.good())break;
list1>>fname1>>fname2>>fname3>>fname4>>fname5>>fname6;
cout<<fname1<<" "<<fname2<<endl;
Double_t macs0=getmacs(fname1);
Double_t macs02=getmacs(fname2);
Double_t macs1=getmacs(fname3);
Double_t macs12=getmacs(fname4);
Double_t macs13=getmacs(fname5);
Double_t macs14=getmacs(fname6);
double swave=macs0+macs02;
double pwave=macs1+macs12+macs13+macs14;
cout<<fname1<<" "<<swave<<" "<<pwave<<" "<<swave+pwave<<endl;
double total=1000*(swave+pwave);
hmacs->Fill(total);
hmacs0->Fill(swave*1000);
hmacs1->Fill(pwave*1000);
}
list1.close();
}
TFile *f=new TFile("macs.root","recreate");
hmacs->Write("macs");
hmacs0->Write("macs0");
hmacs1->Write("macs1");
f->Write();
}
Double_t getmacs(string filename){
int number, numberold;
numberold=0;
float sum =0;
//TString name1;
Double_t energy, spin, ell, wtot, wn,wg,wx, gw, kernel;
//string name1=filename;
ifstream finp;
finp.open(filename);
//cout<<filename<<endl;
for(int i=0;i<28;i++){
string line; getline(finp,line);}
string line;
while(1){
if(!finp)break;
finp>>number>>energy>>spin>>ell>>wtot>>wn>>wg>>wx>>gw>>kernel;
if(number>numberold){
sum=sum+kernel*exp(-energy/kT);
//cout<<number<< " "<<energy<<" "<<sum<<" "<<kernel<<endl;
number=numberold;
}
}
//cout<<mu<<" "<<kb<<" "<<Temp<<" "<<vT<<endl;
double macs=pow(2*pi/mu/kb/Temp,1.5)*hbar*hbar*eVtojoule*sum*1E28/vT;
//cout<<"MACS = "<<macs<<" b"<<endl;
finp.close();
return macs;
}
void test(){
int number, numberold;
numberold=0;
float sum =0;
//TString name1;
Double_t energy, spin, ell, wtot, wn,wg,wx, gw, kernel;
//string name1=filename;
ifstream finp;
finp.open("run1_25m_0l_1");
for(int i=0;i<28;i++){
string line; getline(finp,line);}
string line;
while(1){
if(!finp)break;
finp>>number>>energy>>spin>>ell>>wtot>>wn>>wg>>wx>>gw>>kernel;
if(number>numberold){
sum=sum+kernel*exp(-energy/kT);
cout<<number<< " "<<energy<<" "<<sum<<" "<<kernel<<endl;
number=numberold;
}
}
//cout<<mu<<" "<<kb<<" "<<Temp<<" "<<vT<<endl;
double macs=pow(2*pi/mu/kb/Temp,1.5)*hbar*hbar*eVtojoule*sum*1E28/vT;
cout<<"MACS = "<<macs<<" b"<<endl;
finp.close();
}
|