#include "c_header.h"
#include "root_header.h"
#include "global_variables.h"
#include "functions.h"

// MAIN SUBROUTINE ================================================================================
int main(int argc, char * const argv[]) {
// Various flags
Int_t	calib_flag    = 0;	// calibration run; 0 if false
Int_t	local_flag    = 1;	// local run (prints progress); 0 if false

Int_t	run_in = 0, i;
Char_t	path[200];

TH1F *hamp1, *hamp2, *hamp3, *hamp4, *hamp5, *hamp6;
TH1F *hamp1r, *hamp2r, *hamp3r, *hamp4r, *hamp5r, *hamp6r;
TH1F *harea1, *harea2, *harea3, *harea4, *harea5, *harea6;

if(argc>=1) {
	for(i=1;i<argc;i++)
		if(*argv[i]=='-') {
			switch (argv[i][1]) {
				case 'p' : strcpy(path,argv[i+1]); break;
				case 'r' : sscanf(argv[i+1],"%d",&run_in); break;
				case 'h':
					printf("Use one of the following options\n");
					printf("\n");
					printf("-n path\n");
					printf("-r Run Number\n");
					printf("\n");          
					return 0;
				default: printf("Execute with -h option for help\n");
					return 0;
			}
		}
}
else {
	printf("Execute with -h option for help.\n");
	return 0;
}

// Get ROOT file
Char_t	fname[256];   
//sprintf(fname,"%srun%d_v%d.root",path,run,dstversion);
//gSystem->Exec("cd ../Beam-ON");
sprintf(fname,"run_%d_in.root",run_in);
//sprintf(fname,"run%d.root",run_in);//Correct name for runs

// Check if ROOT file exists, if not skip to next
FileStat_t	stat;
stat.fSize = 0;

if(gSystem -> GetPathInfo(fname,stat) != 0) {
	printf("*** File %s does not exist! ***\n",fname);
	return 0;
}

// Open ROOT file
TFile f(fname);

printf(" Starting run %d\n",run_in);

// Setup TTree for FIMG
TTree *nt_fimg  = (TTree*)f.Get("FIMG");

int RunNumber;  int segment;  int BunchNumber;  int event;  int movie;  int PSpulse;  float PulseIntensity;  int date;  int time; int detn;
double tflash;  double tof;  double peak_tof;  float amp;  float area;  float fwhm;  int pileup1;  int pileup2; int polarity;  int pulseshape;
float risetime;  int satuflag;  int area_0;  float amp_0;  float area2;

float FIMGneutene;

nt_fimg->SetBranchAddress("RunNumber",&RunNumber);
nt_fimg->SetBranchAddress("event",&event);
nt_fimg->SetBranchAddress("BunchNumber",&BunchNumber);
nt_fimg->SetBranchAddress("PSpulse",&PSpulse);
nt_fimg->SetBranchAddress("PulseIntensity",&PulseIntensity);
nt_fimg->SetBranchAddress("date",&date);
nt_fimg->SetBranchAddress("time",&time);
nt_fimg->SetBranchAddress("detn",&detn);
nt_fimg->SetBranchAddress("tof",&tof);
nt_fimg->SetBranchAddress("peak_tof",&peak_tof);
nt_fimg->SetBranchAddress("amp_0",&amp);
nt_fimg->SetBranchAddress("satuflag",&satuflag);
nt_fimg->SetBranchAddress("area",&area);
nt_fimg->SetBranchAddress("fwhm",&fwhm);
nt_fimg->SetBranchAddress("risetime",&risetime);
nt_fimg->SetBranchAddress("tflash",&tflash);
nt_fimg->SetBranchAddress("pileup1",&pileup1);
nt_fimg->SetBranchAddress("pileup2",&pileup2);
nt_fimg->SetBranchAddress("pulseshape",&pulseshape);
nt_fimg->SetBranchAddress("polarity",&polarity);


// Get number of entries
Int_t	entries = nt_fimg->GetEntries();

printf(" Number of entries: %d\n",entries);

// Create neutron energy binning
Int_t		ndec = E_MAX - E_MIN;
Int_t		nbins = (Int_t) ndec*N_BPDEC;
Double_t	ebins[nbins+1];
Double_t	step = (Double_t) ndec / nbins;

for(Int_t i=0; i <= nbins; i++) {
	ebins[i] = (float) pow(10., step * (Double_t) i + E_MIN);
}

// DEFINE HISTOGRAMS

// Amplitude spectra (full range)
hamp1 = new TH1F("hamp1","Amplitude distribution - EMPTY    #1; Amplitude; Counts",255,0,255);
hamp2 = new TH1F("hamp2","Amplitude distribution - ^{238}U  #2; Amplitude; Counts",255,0,255);
hamp3 = new TH1F("hamp3","Amplitude distribution - ^{240}Pu #3; Amplitude; Counts",255,0,255);
hamp4 = new TH1F("hamp4","Amplitude distribution - ^{240}Pu #4; Amplitude; Counts",255,0,255);
hamp5 = new TH1F("hamp5","Amplitude distribution - ^{240}Pu #5; Amplitude; Counts",255,0,255);
hamp6 = new TH1F("hamp6","Amplitude distribution - ^{235}U  #6; Amplitude; Counts",255,0,255);

// Amplitude spectra (restricted range)
hamp1r = new TH1F("hamp1r","Amplitude distribution - EMPTY    #1 - Restricted En; Amplitude; Counts",255,0,255);
hamp2r = new TH1F("hamp2r","Amplitude distribution - ^{238}U  #2 - Restricted En; Amplitude; Counts",255,0,255);
hamp3r = new TH1F("hamp3r","Amplitude distribution - ^{240}Pu #3 - Restricted En; Amplitude; Counts",255,0,255);
hamp4r = new TH1F("hamp4r","Amplitude distribution - ^{240}Pu #4 - Restricted En; Amplitude; Counts",255,0,255);
hamp5r = new TH1F("hamp5r","Amplitude distribution - ^{240}Pu #5 - Restricted En; Amplitude; Counts",255,0,255);
hamp6r = new TH1F("hamp6r","Amplitude distribution - ^{235}U  #6 - Restricted En; Amplitude; Counts",255,0,255);

hamp1->Sumw2(); hamp2->Sumw2(); hamp3->Sumw2();
hamp4->Sumw2(); hamp5->Sumw2(); hamp6->Sumw2();

hamp1r->Sumw2(); hamp2r->Sumw2(); hamp3r->Sumw2();
hamp4r->Sumw2(); hamp5r->Sumw2(); hamp6r->Sumw2();

// Area Spectra
harea1 = new TH1F("harea1","Area distribution - EMPTY    #1; Area; Counts",8000,0,80000);
harea2 = new TH1F("harea2","Area distribution - ^{238}U  #2; Area; Counts",8000,0,80000);
harea3 = new TH1F("harea3","Area distribution - ^{240}Pu #3; Area; Counts",8000,0,80000);
harea4 = new TH1F("harea4","Area distribution - ^{240}Pu #4; Area; Counts",8000,0,80000);
harea5 = new TH1F("harea5","Area distribution - ^{240}Pu #5; Area; Counts",8000,0,80000);
harea6 = new TH1F("harea6","Area distribution - ^{235}U  #6; Area; Counts",8000,0,80000);

harea1->Sumw2(); harea2->Sumw2(); harea3->Sumw2();
harea4->Sumw2(); harea5->Sumw2(); harea6->Sumw2();

// Amplitude vs. FWHM //
TH2F *hfw1 = new TH2F("hfw1","Amplitude vs. FWHM - EMPTY    #1; FWHM; Amplitude; ",80,0,800,255,0,255);
TH2F *hfw2 = new TH2F("hfw2","Amplitude vs. FWHM - ^{238}U  #2; FWHM; Amplitude; ",80,0,800,255,0,255);
TH2F *hfw3 = new TH2F("hfw3","Amplitude vs. FWHM - ^{240}Pu #3; FWHM; Amplitude; ",80,0,800,255,0,255);
TH2F *hfw4 = new TH2F("hfw4","Amplitude vs. FWHM - ^{240}Pu #4; FWHM; Amplitude; ",80,0,800,255,0,255);
TH2F *hfw5 = new TH2F("hfw5","Amplitude vs. FWHM - ^{240}Pu #5; FWHM; Amplitude; ",80,0,800,255,0,255);
TH2F *hfw6 = new TH2F("hfw6","Amplitude vs. FWHM - ^{235}U  #6; FWHM; Amplitude; ",80,0,800,255,0,255);

TH2F *hfw1r = new TH2F("hfw1r","Amplitude vs. FWHM - EMPTY    #1 - Restricted En; FWHM; Amplitude; ",80,0,800,255,0,255);
TH2F *hfw2r = new TH2F("hfw2r","Amplitude vs. FWHM - ^{238}U  #2 - Restricted En; FWHM; Amplitude; ",80,0,800,255,0,255);
TH2F *hfw3r = new TH2F("hfw3r","Amplitude vs. FWHM - ^{240}Pu #3 - Restricted En; FWHM; Amplitude; ",80,0,800,255,0,255);
TH2F *hfw4r = new TH2F("hfw4r","Amplitude vs. FWHM - ^{240}Pu #4 - Restricted En; FWHM; Amplitude; ",80,0,800,255,0,255);
TH2F *hfw5r = new TH2F("hfw5r","Amplitude vs. FWHM - ^{240}Pu #5 - Restricted En; FWHM; Amplitude; ",80,0,800,255,0,255);
TH2F *hfw6r = new TH2F("hfw6r","Amplitude vs. FWHM - ^{235}U  #6 - Restricted En; FWHM; Amplitude; ",80,0,800,255,0,255);

// Amplitude vs. Area //
TH2F *har1 = new TH2F("har1","Amplitude vs. Area - EMPTY    #1; Area; Amplitude; ",8000,0,80000,255,0,255);
TH2F *har2 = new TH2F("har2","Amplitude vs. Area - ^{238}U  #2; Area; Amplitude; ",8000,0,80000,255,0,255);
TH2F *har3 = new TH2F("har3","Amplitude vs. Area - ^{240}Pu #3; Area; Amplitude; ",8000,0,80000,255,0,255);
TH2F *har4 = new TH2F("har4","Amplitude vs. Area - ^{240}Pu #4; Area; Amplitude; ",8000,0,80000,255,0,255);
TH2F *har5 = new TH2F("har5","Amplitude vs. Area - ^{240}Pu #5; Area; Amplitude; ",8000,0,80000,255,0,255);
TH2F *har6 = new TH2F("har6","Amplitude vs. Area - ^{235}U  #6; Area; Amplitude; ",8000,0,80000,255,0,255);

TH2F *har1r = new TH2F("har1r","Amplitude vs. Area - EMPTY    #1 - Restricted En; Area; Amplitude; ",8000,0,80000,255,0,255);
TH2F *har2r = new TH2F("har2r","Amplitude vs. Area - ^{238}U  #2 - Restricted En; Area; Amplitude; ",8000,0,80000,255,0,255);
TH2F *har3r = new TH2F("har3r","Amplitude vs. Area - ^{240}Pu #3 - Restricted En; Area; Amplitude; ",8000,0,80000,255,0,255);
TH2F *har4r = new TH2F("har4r","Amplitude vs. Area - ^{240}Pu #4 - Restricted En; Area; Amplitude; ",8000,0,80000,255,0,255);
TH2F *har5r = new TH2F("har5r","Amplitude vs. Area - ^{240}Pu #5 - Restricted En; Area; Amplitude; ",8000,0,80000,255,0,255);
TH2F *har6r = new TH2F("har6r","Amplitude vs. Area - ^{235}U  #6 - Restricted En; Area; Amplitude; ",8000,0,80000,255,0,255);

// Amplitude vs. En //
TH2F *hampene1 = new TH2F("hampene1","Amplitude vs. NeutEne - EMPTY    #1; Neutron energy (eV); Amplitude; ",nbins,ebins,255,0,255);
TH2F *hampene2 = new TH2F("hampene2","Amplitude vs. NeutEne - ^{238}U  #2; Neutron energy (eV); Amplitude; ",nbins,ebins,255,0,255);
TH2F *hampene3 = new TH2F("hampene3","Amplitude vs. NeutEne - ^{240}Pu #3; Neutron energy (eV); Amplitude; ",nbins,ebins,255,0,255);
TH2F *hampene4 = new TH2F("hampene4","Amplitude vs. NeutEne - ^{240}Pu #4; Neutron energy (eV); Amplitude; ",nbins,ebins,255,0,255);
TH2F *hampene5 = new TH2F("hampene5","Amplitude vs. NeutEne - ^{240}Pu #5; Neutron energy (eV); Amplitude; ",nbins,ebins,255,0,255);
TH2F *hampene6 = new TH2F("hampene6","Amplitude vs. NeutEne - ^{235}U  #6; Neutron energy (eV); Amplitude; ",nbins,ebins,255,0,255);

// Area vs. En //
TH2F *harene1 = new TH2F("harene1","Area vs. NeutEne - EMPTY    #1; Neutron energy (eV); Area; ",nbins,ebins,8000,0,80000);
TH2F *harene2 = new TH2F("harene2","Area vs. NeutEne - ^{238}U  #2; Neutron energy (eV); Area; ",nbins,ebins,8000,0,80000);
TH2F *harene3 = new TH2F("harene3","Area vs. NeutEne - ^{240}Pu #3; Neutron energy (eV); Area; ",nbins,ebins,8000,0,80000);
TH2F *harene4 = new TH2F("harene4","Area vs. NeutEne - ^{240}Pu #4; Neutron energy (eV); Area; ",nbins,ebins,8000,0,80000);
TH2F *harene5 = new TH2F("harene5","Area vs. NeutEne - ^{240}Pu #5; Neutron energy (eV); Area; ",nbins,ebins,8000,0,80000);
TH2F *harene6 = new TH2F("harene6","Area vs. NeutEne - ^{235}U  #6; Neutron energy (eV); Area; ",nbins,ebins,8000,0,80000);

// Yields //
TH1F *hyield1 = new TH1F("hyield1","Counts vs. NeutEne - EMPTY    #1; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyield2 = new TH1F("hyield2","Counts vs. NeutEne - ^{238}U  #2; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyield3 = new TH1F("hyield3","Counts vs. NeutEne - ^{240}Pu #3; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyield4 = new TH1F("hyield4","Counts vs. NeutEne - ^{240}Pu #4; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyield5 = new TH1F("hyield5","Counts vs. NeutEne - ^{240}Pu #5; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyield6 = new TH1F("hyield6","Counts vs. NeutEne - ^{235}U  #6; Neutron energy (eV); Counts (arb. units)",nbins,ebins);

TH1F *hyieldDead1 = new TH1F("hyieldDead1","Dead-Time Corrected Counts vs. NeutEne - EMPTY    #1; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyieldDead2 = new TH1F("hyieldDead2","Dead-Time Corrected Counts vs. NeutEne - ^{238}U  #2; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyieldDead3 = new TH1F("hyieldDead3","Dead-Time Corrected Counts vs. NeutEne - ^{240}Pu #3; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyieldDead4 = new TH1F("hyieldDead4","Dead-Time Corrected Counts vs. NeutEne - ^{240}Pu #4; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyieldDead5 = new TH1F("hyielddead5","Dead-Time Corrected Counts vs. NeutEne - ^{240}Pu #5; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyieldDead6 = new TH1F("hyieldDead6","Dead-Time Corrected Counts vs. NeutEne - ^{235}U  #6; Neutron energy (eV); Counts (arb. units)",nbins,ebins);

TH1F *hDeadCorFac1 = new TH1F("hDeadCorFac1","Dead-Time Correction Factor - EMPTY    #1; Neutron energy (eV); Correction Factor",nbins,ebins);
TH1F *hDeadCorFac2 = new TH1F("hDeadCorFac2","Dead-Time Correction Factor - ^{238}U  #2; Neutron energy (eV); Correction Factor",nbins,ebins);
TH1F *hDeadCorFac3 = new TH1F("hDeadCorFac3","Dead-Time Correction Factor - ^{240}Pu #3; Neutron energy (eV); Correction Factor",nbins,ebins);
TH1F *hDeadCorFac4 = new TH1F("hDeadCorFac4","Dead-Time Correction Factor - ^{240}Pu #4; Neutron energy (eV); Correction Factor",nbins,ebins);
TH1F *hDeadCorFac5 = new TH1F("hDeadCorFac5","Dead-Time Correction Factor - ^{240}Pu #5; Neutron energy (eV); Correction Factor",nbins,ebins);
TH1F *hDeadCorFac6 = new TH1F("hDeadCorFac6","Dead-Time Correction Factor - ^{235}U  #6; Neutron energy (eV); Correction Factor",nbins,ebins);


TH1F *hyieldRatio = new TH1F("hyieldRatio","Yield Ratio ^{238}U/^{235}U; Neutron energy (eV); Ratio",nbins,ebins);
TH1F *hyieldRatioDead = new TH1F("hyieldRatioDead","Dead-Time Corrected Yield Ratio ^{238}U/^{235}U; Neutron energy (eV); Ratio",nbins,ebins);

TH1F *hyieldSumPu1 = new TH1F("hyieldSumPu1","^{240}Pu Yiel; Neutron energy (eV); Counts (arb. units)",nbins,ebins);
TH1F *hyieldSumPu = new TH1F("hyieldSumPu","^{240}Pu Yield; Neutron energy (eV); Counts (arb. units)",nbins,ebins);

hyield1->Sumw2(); hyield2->Sumw2(); hyield3->Sumw2();
hyield4->Sumw2(); hyield5->Sumw2(); hyield6->Sumw2();
hyieldDead1->Sumw2(); hyieldDead2->Sumw2(); hyieldDead3->Sumw2();
hyieldDead4->Sumw2(); hyieldDead5->Sumw2(); hyieldDead6->Sumw2();
hDeadCorFac1->Sumw2(); hDeadCorFac2->Sumw2(); hDeadCorFac3->Sumw2();
hDeadCorFac4->Sumw2(); hDeadCorFac5->Sumw2(); hDeadCorFac6->Sumw2();

hyieldRatio->Sumw2(); hyieldRatioDead->Sumw2(); hyieldSumPu1->Sumw2(); hyieldSumPu->Sumw2();

// Proton count
TH1F *hPROTpulse = new TH1F("hPROTpulse","; Proton pulse type / Total protons; Number of proton pulses / protons",3,0.5,3.5);
TH2F *hPROTinten = new TH2F("hPROTinten","PSpulse vs. PulseIntensity; Pulse type; Pulse intensity;",2,0.5,2.5,800,1.0E12,9.0E12);

// TOF distributions
TH1F *htof1 = new TH1F("htof1","TOF Det. #1 - EMPTY   ; Time-of-flight (ns); Counts",16000000,0,16000000);
TH1F *htof2 = new TH1F("htof2","TOF Det. #2 - ^{238}U ; Time-of-flight (ns); Counts",16000000,0,16000000);
TH1F *htof3 = new TH1F("htof3","TOF Det. #3 - ^{240}Pu; Time-of-flight (ns); Counts",16000000,0,16000000);
TH1F *htof4 = new TH1F("htof4","TOF Det. #4 - ^{240}Pu; Time-of-flight (ns); Counts",16000000,0,16000000);
TH1F *htof5 = new TH1F("htof5","TOF det. #5 - ^{240}Pu; Time-of-flight (ns); Counts",16000000,0,16000000);
TH1F *htof6 = new TH1F("htof6","TOF Det. #6 - ^{235}U ; Time-of-flight (ns); Counts",16000000,0,16000000);

//  Time Differences between consecutive pulses - Dead Time correction
TH1F *hdt1 = new TH1F("hdt1","#Delta t Det. #1 - EMPTY   ; Time (ns); Counts",1500,0,1500);
TH1F *hdt2 = new TH1F("hdt2","#Delta t Det. #2 - ^{238}U ; Time (ns); Counts",1500,0,1500);
TH1F *hdt3 = new TH1F("hdt3","#Delta t Det. #3 - ^{240}Pu; Time (ns); Counts",1500,0,1500);
TH1F *hdt4 = new TH1F("hdt4","#Delta t Det. #4 - ^{240}Pu; Time (ns); Counts",1500,0,1500);
TH1F *hdt5 = new TH1F("hdt5","#Delta t det. #5 - ^{240}Pu; Time (ns); Counts",1500,0,1500);
TH1F *hdt6 = new TH1F("hdt6","#Delta t Det. #6 - ^{235}U ; Time (ns); Counts",1500,0,1500);

// Tflash distribution
TH1F *htflash1 = new TH1F("htflash1","Tflash Det. #1 - EMPTY   ; Time (ns); Counts",1600000,0,16000000);
TH1F *htflash2 = new TH1F("htflash2","Tflash Det. #2 - ^{238}U ; Time (ns); Counts",1600000,0,16000000);
TH1F *htflash3 = new TH1F("htflash3","Tflash Det. #3 - ^{240}Pu; Time (ns); Counts",1600000,0,16000000);
TH1F *htflash4 = new TH1F("htflash4","Tflash Det. #4 - ^{240}Pu; Time (ns); Counts",1600000,0,16000000);
TH1F *htflash5 = new TH1F("htflash5","Tflash det. #5 - ^{240}Pu; Time (ns); Counts",1600000,0,16000000);
TH1F *htflash6 = new TH1F("htflash6","Tflash Det. #6 - ^{235}U ; Time (ns); Counts",1600000,0,16000000);

htof1->Sumw2();htof2->Sumw2();htof3->Sumw2();
htof4->Sumw2();htof5->Sumw2();htof6->Sumw2();

// Progress indicator
if(local_flag > 0) {
	printf("         10        20        30        40        50        60        70        80        90        100 %%\n");
	printf(" ---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|    \n");
}

Float_t	j = entries/100.0;
Int_t	n = 1;

//Pulse Intensity Correction Factor
double PulseIntensity_CorFac = 1.0e-10;
//Proton Correction Factor
double proton_CorFac = 0.96;

// Energy cuts
Float_t	Ehigh = 1.0E7;// High energy cut
Float_t	Elow  = 1.0E5;	//Low energy cut for restricted range

// Amplitude cuts
Int_t	FIMGcut[7];

FIMGcut[0] = 0;
FIMGcut[1] = 50;
FIMGcut[2] = 45;
FIMGcut[3] = 45;
FIMGcut[4] = 45;
FIMGcut[5] = 35;
FIMGcut[6] = 25;

//FIMGcut[1] = 50;
//FIMGcut[2] = 25;
//FIMGcut[3] = 55;
//FIMGcut[4] = 55;
//FIMGcut[5] = 45;
//FIMGcut[6] = 25;

// Read lamda(E) file, 2 columns: E [eV], l [cm]
/**FILE	*flamda;

flamda=fopen("/afs/cern.ch/user/a/atsingan/ANALYSIS/lamda_mean.dat","r");

for (int m=0; ; m++) {
	fscanf(flamda,"%lf",&lamd);
	if (feof(flamda)) {
		break;
	}
	lamda[m] = lamd;
}

fclose(flamda);*/

// Auxilliary variables for proton/event count
Float_t	protons = 0.0; Int_t bunches = 0; Int_t	Event_old = 0;
// delta t
Double_t tof0 = -1.; Int_t event0 = -1; Int_t detn0 = -1;

// Get entries, fill histograms
for (Int_t i = 0; i <= entries; i++) {
	if(local_flag > 0) {
		if(i == 0) {cout << " ";}
		if(i >= n*j) {cout << '|' << flush; n++;}// Progress indicator
	}

	nt_fimg->GetEntry(i);

	if(tflash == 0.) continue;
	if(satuflag == 1) continue;
	if(RunNumber == 201081 || RunNumber == 201082 || RunNumber == 201122) continue;//runs with a few events
	if(PulseIntensity*PulseIntensity_CorFac <= 5.e12) continue;	

	// Fill TOF distributions
	if(detn == 1) {
		htof1->Fill(tof-tflash);
		htflash1->Fill(tflash);
	}
	else if(detn == 2) {
		htof2->Fill(tof-tflash);
		htflash2->Fill(tflash);
	}
	else if(detn == 3) {
		htof3->Fill(tof-tflash);
		htflash3->Fill(tflash);
	}
	else if(detn == 4) {
		htof4->Fill(tof-tflash);
		htflash4->Fill(tflash);
	}
	else if(detn == 5) {
		htof5->Fill(tof-tflash);
		htflash5->Fill(tflash);
	}
	else if(detn == 6) {
		htof6->Fill(tof-tflash);
		htflash6->Fill(tflash);
	}

	// Proton count
	if(BunchNumber != Event_old) {
		hPROTpulse->Fill(PSpulse);
		protons += PulseIntensity*PulseIntensity_CorFac;
		bunches++;
		hPROTpulse->SetBinContent(3,protons*proton_CorFac);
		hPROTinten->Fill(PSpulse,PulseIntensity*PulseIntensity_CorFac);
		Event_old = BunchNumber;
	}

	// Energy calibration (using TOFtoE() function)
	if(tof > t_relcut) {
		FIMGneutene = (Float_t) TOFtoE_classical(tflash, (double) tof);
	}
	else if(tof <= t_relcut) {
		FIMGneutene = (Float_t) TOFtoE_relativistic(tflash, (double) tof);
	}
/**	if(FIMGneutene < E_relcut) {
		FIMGneutene = (Float_t) TOFtoE_classical_lamda(tflash, (double) tof);
	}
	else {
		FIMGneutene = (Float_t) TOFtoE_relativistic_lamda(tflash, (double) tof);
	}*/

	if(calib_flag > 0) FIMGneutene = (Float_t) TOFtoE_relativistic(0.0, (double) tof);	// if calibration run

	// Temporary check for g-flash anomaly: exclude data 8e3 < En < 3e4 eV
	//if(FIMGneutene > 8.E3 && FIMGneutene < 3.E4) continue;

	// Loop over detectors, fill histograms
	if(FIMGneutene <= pow(10.,E_MAX) && FIMGneutene >= pow(10.,E_MIN)) {// General cut
		if(detn == 1) {
			hampene1->Fill(FIMGneutene,amp);
			harene1->Fill(FIMGneutene,area);

			if(FIMGneutene < Ehigh) {
				hamp1->Fill(amp);
				harea1->Fill(area);
				hfw1->Fill(fwhm,amp);
				har1->Fill(area,amp);

				if(amp > FIMGcut[1]){
					hyield1->Fill(FIMGneutene);
				}

			if(FIMGneutene > Elow) {
				hamp1r->Fill(amp);
				hfw1r->Fill(fwhm,amp);
				har1r->Fill(area,amp);
			}
			}
		}
		else if(detn == 2) {
			hampene2->Fill(FIMGneutene,amp);
			harene2->Fill(FIMGneutene,area);

			if(FIMGneutene < Ehigh) {
				hamp2->Fill(amp);
				harea2->Fill(area);
				hfw2->Fill(fwhm,amp);
				har2->Fill(area,amp);

				if(amp > FIMGcut[2]) hyield2->Fill(FIMGneutene);

			if(FIMGneutene > Elow) {
				hamp2r->Fill(amp);
				hfw2r->Fill(fwhm,amp);
				har2r->Fill(area,amp);
			}
			}
		}
		else if(detn == 3) {
			hampene3->Fill(FIMGneutene,amp);
			harene3->Fill(FIMGneutene,area);

			if(FIMGneutene < Ehigh) {
				hamp3->Fill(amp);
				harea3->Fill(area);
				hfw3->Fill(fwhm,amp);
				har3->Fill(area,amp);

				if(amp > FIMGcut[3]) hyield3->Fill(FIMGneutene);

			if(FIMGneutene > Elow) {
				hamp3r->Fill(amp);
				hfw3r->Fill(fwhm,amp);
				har3r->Fill(area,amp);
			}
			}
		}
		else if(detn == 4) {
			hampene4->Fill(FIMGneutene,amp);
			harene4->Fill(FIMGneutene,area);

			if(FIMGneutene < Ehigh) {
				hamp4->Fill(amp);
				harea4->Fill(area);
				hfw4->Fill(fwhm,amp);
				har4->Fill(area,amp);

				if(amp > FIMGcut[4]) hyield4->Fill(FIMGneutene);

			if(FIMGneutene > Elow) {
				hamp4r->Fill(amp);
				hfw4r->Fill(fwhm,amp);
				har4r->Fill(area,amp);
			}
			}
		}
		else if(detn == 5) {
			hampene5->Fill(FIMGneutene,amp);
			harene5->Fill(FIMGneutene,area);

			if(FIMGneutene < Ehigh) {
				hamp5->Fill(amp);
				harea5->Fill(area);
				hfw5->Fill(fwhm,amp);
				har5->Fill(area,amp);

				if(amp > FIMGcut[5]) hyield5->Fill(FIMGneutene);

			if(FIMGneutene > Elow) {
				hamp5r->Fill(amp);
				hfw5r->Fill(fwhm,amp);
				har5r->Fill(area,amp);
			}
			}
		}
		else if(detn == 6) {
			hampene6->Fill(FIMGneutene,amp);
			harene6->Fill(FIMGneutene,area);

			if(FIMGneutene < Ehigh) {
				hamp6->Fill(amp);
				harea6->Fill(area);
				hfw6->Fill(fwhm,amp);
				har6->Fill(area,amp);

				if(amp > FIMGcut[6]){
					hyield6->Fill(FIMGneutene);
				}

			if(FIMGneutene > Elow) {
				hamp6r->Fill(amp);
				hfw6r->Fill(fwhm,amp);
				har6r->Fill(area,amp);
			}
			}
		}

	}// end of quality check loop
	
	if (BunchNumber == event0 && detn == detn0 && i>1 && (tof-tflash)>1e5){
		if (detn ==1){
			hdt1->Fill(tof-tof0);
			tof0 = tof;
		}
		else if (detn ==2){
			hdt2->Fill(tof-tof0);
			tof0 = tof;
		}
		else if (detn ==3){
			hdt3->Fill(tof-tof0);
			tof0 = tof;
		}
		else if (detn ==4){
			hdt4->Fill(tof-tof0);
			tof0 = tof;
		}
		else if (detn ==5){
			hdt5->Fill(tof-tof0);
			tof0 = tof;
		}
		else if (detn ==6){
			hdt6->Fill(tof-tof0);
			tof0 = tof;
		}
	}
	else {
		event0 = BunchNumber;
		detn0 = detn;
		tof0 = tof;
	}
}// end of loop over entries
for (int i=0; i<hyield1->GetNbinsX(); i++){
		hDeadCorFac1->SetBinContent(i+1,1./(1.-(DeadTime[1]*hyield1->GetBinContent(i+1)/bunches/abs(Energy_to_Time_relativistic(hyield1->GetBinLowEdge(i+1))-Energy_to_Time_relativistic(hyield1->GetBinLowEdge(i+1)+hyield1->GetBinWidth(i+1))))));
		hyieldDead1->SetBinContent(i+1,hyield1->GetBinContent(i+1)*hDeadCorFac1->GetBinContent(i+1));
}
for (int i=0; i<hyield2->GetNbinsX(); i++){
		hDeadCorFac2->SetBinContent(i+1,1./(1.-(DeadTime[2]*hyield2->GetBinContent(i+1)/bunches/abs(Energy_to_Time_relativistic(hyield2->GetBinLowEdge(i+1))-Energy_to_Time_relativistic(hyield2->GetBinLowEdge(i+1)+hyield2->GetBinWidth(i+1))))));
		hyieldDead2->SetBinContent(i+1,hyield2->GetBinContent(i+1)*hDeadCorFac2->GetBinContent(i+1));
}
for (int i=0; i<hyield3->GetNbinsX(); i++){
		hDeadCorFac3->SetBinContent(i+1,1./(1.-(DeadTime[3]*hyield3->GetBinContent(i+1)/bunches/abs(Energy_to_Time_relativistic(hyield3->GetBinLowEdge(i+1))-Energy_to_Time_relativistic(hyield3->GetBinLowEdge(i+1)+hyield3->GetBinWidth(i+1))))));
		hyieldDead3->SetBinContent(i+1,hyield3->GetBinContent(i+1)*hDeadCorFac3->GetBinContent(i+1));
}
for (int i=0; i<hyield4->GetNbinsX(); i++){
		hDeadCorFac4->SetBinContent(i+1,1./(1.-(DeadTime[4]*hyield4->GetBinContent(i+1)/bunches/abs(Energy_to_Time_relativistic(hyield4->GetBinLowEdge(i+1))-Energy_to_Time_relativistic(hyield4->GetBinLowEdge(i+1)+hyield4->GetBinWidth(i+1))))));
		hyieldDead4->SetBinContent(i+1,hyield4->GetBinContent(i+1)*hDeadCorFac4->GetBinContent(i+1));
}
for (int i=0; i<hyield5->GetNbinsX(); i++){
		hDeadCorFac5->SetBinContent(i+1,1./(1.-(DeadTime[5]*hyield5->GetBinContent(i+1)/bunches/abs(Energy_to_Time_relativistic(hyield5->GetBinLowEdge(i+1))-Energy_to_Time_relativistic(hyield5->GetBinLowEdge(i+1)+hyield5->GetBinWidth(i+1))))));
		hyieldDead5->SetBinContent(i+1,hyield5->GetBinContent(i+1)*hDeadCorFac5->GetBinContent(i+1));
}
for (int i=0; i<hyield6->GetNbinsX(); i++){
		hDeadCorFac6->SetBinContent(i+1,1./(1.-(DeadTime[6]*hyield6->GetBinContent(i+1)/bunches/abs(Energy_to_Time_relativistic(hyield6->GetBinLowEdge(i+1))-Energy_to_Time_relativistic(hyield6->GetBinLowEdge(i+1)+hyield6->GetBinWidth(i+1))))));
		hyieldDead6->SetBinContent(i+1,hyield6->GetBinContent(i+1)*hDeadCorFac6->GetBinContent(i+1));
}

// Save number of events for calibration runs
TH1F *hcalibevts = new TH1F("hcalibevts","; ; Number of events",1,0.5,1.5);
if(calib_flag > 0) {
	nt_fimg->GetEntry(entries);
	hcalibevts->SetBinContent(1,event);
}

hyieldRatio->Divide(hyield2,hyield6,mass[6],mass[2]);
hyieldRatioDead->Divide(hyieldDead2,hyieldDead6,mass[6],mass[2]);
hyieldSumPu->Add(hyield3);
hyieldSumPu->Add(hyield4);
hyieldSumPu->Add(hyield5);

cout << endl;

cout << "Total Number of Bunches found  = " << bunches << endl;

// Write output
//sprintf(fname,"%srunv%d_%d_histos_v4_L188.1_2000bpd.root",path,dstversion,run);
//gSystem->Exec("cd Outputs");
sprintf(fname,"run%d_out.root",run_in);
//if(calib_flag > 0) sprintf(fname,"runv%d_%d_histos_beamoff_NEW.root",dstversion,run_in);
cout << " Saving file: " << fname << endl;
TFile fout(fname,"RECREATE");

hamp1->Write();hamp2->Write();hamp3->Write();
hamp4->Write();hamp5->Write();hamp6->Write();

harea1->Write();harea2->Write();harea3->Write();
harea4->Write();harea5->Write();harea6->Write();


hamp1r->Write();hamp2r->Write();hamp3r->Write();
hamp4r->Write();hamp5r->Write();hamp6r->Write();

hfw1->Write();hfw2->Write();hfw3->Write();
hfw4->Write();hfw5->Write();hfw6->Write();

hfw1r->Write();hfw2r->Write();hfw3r->Write();
hfw4r->Write();hfw5r->Write();hfw6r->Write();

har1->Write();har2->Write();har3->Write();
har4->Write();har5->Write();har6->Write();

har1r->Write();har2r->Write();har3r->Write();
har4r->Write();har5r->Write();har6r->Write();

hampene1->Write();hampene2->Write();hampene3->Write();
hampene4->Write();hampene5->Write();hampene6->Write();

harene1->Write();harene2->Write();harene3->Write();
harene4->Write();harene5->Write();harene6->Write();

hyield1->Write();hyield2->Write();hyield3->Write();
hyield4->Write();hyield5->Write();hyield6->Write();
hDeadCorFac1->Write();hDeadCorFac2->Write();hDeadCorFac3->Write();
hDeadCorFac4->Write();hDeadCorFac5->Write();hDeadCorFac6->Write();
hyieldDead1->Write();hyieldDead2->Write();hyieldDead3->Write();
hyieldDead4->Write();hyieldDead5->Write();hyieldDead6->Write();
hyieldRatio->Write(); hyieldRatioDead->Write(); hyieldSumPu->Write(); hyieldSumPu->Write();	

hPROTpulse->Write();
hPROTinten->Write();

htof1->Write();htof2->Write();htof3->Write();
htof4->Write();htof5->Write();htof6->Write();

htflash1->Write();htflash2->Write();htflash3->Write();
htflash4->Write();htflash5->Write();htflash6->Write();

hdt1->Write();hdt2->Write();hdt3->Write();
hdt4->Write();hdt5->Write();hdt6->Write();

hcalibevts->Write();

fout.Close();
return 0;

}