#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();

}