//ROOT macro for calculation of neutron rates in the BRIKEN detector
//	A. Tarifeno-Saldivia && A. Tolosa, Date: 12Oct16, update 12Oct16

TH1F* h_sum = NULL;

void BRIKEN_rates_v3(){
  gStyle->SetOptStat("nei");
//****************  SETUP  ****************
int Ndet_i=1;
int Ndet_f=140;
Double_t deltaPulser=30.0;
Double_t PulserLowLimit=1350.0;

Double_t Emin = 160.0; // in keV
Double_t Emax = 900.0; // in keV

Double_t PeakAmplitude_threshold=0.005;	// 0.005 should work. 
					// Higher values than 0.08 may affect finding pulser position  

//*****************************************
TString input;//="0cm_2.root";
cout<<"Existing root files in current directory: \n";

 gSystem->Exec("ls *.root");

cout<<"Input root file: \n";
cin >> input;


TFile * file = TFile::Open(input);

ofstream outFile(input + "_data");

TH1F * inter = (TH1F*) file->Get("He001");

h_sum = (TH1F*) inter->Clone();
h_sum -> Reset();
h_sum->SetName("Sum histogram");
h_sum->SetTitle("BRIKEN full histogram: " + input);

TString histName;

TCanvas *C1= new TCanvas("C1", "C1", 1); 
C1->cd(); 

double Counter_nRates[141];
for(int i=0 ; i < 141 ; i++){Counter_nRates[i]=0.0;}

Double_t TotalnRate=0.0;

outFile<<"# File:\t"<<input<<endl;
outFile<<"# Histograms size:\t"<<h_sum->GetSize()<<endl;
outFile<<"# Emin:\t"<<Emin<<" keV\tEmax:\t"<<Emax<<" keV\n";
outFile<<"# deltaPulser:\t"<<deltaPulser<<" keV\tPulserLowLimit:\t"<<PulserLowLimit<<" keV\n";
outFile<<"#####################################################################################################\n";
outFile<<"# Detector\tnRate\tNeutronCounts\tPulserCounts\tpulserPos\tPeakBINo\tPeakBINf\n";
outFile<<"#\t\t(counts\t(counts)\t  (counts)\t(keV)\t \t(bin)\t \t(bin)\t\n";
outFile<<"#\t\t/pulser)\n";
outFile<<"#####################################################################################################\n";


for(int i=Ndet_i; i<=Ndet_f ; i++){
  if(i>0 && i<10) histName= Form("He00%d", i);
  if(i>=10 && i<100) histName= Form("He0%d", i);
  if(i>=100 && i<=140) histName= Form("He%d", i);

   TH1F* hi =(TH1F*)file -> Get(histName);
  
  // Find peak for pulser
   TSpectrum *s = new TSpectrum(100);
   Int_t nfound = s->Search(hi,10,"",PeakAmplitude_threshold);
   TList *functions = hi->GetListOfFunctions();
   TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
   Double_t* pmarrayX=pm->GetX() ; // Find polymarkers position

   Double_t pulserPos, PulserCounts, NeutronCounts, nRate; 
   bool found=false;
   int peakCounter=0;
   
   int BinEmin= hi->GetXaxis()->FindBin(Emin); 
   int BinEmax= hi->GetXaxis()->FindBin(Emax); 

 for (int j = 0; j < nfound; j++) {
      pulserPos=pmarrayX[j];
   if(pulserPos>=PulserLowLimit){
     peakCounter++;
     found=true;
   TF1 *f1 = new TF1("f1","gaus",-.5,.5);
   f1->SetParameters(1000.0, pulserPos, deltaPulser); 
   hi -> Fit("f1", "Q", "",pulserPos-deltaPulser,pulserPos+deltaPulser); 
   Double_t mean = f1->GetParameter(1);
   Double_t sd = f1->GetParameter(2);
   Double_t Peak_o=mean-3.5*sd;
   Double_t Peak_f=mean+3.5*sd;

     int PeakBINo= hi->GetXaxis()->FindBin(Peak_o); 
     int PeakBINf= hi->GetXaxis()->FindBin(Peak_f); 

       PulserCounts = hi->Integral(PeakBINo,PeakBINf);
       NeutronCounts = hi->Integral(BinEmin,BinEmax);
       nRate=NeutronCounts/PulserCounts;
       Counter_nRates[i]=nRate;
       //TotalnRate+=nRate;
       
     outFile<<histName<<"\t"<<nRate<<"\t"<<NeutronCounts<<"\t"<<PulserCounts<<"\t";
     outFile<<pulserPos<<"\t"<<Peak_o<<"\t"<<Peak_f<<"\t"<<endl;
     
    // cout<<histName<<"\t"<<nRate<<"\t"<<NeutronCounts<<"\t"<<PulserCounts<<"\t";
    // cout<<pulserPos<<"\t"<<PeakBINo<<"\t"<<PeakBINf<<"\t"<<hi->GetSize()<<endl;
     

   }
   if(j==nfound-1 && found==false){
     outFile<<histName<<"\t ------->    ERROR FINDING THE PULSER, CHECK PeakAmplitude_threshold !!!!!\n";
     cout<<histName<<"\t ------->    ERROR FINDING THE PULSER, CHECK PeakAmplitude_threshold !!!!!\n";

  }
  if(j==nfound-1 && found==true && peakCounter>1){
    outFile<<histName<<"\t ------->    MORE THAN ONE PULSER POSITION, CHECK PeakAmplitude_threshold !!!!!\n";
    cout<<histName<<"\t ------->    MORE THAN ONE PULSER POSITION, CHECK PeakAmplitude_threshold !!!!!\n";
   
  }
}
h_sum -> Add(hi);

}

outFile<<"\n#######################################################################\n";
outFile<<"#			NEUTRON RATES (counts/pulser)\n";
outFile<<"#######################################################################\n";

//**************** Rates Calculation:
// ADC arrays
int V1A[6][16]={{1,2,3,4,5,6,7,8,9,10,11,12,21,22,23,24},
	    {13,14,15,16,17,18,19,20,31,32,33,34,35,36,37,38},
	    {25,26,27,28,29,30,43,44,45,46,47,48,49,50,51,52},
	    {39,40,41,42,53,54,55,56,57,58,59,60,61,62,63,64},
	    {65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80},
	    {81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,0}};

int V2A[6][8]=	{{96,97,98,99,100,101,102,103},
		{104,105,106,107,108,109,110,0},
		{111,112,113,114,115,116,117,118},
		{119,120,121,122,123,124,125,0},
		{126,127,128,129,130,131,132,133},
		{134,135,136,137,138,139,140,0}};


int Ri[7][28]={{1,2,3,7,8,9,21,13,14,15,19,20,26,27,28,29,30,46,47,48,49,50,0,0,0,0,0,0},
{4,5,10,11,22,23,16,17,33,34,35,36,53,54,55,56,0,0,0,0,0,0,0,0,0,0,0,0},
{6,12,24,18,31,32,37,38,25,43,44,45,51,52,39,40,41,42,57,58,59,60,61,62,63,64,0,0},
{66,67,74,75,81,82,83,84,85,111,112,113,114,115,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{65,68,73,76,86,87,88,89,90,91,116,117,118,119,120,121,0,0,0,0,0,0,0,0,0,0,0,0},
{92,93,94,95,96,97,98,99,100,122,123,124,125,126,127,128,129,130,0,0,0,0,0,0,0,0,0,0},
{69,70,71,72,77,78,79,80,101,102,103,104,105,106,107,108,109,110,131,132,133,134,135,136,137,138,139,140}};


//total
for(int i=0 ; i < 141 ; i++){TotalnRate+=Counter_nRates[i];}
cout<<"#########################################\n";
cout<<"# Total neutron rate (counts/pulser):\t"<<TotalnRate<<endl;
cout<<"#########################################\n";
cout<<"	RATE BY RING:"<<endl;
cout<<"---------------"<<endl;


outFile<<"Total neutron rate:\t"<<TotalnRate<<endl;
outFile<<"#########################################\n";
outFile<<"	BY RING:"<<endl;
outFile<<"#########################################\n";
for(int i=0 ; i < 7 ; i++){
  double ringRATE=0.0;
  for(int j=0 ; j < 28 ; j++){
      ringRATE+=Counter_nRates[Ri[i][j]];
  }
cout<<"R"<<i+1<<"\t"<<ringRATE<<endl;
outFile<<"R"<<i+1<<"\t"<<ringRATE<<endl;
}

outFile<<"#########################################\n";
outFile<<"	BY ADC:"<<endl;
outFile<<"#########################################\n";
cout<<"#########################################\n";
cout<<"	RATE BY ADC:"<<endl;
cout<<"---------------"<<endl;

//By ADC
  //VME1
for(int i=0 ; i < 6 ; i++){
  double vmeRATE=0.0;
  for(int j=0 ; j < 16 ; j++){
      vmeRATE+=Counter_nRates[V1A[i][j]];
  }
cout<<"V1A"<<i+1<<"\t"<<vmeRATE<<endl;
outFile<<"V1A"<<i+1<<"\t"<<vmeRATE<<endl;
}
  //VME2
for(int i=0 ; i < 6 ; i++){
  double vmeRATE=0.0;
  for(int j=0 ; j < 8 ; j++){
      vmeRATE+=Counter_nRates[V2A[i][j]];
  }
cout<<"V2A"<<i+1<<"\t"<<vmeRATE<<endl;
outFile<<"V1A"<<i+1<<"\t"<<vmeRATE<<endl;
}



h_sum -> Draw();

gPad->SetLogy();
outFile.close();
//gApplication->Terminate();

}