diff --git a/.gitignore b/.gitignore
index 04b07e5..9a17e9a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,14 @@
 *.root
 *.evt
+*.evt.to
 *.ev2
 *.log
+
 xia2root
 xia2ev2*
+pixie2root
+scan
+
+*.so
+*.d
+*.pcm
diff --git a/Analyzer.C b/Analyzer.C
new file mode 100644
index 0000000..d9b3caa
--- /dev/null
+++ b/Analyzer.C
@@ -0,0 +1,189 @@
+#define Analyzer_cxx
+
+#include "Analyzer.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TH1.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <vector>
+#include <TStopwatch.h>
+
+//############################################ User setting
+
+int rawEnergyRange[2] = {100,  2000};
+
+double BGO_threshold = 100;
+
+//############################################ end of user setting
+
+ULong64_t NumEntries = 0;
+ULong64_t ProcessedEntries = 0;
+Float_t Frac = 0.1; ///Progress bar
+TStopwatch StpWatch;
+
+//############################################ histogram declaration
+
+TH2F * heVID;
+TH1F * he[NCLOVER];
+
+TH2F * hgg[NCLOVER][NCLOVER];
+
+TH2F * hcoin;
+
+///----- after calibration
+TH2F * heCalVID;
+TH1F * heCal[NCLOVER]; // BGO veto
+
+
+void Analyzer::Begin(TTree * tree){
+
+   TString option = GetOption();
+
+   NumEntries = tree->GetEntries();
+
+   printf("======================== histogram declaration\n");
+   heVID    = new TH2F("heVID",                                              "e vs ID; det ID; e [ch]", NCLOVER, 0, NCLOVER, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
+   heCalVID = new TH2F("heCalVID", Form("eCal vs ID (BGO veto > %.1f); det ID; e [ch]", BGO_threshold), NCLOVER, 0, NCLOVER, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
+   for( int i = 0; i < NCLOVER; i ++){
+      he[i]    = new TH1F(   Form("he%02d", i),                                  Form("e -%02d", i), rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
+      heCal[i] = new TH1F(Form("heCal%02d", i), Form("e -%02d (BGO veto > %.1f)", i, BGO_threshold), rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
+   }
+   
+   for( int i = 0; i < NCLOVER; i++){
+      for( int j = i; j < NCLOVER; j++){
+         hgg[i][j] = new TH2F(Form("hg%02dg%02d", i, j), Form(" e%02d vs e%02d; e%02d; e%02d", i, j, i, j), 
+               rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1], 
+               rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
+      }
+   }
+   
+   hcoin = new TH2F("hcoin", "detector coin.; det ID; det ID", NCLOVER, 0, NCLOVER, NCLOVER, 0 , NCLOVER); 
+   
+   printf("======================== End of histograms Declaration\n");
+   StpWatch.Start();
+
+}
+
+Bool_t Analyzer::Process(Long64_t entry){
+
+   ProcessedEntries++;
+   
+   /*********** Progress Bar ******************************************/ 
+   if (ProcessedEntries>NumEntries*Frac-1) {
+      TString msg; msg.Form("%llu", NumEntries/1000);
+      int len = msg.Sizeof();
+      printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n",
+               Frac*100, len, ProcessedEntries/1000,NumEntries/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac);
+      StpWatch.Start(kFALSE);
+      Frac+=0.1;
+   }
+
+   b_energy->GetEntry(entry);
+   b_time->GetEntry(entry);
+   b_pileup->GetEntry(entry);
+   b_bgo->GetEntry(entry);
+   b_other->GetEntry(entry);
+   b_multiplicity->GetEntry(entry);
+   
+   if( multi == 0 ) return kTRUE;
+
+   for( int detID = 0; detID < NCLOVER ; detID ++){
+      
+      //======== baics gate when no energy or pileup 
+      if( TMath::IsNaN(e[detID])) continue;
+      //if( pileup[detID] == 1 ) continue;
+      
+      //======== Fill raw data
+      heVID->Fill( detID, e[detID]);
+      he[detID]->Fill(e[detID]);
+      
+      for( int detJ = detID; detJ < NCLOVER; detJ++) {
+         if( TMath::IsNaN(e[detJ])) continue;
+         hgg[detID][detJ]->Fill(e[detID], e[detJ]); // x then y
+         hcoin->Fill(detID, detJ); 
+      }
+      
+      //======== BGO veto
+      for( int kk = 0; kk < NBGO ; kk++){
+         if( TMath::IsNaN(bgo[kk]) ) continue;
+         if( bgo[kk] > BGO_threshold ) {
+            return kTRUE;
+         }
+      }
+      
+      //========= apply correction
+      double eCal = e[detID];
+      
+      heCalVID->Fill( detID, eCal);
+      heCal[detID]->Fill(eCal);
+      
+   }
+   
+   return kTRUE;
+}
+
+
+void listDraws(void) {
+  printf("------------------- List of Plots -------------------\n");
+  printf("      rawID() - Raw e vs ID\n");
+  printf("      drawE() - Raw e for all %d detectors\n", NCLOVER);
+  printf("-----------------------------------------------------\n");
+}
+
+void rawID(){
+  TCanvas * cRawID = (TCanvas *) gROOT->FindObjectAny("cRawID");
+  if( cRawID == NULL ) cRawID = new TCanvas("cRawID", "raw ID", 1000, 800);
+  cRawID->cd(1)->SetGrid();
+  heVID->Draw("colz");
+}
+
+void drawE(bool isLogy = false, bool cali = false){
+
+   TCanvas *cRawE = (TCanvas *) gROOT->FindObjectAny("cRawE");
+   if( cRawE == NULL ) cRawE = new TCanvas("cRawE", "raw e", 800, 1500);
+   cRawE->Clear();cRawE->Divide(4,9);
+   for (Int_t i = 0; i < 36; i++) {
+      cRawE->cd(i+1); 
+      cRawE->cd(i+1)->SetGrid();
+      if( isLogy ) cRawE->cd(i+1)->SetLogy();
+      if( cali ) {
+         heCal[i]->Draw("");
+      }else{
+         he[i]->Draw("");
+      }
+   }
+
+}
+
+
+
+
+void Analyzer::Terminate(){
+
+   printf("============================== finishing.\n");
+   gROOT->cd();
+
+   int canvasXY[2] = {1200 , 1600} ;// x, y
+   int canvasDiv[2] = {1,2};
+   TCanvas *cCanvas  = new TCanvas("cCanvas", "" ,canvasXY[0],canvasXY[1]);
+   cCanvas->Modified(); cCanvas->Update();
+   cCanvas->cd(); cCanvas->Divide(canvasDiv[0],canvasDiv[1]);
+
+   gStyle->SetOptStat("neiou");
+
+   cCanvas->cd(1);
+   cCanvas->cd(1)->SetLogz(1);
+   heVID->Draw("colz");
+
+   cCanvas->cd(2);
+   cCanvas->cd(2)->SetLogz(1);
+   heCalVID->Draw("colz");
+   
+   listDraws();
+   
+   gROOT->ProcessLine(".L AutoFit.C");
+   printf("=============== loaded AutoFit.C\n");
+   
+
+}
diff --git a/Analyzer.h b/Analyzer.h
new file mode 100644
index 0000000..26c1c81
--- /dev/null
+++ b/Analyzer.h
@@ -0,0 +1,115 @@
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Mon Dec 13 18:37:46 2021 by ROOT version 6.24/06
+// from TTree tree/tree
+// found on file: efEu152b-000.root
+//////////////////////////////////////////////////////////
+
+#ifndef Analyzer_h
+#define Analyzer_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TSelector.h>
+
+#include "mapping.h"
+
+// Header file for the classes stored in the TTree if any.
+
+class Analyzer : public TSelector {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   ULong64_t       evID;
+   Double_t        e[NCLOVER];
+   ULong64_t       t[NCLOVER];
+   UShort_t        p[NCLOVER];
+   Double_t        bgo[NBGO];
+   Double_t        other[NOTHER];
+   Int_t           multi;
+
+   // List of branches
+   TBranch        *b_event_ID;   //!
+   TBranch        *b_energy;   //!
+   TBranch        *b_time;   //!
+   TBranch        *b_pileup;   //!
+   TBranch        *b_bgo;   //!
+   TBranch        *b_other;   //!
+   TBranch        *b_multiplicity;   //!
+
+   Analyzer(TTree * /*tree*/ =0) : fChain(0) { }
+   virtual ~Analyzer() { }
+   virtual Int_t   Version() const { return 2; }
+   virtual void    Begin(TTree *tree);
+   virtual void    SlaveBegin(TTree *tree);
+   virtual void    Init(TTree *tree);
+   virtual Bool_t  Notify();
+   virtual Bool_t  Process(Long64_t entry);
+   virtual Int_t   GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
+   virtual void    SetOption(const char *option) { fOption = option; }
+   virtual void    SetObject(TObject *obj) { fObject = obj; }
+   virtual void    SetInputList(TList *input) { fInput = input; }
+   virtual TList  *GetOutputList() const { return fOutput; }
+   virtual void    SlaveTerminate();
+   virtual void    Terminate();
+
+   ClassDef(Analyzer,0);
+};
+
+#endif
+
+#ifdef Analyzer_cxx
+void Analyzer::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the branch addresses and branch
+   // pointers of the tree will be set.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("evID",   &evID, &b_event_ID);
+   fChain->SetBranchAddress("e",          e, &b_energy);
+   fChain->SetBranchAddress("t",          t, &b_time);
+   fChain->SetBranchAddress("p",          p, &b_pileup);
+   fChain->SetBranchAddress("bgo",      bgo, &b_bgo);
+   fChain->SetBranchAddress("other",  other, &b_other);
+   fChain->SetBranchAddress("multi", &multi, &b_multiplicity);
+}
+
+Bool_t Analyzer::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+
+
+void Analyzer::SlaveBegin(TTree * /*tree*/){
+
+   TString option = GetOption();
+
+}
+
+
+void Analyzer::SlaveTerminate(){
+
+}
+
+
+#endif // #ifdef Analyzer_cxx
diff --git a/AutoFit.C b/AutoFit.C
new file mode 100644
index 0000000..0773037
--- /dev/null
+++ b/AutoFit.C
@@ -0,0 +1,2701 @@
+#ifndef AutoFit_C
+#define AutoFit_C
+
+#include <TF1.h>
+#include <TGraph.h>
+#include <TColor.h>
+#include <TSpectrum.h>
+#include <TMath.h>
+#include <vector>
+
+void ShowFitMethod(){
+  printf("----------------------- Method of Fitting ---------------\n");
+  printf("---------------------------------------------------------\n");
+  printf("         fitAuto() - estimate BG, find peak, and fit n-Gauss \n");
+  printf("     fitGaussPol() - fit 1 Gauss + pol-n BG\n");
+  printf("     fit2GaussP1() - fit 2 Gauss + pol-1 BG \n");
+  printf("       fitGF3Pol() - fit GF3 + pol-n BG \n");
+  //printf("         fitNGF3() - fit n-GF3, estimated BG \n");
+  printf("       fitNGauss() - fit n-Gauss, estimated BG, need input\n");
+  printf("    fitNGaussSub() - fit estimated BG with Pol, subtract, fit n-Guass\n");
+  printf("    fitNGaussPol() - fit n-Gauss + pol-n BG \n");
+  printf(" fitNGaussPolSub() - subtract Pol-n BG, fit n-Gauss \n");
+  printf("\n");
+  printf("------- Mouse click Fit : \n");
+  printf("    clickFitNGaussPol() - fit n-Gauss + pol-n BG \n");
+  printf(" clickFitNGaussPolSub() - Fit Pol-n BG, subtract, fit n-Gauss\n");
+  printf("          SaveFitPara() - Save the initial guess parameters.\n");
+  printf("---------------------------------------------------------\n");
+}
+
+std::vector<std::string> SplitStrAF(std::string tempLine, std::string splitter, int shift = 0){
+
+  std::vector<std::string> output;
+
+  size_t pos;
+  do{
+    pos = tempLine.find(splitter); // fine splitter
+    if( pos == 0 ){ //check if it is splitter again
+      tempLine = tempLine.substr(pos+1);
+      continue;
+    }
+
+    std::string secStr;
+    if( pos == std::string::npos ){
+      secStr = tempLine;
+    }else{
+      secStr = tempLine.substr(0, pos+shift);
+      tempLine = tempLine.substr(pos+shift);
+    }
+
+    //check if secStr is begin with space
+    if( secStr.substr(0, 1) == " "){
+      secStr = secStr.substr(1);
+    }
+
+    output.push_back(secStr);
+    //printf(" |%s---\n", secStr.c_str());
+    
+  }while(pos != std::string::npos );
+
+  return output;
+}
+
+TColor RGBWheel(double ang){
+
+  ang = ang * TMath::DegToRad();
+  
+  double r = max(0., (1+2*cos(ang))/3.);
+  double g = max(0., (1 - cos(ang) + sqrt(3)* sin(ang))/3.);
+  double b = max(0., (1 - cos(ang) - sqrt(3)* sin(ang))/3.);
+
+  TColor col(r,g,b);
+
+  return col;
+
+}
+
+int nPeaks = 16;
+Double_t nGauss(Double_t *x, Double_t *par) {
+   Double_t result = 0;
+   for (Int_t p=0;p<nPeaks;p++) {
+      Double_t norm  = par[3*p+0];
+      Double_t mean  = par[3*p+1];
+      Double_t sigma = par[3*p+2];
+      result += norm * TMath::Gaus(x[0],mean,sigma, 1); // normalized Gaussian
+   }
+   return result;
+}
+
+int nPols = 1;
+
+Double_t nPolFunc(Double_t *x, Double_t *par) {
+  Double_t result = 0;
+  
+  for( int p = 0; p < nPols+1; p++){
+    result += par[p]*TMath::Power(x[0], p);
+  }
+  return result;
+}
+
+Double_t nGaussPol(Double_t *x, Double_t *par) {
+
+  Double_t result = 0;
+  for (Int_t p=0;p<nPeaks;p++) {
+    Double_t norm  = par[3*p+0];
+    Double_t mean  = par[3*p+1];
+    Double_t sigma = par[3*p+2];
+    result += norm * TMath::Gaus(x[0],mean,sigma, 1); // normalized Gaussian
+  }
+
+  for( int p = 0; p < nPols+1; p++){
+    result += par[3*nPeaks + p]*TMath::Power(x[0], p);
+  }
+
+  return result;
+}
+
+Double_t nGF3(Double_t *x, Double_t *par){
+
+  /** this is the fitting function for gamma peak from gf3, RadWare **/
+  Double_t result = 0;
+  for( Int_t p = 0; p < nPeaks; p++){
+    
+    Double_t norm  = par[6*p+0];
+    Double_t mean  = par[6*p+1];
+    Double_t sigma = par[6*p+2];
+    Double_t ratio = par[6*p+3];
+    Double_t beta  = par[6*p+4]; // skewness
+    Double_t step  = par[6*p+5];
+    
+    result += norm * (1-ratio)* TMath::Gaus(x[0], mean, sigma, 1) ;
+    result += norm * ratio * exp( sigma * sigma/2/beta/beta) / (2* beta )* exp((x[0]-mean)/beta) * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) + sigma/beta/sqrt(2)) ;
+    result += norm * step * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) );
+  }  
+  return result;
+}
+
+
+Double_t nGF3Pol(Double_t *x, Double_t *par){
+
+  /** this is the fitting function for gamma peak from gf3, RadWare **/
+  Double_t result = 0;
+  for( Int_t p = 0; p < nPeaks; p++){
+    Double_t norm  = par[6*p+0];
+    Double_t mean  = par[6*p+1];
+    Double_t sigma = par[6*p+2];
+    Double_t ratio = par[6*p+3];
+    Double_t beta  = par[6*p+4]; // skewness
+    Double_t step  = par[6*p+5];
+
+    result += norm * (1.0-ratio)* TMath::Gaus(x[0], mean, sigma, 1) ;
+    result += norm * ratio * exp( sigma * sigma/2/beta/beta) / (2* beta )* exp((x[0]-mean)/beta) * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) + sigma/beta/sqrt(2)) ;
+    result += norm * step * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) );
+  }
+
+  for( int p = 0; p < nPols + 1; p++){
+    result += par[6*nPeaks + p]*TMath::Power(x[0], p);
+  }
+  
+  return result;
+}
+
+void PrintPar(TF1 * fit, int numParPerPeak){
+   int totPar = fit->GetNpar();
+   int count = totPar/numParPerPeak;
+   printf("ID  ");
+   for( int i = 0; i < numParPerPeak; i++){
+     printf("par%d       ", i);
+   }
+   printf("\n");
+   
+   for( int i = 0; i < count ; i++){
+     printf("%3d ", i);
+     for( int j = 0; j < numParPerPeak; j++){
+       int parID = numParPerPeak * i + j;
+       printf("%.3f(%.3f) ", fit->GetParameter(parID), fit->GetParError(parID));
+     }
+     printf("\n");
+   }
+}
+
+void GoodnessofFit(TH1F * hist, TF1 * fit){
+  
+   int nBin = hist->GetNbinsX();
+   int effBin = 0;   
+   double mean = 0;
+   double ysq = 0;
+   double SSR = 0;
+   double chisq = 0; //with estimated error be sqrt(y)
+   double Xsq = 0; // for Pearson's chi-sq test
+   for( int i = 1; i <= nBin; i++){
+
+      double e = hist->GetBinError(i);
+      if( e > 0  ) {
+         effBin ++;
+         double y = hist->GetBinContent(i);
+         double x = hist->GetBinCenter(i);
+         double ybar = fit->Eval(x);
+         ysq += y*y;
+         mean + y;
+         SSR += (y - ybar)*(y-ybar);
+         chisq += (y - ybar)*(y-ybar)/e/e;
+         
+         
+         if( ybar > e ) {
+           Xsq += (y - ybar)*(y-ybar)/ybar;
+         }else{
+           Xsq += (y - ybar)*(y-ybar)/e;
+         }
+         //printf(" %d | x : %f, y : %f, ybar : %f , X-sq : %f\n", i, x, y, ybar, Xsq);
+      }
+   }
+   mean = mean / nBin;
+   double SSTotal = ysq + mean*mean;
+   
+   int npar = fit->GetNpar();
+   int ndf = effBin - npar;
+   printf("#################################################\n");
+   printf("##               Goodness of Fit.              ##\n");
+   printf("#################################################\n");
+   printf(" eff. Bin(%d) - numPar(%d) = ndf = %d \n", effBin, npar, ndf);
+
+   printf("============== Regression analysis\n");
+   printf("----------------- R-sq \n");
+   printf("    SSTotal = %f \n", SSTotal);
+   printf("        SSR = %f \n", SSR);
+   printf("        MSR = %f  <-- is it similar to sample variance?\n", SSR/ndf);
+   double Rsq = 1 - SSR/SSTotal;
+   printf("       R-sq = %f \n", Rsq );
+   
+   printf("----------------- Chi-sq \n");
+   printf("     Chi-sq = %f \n", chisq);
+   printf(" rd. Chi-sq = %f \n", chisq/ndf);
+   printf("ROOT Chi-Sq = %f , NDF = %d \n", fit->GetChisquare(), fit->GetNDF());
+   //================ chi-sq test
+   printf("============== Hypothesis testing\n");
+   printf(" Null Hypothesis : the fitting model is truth. \n");
+   printf(" * p-value = prob. that Null Hypo. is truth. \n");
+   printf(" * the Pearson's test in here only for background free data \n");
+   printf("     Pearson's X-sq = %.2f \n", Xsq);
+   double p = 1 - TMath::Prob(Xsq, ndf);
+   printf("  Pearson's p-value = %.2f %s 0.05 | %s\n", p, p < 0.05 ? "<": ">", p < 0.05 ? "reject" : "cannot reject");
+   double pchi = 1 - TMath::Prob(chisq, ndf);
+   printf("     Chi-sq p-value = %.2f %s 0.05 | %s\n", pchi, pchi < 0.05 ? "<": ">", pchi < 0.05 ? "reject" : "cannot reject");
+   double pRoot = 1- fit->GetProb();
+   printf("ROOT Chi-sq p-value = %.2f %s 0.05 | %s\n", pRoot, pRoot < 0.05 ? "<": ">", pRoot < 0.05 ? "reject" : "cannot reject");
+   printf("################################################\n");  
+}
+
+vector<double> energy, height, sigma, lowE, highE ;
+vector<int> energyFlag, sigmaFlag;
+
+bool loadFitParameters(TString fitParaFile){
+  
+  energy.clear(); energyFlag.clear();
+  sigma.clear();  sigmaFlag.clear();
+  
+  lowE.clear(); highE.clear();
+  
+  height.clear();
+  
+  bool paraRead = false;
+
+  printf("====================================================================== \n");
+  printf("----- loading fit parameters from : %s", fitParaFile.Data());
+  ifstream file;
+  file.open(fitParaFile.Data());
+
+  if( !file){
+    printf("\ncannot read file : %s \n", fitParaFile.Data());
+    return 0;
+  }
+
+  while( file.good()) {
+    string tempLine;
+    getline(file, tempLine);
+
+    if( tempLine.substr(0, 1) == "#" ) continue;
+    if( tempLine.substr(0, 2) == "//" ) continue;
+    if( tempLine.size() < 5 ) continue;
+    
+    ///printf("%s\n", tempLine.c_str());
+
+    vector<string> temp = SplitStrAF(tempLine, " ");
+
+    if( temp.size() < 7 ) continue;
+
+    energy.push_back(    atof(temp[0].c_str()));
+    lowE.push_back(      atof(temp[1].c_str()));
+    highE.push_back(     atof(temp[2].c_str()));
+    energyFlag.push_back(atoi(temp[3].c_str()));         
+    sigma.push_back(     atof(temp[4].c_str()));
+    sigmaFlag.push_back( atoi(temp[5].c_str()));
+    height.push_back(    atof(temp[6].c_str()));
+
+  }
+
+  printf("... done.\n");
+    
+  int n = energy.size();
+  TString limStr = "(fixed)";
+  printf("%2s| %34s | %10s \n", "ID", "Peak [MeV]", "Sigma [MeV]");
+  for( int j = 0; j < n; j ++){    
+    if( energyFlag[j] == 0 ) limStr.Form("(limited, %6.3f - %6.3f)", lowE[j], highE[j]);
+    printf("%2d| %7.4f %-26s | %7.4f (%5s) \n", j, energy[j], limStr.Data(), sigma[j], sigmaFlag[j] == 1 ? "fixed" : "free"); 
+  }
+    
+  paraRead = true;
+
+  printf("====================================================================== \n");
+
+  return  paraRead;
+  
+}
+
+TCanvas * NewCanvas(TString name, TString title, int divX, int divY, int padSizeX, int padSizeY){
+  TCanvas * output = NULL;
+  if( gROOT->FindObjectAny(name) == NULL ){
+    output = new TCanvas(name, title, divX * padSizeX, divY * padSizeY);
+  }else{
+    output = (TCanvas *) gROOT->FindObjectAny(name) ;
+    output->Clear();
+  }
+  output->Divide(divX, divY);
+  return output;
+}
+
+void ScaleAndDrawHist(TH1F * hist, double xMin, double xMax){
+
+  if ( xMin != xMax ) hist->GetXaxis()->SetRangeUser(xMin, xMax);
+  int maxBin = hist->GetMaximumBin();
+  double ymax = hist->GetBinContent(maxBin);
+  hist->GetYaxis()->SetRangeUser(0, 1.1 * ymax);
+  hist->Draw();
+  
+}
+
+void PlotResidual(TH1F * hist, TF1 * fit){
+
+  TH1F * hRes = (TH1F*) hist->Clone();
+  hRes->GetListOfFunctions()->Clear();
+  hRes->SetTitle("Residual");
+  hRes->SetName("hRes");
+  hRes->SetYTitle("Hist - fit");
+  hRes->Sumw2(0);
+  hRes->Sumw2(1);
+  hRes->Add(fit, -1);
+  hRes->Draw();
+  
+}
+
+
+//########################################
+//### Fit a Gauss + Pol-n
+//########################################
+void fitGaussPol(TH1F * hist, double mean, double sigmaMax, int degPol, double xMin, double xMax, TString optStat = ""){
+
+  printf("=========================================================\n");
+  printf("================ fit 1-Gauss + Pol-%d BG ================\n", degPol);
+  printf(" * mean Range  +- 5 sigmaMax \n");
+  printf(" * inital parameters of the polynomial is random/pow(10, i) \n"); 
+  printf("==========================================================\n");
+  
+  gStyle->SetOptStat(optStat);
+  TCanvas * cFitGaussPol = NewCanvas("cFitGaussPol", Form("fit Gauss & Pol-%d | fitGaussPol", degPol), 1, 2, 800, 350);  
+  cFitGaussPol->cd(1);
+  
+  ScaleAndDrawHist(hist, xMin, xMax);
+
+  const int nPar = 3 + degPol + 1;
+
+  TString funcExp = "[0] * TMath::Gaus(x, [1], [2], 1)";
+  for( int i = 0; i < degPol+1 ; i++){
+    funcExp += Form(" + [%d]*TMath::Power(x,%d)", i+3 , i);
+  }  
+  TF1 * fit = new TF1("fit", funcExp.Data(), xMin, xMax);
+
+  double * para = new double[nPar]; 
+  para[0] = 100 * 0.05 * TMath::Sqrt(TMath::TwoPi());
+  para[1] = mean;
+  para[2] = sigmaMax/2.;
+  for( int i = 0 ; i < degPol+1; i++){
+    para[3 + i ] = gRandom->Rndm()/TMath::Power(10, i);
+  }
+
+  fit->SetLineWidth(2);
+  fit->SetLineColor(1);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+
+  fit->SetParLimits(0, 0, 1e9);
+  fit->SetParLimits(1, mean - 5*sigmaMax, mean + 5 * sigmaMax);
+  fit->SetParLimits(2, 0, sigmaMax);
+
+  hist->Fit("fit", "Rq");
+
+  const Double_t* paraE = fit->GetParErrors();
+  const Double_t* paraA = fit->GetParameters();
+  
+  double bw = hist->GetBinWidth(1);
+
+  printf("histogram name : %s \n====== Gaussian:\ncount: %8.0f(%3.0f)\nmean : %8.4f(%8.4f)\nsigma: %8.4f(%8.4f) \n", 
+            hist->GetName(),
+            paraA[0] / bw,   paraE[0] /bw, 
+            paraA[1], paraE[1],
+            paraA[2], paraE[2]);
+
+  printf("------- the polnomail BG:\n");
+  for(int i = 0 ; i < degPol+1; i++){
+    printf("%2d | %8.4f(%8.4f) \n", i, paraA[3+i], paraE[3+i]);
+  }
+                        
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.04);
+  
+  double chi2 = fit->GetChisquare();
+  int ndf = fit->GetNDF();
+  text.DrawLatex(0.12, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+  text.DrawLatex(0.12, 0.75,Form("count: %4.0f(%3.0f)",  paraA[0] / bw,   paraE[0] /bw));
+  text.DrawLatex(0.12, 0.70,Form("E_{x}: %6.3f(%5.3f) MeV",  paraA[1], paraE[1]));
+  text.DrawLatex(0.12, 0.65,Form("#sigma: %3.0f(%3.0f) keV",  paraA[2] * 1000., paraE[2] * 1000.));
+
+  for( int i = 0; i < degPol + 1; i++){
+    text.DrawLatex(0.60, 0.85 - 0.05*i ,Form("%3s: %8.3f(%8.3f)", Form("p%d", i),  paraA[3+i], paraE[3+i]));
+  }
+
+  TString expression = "[0] ";
+  for( int j = 1; j < degPol + 1; j++){
+    expression += Form(" + [%d]*TMath::Power(x, %d)", j, j);
+  }
+  TF1 * g0 = new TF1("g0", expression.Data(), xMin, xMax);
+  for( int j = 0; j < degPol + 1 ; j++){
+    g0->SetParameter(j, paraA[3+j]);
+  }
+  g0->SetLineColor(4);
+  g0->Draw("same");
+
+  TF1* f0 = new TF1("f0", "[0] * TMath::Gaus(x, [1], [2], 1)", xMin, xMax);
+  f0->SetParameter(0, paraA[0]);
+  f0->SetParameter(1, paraA[1]);
+  f0->SetParameter(2, paraA[2]);
+  f0->SetLineColor(2);
+  f0->SetNpx(1000);
+  f0->Draw("same");
+
+//  GoodnessofFit(hist, fit);
+
+  cFitGaussPol->cd(2);
+  PlotResidual(hist, fit);
+
+}
+
+//########################################
+//#### fit 2 gauss + pol-1  // not updated
+//########################################
+vector<double> fit2GaussP1(TH1F * hist, double mean1, double sigma1, 
+                                       double mean2, double sigma2, 
+                           double xMin, double xMax, TString optStat = "", bool newCanvas = false){
+
+
+  printf("=========================================================\n");
+  printf("================ fit 2-Gauss + Pol-1 BG ================\n" );
+  printf(" NOT updated. It works, but the code is old \n"); 
+  printf("==========================================================\n");
+  
+  
+  vector<double> output;
+  output.clear();
+
+  gStyle->SetOptStat(optStat);  
+  TCanvas * cFit2GaussP1 = NewCanvas("cFit2GaussP1", "fit Gauss & P1  | fit2GaussP1", 1, 1, 800, 350);
+  cFit2GaussP1->cd(1);
+
+  ScaleAndDrawHist(hist, xMin, xMax);
+  
+  TF1 * fit = new TF1("fit", "[0] * TMath::Gaus(x, [1], [2], 1) + [3] * TMath::Gaus(x, [4], [5], 1) + [6] + [7]*x", xMin, xMax);
+  
+  double * para = new double[8]; 
+  para[0] = 20 * 0.05 * TMath::Sqrt(TMath::TwoPi());
+  para[1] = mean1;
+  para[2] = sigma1;
+  para[3] = 100 * 0.05 * TMath::Sqrt(TMath::TwoPi());
+  para[4] = mean2;
+  para[5] = sigma2;
+  para[6] = 1;
+  para[7] = 0;
+
+  fit->SetLineWidth(2);
+  fit->SetLineColor(2);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+  
+  hist->Fit("fit", "Rq");
+  
+  const Double_t* paraE = fit->GetParErrors();
+  const Double_t* paraA = fit->GetParameters();
+  
+  double bw = hist->GetBinWidth(1);
+
+  printf("%7s ====== count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", 
+            hist->GetName(),
+            paraA[0] / bw,   paraE[0] /bw, 
+            paraA[1], paraE[1],
+            paraA[2], paraE[2]);
+  printf("%7s ====== count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", 
+            "",
+            paraA[3] / bw,   paraE[3] /bw, 
+            paraA[4], paraE[4],
+            paraA[5], paraE[5]);
+            
+  output.push_back( paraA[0]/bw);
+  output.push_back( paraE[0]/bw);
+  output.push_back( paraA[1]);
+  output.push_back( paraE[1]);
+  output.push_back( paraA[2]);
+  output.push_back( paraE[2]);
+
+  output.push_back( paraA[3]/bw);
+  output.push_back( paraE[3]/bw);
+  output.push_back( paraA[4]);
+  output.push_back( paraE[4]);
+  output.push_back( paraA[5]);
+  output.push_back( paraE[5]);
+            
+            
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.04);
+  
+  double chi2 = fit->GetChisquare();
+  int ndf = fit->GetNDF();
+  text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+
+  text.DrawLatex(0.15, 0.75,Form("count: %4.0f(%3.0f), E_{x}: %6.3f(%5.3f) MeV, #sigma: %3.0f(%3.0f) keV ", 
+                                    paraA[0] / bw,   paraE[0] /bw,
+                                    paraA[1], paraE[1],
+                                    paraA[2] * 1000., paraE[2] * 1000.));
+  text.DrawLatex(0.15, 0.7, Form("count: %4.0f(%3.0f), E_{x}: %6.3f(%5.3f) MeV, #sigma: %3.0f(%3.0f) keV  ", 
+                                    paraA[3] / bw,   paraE[3] /bw,
+                                    paraA[4], paraE[4],
+                                    paraA[5] * 1000., paraE[5] * 1000.));
+                                    
+  text.DrawLatex(0.15, 0.6, Form("Line : %6.3f(%5.3f) + %6.3f(%5.3f)x ",
+                                    paraA[6], paraE[6],
+                                    paraA[7], paraE[7]));
+                                    
+  GoodnessofFit(hist, fit);
+  return output;
+}
+
+
+//########################################
+//#### fit for gamma + pol-n BG
+//########################################
+void fitGF3Pol(TH1F * hist, double mean, double sigmaMax, double ratio, double beta, double step, int degPol, double xMin, double xMax, TString optStat = ""){
+
+  printf("=========================================================\n");
+  printf("================ fit GF1 + Pol-%d BG ================\n", degPol);
+  printf(" * mean Range  = xMin, xMax \n");
+  printf(" * inital parameters of the polynomial is random/pow(10, i) \n"); 
+  printf("==========================================================\n");
+  
+
+  gStyle->SetOptStat(optStat);
+
+  gStyle->SetOptStat(optStat);  
+  TCanvas * cFitGF3Pol = NewCanvas("cFitGF3Pol", Form("fit GF3 + pol-%d | fitGF3Pol", degPol), 1, 2, 800, 350);
+  cFitGF3Pol->cd(1);
+
+  ScaleAndDrawHist(hist, xMin, xMax);
+  
+  nPeaks = 1;
+  nPols = degPol;
+  int nPar = 6*nPeaks + degPol + 1;
+
+  TF1 * fit = new TF1("fit", nGF3Pol, xMin, xMax, nPar);
+
+  fit->Print();
+
+  double * para = new double[nPar]; 
+  para[0] = hist->GetMaximum() *4;
+  para[1] = mean;
+  para[2] = sigmaMax/2.;
+  para[3] = ratio ;
+  para[4] = beta ;
+  para[5] = step ;
+  for( int i = 0 ; i < degPol + 1; i++){
+    para[6+i] = gRandom->Rndm()/TMath::Power(10, i);
+  }
+  
+  fit->SetLineWidth(2);
+  fit->SetLineColor(1);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+
+  fit->SetParLimits(0, 0, 1e9);
+  fit->SetParLimits(1, xMin, xMax);
+  fit->SetParLimits(2, 0.00000001, sigmaMax);
+  fit->SetParLimits(3, 0, 0.5);
+  fit->SetParLimits(4, 1, 400);
+  fit->SetParLimits(5, 0, 0.5);
+  
+  hist->Fit("fit", "Rq");
+  
+  const Double_t* paraE = fit->GetParErrors();
+  const Double_t* paraA = fit->GetParameters();
+  
+  double chisquare = fit->GetChisquare();
+  int ndf = fit->GetNDF();
+  double bw = hist->GetBinWidth(1);
+
+  printf("histogram : %s \n", hist->GetName());
+  printf("========= The Gaussian \n");
+  printf("count: %8.0f(%3.0f)\n",  paraA[0] / bw,   paraE[0] /bw);
+  printf("mean : %8.4f(%8.4f)\n",  paraA[1], paraE[1]);
+  printf("sigma: %8.4f(%8.4f)\n",  paraA[2], paraE[2]);
+            
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.04);
+  text.SetTextColor(1);
+  
+  text.DrawLatex(0.12, 0.65, Form("count : %5.0f(%5.0f)", paraA[0]/bw, paraE[0]/bw));
+  text.DrawLatex(0.12, 0.60, Form(" mean : %5.3f(%5.3f) keV", paraA[1], paraE[1]));
+  text.DrawLatex(0.12, 0.55, Form("sigma : %5.3f(%5.3f) keV", paraA[2], paraE[2]));
+  text.DrawLatex(0.12, 0.50, Form(" FWHM : %5.3f(%5.3f) keV", paraA[2] *2.355, paraE[2]*2.355));
+   
+  text.DrawLatex(0.12, 0.40, Form("#chi^2/ndf : %5.3f", chisquare/ndf));
+   
+   //GoodnessofFit(hist, fit);
+
+  ///                  0       1        2        3       4       5   
+  string label[8] = {"Area", "mean", "sigma", "ratio", "beta", "step"};
+  printf("---------- The detail\n");
+  for(int i = 0 ; i < 6 ; i++){
+    printf("%d | %8s | %f (%f) \n", i, label[i].c_str(), paraA[i], paraE[i]);
+    text.DrawLatex(0.65, 0.85-0.05*i, Form("%6s: %5.3f(%5.3f)", label[i].c_str(), paraA[i], paraE[i]));
+  }
+  for(int i = 6 ; i < nPar; i++){
+    printf("%d | %8s | %f (%f) \n", i, Form("p%d", i-6), paraA[i], paraE[i]);
+    text.DrawLatex(0.65, 0.85-0.05*i, Form("%6s: %5.3f (%5.3f) \n", Form("p%d", i-6), paraA[i], paraE[i]));
+  }
+
+  /// norm * (1-ratio)* TMath::Gaus(x[0], mean, sigma, 1)
+  TF1 * g0 = new TF1("g0", "[0] * (1.0-[3]) * TMath::Gaus(x, [1], [2], 1)", xMin, xMax);
+  g0->SetParameter(0, paraA[0]);
+  g0->SetParameter(1, paraA[1]);
+  g0->SetParameter(2, paraA[2]);
+  g0->SetParameter(3, paraA[3]);
+  g0->SetNpx(1000);
+  g0->SetLineColor(kRed);
+
+  ///                     norm * ratio * exp( sigma * sigma/2/beta/beta)* exp((x[0]-mean)/beta) * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) + sigma/beta/sqrt(2)) ;
+  TF1 * g1 = new TF1("g1", "[0] * [3] * exp( [2] * [2]/2/[4]/[4]) / (2* [4])* exp((x-[1])/[4]) * TMath::Erfc( (x-[1])/([2] * sqrt(2)) + [2]/[4]/sqrt(2)) ", xMin, xMax);
+  g1->SetParameter(0, paraA[0]);
+  g1->SetParameter(1, paraA[1]);
+  g1->SetParameter(2, paraA[2]);
+  g1->SetParameter(3, paraA[3]);
+  g1->SetParameter(4, paraA[4]);
+  g1->SetNpx(1000);
+  g1->SetLineColor(kGreen +3);
+
+  /// norm * step * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) );
+  TF1 * g2 = new TF1("g2", "[0] * [3] * TMath::Erfc( (x-[1])/([2] * sqrt(2)) );", xMin, xMax);
+  g2->SetParameter(0, paraA[0]);
+  g2->SetParameter(1, paraA[1]);
+  g2->SetParameter(2, paraA[2]);
+  g2->SetParameter(3, paraA[5]);
+  g2->SetNpx(1000);
+  g2->SetLineColor(kViolet);
+
+
+  TString expression = "[0] ";
+  for( int j = 1; j < degPol + 1; j++){
+    expression += Form(" + [%d]*TMath::Power(x, %d)", j, j);
+  }
+  TF1 * g3 = new TF1("g3", expression.Data(), xMin, xMax);
+  for( int j = 0; j < degPol + 1 ; j++){
+    g3->SetParameter(j, paraA[6+j]);
+  }
+  g3->SetLineColor(kBlue);
+  g3->Draw("same");
+
+
+  g0->Draw("same");
+  g1->Draw("same");
+  g2->Draw("same");
+  g3->Draw("same");
+
+  cFitGF3Pol->cd(2);
+  PlotResidual(hist, fit);
+   
+}
+
+//##############################################
+//##### Auto Fit n-Gauss with estimated BG 
+//##############################################
+vector<double> fitAuto(TH1F * hist, int bgEst = 10, 
+                       double peakThreshold = 0.05, 
+                       double sigmaMax = 0, 
+                       int peakDensity = 4, 
+                       TString optStat = ""){
+
+  printf("================================================================\n");
+  printf("========== Auto Fit n-Gauss with estimated BG ==================\n");
+  printf(" * bgEst = parameter of BG estimation, larger BG, more linear \n");
+  printf(" * peakThreshold = precentage of the highest peak that will count \n");
+  printf(" * sigmaMax  =  maximum sigma, if -1, program try to find the sigma \n");
+  printf(" * peakDensity = peak will closer when the number is larger ");
+  printf(" \n");
+  printf(" after peaks found, the i-th peaks will be limited by the mid-point\n");
+  printf("   by the (i-1)-th peak and the i-th peak, and the mid-point of the\n");
+  printf("    i-th peak and (i+1)-th peak \n");
+  printf("   i.e. [peak(i-1)+peak(i)]/2  < limit of peak(i) < [peak(i)+peak(i+1)]/2 \n");
+  printf("================================================================\n");
+
+  gStyle->SetOptStat(optStat);
+  TCanvas *cFitAuto = NewCanvas("cFitAuto","Auto Fitting | fitAuto", 1, 4, 800, 300);
+  cFitAuto->cd(1);
+
+  ScaleAndDrawHist(hist, 0, 0);
+  
+  TH1F * specS = (TH1F*) hist->Clone();
+  double xMin = hist->GetXaxis()->GetXmin();
+  double xMax = hist->GetXaxis()->GetXmax();
+  int xBin = hist->GetXaxis()->GetNbins();
+  
+  TString titleH;
+  titleH.Form("fitted spectrum (BG=%d); Ex [MeV]; Count / %4.0f keV", bgEst, (xMax-xMin)*1000./xBin );
+  specS->SetTitle(titleH);   
+  specS->SetName("specS");
+  ///specS->GetXaxis()->SetTitleSize(0.06);
+  ///specS->GetYaxis()->SetTitleSize(0.06);
+  ///specS->GetXaxis()->SetTitleOffset(0.7);
+  ///specS->GetYaxis()->SetTitleOffset(0.6);
+  
+  //=================== find peak and fit
+  gStyle->SetOptFit(0);
+  TSpectrum * peak = new TSpectrum(50);
+  nPeaks = peak->Search(hist, peakDensity, "", peakThreshold); 
+  
+  if( bgEst > 0 ) {
+    printf("============= estimating background...\n");
+    TH1 * h1 = peak->Background(hist, bgEst);
+    h1->Draw("same");
+    printf("============= substracting the linear background...\n");
+    specS->Sumw2();
+    specS->Add(h1, -1.);
+  }
+  
+  cFitAuto->cd(2)->SetGrid();
+  cFitAuto->cd(2);
+  specS->Draw("hist");
+
+  //========== Fitting 
+  printf("============= Fitting.....");
+  printf(" found %d peaks \n", nPeaks);
+  
+  double * xpos = peak->GetPositionX();
+  double * ypos = peak->GetPositionY();
+
+  int * inX = new int[nPeaks];
+  TMath::Sort(nPeaks, xpos, inX, 0 );
+  vector<double> energy, height;
+  for( int j = 0; j < nPeaks; j++){
+    energy.push_back(xpos[inX[j]]);
+    height.push_back(ypos[inX[j]]);
+  }
+  for( int j = 0; j < nPeaks; j++){
+    printf(" energy : %f , %f \n", energy[j], height[j]);
+  }
+  
+  
+  if( sigmaMax == 0 ){
+    printf("------------- Estimate sigma for each peak \n");
+    sigma.clear();
+    int binMin = hist->FindBin(xMin);
+    int binMax = hist->FindBin(xMax);
+    for( int i = 0; i < nPeaks ; i++){
+      int b0 = hist->FindBin(energy[i]);
+      double sMin = (xMax-xMin)/5., sMax = (xMax-xMin)/5.;
+      //---- backward search, stop when 
+      for( int b = b0-1 ; b > binMin ; b-- ){
+        double y = hist->GetBinContent(b);
+        double x = hist->GetBinCenter(b);
+        if( y  < (height[i])/2. ) {
+          sMin = energy[i] - hist->GetBinCenter(b);
+          break;
+        }
+      }
+      //---- forward search, stop when 
+      for( int b = b0+1 ; b < binMax  ; b++ ){
+        double y = hist->GetBinContent(b);
+        double x = hist->GetBinCenter(b);
+        if( y < (height[i])/2. ) {
+          sMax = hist->GetBinCenter(b) - energy[i];
+          break;
+        }
+      }
+
+      double temp = TMath::Min(sMin, sMax);
+      /// When there are multiple peaks closely packed : 
+      if( i > 0 && temp > 2.5 * sigma.back() ) temp = sigma.back();
+      sigma.push_back(temp);
+      
+      printf("%2d | x : %8.2f | sigma(est) %f \n", i, energy[i], sigma[i]);
+    }
+  }else if( sigmaMax < 0 ){
+    printf("========== use user input sigma : %f (fixed)\n", abs(sigmaMax));
+    sigma.clear();
+    for( int i = 0; i < nPeaks ; i++) sigma.push_back(abs(sigmaMax));
+  }else if( sigmaMax > 0 ){
+    printf("========== use user input sigma : %f/2. \n", sigmaMax/2.);
+    sigma.clear();
+    for( int i = 0; i < nPeaks ; i++) sigma.push_back(sigmaMax/2.);
+  }
+    
+  
+  
+  int numParPerPeak = 3;  
+  const int  n = numParPerPeak * nPeaks;
+  double * para = new double[n]; 
+  for(int i = 0; i < nPeaks ; i++){
+    para[numParPerPeak*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi());
+    para[numParPerPeak*i+1] = energy[i];
+    if( sigmaMax == 0 ){
+      para[numParPerPeak*i+2] = sigma[i];
+    }else if(sigmaMax < 0 ){
+      para[numParPerPeak*i+2] = abs(sigmaMax);
+    }else if(sigmaMax > 0 ){
+      para[numParPerPeak*i+2] = sigmaMax/2.;
+    }
+  }
+
+  TF1 * fit = new TF1("fit", nGauss, xMin, xMax, 3 * nPeaks );
+  fit->SetLineWidth(2);
+  fit->SetLineColor(2);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+  
+  if( nPeaks > 1 ){
+    for( int i = 0; i < nPeaks; i++){
+      fit->SetParLimits(numParPerPeak*i+0, 0, 1e+9);
+      double de1 = 1, de2 = 1;
+      if( i == 0 ){ 
+        de2 = (energy[i+1] - energy[i])/2.;
+        de1 = de2;
+      }else if( i < nPeaks -1 ){
+        de1 = (energy[i] - energy[i-1])/2.;
+        de2 = (energy[i+1] - energy[i])/2.;
+      }else{
+        de1 = (energy[i] - energy[i-1])/2.;
+        de2 = de1;
+      }
+      
+      fit->SetParLimits(numParPerPeak*i+1, energy[i] - de1 , energy[i] + de2);
+      if( sigmaMax== 0 ) fit->SetParLimits(numParPerPeak*i+2, 0, 1.5*sigma[i]); // add 50% margin of sigma
+      if( sigmaMax < 0 ) fit->FixParameter(numParPerPeak*i+2, abs(sigmaMax));
+      if( sigmaMax > 0 ) fit->SetParLimits(numParPerPeak*i+2, 0, sigmaMax);
+    }
+  }else{
+    fit->SetParLimits(0, 0, 1e+9);
+    fit->SetParLimits(2, 0, sigmaMax); 
+  } 
+  
+  specS->Fit("fit", "q");
+  
+  
+  const Double_t* paraE = fit->GetParErrors();
+  const Double_t* paraA = fit->GetParameters();
+  
+  //======== calculate reduce chi-squared
+  //GoodnessofFit(specS, fit);
+  
+  double bw = specS->GetBinWidth(1);
+
+  vector<double> exPos;
+  for(int i = 0; i < nPeaks ; i++){
+    exPos.push_back(paraA[numParPerPeak*i+1]);
+    printf("%2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", 
+            i, 
+            paraA[numParPerPeak*i] / bw,   paraE[numParPerPeak*i] /bw, 
+            paraA[numParPerPeak*i+1], paraE[numParPerPeak*i+1],
+            paraA[numParPerPeak*i+2], paraE[numParPerPeak*i+2]);
+    
+  }
+  cFitAuto->Update();
+
+  //draw the indivual fit
+  fit->Draw("same");
+  
+  const int nn = nPeaks;
+  TF1 ** gFit = new TF1 *[nn];
+  for( int i = 0; i < nPeaks; i++){
+      gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax);
+      gFit[i]->SetParameter(0, paraA[numParPerPeak*i]);
+      gFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]);
+      gFit[i]->SetParameter(2, paraA[numParPerPeak*i+2]);
+      gFit[i]->SetLineColor(i+1);
+      gFit[i]->SetNpx(1000);
+      gFit[i]->SetLineWidth(1);
+      gFit[i]->Draw("same");
+  }
+  
+  specS->Draw("hist same");
+
+  //======== print text on plot
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.04);
+
+  double chi2 = fit->GetChisquare();
+  int ndf = fit->GetNDF();
+  text.SetTextSize(0.06);
+  text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+
+  cFitAuto->cd(3);
+  PlotResidual(specS, fit);
+
+  cFitAuto->cd(4);
+  text.SetTextSize(0.05);
+  text.SetTextColor(2);
+  
+  text.DrawLatex(0.1, 0.9, Form("      %13s, %18s, %18s", "count", "mean", "sigma"));
+  
+  for( int i = 0; i < nPeaks; i++){
+  text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]));
+  }
+  
+  return exPos;
+  
+}
+
+
+//########################################
+//###### NOT tested
+//########################################
+vector<double> fitNGF3(TH1 * hist, int bgEst = 10, 
+                       double peakThreshold = 0.1, 
+                       double sigmaMax = 20, 
+                       int peakDensity = 4, 
+                       TString optStat = "", bool newPlot = true){
+
+  TCanvas *cFitAuto = NULL;
+  if( newPlot ){
+    cFitAuto = new TCanvas("cFitAuto","Auto Fitting", 100, 100, 800,800);
+    cFitAuto->Divide(1,2);
+    
+    gStyle->SetOptStat(optStat);
+    cFitAuto->cd(1);
+    hist->Draw();
+  }
+  
+  TH1F * specS = (TH1F*) hist->Clone();
+  double xMin = hist->GetXaxis()->GetXmin();
+  double xMax = hist->GetXaxis()->GetXmax();
+  int xBin = hist->GetXaxis()->GetNbins();
+  
+  TString titleH;
+  titleH.Form("fitted spectrum (BG=%d); Ex [MeV]; Count / %4.0f keV", bgEst, (xMax-xMin)*1000./xBin );
+  specS->SetTitle(titleH);   
+  specS->SetName("specS");
+  ///specS->GetXaxis()->SetTitleSize(0.06);
+  ///specS->GetYaxis()->SetTitleSize(0.06);
+  ///specS->GetXaxis()->SetTitleOffset(0.7);
+  ///specS->GetYaxis()->SetTitleOffset(0.6);
+  
+  //=================== find peak and fit
+  gStyle->SetOptFit(0);
+  TSpectrum * peak = new TSpectrum(50);
+  nPeaks = peak->Search(hist, peakDensity, "", peakThreshold); 
+  
+  if( bgEst > 0 ) {
+    printf("============= estimating background...\n");
+    TH1 * h1 = peak->Background(hist, bgEst);
+    h1->Draw("same");
+    printf("============= substracting the linear background...\n");
+    specS->Sumw2();
+    specS->Add(h1, -1.);
+  }
+  
+  if( newPlot ){
+    cFitAuto->cd(2)->SetGrid();
+    cFitAuto->cd(2);
+  }
+  specS->Draw("hist");
+  
+
+  //========== Fitting 
+  if( newPlot ){
+    printf("============= Fitting.....");
+    printf(" found %d peaks \n", nPeaks);
+  }
+  double * xpos = peak->GetPositionX();
+  double * ypos = peak->GetPositionY();
+
+  int * inX = new int[nPeaks];
+  TMath::Sort(nPeaks, xpos, inX, 0 );
+  vector<double> energy, height;
+  for( int j = 0; j < nPeaks; j++){
+    energy.push_back(xpos[inX[j]]);
+    height.push_back(ypos[inX[j]]);
+  }
+  if( newPlot ){
+    for( int j = 0; j < nPeaks; j++){
+      printf(" energy : %f , %f \n", energy[j], height[j]);
+    }
+  }
+  
+  int numParPerPeak = 6;  
+  const int  n = numParPerPeak * nPeaks;
+  double * para = new double[n]; 
+  for(int i = 0; i < nPeaks ; i++){
+    para[numParPerPeak*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi());
+    para[numParPerPeak*i+1] = energy[i];
+    para[numParPerPeak*i+2] = sigmaMax;
+    para[numParPerPeak*i+3] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()) * 0.1;
+    para[numParPerPeak*i+4] = sigmaMax;
+    para[numParPerPeak*i+5] = -4;
+  }
+
+  TF1 * fit = new TF1("fit", nGF3, xMin, xMax, numParPerPeak * nPeaks );
+  fit->SetLineWidth(2);
+  fit->SetLineColor(2);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+  
+  if( nPeaks > 1 ){
+    for( int i = 0; i < nPeaks; i++){
+      fit->SetParLimits(numParPerPeak*i+0, 0, 1e+9);
+      double de1 = 1, de2 = 1;
+      if( i == 0 ){ 
+        de2 = (energy[i+1] - energy[i])/2.;
+        de1 = de2;
+      }else if( i < nPeaks -1 ){
+        de1 = (energy[i] - energy[i-1])/2.;
+        de2 = (energy[i+1] - energy[i])/2.;
+      }else{
+        de1 = (energy[i] - energy[i-1])/2.;
+        de2 = de1;
+      }
+      
+      fit->SetParLimits(numParPerPeak*i+1, energy[i] - de1 , energy[i] + de2); 
+      fit->SetParLimits(numParPerPeak*i+2, 0, sigmaMax * 5); 
+      
+      fit->SetParLimits(numParPerPeak*i+3, 0, 1e+9);
+      fit->SetParLimits(numParPerPeak*i+4, 0, sigmaMax); 
+      fit->SetParLimits(numParPerPeak*i+5, -10, -2); 
+          
+    }
+  }else{
+    fit->SetParLimits(0, 0, 1e+9);
+    fit->SetParLimits(2, 0, sigmaMax); 
+    fit->SetParLimits(3, 0, 1e+9);
+    fit->SetParLimits(4, 0, sigmaMax); 
+    fit->SetParLimits(5, -10, -2);
+  }
+  specS->Fit("fit", "q");
+  
+  
+  const Double_t* paraE = fit->GetParErrors();
+  const Double_t* paraA = fit->GetParameters();
+  
+  //======== calculate reduce chi-squared
+  if( newPlot ) GoodnessofFit(specS, fit);
+  
+  double bw = specS->GetBinWidth(1);
+
+  vector<double> exPos;
+  for(int i = 0; i < nPeaks ; i++){
+    exPos.push_back(paraA[numParPerPeak*i+1]);
+    double totCount = paraA[numParPerPeak*i] + paraA[numParPerPeak*i+3];
+    double totCountErr = sqrt(paraE[numParPerPeak*i]*paraE[numParPerPeak*i] + paraE[numParPerPeak*i+3]*paraE[numParPerPeak*i+3]);
+    printf("%2d , count: %8.0f(%3.0f)+%8.0f(%3.0f)=%8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f), skewneww: %4.1f(%4.1f) \n", 
+            i,
+            paraA[numParPerPeak*i] / bw,   paraE[numParPerPeak*i] /bw, 
+            paraA[numParPerPeak*i+3] / bw,   paraE[numParPerPeak*i+3] /bw, 
+            totCount / bw,   totCountErr /bw, 
+            paraA[numParPerPeak*i+1], paraE[numParPerPeak*i+1],
+            paraA[numParPerPeak*i+2], paraE[numParPerPeak*i+2],
+            paraA[numParPerPeak*i+5], paraE[numParPerPeak*i+5]);
+    
+    //PrintPar(fit, numParPerPeak);
+  }
+  if( newPlot ) cFitAuto->Update();
+
+  //draw the indivual fit
+  fit->Draw("same");
+  
+  const int nn = nPeaks;
+  TF1 ** gFit = new TF1 *[nn];
+  TF1 ** kFit = new TF1 *[nn];
+  TF1 ** zFit = new TF1 *[nn];
+  for( int i = 0; i < nPeaks; i++){
+      gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1) + [3] * TMath::Gaus(x,[1],[4], 1) * ( 1 + TMath::Erf( [5]*(x-[1])/sqrt(2)/[4] ))", xMin, xMax);
+      gFit[i]->SetParameter(0, paraA[numParPerPeak*i]);
+      gFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]);
+      gFit[i]->SetParameter(2, paraA[numParPerPeak*i+2]);
+      gFit[i]->SetParameter(3, paraA[numParPerPeak*i+3]);
+      gFit[i]->SetParameter(4, paraA[numParPerPeak*i+4]);
+      gFit[i]->SetParameter(5, paraA[numParPerPeak*i+5]);
+      gFit[i]->SetLineColor(i+1);
+      gFit[i]->SetNpx(1000);
+      gFit[i]->SetLineWidth(1);
+      gFit[i]->Draw("same");
+      
+      kFit[i] = new TF1(Form("kFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1) * ( 1 + TMath::Erf( [3]*(x-[1])/sqrt(2)/[2] ))", xMin, xMax);
+      kFit[i]->SetParameter(0, paraA[numParPerPeak*i+3]);
+      kFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]);
+      kFit[i]->SetParameter(2, paraA[numParPerPeak*i+4]);
+      kFit[i]->SetParameter(3, paraA[numParPerPeak*i+5]);
+      kFit[i]->SetLineColor(i+1);
+      kFit[i]->SetNpx(1000);
+      kFit[i]->SetLineWidth(1);
+      kFit[i]->Draw("same");
+      
+      zFit[i] = new TF1(Form("zFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax);
+      zFit[i]->SetParameter(0, paraA[numParPerPeak*i]);
+      zFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]);
+      zFit[i]->SetParameter(2, paraA[numParPerPeak*i+2]);
+      zFit[i]->SetLineColor(i+1);
+      zFit[i]->SetNpx(1000);
+      zFit[i]->SetLineWidth(1);
+      zFit[i]->Draw("same");
+      
+  }
+  
+  specS->Draw("hist same");
+  
+  return exPos;
+  
+}
+
+//########################################
+//###### fir N Gauss with estimated BG
+//########################################
+void fitNGauss(TH1F * hist, int bgEst = 10, TString fitFile = "AutoFit_para.txt", TString optStat = ""){
+
+  printf("================================================================\n");
+  printf("================ fit N-Gauss with estimated  BG ================\n");
+  printf(" * bgEst = larger of bgEst, more linear the estimated BG \n");
+  printf(" * need the file input \n");
+  printf(" \n");
+  printf(" 1) The histogram will be subtracted by the estimated BG. \n");
+  printf(" 2) n-Gauss will then be fitted the BG subtracted histogram \n");
+  printf("================================================================\n");
+  
+
+  bool isParaRead = loadFitParameters(fitFile);
+  if( !isParaRead ) {
+    printf("Please provide a valid input file\n");
+    return;
+  }
+  
+  nPeaks = energy.size();
+  
+  gStyle->SetOptStat(optStat);
+  TCanvas *cFitNGauss = NewCanvas("cFitNGauss","Fit n-Gauss | fitNGauss", 1,4, 800, 300);;
+  cFitNGauss->cd(1);
+
+  ScaleAndDrawHist(hist, 0, 0);
+  
+  TH1F * specS = (TH1F*) hist->Clone();
+  double xMin = hist->GetXaxis()->GetXmin();
+  double xMax = hist->GetXaxis()->GetXmax();
+  int xBin = hist->GetXaxis()->GetNbins();
+  
+  TString titleH;
+  titleH.Form("fitNGauss (BG = %2d); Ex [MeV]; Count / %4.0f keV", bgEst, (xMax-xMin)*1000./xBin );
+  specS->SetTitle(titleH);   
+  specS->SetName("specS");
+  
+  //=================== find peak and fi
+  
+  gStyle->SetOptFit(0);
+  //cFitNGauss->cd(2)->SetGrid();
+  cFitNGauss->cd(2);
+  
+  if( bgEst > 0 ) {
+    printf("============= estimating background...\n");
+    TSpectrum * peak = new TSpectrum(50);
+    TH1 * h1 = peak->Background(hist, bgEst);
+    cFitNGauss->cd(1);
+    h1->Draw("same");
+    cFitNGauss->cd(2);
+    printf("============= substracting the estimated background...\n");
+    specS->Sumw2();
+    specS->Add(h1, -1.);
+  }
+  
+  specS->Draw("hist");
+
+  //========== Fitting 
+  printf("============= Fitting %d-Gauss..... \n", nPeaks);
+  
+  const int  n = 3 * nPeaks;
+  double * para = new double[n]; 
+  for(int i = 0; i < nPeaks ; i++){
+    para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi());
+    para[3*i+1] = energy[i];
+    para[3*i+2] = sigma[i]/2.;
+  }
+
+  TF1 * fit = new TF1("fit", nGauss, xMin, xMax, 3* nPeaks );
+  fit->SetLineWidth(3);
+  fit->SetLineColor(1);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+  
+  //fixing parameters
+  for( int i = 0; i < nPeaks; i++){
+    fit->SetParLimits(3*i  , 0, 1e9); 
+    
+    if( energyFlag[i] == 1 ) {
+      fit->FixParameter(3*i+1, energy[i]);
+    }else{
+      fit->SetParLimits(3*i+1, lowE[i], highE[i]); 
+    }
+    if( sigmaFlag[i]  == 1 ) {
+      fit->FixParameter(3*i+2, sigma[i]);
+    }else{
+      fit->SetParLimits(3*i+2, 0, sigma[i]);
+    }
+  }
+  
+  specS->Fit("fit", "q");
+  
+  const Double_t* paraE = fit->GetParErrors();
+  const Double_t* paraA = fit->GetParameters();
+  
+  //======== calculate reduce chi-squared
+  //GoodnessofFit(specS, fit);
+  
+  double bw = specS->GetBinWidth(1);
+
+  for(int i = 0; i < nPeaks ; i++){
+    printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]);
+  }
+  printf("\n");
+
+  //draw the indivual fit
+  specS->Draw("hist");
+  fit->Draw("same");
+  
+  const int nn = nPeaks;
+  TF1 ** gFit = new TF1 *[nn];
+  for( int i = 0; i < nPeaks; i++){
+      gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax);
+      gFit[i]->SetParameter(0, paraA[3*i]);
+      gFit[i]->SetParameter(1, paraA[3*i+1]);
+      gFit[i]->SetParameter(2, paraA[3*i+2]);
+      gFit[i]->SetLineColor(i+1);
+      gFit[i]->SetNpx(1000);
+      gFit[i]->SetLineWidth(1);
+      gFit[i]->Draw("same");
+  }
+  
+  specS->Draw("hist same");
+  //specS->Draw("E same");
+
+  
+  //======== print text on plot
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.06);
+
+  double chi2 = fit->GetChisquare();
+  int ndf = fit->GetNDF();
+  text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+
+  cFitNGauss->cd(3);
+  PlotResidual(specS, fit);
+  
+  cFitNGauss->cd(4);
+  
+  text.SetTextSize(0.05);
+  text.SetTextColor(2);
+  
+  text.DrawLatex(0.1, 0.9, Form("      %13s, %18s, %18s", "count", "mean", "sigma"));
+  
+  for( int i = 0; i < nPeaks; i++){
+  text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]));
+  }
+  
+  cFitNGauss->Update();
+}
+
+
+//########################################
+//#### not updated
+//########################################
+void fitNGaussSub(TH1F * hist, int bgEst = 10, int degPol = 1, TString fitFile = "AutoFit_para.txt", TString optStat = ""){
+
+
+  printf("==================================================================\n");
+  printf("======== fit N-Gauss with estimated  BG (method-2) ================\n");
+  printf(" * bgEst = larger of bgEst, more linear the estimated BG \n");
+  printf(" * degPol = degree of polynomial \n");
+  printf(" * need the file input \n");
+  printf(" \n");
+  printf(" 1) A BG is estimated, and then the BG is fitted by a polynomial. \n");
+  printf(" 2) The histogram will be subtracted by the polynomial. \n"); 
+  printf(" 3) n-Gauss will then be fitted the subtracted histogram \n");
+  printf("================================================================\n");
+  
+
+  bool isParaRead = loadFitParameters(fitFile);
+  if( !isParaRead ) {
+    printf("Please provide a valid input file\n");
+    return;
+  }
+
+  nPeaks = energy.size();
+  nPols = degPol;
+
+  gStyle->SetOptStat(optStat);
+  TCanvas *cFitNGaussSub = NewCanvas("cFitNGaussSub","Fit n-Gauss, replace estimated BG with Pol-n  | fitNGauss2", 1, 4, 800, 300 );  
+  //if(! cFitNGaussSub->GetShowEventStatus() ) cFitNGaussSub->ToggleEventStatus();
+  
+  cFitNGaussSub->cd(1);
+  ScaleAndDrawHist(hist, 0, 0);
+  
+  TH1F * specS = (TH1F*) hist->Clone();
+  double xMin = hist->GetXaxis()->GetXmin();
+  double xMax = hist->GetXaxis()->GetXmax();
+  int xBin = hist->GetXaxis()->GetNbins();
+  
+  TString titleH;
+  titleH.Form("fitNGauss2 (replace Est. BG with Pol-%d) (BG = %2d); Ex [MeV]; Count / %4.0f keV", degPol, bgEst, (xMax-xMin)*1000./xBin );
+  specS->SetTitle(titleH);   
+  specS->SetName("specS");
+  
+  printf("============= estimating background...\n");
+  TSpectrum * peak = new TSpectrum(50);
+  TH1 * h1 = peak->Background(hist, bgEst);
+  
+  printf("============= fit the est-background with a polynomial function...\n");
+
+  TString polExp = "[0]";
+  for( int i = 1; i < degPol + 1; i++){
+    polExp += Form("+[%d]*TMath::Power(x,%d)", i, i );
+  }
+  TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax);
+  
+  bg->SetParameter(0, 50);
+  bg->SetParameter(0, 0);
+  bg->SetLineColor(2);
+  bg->SetNpx(1000);
+  h1->Fit("bg", "R", "");
+  
+  hist->Draw();
+  bg->Draw("same");
+
+
+  //======== print text on plot
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.04);
+
+  const Double_t * paraAt = bg->GetParameters();
+  const Double_t * paraEt = bg->GetParErrors();
+
+  for( int i = 0; i < degPol + 1; i++){
+    text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraAt[i], paraEt[i]));
+  }
+  
+  gStyle->SetOptFit(0);
+//  cFitNGaussSub->cd(2)->SetGrid();
+  cFitNGaussSub->cd(2);
+
+  printf("============= substracting the polynomial background...\n");
+  specS->Sumw2();
+  specS->Add(bg, -1.);
+  specS->Draw("hist");
+
+  //========== Fitting 
+  printf("============= Fitting..... \n");
+  
+  const int  n = 3 * nPeaks;
+  double * para = new double[n]; 
+  for(int i = 0; i < nPeaks ; i++){
+    para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi());
+    para[3*i+1] = energy[i];
+    para[3*i+2] = sigma[i]/2.;
+  }
+
+  TF1 * fit = new TF1("fit", nGauss, xMin, xMax, 3* nPeaks );
+  fit->SetLineWidth(3);
+  fit->SetLineColor(1);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+  
+  //fixing parameters
+  for( int i = 0; i < nPeaks; i++){
+    fit->SetParLimits(3*i  ,       0, 1e9); 
+    
+    if( energyFlag[i] == 1 ) {
+      fit->FixParameter(3*i+1, energy[i]);
+    }else{
+      fit->SetParLimits(3*i+1, lowE[i], highE[i]); 
+    }
+    if( sigmaFlag[i]  == 1 ) {
+      fit->FixParameter(3*i+2, sigma[i]);
+    }else{
+      fit->SetParLimits(3*i+2, 0, sigma[i]); 
+    }
+  }
+  
+  specS->Fit("fit", "q");
+  
+  const Double_t* paraE = fit->GetParErrors();
+  const Double_t* paraA = fit->GetParameters();
+  
+  GoodnessofFit(specS, fit);  
+  
+  double bw = specS->GetBinWidth(1);
+   
+  for(int i = 0; i < nPeaks ; i++){
+    printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]);
+  }
+  printf("\n");
+  
+
+  //draw the indivual fit
+  specS->Draw("hist");
+  fit->Draw("same");
+  
+  const int nn = nPeaks;
+  TF1 ** gFit = new TF1 *[nn];
+  for( int i = 0; i < nPeaks; i++){
+      gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax);
+      gFit[i]->SetParameter(0, paraA[3*i]);
+      gFit[i]->SetParameter(1, paraA[3*i+1]);
+      gFit[i]->SetParameter(2, paraA[3*i+2]);
+      gFit[i]->SetLineColor(i+1);
+      gFit[i]->SetNpx(1000);
+      gFit[i]->SetLineWidth(1);
+      gFit[i]->Draw("same");
+  }
+  
+  specS->Draw("hist same");
+  //specS->Draw("E same");
+
+  double chi2 = fit->GetChisquare();
+  int ndf = fit->GetNDF();
+  text.SetTextSize(0.06);
+  text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+
+
+  cFitNGaussSub->cd(3);
+  PlotResidual(specS, fit);
+
+  cFitNGaussSub->cd(4);
+
+  text.SetTextSize(0.05);
+  text.SetTextColor(2);
+  
+  text.DrawLatex(0.1, 0.9, Form("      %13s, %18s, %18s", "count", "mean", "sigma"));
+  
+  for( int i = 0; i < nPeaks; i++){
+  text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]));
+  }
+  
+}
+
+
+//########################################
+//#### fit N-Gauss with pol-n BG
+//########################################
+void fitNGaussPol(TH1F * hist, int degPol,  TString fitFile = "AutoFit_para.txt",double xMin = 0, double xMax = 0, TString optStat = ""){
+
+  printf("================================================================\n");
+  printf("================ fit N-Gauss with Pol-%1d BG ==================\n", degPol);
+  printf(" * degPol = degree of polynomial \n");
+  printf(" * need the file input \n");
+  printf(" * xMin, xMax = if left empty, full range will be used\n");
+  printf(" \n");
+  printf(" 1) The histogram will be fitted by n-Gauss + Pol \n");
+  printf("================================================================\n");
+  
+
+  bool isParaRead = loadFitParameters(fitFile);
+  if( !isParaRead ) {
+    printf("Please provide a valid input file\n");
+    return;
+  }
+
+  gStyle->SetOptStat(optStat);
+  nPeaks = energy.size();
+  nPols = degPol;
+
+  TCanvas * cFitNGaussPol = NewCanvas("cFitNGaussPol", Form("Fitting with n-Gauss + pol-%d  | fitNGaussPol", degPol), 1, 3, 800, 300);
+  //if(! cFitNGaussPol->GetShowEventStatus() ) cFitNGaussPol->ToggleEventStatus();
+  cFitNGaussPol->cd(1);
+
+  ScaleAndDrawHist(hist, xMin, xMax);
+  
+  if( xMin == xMax){ 
+    xMin = hist->GetXaxis()->GetXmin();
+    xMax = hist->GetXaxis()->GetXmax();
+  }
+  int xBin = hist->GetXaxis()->GetNbins();
+
+  
+  printf("============= find the polynomial background ..... \n");
+  
+  int  nPar = 3 * nPeaks + nPols + 1;
+  double * para = new double[nPar]; 
+  for(int i = 0; i < nPeaks ; i++){
+    para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi());
+    para[3*i+1] = energy[i];
+    para[3*i+2] = sigma[i]/2.;
+  }
+
+  for(int i = 0; i < nPols + 1; i++){
+    //if( ggg == NULL ){
+      para[3*nPeaks+i] = gRandom->Rndm()/TMath::Power(10, 3*i);
+    //}else{
+    //  para[3*nPeaks+i] = gPara[i];
+    //}
+  }
+  
+  TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar );
+  fit->SetLineWidth(3);
+  fit->SetLineColor(1);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+
+  //fixing parameters
+  for( int i = 0; i < nPeaks; i++){
+    fit->SetParLimits(3*i  ,       0, 1e9); 
+    if( energyFlag[i] == 1 ) {
+      fit->FixParameter(3*i+1, energy[i]);
+    }else{
+      fit->SetParLimits(3*i+1, lowE[i], highE[i]); 
+    }
+    if( sigmaFlag[i]  == 1 ) {
+      fit->FixParameter(3*i+2, sigma[i]);
+    }else{
+      fit->SetParLimits(3*i+2, 0, sigma[i]);
+    }
+  }
+  
+  hist->Fit("fit", "Rq");
+  
+  //=========== get the polynomial part 
+  const Double_t* paraA = fit->GetParameters();
+  const Double_t* paraE = fit->GetParErrors();
+
+  TString polExp = "[0]";
+  for( int i = 1; i < degPol + 1; i++){
+    polExp += Form("+[%d]*TMath::Power(x,%d)", i, i );
+  }
+
+  TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax);
+  for( int i = 0; i < degPol + 1; i++){
+    bg->SetParameter(i, paraA[3*nPeaks+i]);
+  }
+  bg->SetNpx(1000);
+
+  for( int i = 0; i < degPol + 1; i++){
+    printf("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i]);
+  }
+  printf("====================================\n");
+  
+  cFitNGaussPol->cd(1);
+  bg->Draw("same");
+  
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.04);
+  text.SetTextColor(1);
+
+  for( int i = 0; i < degPol + 1; i++){
+    text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i]));
+  }
+
+  double chi2 = fit->GetChisquare();
+  int ndf = fit->GetNDF();
+  text.SetTextSize(0.06);
+  text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+
+  
+  //GoodnessofFit(specS, fit); 
+  
+  double bw = hist->GetBinWidth(1);
+
+  for(int i = 0; i < nPeaks ; i++){
+    printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]);
+  }
+  printf("\n");
+  
+  
+  const int nn = nPeaks;
+  TF1 ** gFit = new TF1 *[nn];
+  for( int i = 0; i < nPeaks; i++){
+      gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax);
+      gFit[i]->SetParameter(0, paraA[3*i]);
+      gFit[i]->SetParameter(1, paraA[3*i+1]);
+      gFit[i]->SetParameter(2, paraA[3*i+2]);
+      gFit[i]->SetLineColor(i+1);
+      gFit[i]->SetNpx(1000);
+      gFit[i]->SetLineWidth(1);
+      gFit[i]->Draw("same");
+  }
+  
+  cFitNGaussPol->Update();
+
+  cFitNGaussPol->cd(2);
+  PlotResidual(hist, fit);
+
+  cFitNGaussPol->cd(3);
+  
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.05);
+  text.SetTextColor(2);
+  
+  text.DrawLatex(0.1, 0.9, Form("      %13s, %18s, %18s", "count", "mean", "sigma"));
+  
+  for( int i = 0; i < nPeaks; i++){
+  text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]));
+  }
+  
+}
+
+
+
+
+//########################################
+//#### fit N-Gauss with pol-n BG
+//########################################
+void fitNGaussPolSub(TH1F * hist, int degPol,  TString fitFile = "AutoFit_para.txt",double xMin = 0, double xMax = 0, TString optStat = ""){
+
+  printf("================================================================\n");
+  printf("========= fit N-Gauss with Pol-%1d BG / 2nd method ============\n", degPol);
+  printf(" * degPol = degree of polynomial \n");
+  printf(" * need the file input \n");
+  printf(" * xMin, xMax = if left empty, full range will be used\n");
+  printf(" \n");
+  printf(" 1) The histogram will be fitted by n-Gauss + Pol -> to get the estimated BG \n");
+  printf(" 2) The histogram will be subtracted by the Pol BG. \n");
+  printf(" 3) n-Gauss will then be fitted the BG subtracted histogram \n");
+  printf("================================================================\n");
+  
+
+  bool isParaRead = loadFitParameters(fitFile);
+  if( !isParaRead ) {
+    printf("Please provide a valid input file\n");
+    return;
+  }
+
+  gStyle->SetOptStat(optStat);
+  nPeaks = energy.size();
+  nPols = degPol;
+
+  TCanvas * cFitNGaussPolSub = NewCanvas("cFitNGaussPolSub", Form("Fitting with n-Gauss + pol-%d (method-2)  | fitGaussPol2", degPol), 1, 4, 800, 300);
+  //if(! cFitNGaussPol->GetShowEventStatus() ) cFitNGaussPol->ToggleEventStatus();
+  cFitNGaussPolSub->cd(1);
+
+  ScaleAndDrawHist(hist, xMin, xMax);
+  
+  if( xMin == xMax){ 
+    xMin = hist->GetXaxis()->GetXmin();
+    xMax = hist->GetXaxis()->GetXmax();
+  }
+  int xBin = hist->GetXaxis()->GetNbins();
+  
+  printf("============= find the polynomial background ..... \n");
+  
+  int  nPar = 3 * nPeaks + nPols + 1;
+  double * para = new double[nPar]; 
+  for(int i = 0; i < nPeaks ; i++){
+    para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi());
+    para[3*i+1] = energy[i];
+    para[3*i+2] = sigma[i]/2.;
+  }
+
+  for(int i = 0; i < nPols + 1; i++){
+    para[3*nPeaks+i] = gRandom->Rndm()/TMath::Power(10, 3*i);
+  }
+  
+  TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar );
+  fit->SetLineWidth(3);
+  fit->SetLineColor(1);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+  
+  //fixing parameters
+  for( int i = 0; i < nPeaks; i++){
+    fit->SetParLimits(3*i  ,       0, 1e9); 
+    if( energyFlag[i] == 1 ) {
+      fit->FixParameter(3*i+1, energy[i]);
+    }else{
+      fit->SetParLimits(3*i+1, lowE[i], highE[i]); 
+    }
+    if( sigmaFlag[i]  == 1 ) {
+      fit->FixParameter(3*i+2, sigma[i]);
+    }else{
+      fit->SetParLimits(3*i+2, 0, sigma[i]);
+    }
+  }
+  
+  hist->Fit("fit", "nq");
+  
+  //=========== get the polynomial part and substract
+  const Double_t* paraAt = fit->GetParameters();
+  const Double_t* paraEt = fit->GetParErrors();
+
+  TString polExp = "[0]";
+  for( int i = 1; i < degPol + 1; i++){
+    polExp += Form("+[%d]*TMath::Power(x,%d)", i, i );
+  }
+
+  TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax);
+  for( int i = 0; i < degPol + 1; i++){
+    bg->SetParameter(i, paraAt[3*nPeaks+i]);
+  }
+  bg->SetNpx(1000);
+
+  for( int i = 0; i < degPol + 1; i++){
+    printf("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraAt[3*nPeaks+i], paraEt[3*nPeaks+i]);
+  }
+  printf("====================================\n");
+  
+  cFitNGaussPolSub->cd(1);
+  bg->Draw("same");
+  
+  TH1F * specS = (TH1F*) hist->Clone();
+  TString titleH;
+  titleH.Form("pol-%d BG Subtracted spectrum (fitNGaussPol-2); Ex [MeV]; Count / %4.0f keV", degPol, (xMax-xMin)*1000./xBin );
+  specS->SetTitle(titleH);   
+  specS->SetName("specS");
+  
+  //=================== find peak and fit
+  
+  gStyle->SetOptFit(0);
+  ///cFitNGaussPol->cd(2)->SetGrid();
+  cFitNGaussPolSub->cd(2);
+  
+  printf("============= substracting the polynomial background...\n");
+  specS->Sumw2();
+  specS->Add(bg, -1.);
+  specS->Draw("hist");
+
+  //======= fit again
+  printf("============= fitting the subtracted spectrum.... \n");
+  
+  nPar = 3* nPeaks;
+  TF1 * fita = new TF1("fita", nGauss, xMin, xMax, nPar );
+  fita->SetLineWidth(3);
+  fita->SetLineColor(1);
+  fita->SetNpx(1000);
+  fita->SetParameters(para);
+  
+  //fixing parameters
+  for( int i = 0; i < nPeaks; i++){
+    fita->SetParLimits(3*i  ,       0, 1e9); 
+    
+    if( energyFlag[i] == 1 ) {
+      fita->FixParameter(3*i+1, energy[i]);
+    }else{
+      fita->SetParLimits(3*i+1, lowE[i], highE[i]); 
+    }
+    if( sigmaFlag[i]  == 1 ) {
+      fita->FixParameter(3*i+2, sigma[i]);
+    }else{
+      fita->SetParLimits(3*i+2, 0, sigma[i]);
+    }
+  }
+  
+  specS->Fit("fita", "q");
+  
+  const Double_t* paraE = fita->GetParErrors();
+  const Double_t* paraA = fita->GetParameters();
+  
+  //GoodnessofFit(specS, fit); 
+  
+  double bw = specS->GetBinWidth(1);
+
+  for(int i = 0; i < nPeaks ; i++){
+    printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]);
+  }
+  printf("\n");
+  
+
+  //draw the indivual fit
+  specS->Draw("hist");
+  fita->Draw("same");
+  
+  const int nn = nPeaks;
+  TF1 ** gFit = new TF1 *[nn];
+  for( int i = 0; i < nPeaks; i++){
+      gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax);
+      gFit[i]->SetParameter(0, paraA[3*i]);
+      gFit[i]->SetParameter(1, paraA[3*i+1]);
+      gFit[i]->SetParameter(2, paraA[3*i+2]);
+      gFit[i]->SetLineColor(i+1);
+      gFit[i]->SetNpx(1000);
+      gFit[i]->SetLineWidth(1);
+      gFit[i]->Draw("same");
+  }
+  
+  specS->Draw("hist same");
+  //specS->Draw("E same");
+
+
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.04);
+  text.SetTextColor(1);
+
+  for( int i = 0; i < degPol + 1; i++){
+    text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i]));
+  }
+
+  double chi2 = fita->GetChisquare();
+  int ndf = fita->GetNDF();
+  text.SetTextSize(0.06);
+  text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+
+
+  cFitNGaussPolSub->cd(3);
+  PlotResidual(specS, fita);
+
+  cFitNGaussPolSub->cd(4);
+  
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.05);
+  text.SetTextColor(2);
+  
+  text.DrawLatex(0.1, 0.9, Form("      %13s, %18s, %18s", "count", "mean", "sigma"));
+  
+  for( int i = 0; i < nPeaks; i++){
+  text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]));
+  }
+  
+}
+
+//#################################################
+//#### fit N-Gauss with pol-n BG using mouse click
+//#################################################
+
+int nClick = 0;
+bool peakFlag = 1;
+vector<double> xPeakList;
+vector<double> yPeakList;
+vector<double> xBGList;
+vector<double> yBGList;
+
+TH1F * tempHist;
+int markerStyle = 23;
+int markerColor = 2;
+int markerBGColor = 4;
+
+void Clicked() {
+  int event = gPad->GetEvent();
+  if (event != 11) return;
+  TObject *select = gPad->GetSelected();
+  if (!select) return;
+  
+  TH1F *h = (TH1F*)select;
+  int px = gPad->GetEventX();
+  double xd = gPad->AbsPixeltoX(px);
+  float x = gPad->PadtoX(xd);
+
+  if( peakFlag ) {
+    xPeakList.push_back(x);
+  }else{
+    xBGList.push_back(x);
+  }
+  int b = tempHist->FindBin(x);
+  double y = tempHist->GetBinContent(b);
+  if( peakFlag ){
+    yPeakList.push_back(y);
+  }else{
+    yBGList.push_back(y);
+  }
+  // add marker in the histogram
+  TMarker * mark = new TMarker(x,y, markerStyle);
+  mark->SetMarkerColor(markerColor);
+  tempHist->GetListOfFunctions()->Add(mark);
+
+  printf("%2d | x : %8.2f , y : %.0f \n", nClick, x, y);
+
+  nClick ++;
+
+}
+
+
+void SaveFitPara(TString fileName = "AutoFit_para.txt"){
+  printf("Save to : %s \n", fileName.Data()); 
+  FILE * file_out;
+  file_out = fopen (fileName.Data(), "w+");
+
+  fprintf(file_out, "# for n-Gauss fit, can use \"#\", or \"//\" to comment out whole line\n");
+  fprintf(file_out, "# peak    low    high   fixed?  sigma_Max    fixed?  hight\n");
+
+  for( int i = 0 ; i < xPeakList.size() ; i++){
+    fprintf(file_out, "%.3f    %.3f    %.3f    0    %.3f     0    %.0f\n",
+                      xPeakList[i],
+                      xPeakList[i] - 5*sigma[i],
+                      xPeakList[i] + 5*sigma[i],
+                      sigma[i],
+                      yPeakList[i]);
+  }
+  fclose(file_out);
+}
+
+void clickFitNGaussPol(TH1F * hist, int degPol, double sigmaMax = 0){
+
+  printf("=========================================================\n");
+  printf("======== fit n-Gauss + Pol-%d BG using mouse click =====\n", degPol );
+  printf("==========================================================\n");
+
+  gStyle->SetOptStat("");
+  gStyle->SetOptFit(0);
+  
+  TCanvas * cClickFitNGaussPol = NULL;
+  if( gROOT->FindObjectAny("cClickFitGaussPol") == NULL ){
+    cClickFitNGaussPol = new TCanvas("cClickFitNGaussPol", Form("fit Gauss & Pol-%d by mouse click | clickFitGaussPol", degPol), 1400, 1200);
+  }else{
+    delete gROOT->FindObjectAny("cClickFitNGaussPol") ;
+    cClickFitNGaussPol = new TCanvas("cClickFitNGaussPol", Form("fit Gauss & Pol-%d by mouse click | clickFitGaussPol", degPol), 1400, 1200);
+  }
+  cClickFitNGaussPol->Divide(1, 4);
+  for(int i = 1; i <= 4 ; i++){
+    cClickFitNGaussPol->cd(i)->SetGrid(0,0);
+  }
+  
+  if(! cClickFitNGaussPol->GetShowEventStatus() ) cClickFitNGaussPol->ToggleEventStatus();
+  cClickFitNGaussPol->cd(1);
+
+
+  hist->GetListOfFunctions()->Clear();
+  ScaleAndDrawHist(hist, 0, 0);
+  TH1F* hspec = (TH1F*) hist->Clone();
+  hspec->Sumw2();
+  cClickFitNGaussPol->Update();
+  cClickFitNGaussPol->Draw();
+
+  TLatex helpMsg;
+  helpMsg.SetNDC();
+  helpMsg.SetTextFont(82);
+  helpMsg.SetTextSize(0.06);
+  helpMsg.SetTextColor(kRed);
+  helpMsg.DrawLatex(0.15, 0.8, "Click for peaks: (Double-click / x / Ctrl to end) ");
+
+  printf("--------- Click the histogram for peaks: (Double-click / x / Ctrl to end) \n");
+  nClick = 0;
+  xPeakList.clear();
+  yPeakList.clear();
+  markerColor = 2;
+  markerStyle = 23;
+  peakFlag = 1;
+  cClickFitNGaussPol->cd(1)->SetCrosshair(1);
+  cClickFitNGaussPol->cd(1)->AddExec("ex", "Clicked()");
+  tempHist = hist;
+  TObject * obj ;
+  do{
+    obj = gPad->WaitPrimitive();
+    if( obj == NULL  ) break;
+  }while( obj != NULL);
+
+  if( degPol >= 0 ){
+    printf("--------- Click the histogram for Background: (Double-click / x / Ctrl to end) \n");
+    printf("                         * when no input, program will estimate  \n");
+
+    cClickFitNGaussPol->cd(1)->Clear();
+    hist->Draw();
+    helpMsg.SetTextColor(markerBGColor);
+    helpMsg.DrawLatex(0.15, 0.8, "Click for BG: (Double-click / x / Ctrl to end) ");
+    helpMsg.DrawLatex(0.15, 0.75, "* when no input, program will estimate");
+    cClickFitNGaussPol->Update();
+    cClickFitNGaussPol->Draw();
+
+    nClick = 0;
+    xBGList.clear();
+    yBGList.clear();
+    markerColor = markerBGColor;
+    markerStyle = 33;
+    peakFlag = 0;
+    cClickFitNGaussPol->cd(1)->AddExec("ex", "Clicked()");
+    do{
+      obj = gPad->WaitPrimitive();
+      if( obj == NULL  ) break;
+    }while( obj != NULL);
+
+  }
+  
+  cClickFitNGaussPol->cd(1)->DeleteExec("ex");
+  cClickFitNGaussPol->cd(1)->SetCrosshair(0);
+  cClickFitNGaussPol->cd(1)->Clear();
+  hist->Draw();
+
+  tempHist = NULL;
+
+  cClickFitNGaussPol->Update();
+  cClickFitNGaussPol->Draw();
+
+  nPols = degPol;
+  double xMin = hspec->GetXaxis()->GetXmin();
+  double xMax = hspec->GetXaxis()->GetXmax();
+
+  TString polExp = "[0]";
+  for( int i = 1; i < degPol + 1; i++) polExp += Form("+[%d]*TMath::Power(x,%d)", i, i );
+  TF1 *bg = new TF1("bg", polExp.Data(), xMin, xMax);
+  bg->SetNpx(1000);
+  bg->SetLineColor(markerBGColor);
+  bg->SetLineWidth(1);
+  if( xBGList.size() > 0 ) {
+    printf("---------------- fit the BG with Pol-%d \n", nPols);
+    TGraph * gBG = new TGraph((int) xBGList.size(), &xBGList[0], &yBGList[0]);
+    for( int i = 0; i < degPol + 1; i++) bg->SetParameter(i, gRandom->Rndm()/TMath::Power(10, 3*i));
+    gBG->Fit("bg", "Rq");
+    bg->Draw("same");
+    //printf("--------------- Subtracting the BG \n");
+    //hspec->Add(bg, -1);
+  }else{
+    for( int i = 0; i < nPols+1; i++) bg->SetParameter(i, 0.);
+  }
+
+  nPeaks = (int) xPeakList.size();
+  if( sigmaMax == 0 ){
+    printf("------------- Estimate sigma for each peak \n");
+    sigma.clear();
+    int binMin = hist->FindBin(xMin);
+    int binMax = hist->FindBin(xMax);
+    for( int i = 0; i < nPeaks ; i++){
+      int b0 = hist->FindBin(xPeakList[i]);
+      double estBG = bg->Eval(xPeakList[i]);    
+      double sMin = (xMax-xMin)/5., sMax = (xMax-xMin)/5.;
+      //---- backward search, stop when 
+      for( int b = b0-1 ; b > binMin ; b-- ){
+        double y = hist->GetBinContent(b);
+        double x = hist->GetBinCenter(b);
+        if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) {
+          sMin = xPeakList[i] - hist->GetBinCenter(b);
+          break;
+        }
+      }
+      //---- forward search, stop when 
+      for( int b = b0+1 ; b < binMax  ; b++ ){
+        double y = hist->GetBinContent(b);
+        double x = hist->GetBinCenter(b);
+        if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) {
+          sMax = hist->GetBinCenter(b) - xPeakList[i];
+          break;
+        }
+      }
+
+      double temp = TMath::Min(sMin, sMax);
+      /// When there are multiple peaks closely packed : 
+      if( i > 0 && temp > 2.5 * sigma.back() ) temp = sigma.back();
+      sigma.push_back(temp);
+      
+      printf("%2d | x : %8.2f | sigma(est) %f \n", i, xPeakList[i], sigma[i]);
+
+    }
+
+    //---- use the mean of the sigma
+    double sigma0 = 0;
+    for( int i = 0; i < nPeaks ; i++) sigma0 += sigma[i];
+    sigma0 = sigma0/(nPeaks+1);
+    for( int i = 0; i < nPeaks ; i++) sigma[i] = sigma0;
+    printf("========== use the mean sigma : %f \n", sigma0);
+  }else if( sigmaMax < 0 ){
+    printf("========== use user input sigma : %f (fixed)\n", abs(sigmaMax));
+    sigma.clear();
+    for( int i = 0; i < nPeaks ; i++) sigma.push_back(abs(sigmaMax));
+  }else if( sigmaMax > 0 ){
+    printf("========== use user input sigma : %f/2. \n", sigmaMax/2.);
+    sigma.clear();
+    for( int i = 0; i < nPeaks ; i++) sigma.push_back(sigmaMax/2.);
+  }
+
+  printf("-------------- Fit the spectrum with %d-Gauss + Pol-%d\n", nPeaks, nPols );
+  cClickFitNGaussPol->cd(2);
+  hspec->Draw("hist");
+
+  int  nPar = 3 * nPeaks + nPols + 1;
+  double * para = new double[nPar]; 
+  for(int i = 0; i < nPeaks ; i++){
+    para[3*i+0] = yPeakList[i] * 0.05 * TMath::Sqrt(TMath::TwoPi());
+    para[3*i+1] = xPeakList[i];
+    if( sigmaMax == 0){
+      para[3*i+2] = sigma[i];  
+    }else if(sigmaMax < 0 ){
+      para[3*i+2] = abs(sigmaMax);
+    }else if(sigmaMax > 0 ){
+      para[3*i+2] = sigmaMax/2.;
+    }
+  }
+  for(int i = 0; i < nPols+1 ; i++){
+    para[3*nPeaks+i] = bg->GetParameter(i);
+  }
+
+  TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar );
+  fit->SetLineWidth(3);
+  fit->SetLineColor(1);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+
+  //limit parameters
+  for( int i = 0; i < nPeaks; i++){
+    fit->SetParLimits(3*i  ,       0, 1e9);
+    
+    double dE1, dE2;
+    if( i == 0 ){
+      dE2 = xPeakList[i+1] - xPeakList[i];
+      dE1 = dE2;
+    }else if ( i == nPeaks-1 ){
+      dE1 = xPeakList[i] - xPeakList[i-1];
+      dE2 = dE1;
+    }else{
+      dE1 = xPeakList[i] - xPeakList[i-1];
+      dE2 = xPeakList[i+1] - xPeakList[i];
+    }
+    
+    fit->SetParLimits(3*i+1, xPeakList[i] - dE1 , xPeakList[i] + dE2 );
+    if( sigmaMax== 0 ) fit->SetParLimits(3*i+2, 0, 1.5*sigma[i]); // add 50% margin of sigma
+    if( sigmaMax < 0 ) fit->FixParameter(3*i+2, abs(sigmaMax));
+    if( sigmaMax > 0 ) fit->SetParLimits(3*i+2, 0, sigmaMax);
+  }
+  
+  hspec->Fit("fit", "Rq");
+  fit->Draw("same");
+
+  //=========== get the polynomial part 
+  const Double_t* paraA = fit->GetParameters();
+  const Double_t* paraE = fit->GetParErrors();
+
+  double chi2 = fit->GetChisquare();
+  int ndf = fit->GetNDF();
+  double bw = hspec->GetBinWidth(1);
+
+  for(int i = 0; i < nPeaks ; i++){
+    printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]);
+  }
+  printf("\n");
+
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.06);
+  text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+
+  for( int i = 0; i < degPol + 1; i++){
+    text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i]));
+  }
+
+  TF1 ** gFit = new TF1 *[nPeaks];
+  for( int i = 0; i < nPeaks; i++){
+      gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax);
+      gFit[i]->SetParameter(0, paraA[3*i]);
+      gFit[i]->SetParameter(1, paraA[3*i+1]);
+      gFit[i]->SetParameter(2, paraA[3*i+2]);
+      gFit[i]->SetLineColor(i+1);
+      gFit[i]->SetNpx(1000);
+      gFit[i]->SetLineWidth(1);
+      gFit[i]->Draw("same");
+  }
+
+  if( degPol >= 0 ){
+    TF1 *bg2 = new TF1("bg", polExp.Data(), xMin, xMax);
+    bg2->SetNpx(1000);
+    bg2->SetLineColor(markerBGColor);
+    bg2->SetLineWidth(1);
+    for( int i = 0; i < nPols + 1; i++){
+      bg2->SetParameter(i, paraA[3*nPeaks+i]);
+    }
+    bg2->Draw("same");
+  }
+
+  cClickFitNGaussPol->cd(3);
+  PlotResidual(hspec, fit);
+  
+  cClickFitNGaussPol->cd(4);
+
+  text.SetTextSize(0.05);
+  text.SetTextColor(2);
+  text.DrawLatex(0.1, 0.9, Form("      %13s, %18s, %18s", "count", "mean", "sigma"));
+  
+  for( int i = 0; i < nPeaks; i++){
+  text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]));
+  }
+}
+
+
+
+void clickFitNGaussPolSub(TH1F * hist, int degPol, double sigmaMax = 0){
+
+  printf("=========================================================\n");
+  printf("= fit n-Gauss + Pol-%d BG using mouse click (method-2) =\n", degPol );
+  printf("==========================================================\n");
+
+  gStyle->SetOptStat("");
+  gStyle->SetOptFit(0);
+  
+  TCanvas * cClickFitNGaussPolsub = NULL;
+  if( gROOT->FindObjectAny("cClickFitGaussPolsub") == NULL ){
+    cClickFitNGaussPolsub = new TCanvas("cClickFitNGaussPol", Form("fit Gauss & Pol-%d by mouse click (sub) | clickFitGaussPolsub", degPol), 1400, 1200);
+  }else{
+    delete gROOT->FindObjectAny("cClickFitNGaussPolsub") ;
+    cClickFitNGaussPolsub = new TCanvas("cClickFitNGaussPolsub", Form("fit Gauss & Pol-%d by mouse click (sub) | clickFitGaussPolsub", degPol), 1400, 1200);
+  }
+  cClickFitNGaussPolsub->Divide(1, 4);
+  for(int i = 1; i <= 4 ; i++){
+    cClickFitNGaussPolsub->cd(i)->SetGrid(0,0);
+  }
+  
+  if(! cClickFitNGaussPolsub->GetShowEventStatus() ) cClickFitNGaussPolsub->ToggleEventStatus();
+  cClickFitNGaussPolsub->cd(1);
+
+
+  hist->GetListOfFunctions()->Clear();
+  ScaleAndDrawHist(hist, 0, 0);
+  TH1F* hspec = (TH1F*) hist->Clone();
+  hspec->Sumw2();
+  cClickFitNGaussPolsub->Update();
+  cClickFitNGaussPolsub->Draw();
+
+  TLatex helpMsg;
+  helpMsg.SetNDC();
+  helpMsg.SetTextFont(82);
+  helpMsg.SetTextSize(0.06);
+  helpMsg.SetTextColor(kRed);
+  helpMsg.DrawLatex(0.15, 0.8, "Click for peaks: (Double-click / x / Ctrl to end) ");
+
+
+  printf("--------- Click the histogram for peaks: (Double-click / x / Ctrl to end) \n");
+  nClick = 0;
+  xPeakList.clear();
+  yPeakList.clear();
+  markerColor = 2;
+  markerStyle = 23;
+  peakFlag = 1;
+  cClickFitNGaussPolsub->cd(1)->SetCrosshair(1);
+  cClickFitNGaussPolsub->cd(1)->AddExec("ex", "Clicked()");
+  tempHist = hist;
+  TObject * obj ;
+  do{
+    obj = gPad->WaitPrimitive();
+    if( obj == NULL  ) break;
+  }while( obj != NULL);
+
+  if( degPol >= 0 ){
+    printf("--------- Click the histogram for Background: (Double-click / x / Ctrl to end) \n");
+
+    cClickFitNGaussPolsub->cd(1)->Clear();
+    hist->Draw();
+    helpMsg.SetTextColor(markerBGColor);
+    helpMsg.DrawLatex(0.15, 0.8, "Click for BG: (Double-click / x / Ctrl to end) ");
+    cClickFitNGaussPolsub->Update();
+    cClickFitNGaussPolsub->Draw();
+
+    nClick = 0;
+    xBGList.clear();
+    yBGList.clear();
+    markerColor = markerBGColor;
+    markerStyle = 33;
+    peakFlag = 0;  
+    do{
+      obj = gPad->WaitPrimitive();
+      if( obj == NULL  ) break;
+    }while( obj != NULL);
+  }
+  cClickFitNGaussPolsub->cd(1)->DeleteExec("ex");
+  cClickFitNGaussPolsub->cd(1)->SetCrosshair(0);
+  cClickFitNGaussPolsub->cd(1)->Clear();
+  hist->Draw();
+
+  tempHist = NULL;
+
+  if( xBGList.size() == 0 ) helpMsg.DrawLatex(0.15, 0.75, " No BG defined ! fitting could be problematics. ");
+
+  cClickFitNGaussPolsub->Update();
+  cClickFitNGaussPolsub->Draw();
+
+  nPols = degPol;
+  double xMin = hspec->GetXaxis()->GetXmin();
+  double xMax = hspec->GetXaxis()->GetXmax();
+
+  TString polExp = "[0]";
+  for( int i = 1; i < degPol + 1; i++) polExp += Form("+[%d]*TMath::Power(x,%d)", i, i );
+  TF1 *bg = new TF1("bg", polExp.Data(), xMin, xMax);
+  bg->SetNpx(1000);
+  bg->SetLineColor(markerBGColor);
+  bg->SetLineWidth(1);
+  if( xBGList.size() > 0 ) {
+    printf("---------------- fit the BG with Pol-%d \n", nPols);
+    
+    TGraph * gBG = new TGraph((int) xBGList.size(), &xBGList[0], &yBGList[0]);
+
+    for( int i = 0; i < degPol + 1; i++) bg->SetParameter(i, gRandom->Rndm()/TMath::Power(10, 3*i));
+
+    gBG->Fit("bg", "Rq");
+    bg->Draw("same");
+
+    printf("--------------- Subtracting the BG \n");
+    hspec->Add(bg, -1);
+  }else{
+    for( int i = 0; i < nPols+1; i++) bg->SetParameter(i, 0.);
+  }
+  
+  
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.06);
+  
+  double chi2BG = bg->GetChisquare();
+  int ndfBG = bg->GetNDF();
+  
+  text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2BG/ndfBG));
+
+  //=========== get the polynomial BG
+  const Double_t* paraAt = bg->GetParameters();
+  const Double_t* paraEt = bg->GetParErrors();
+  
+  for( int i = 0; i < degPol + 1; i++){
+    text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.2f(%8.2f)\n", Form("p%d", i), paraAt[i], paraEt[i]));
+  }
+
+  nPeaks = (int) xPeakList.size();
+  if( sigmaMax == 0 ){
+    printf("------------- Estimate sigma for each peak \n");
+    sigma.clear();
+    int binMin = hist->FindBin(xMin);
+    int binMax = hist->FindBin(xMax);
+    
+    for( int i = 0; i < nPeaks ; i++){
+      int b0 = hist->FindBin(xPeakList[i]);
+      double estBG = bg->Eval(xPeakList[i]);    
+      double sMin = (xMax-xMin)/5., sMax = (xMax-xMin)/5.;
+      //---- backward search, stop when 
+      for( int b = b0-1 ; b > binMin ; b-- ){
+        double y = hist->GetBinContent(b);
+        double x = hist->GetBinCenter(b);
+        if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) {
+          sMin = xPeakList[i] - hist->GetBinCenter(b);
+          break;
+        }
+      }
+      //---- forward search, stop when 
+      for( int b = b0+1 ; b < binMax  ; b++ ){
+        double y = hist->GetBinContent(b);
+        double x = hist->GetBinCenter(b);
+        if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) {
+          sMax = hist->GetBinCenter(b) - xPeakList[i];
+          break;
+        }
+      }
+
+      double temp = TMath::Min(sMin, sMax);
+      /// When there are multiple peaks closely packed : 
+      if( i > 0 && temp > 2.5 * sigma.back() ) temp = sigma.back();
+      sigma.push_back(temp);
+
+      printf("%2d | x : %8.2f | sigma(est) %f \n", i, xPeakList[i], sigma[i]);
+
+    }  
+  }else if( sigmaMax < 0 ){
+    printf("========== use user input sigma : %f (fixed)\n", abs(sigmaMax));
+    sigma.clear();
+    for( int i = 0; i < nPeaks ; i++) sigma.push_back(abs(sigmaMax));
+  }else if( sigmaMax > 0 ){
+    printf("========== use user input sigma : %f/2. \n", sigmaMax/2.);
+    sigma.clear();
+    for( int i = 0; i < nPeaks ; i++) sigma.push_back(sigmaMax/2.);
+  }
+    
+  printf("-------------- Fit the subtracted spectrum with %d-Gauss\n", nPeaks );
+  cClickFitNGaussPolsub->cd(2);
+  hspec->Draw("hist");
+
+  int  nPar = 3 * nPeaks;
+  double * para = new double[nPar]; 
+  for(int i = 0; i < nPeaks ; i++){
+    para[3*i+0] = yPeakList[i] * 0.05 * TMath::Sqrt(TMath::TwoPi());
+    para[3*i+1] = xPeakList[i];
+    if( sigmaMax == 0){
+      para[3*i+2] = sigma[i];  
+    }else if(sigmaMax < 0 ){
+      para[3*i+2] = abs(sigmaMax);
+    }else if(sigmaMax > 0 ){
+      para[3*i+2] = sigmaMax/2.;
+    }
+  }
+
+  TF1 * fit = new TF1("fit", nGauss, xMin, xMax, nPar );
+  fit->SetLineWidth(3);
+  fit->SetLineColor(1);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+
+  //limit parameters
+  for( int i = 0; i < nPeaks; i++){
+    fit->SetParLimits(3*i  ,       0, 1e9);
+
+    double dE1, dE2;
+    if( i == 0 ){
+      dE2 = xPeakList[i+1] - xPeakList[i];
+      dE1 = dE2;
+    }else if ( i == nPeaks-1 ){
+      dE1 = xPeakList[i] - xPeakList[i-1];
+      dE2 = dE1;
+    }else{
+      dE1 = xPeakList[i] - xPeakList[i-1];
+      dE2 = xPeakList[i+1] - xPeakList[i];
+    }
+    
+    fit->SetParLimits(3*i+1, xPeakList[i] - dE1 , xPeakList[i] + dE2 ); 
+    if( sigmaMax== 0 ) fit->SetParLimits(3*i+2, 0, 1.5*sigma[i]); // add 50% margin of sigma
+    if( sigmaMax < 0 ) fit->FixParameter(3*i+2, abs(sigmaMax));
+    if( sigmaMax > 0 ) fit->SetParLimits(3*i+2, 0, sigmaMax);
+  }
+  
+  hspec->Fit("fit", "Rq");
+  fit->Draw("same");
+
+  //=========== get the fit parameters 
+  const Double_t* paraA = fit->GetParameters();
+  const Double_t* paraE = fit->GetParErrors();
+
+  double chi2 = fit->GetChisquare();
+  int ndf = fit->GetNDF();
+  double bw = hspec->GetBinWidth(1);
+  
+  text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+
+  for(int i = 0; i < nPeaks ; i++){
+    printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]);
+  }
+  printf("\n");
+
+  TF1 ** gFit = new TF1 *[nPeaks];
+  for( int i = 0; i < nPeaks; i++){
+      gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax);
+      gFit[i]->SetParameter(0, paraA[3*i]);
+      gFit[i]->SetParameter(1, paraA[3*i+1]);
+      gFit[i]->SetParameter(2, paraA[3*i+2]);
+      gFit[i]->SetLineColor(i+1);
+      gFit[i]->SetNpx(1000);
+      gFit[i]->SetLineWidth(1);
+      gFit[i]->Draw("same");
+  }
+
+  cClickFitNGaussPolsub->cd(3);
+  PlotResidual(hspec, fit);
+  
+  cClickFitNGaussPolsub->cd(4);
+
+  text.SetTextSize(0.05);
+  text.SetTextColor(2);
+  text.DrawLatex(0.1, 0.9, Form("      %13s, %18s, %18s", "count", "mean", "sigma"));
+  
+  for( int i = 0; i < nPeaks; i++){
+  text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]));
+  }
+}
+
+
+//########################################
+//#### fit N-Gauss with pol-1 BG at fixed for < 0
+//########################################
+void fitSpecial(TH1F * hist, TString fitFile = "AutoFit_para.txt"){
+
+  printf("================================================================\n");
+  printf("================ Special fit for h074_82Kr =====================\n");
+  printf(" * need the file input \n");
+  printf("================================================================\n");
+  
+
+  bool isParaRead = loadFitParameters(fitFile);
+  if( !isParaRead ) {
+    printf("Please provide a valid input file\n");
+    return;
+  }
+
+  gStyle->SetOptStat("");
+  gStyle->SetOptFit(0);
+  nPeaks = energy.size();
+  nPols = 1;
+
+  TCanvas * cFitSpecial = NewCanvas("cFitSpecial", "Fitting for h074_82Kr", 1, 3, 800, 300);
+  //if(! cFitSpecial->GetShowEventStatus() ) cFitSpecial->ToggleEventStatus();
+  cFitSpecial->cd(1);
+
+  ScaleAndDrawHist(hist, 0, 0);
+  
+  double  xMin = hist->GetXaxis()->GetXmin();
+  double  xMax = hist->GetXaxis()->GetXmax();
+  int     xBin = hist->GetXaxis()->GetNbins();
+
+  printf("============= find the pol-1 background ..... \n");
+  
+  TF1 * fitpol1 = new TF1("fitpol1", "pol1", xMin, -0.3);
+  fitpol1->SetParameter(0, gRandom->Rndm());
+  fitpol1->SetParameter(1, gRandom->Rndm());
+  
+  hist->Fit("fitpol1", "Rq");
+  
+  double x0 = fitpol1->GetParameter(0);
+  double x1 = fitpol1->GetParameter(1);
+  
+  
+  int  nPar = 3 * nPeaks + nPols + 1;
+  double * para = new double[nPar]; 
+  for(int i = 0; i < nPeaks ; i++){
+    para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi());
+    para[3*i+1] = energy[i];
+    para[3*i+2] = sigma[i]/2.;
+  }
+
+  para[3*nPeaks+0] = x0;
+  para[3*nPeaks+1] = x1;
+  
+  
+  TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar );
+  fit->SetLineWidth(3);
+  fit->SetLineColor(1);
+  fit->SetNpx(1000);
+  fit->SetParameters(para);
+
+  //fixing parameters
+  for( int i = 0; i < nPeaks; i++){
+    fit->SetParLimits(3*i  ,       0, 1e9); 
+    if( energyFlag[i] == 1 ) {
+      fit->FixParameter(3*i+1, energy[i]);
+    }else{
+      fit->SetParLimits(3*i+1, lowE[i], highE[i]); 
+    }
+    if( sigmaFlag[i]  == 1 ) {
+      fit->FixParameter(3*i+2, sigma[i]);
+    }else{
+      fit->SetParLimits(3*i+2, 0, sigma[i]);
+    }
+  }
+  
+  fit->FixParameter(3*nPeaks + 0, x0);
+  fit->FixParameter(3*nPeaks + 1, x1);
+  
+  hist->Fit("fit", "Rq");
+  
+  //=========== get the polynomial part 
+  const Double_t* paraA = fit->GetParameters();
+  const Double_t* paraE = fit->GetParErrors();
+
+  TString polExp = "[0]";
+  for( int i = 1; i < nPols + 1; i++){
+    polExp += Form("+[%d]*TMath::Power(x,%d)", i, i );
+  }
+
+  TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax);
+  for( int i = 0; i < nPols + 1; i++){
+    bg->SetParameter(i, paraA[3*nPeaks+i]);
+  }
+  bg->SetNpx(1000);
+
+  for( int i = 0; i < nPols + 1; i++){
+    printf("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i]);
+  }
+  printf("====================================\n");
+  
+  cFitSpecial->cd(1);
+  bg->Draw("same");
+  
+  TLatex text;
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.04);
+  text.SetTextColor(1);
+
+  for( int i = 0; i < nPols + 1; i++){
+    text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i]));
+  }
+
+  double chi2 = fit->GetChisquare();
+  int ndf = fit->GetNDF();
+  text.SetTextSize(0.06);
+  text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
+
+  
+  //GoodnessofFit(specS, fit); 
+  
+  double bw = hist->GetBinWidth(1);
+
+  for(int i = 0; i < nPeaks ; i++){
+    printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]);
+  }
+  printf("\n");
+  
+  
+  const int nn = nPeaks;
+  TF1 ** gFit = new TF1 *[nn];
+  for( int i = 0; i < nPeaks; i++){
+      gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax);
+      gFit[i]->SetParameter(0, paraA[3*i]);
+      gFit[i]->SetParameter(1, paraA[3*i+1]);
+      gFit[i]->SetParameter(2, paraA[3*i+2]);
+      gFit[i]->SetLineColor(i+1);
+      gFit[i]->SetNpx(1000);
+      gFit[i]->SetLineWidth(1);
+      gFit[i]->Draw("same");
+  }
+  
+  cFitSpecial->Update();
+
+  cFitSpecial->cd(2);
+  PlotResidual(hist, fit);
+
+  cFitSpecial->cd(3);
+  
+  text.SetNDC();
+  text.SetTextFont(82);
+  text.SetTextSize(0.05);
+  text.SetTextColor(2);
+  
+  text.DrawLatex(0.1, 0.9, Form("      %13s, %18s, %18s", "count", "mean", "sigma"));
+  
+  for( int i = 0; i < nPeaks; i++){
+  text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", 
+            i, 
+            paraA[3*i] / bw,   paraE[3*i] /bw, 
+            paraA[3*i+1], paraE[3*i+1],
+            paraA[3*i+2], paraE[3*i+2]));
+  }
+  
+}
+
+
+
+
+#endif 
diff --git a/makefile b/makefile
index 3e260d4..8a44af3 100644
--- a/makefile
+++ b/makefile
@@ -1,10 +1,15 @@
 CC=g++
 
-all: xia2root
-#all: xia2root  xia2ev2
+all: xia2root  xia2ev2_nopart pixie2root scan
 
 xia2root: xia2root.cpp
 	$(CC) xia2root.cpp -o xia2root `root-config --cflags --glibs`
 
-xia2ev2: xia2ev2_nopart.cpp
+xia2ev2_nopart: xia2ev2_nopart.cpp
 	$(CC) xia2ev2_nopart.cpp -o xia2ev2_nopart
+
+pixie2root: pixie2root.cpp
+	$(CC) pixie2root.cpp -o pixie2root `root-config --cflags --glibs`
+
+scan: scan.c
+	$(CC) scan.c -o scan
diff --git a/mapping.h b/mapping.h
new file mode 100644
index 0000000..9dfae75
--- /dev/null
+++ b/mapping.h
@@ -0,0 +1,28 @@
+/************************************
+Clover :   0 -  99
+BGO    : 100 - 199 
+Other  : 200 - 299
+
+ * *********************************/
+
+//==================== mapping
+
+#define NCLOVER 36
+#define NBGO     9
+#define NOTHER  62
+
+//                 0    1    2    3    4    5    6    7    8    9
+int map[130] = {   0,   1,   2,   3, 100,   4,   5,   6,   7, 101,   //   0
+                   8,   9,  10,  11, 102,  -1,  12,  13,  14,  15,   //  10
+                 103,  16,  17,  18,  19, 104,  -1,  -1,  -1,  -1,   //  20
+                  -1,  -1,  20,  21,  22,  23, 105,  24,  25,  26,   //  30
+                  27, 106,  28,  29,  30,  31, 107,  -1,  -1,  -1,   //  40
+                  -1,  -1,  -1,  32,  33,  34,  35, 108,  -1,  -1,   //  50
+                 200, 201, 202, 203, 204, 205, 206, 207, 208, 209,   //  60
+                 210, 211, 212, 213, 214, 215, 216, 217, 218, 219,   //  70
+                 220, 221, 222, 223, 224, 225, 226, 227, 228, 229,   //  80
+                 230, 231, 232, 233, 234, 235, 236, 237, 238, 239,   //  90
+                 240, 241, 242, 243, 244, 245, 246, 247, 248, 249,   // 100
+                 250, 251, 252, 253, 254, 255, 256, 257, 258, 259,   // 110
+                 260, 261,  -1,  -1,  -1,  -1,  -1,  -1,  -1,  -1};   // 120
+                 
diff --git a/pixie2root.cpp b/pixie2root.cpp
new file mode 100644
index 0000000..7f618e1
--- /dev/null
+++ b/pixie2root.cpp
@@ -0,0 +1,474 @@
+/**********************************************************/
+/* PXI SCAN CODE -- J.M. Allmond (ORNL) -- July 2016      */
+/*                                                        */
+/* !unpak data from Pixie-16 digitizers, event build,     */
+/* !and create detctors and user defined spectra          */
+/*                                                        */
+/* gcc -o pxi-scan pxi-scan.c                             */
+/* ./pxi-scan -op datafile calibrationfile mapfile        */
+/*                                                        */
+/* ..... calibration file optional                        */
+/* ..... map file optional                                */
+/* ..... u for update spectra                             */
+/* ..... o for overwrite spectra                          */
+/* ..... p for print realtime stats                       */
+/**********************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <stdbool.h>
+
+#include "TFile.h"
+#include "TTree.h"
+#include "TMath.h"
+#include "TBenchmark.h"
+
+#define PRINT_CAL 1
+#define PRINT_MAP 1
+
+#define RAND ((float) rand() / ((unsigned int) RAND_MAX + 1))   // random number in interval (0,1)
+#define TRUE  1
+#define FALSE 0
+
+#define LINE_LENGTH 120
+
+#define MAX_CRATES 2
+#define MAX_BOARDS_PER_CRATE 13
+#define MAX_CHANNELS_PER_BOARD 16
+#define BOARD_START 2
+
+#define MAX_ID MAX_CRATES*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD
+
+#define HEADER_LENGTH 4     //unit = words with 4 bytes per word
+#define MAX_SUB_LENGTH 2016 //unit = words with 4 bytes per word ; 2004 --> 40 micro-second trace + 4 word header
+
+#define EVENT_BUILD_TIME 109 // 100 = 1 micro-second ; should be < L + G ~ 5.04 us (note 0.08 us scale factor in set file)
+
+#define RAWE_REBIN_FACTOR 2.0 // Rebin 32k pixie16 spectra to something smaller to fit better into 8k.
+
+#include "mapping.h"  
+
+/////////////////////
+// RAW EVENT TYPES //
+/////////////////////
+struct subevent
+{
+    int chn; 
+    int sln;
+    int crn;
+    int id;
+    int hlen;
+    int elen;
+    int trlen;          //number of samples
+    int trwlen;         //number of words (two samples per word)
+    int fcode;          //pileup flag
+    long long int time;
+    int ctime;
+    int ctimef;
+    int energy;
+    int extra;
+    short tr[4096];
+    int esum[4];
+    int qsum[8];        
+}; 
+struct subevent subevt[MAX_ID]={0};
+int sevtmult=0;
+unsigned long long int sevtcount=0;
+unsigned long long int pileupcount=0;
+unsigned long long int evtcount=0;
+
+int mult[1][4096]={0};
+
+int tdifid[MAX_ID][8192]={0};
+
+///////////////////////
+// Write 2-byte data //
+///////////////////////
+void write_data2(char *filename, short *data, int xdim, int ydim, int overwrite) { //2byte per channel Write / Add to previous
+    
+    FILE *FP;
+    int i;
+    short *previous;
+    
+    if(!overwrite) {
+        //allocate memory for 1d-array for reading in rows of 2d Radware matrix
+        if ( ( previous = (short *)malloc(xdim * ydim * sizeof(short)) ) == NULL ) {
+            printf("\nError, memory not allocated.\n");
+            exit(1);
+        }
+      
+        //open previous spectra file  
+        if( (FP=fopen(filename, "r")) != NULL ){
+            fread(previous, sizeof(short)*xdim*ydim, 1, FP);        
+            fclose(FP);
+            //update spectra
+            for (i=0; i<xdim*ydim; i++) {
+                if(previous[i] < (powf(2,sizeof(short)*8.0)-2)) 
+                    data[i] = data[i] + previous[i];
+                }   
+            }
+        else{
+            printf("%s did not previously exist, creating ...\n", filename);       
+        }     
+   
+        //Deallocate previous data
+        free(previous);
+    }
+  
+    FP=fopen(filename, "w");
+    fwrite(data, sizeof(short)*xdim, ydim, FP);
+    fclose(FP);            
+}    
+
+
+///////////////////////
+// Write 4-byte data //
+///////////////////////
+void write_data4(char *filename, int *data, int xdim, int ydim, int overwrite) { //4byte per channel Write / Add to previous
+    
+    FILE *FP;
+    int i;
+    int *previous;
+   
+    if(!overwrite) {
+        //allocate memory for 1d-array for reading in rows of 2d Radware matrix
+        if ( ( previous = (int *)malloc(xdim * ydim * sizeof(int)) ) == NULL ) {
+            printf("\nError, memory not allocated.\n");
+            exit(1);
+        }
+      
+        //open previous spectra file  
+        if( (FP=fopen(filename, "r")) != NULL ){
+            fread(previous, sizeof(int)*xdim*ydim, 1, FP);        
+            fclose(FP);
+            //update spectra
+            for (i=0; i<xdim*ydim; i++) {
+                if(previous[i] < (powf(2,sizeof(int)*8.0)-2)) 
+                    data[i] = data[i] + previous[i];
+            }   
+        } 
+        else{
+            printf("%s did not previously exist, creating ...\n", filename);       
+        }   
+       
+        //Deallocate previous data
+        free(previous);
+    }
+  
+    FP=fopen(filename, "w");
+    fwrite(data, sizeof(int)*xdim, ydim, FP);
+    fclose(FP);            
+}    
+
+
+///////////////////////////////////
+// START OF MAIN FUNCTION        //
+///////////////////////////////////
+int main(int argc, char **argv) {
+  
+  int i=0, j=0, k=0;
+  float tempf=0;
+  int max1=0, min1=0;
+  int max2=0, min2=0;
+  int maxid1=-1, minid1=-1;
+  int maxid2=-1, minid2=-1;
+  div_t e_div;
+  lldiv_t lle_div;
+
+  int overwrite = 1;    
+
+  double etrace0,etrace1,btrace0,btrace1;
+  double ptrace0,ptrace1,ttrace0,ttrace1,tautrace0,tautrace1;
+  int dbcount = 0;
+  long long int strace[500];
+  memset(strace, 0, sizeof(strace));
+
+  //temp buffer for each sub event
+  unsigned int sub[MAX_SUB_LENGTH];
+  memset(sub, 0, sizeof(sub));
+  
+  //Reference time and difference for event building
+  long long int etime, tdif, idtime[MAX_ID]={0}, temptime;
+  
+  // Check that the corrent number of arguments were provided.
+  if (argc != 2 && argc != 3 )    {
+    printf("Incorrect number of arguments:\n");
+    printf("%s datafile <outFile>\n", argv[0]);
+    return 1;
+  }
+
+  printf("=====================================\n");
+  printf("===       evt --> root            ===\n");
+  printf("=====================================\n");
+  
+  //CERN ROOT things
+  TString inFileName = argv[1];
+  TString outFileName = inFileName;
+  
+  if( argc >= 3 ){
+    outFileName = argv[2];
+  }else{
+    outFileName.Remove(inFileName.First('.'));
+    outFileName.Append(".root");
+  }
+  printf("  in file : %s \n", inFileName.Data());
+  printf(" our file : %s \n", outFileName.Data());
+
+  printf(" number of detector channal: %d \n", MAX_ID);
+
+  TFile * outRootFile = new TFile(outFileName, "recreate");
+  outRootFile->cd();
+  TTree * tree = new TTree("tree", "tree");
+
+  unsigned long long evID = -1;
+  double  energy[NCLOVER];  
+  unsigned long long etimestamp[NCLOVER];
+  double  bgo[NBGO]; 
+  double  other[NOTHER];
+  unsigned short pileup[NCLOVER];
+  
+  //const int maxMulti = 40;
+  //double energy[maxMulti];  
+  //unsigned timestamp[maxMulti];  
+  //short detID[maxMulti];
+  int multi;
+
+  tree->Branch("evID", &evID, "event_ID/l"); 
+  ///tree->Branch("detID", detID, Form("det ID[%d]/B", NCLOVER));
+  tree->Branch("e", energy, Form("energy[%d]/D", NCLOVER));
+  tree->Branch("t", etimestamp, Form("energy_time_stamp[%d]/l", NCLOVER));
+  tree->Branch("p", pileup, Form("pile_up_flag[%d]/s", NCLOVER));
+  
+  tree->Branch("bgo", bgo, Form("BGO_energy[%d]/D", NBGO));
+  tree->Branch("other", other, Form("other_energy[%d]/D", NOTHER));
+  
+  tree->Branch("multi", &multi, "multiplicity/I");
+  
+  //open list-mode data file from PXI digitizer  
+  FILE *fpr;
+  long int fprsize,fprpos;
+  if ((fpr = fopen(argv[1], "r")) == NULL) {
+    fprintf(stderr, "Error, cannot open input file %s\n", argv[2]);
+    return 1;
+  }
+ 
+  //get file size
+  fseek(fpr, 0L, SEEK_END);
+  fprsize = ftell(fpr);
+  rewind(fpr);
+  
+  TBenchmark gClock;
+  gClock.Reset();
+  gClock.Start("timer");
+  
+  /////////////////////
+  // MAIN WHILE LOOP //
+  /////////////////////
+  while (1) { //main while loop 
+
+      /////////////////////////////////
+      // UNPACK DATA AND EVENT BUILD //
+      /////////////////////////////////
+      
+      //CERN data clear
+      for( int haha = 0; haha < NCLOVER; haha++){
+         energy[haha] = TMath::QuietNaN();
+         etimestamp[haha] = 0;
+         pileup[haha] = 0;
+      }
+      for( int haha = 0; haha < NBGO; haha++) bgo[haha] = TMath::QuietNaN();
+      for( int haha = 0; haha < NOTHER; haha++) other[haha] = TMath::QuietNaN();
+      multi = 0;
+      evID++;	
+      
+      etime=-1; tdif=-1; sevtmult=0;  
+      //memset(&subevt, 0, sizeof(subevt)); //not needed since everything is redefined (except maybe trace on pileup evts)
+      while (1) { //get subevents and event build for one "event" 
+        
+       // memset(&subevt[sevtmult], 0, sizeof(subevt[sevtmult])); //not needed since everything is redefined (except maybe trace on pileup evts)
+        
+        //read 4-byte header
+        if (fread(sub, sizeof(int)*HEADER_LENGTH, 1, fpr) != 1) break;
+        subevt[sevtmult].chn = sub[0] & 0xF;     /// channel in digitizer
+        subevt[sevtmult].sln = (sub[0] & 0xF0) >> 4; /// digitizer ID
+        subevt[sevtmult].crn = (sub[0] & 0xF00) >> 8;  /// crate
+        subevt[sevtmult].id = subevt[sevtmult].crn*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (subevt[sevtmult].sln-BOARD_START)*MAX_CHANNELS_PER_BOARD + subevt[sevtmult].chn;   
+        subevt[sevtmult].hlen = (sub[0] & 0x1F000) >> 12;
+        subevt[sevtmult].elen = (sub[0] & 0x7FFE0000) >> 17;
+        subevt[sevtmult].fcode = (sub[0] & 0x80000000) >> 31;
+        subevt[sevtmult].time = ( (long long int)(sub[2] & 0xFFFF) << 32) + sub[1];
+        subevt[sevtmult].ctime = (sub[2] & 0x7FFF0000) >> 16;
+        subevt[sevtmult].ctimef = (sub[2] & 0x80000000) >> 31;
+        subevt[sevtmult].energy = (sub[3] & 0xFFFF);
+        subevt[sevtmult].trlen = (sub[3] & 0x7FFF0000) >> 16;
+        subevt[sevtmult].trwlen = subevt[sevtmult].trlen / 2;
+        subevt[sevtmult].extra = (sub[3] & 0x80000000) >> 31;       
+ 
+        //rebin raw trap energy from 32k to ....            
+        tempf = (float)subevt[sevtmult].energy/RAWE_REBIN_FACTOR;// + RAND;       
+        subevt[sevtmult].energy = (int)tempf;                        
+ 
+        //check lengths (sometimes all of the bits for trace length are turned on ...)
+       /* if (subevt[sevtmult].elen - subevt[sevtmult].hlen != subevt[sevtmult].trwlen) {
+            printf("SEVERE ERROR: event, header, and trace length inconsistencies found\n");
+            printf("event length = %d\n", subevt[sevtmult].elen);
+            printf("header length = %d\n", subevt[sevtmult].hlen);
+            printf("trace length = %d\n", subevt[sevtmult].trwlen);  
+            printf("Extra = %d\n", subevt[sevtmult].extra); 
+            printf("fcode = %d\n", subevt[sevtmult].fcode);              
+            //sleep(1);          
+            //return 0;
+        } */ 
+        
+        
+        //CERN fill tree
+        
+        ///========== need a mapping, can reduce the array size, speed up.
+        
+        int ch = map[subevt[sevtmult].id];
+        if ( 0 <= ch && ch < NCLOVER ){
+            energy[ch] = subevt[sevtmult].energy;
+            etimestamp[ch] = subevt[sevtmult].time;
+            pileup[ch] = subevt[sevtmult].fcode;
+            multi++;
+        }
+        if ( 100 <= ch && ch < 100 + NBGO ){
+            bgo[ch-100] = subevt[sevtmult].energy;
+        }
+        if ( 200 <= ch && ch < 200 + NOTHER){
+            other[ch-200] = subevt[sevtmult].energy;
+        }
+       
+        //Set reference time for event building
+        if (etime == -1) {
+            etime = subevt[sevtmult].time;
+            tdif = 0;
+        }
+        else {
+            tdif = subevt[sevtmult].time - etime;
+            if (tdif < 0) {
+                printf("SEVERE ERROR: tdiff < 0, file must be time sorted\n");
+                printf("etime = %lld, time = %lld, and tdif = %lld\n", etime, subevt[sevtmult].time, tdif);                
+                return 0;   
+            }    
+        }    
+      
+        //Check for end of event, rewind, and break out of while loop
+        if (tdif > EVENT_BUILD_TIME) {
+            fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR); //fwrite/fread is buffered by system ; storing this in local buffer is no faster!
+            break;           
+        }    
+               
+        
+        //time between sequential events for a single channel ; useful for determining optimal event building time
+        temptime = (subevt[sevtmult].time - idtime[subevt[sevtmult].id])/100; //rebin to 1 micro-second
+        if ( temptime >= 0 && temptime < 8192) {
+            tdifid[subevt[sevtmult].id][temptime]++;
+        }    
+        idtime[subevt[sevtmult].id]=subevt[sevtmult].time; //store time for next subevent of channel    
+    
+        // total pileups
+        if (subevt[sevtmult].fcode==1) {
+            pileupcount++;
+        }
+        
+        
+        //more data than just the header; read entire sub event
+        fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR);
+        if (fread(sub, sizeof(int)*subevt[sevtmult].elen, 1, fpr) != 1) break;
+                              
+        /*                      
+        //trace
+        k=0;
+        for (i = subevt[sevtmult].hlen; i < subevt[sevtmult].elen; i++) {      
+            subevt[sevtmult].tr[i - subevt[sevtmult].hlen + k] = sub[i] & 0x3FFF;
+            subevt[sevtmult].tr[i - subevt[sevtmult].hlen + k + 1] = (sub[i]>>16) & 0x3FFF;
+            k=k+1;
+        } 
+        
+     // if (subevt[sevtmult].id == 4 && subevt[sevtmult].fcode == 1) DB(subevt[sevtmult].tr);
+            
+        //continue if no esum or qsum   
+        if (subevt[sevtmult].hlen==HEADER_LENGTH) {
+            sevtmult++;        
+            continue;
+        }
+        
+        //esum
+        if (subevt[sevtmult].hlen==8 || subevt[sevtmult].hlen==16) { 
+            for (i=4; i < 8; i++) {
+                subevt[sevtmult].esum[i-4] = sub[i];
+            }
+        }
+    
+        //qsum
+        if (subevt[sevtmult].hlen==12) { 
+            for (i=4; i < 12; i++) {
+                subevt[sevtmult].qsum[i-4] = sub[i];
+            }
+        }
+    
+        //qsum
+        if (subevt[sevtmult].hlen==16) { 
+            for (i=8; i < 16; i++) {
+                subevt[sevtmult].qsum[i-8] = sub[i];
+            }
+        }    
+        */
+        sevtmult++;
+     
+      } //end while loop for unpacking sub events and event building for one "event"
+      if (sevtmult==0) break; //end main WHILE LOOP when out of events 
+      mult[0][sevtmult]++; //Histogram raw sub event multiplicity 
+      sevtcount += sevtmult;
+      evtcount++; //event-built number
+      /////////////////////////////////////
+      // END UNPACK DATA AND EVENT BUILD //
+      /////////////////////////////////////
+
+      //event stats, print status every 10000 events
+      lle_div=lldiv(evtcount,10000);
+      if ( lle_div.rem == 0 ) {
+        fprpos = ftell(fpr);
+        tempf = (float)fprsize/(1024.*1024.*1024.);
+        gClock.Stop("timer");
+        double time = gClock.GetRealTime("timer");
+        gClock.Start("timer");
+        printf("Total SubEvents: \x1B[32m%llu \x1B[31m(%d%% pileup)\x1B[0m\nTotal Events: \x1B[32m%llu (%.1f <mult>)\x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[3A\r", 
+                   sevtcount, (int)((100*pileupcount)/sevtcount), evtcount, (float)sevtcount/(float)evtcount, (100*fprpos/fprsize), tempf,  TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
+      }      
+
+
+      //cern fill tree
+      outRootFile->cd();
+      tree->Fill();
+
+          
+  } // end main while loop 
+  /////////////////////////
+  // END MAIN WHILE LOOP //
+  /////////////////////////
+  fprpos = ftell(fpr);
+  tempf = (float)fprsize/(1024.*1024.*1024.);
+  printf("Total SubEvents: \x1B[32m%llu \x1B[31m(%d%% pileup)\x1B[0m\nTotal Events: \x1B[32m%llu (%.1f <mult>)\x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\n\033[3A\n", 
+  sevtcount, (int)((100*pileupcount)/sevtcount), evtcount, (float)sevtcount/(float)evtcount, (100*fprpos/fprsize), tempf);
+           
+  //cern save root
+  outRootFile->cd();
+  tree->Write();
+  outRootFile->Close();
+  
+  fclose(fpr);
+  
+  gClock.Stop("timer");
+  double time = gClock.GetRealTime("timer");
+  printf("\n==================== finished.\r\n");
+  printf("Total time spend : %3.0f min %5.2f sec\n", TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
+  
+  
+  return 0;
+}
+
+
diff --git a/scan.c b/scan.c
new file mode 100644
index 0000000..5e2532a
--- /dev/null
+++ b/scan.c
@@ -0,0 +1,1208 @@
+/**********************************************************/
+/* PXI SCAN CODE -- J.M. Allmond (ORNL) -- July 2016      */
+/*                                                        */
+/* !unpak data from Pixie-16 digitizers, event build,     */
+/* !and create detctors and user defined spectra          */
+/*                                                        */
+/* gcc -o pxi-scan pxi-scan.c                             */
+/* ./pxi-scan -op datafile calibrationfile mapfile        */
+/*                                                        */
+/* ..... calibration file optional                        */
+/* ..... map file optional                                */
+/* ..... u for update spectra                             */
+/* ..... o for overwrite spectra                          */
+/* ..... p for print realtime stats                       */
+/**********************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <stdbool.h>
+
+#define PRINT_CAL 1
+#define PRINT_MAP 1
+
+#define DB(x) fwrite(x, sizeof(x), 1, debugfile);
+#define DEBUGFN "debug.mat"
+
+#define RAND ((float) rand() / ((unsigned int) RAND_MAX + 1))   // random number in interval (0,1)
+#define TRUE  1
+#define FALSE 0
+
+#define LINE_LENGTH 120
+
+#define MAX_CRATES 2
+#define MAX_BOARDS_PER_CRATE 13
+#define MAX_CHANNELS_PER_BOARD 16
+#define BOARD_START 2
+
+#define MAX_ID MAX_CRATES*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD
+
+#define HEADER_LENGTH 4     //unit = words with 4 bytes per word
+#define MAX_SUB_LENGTH 2016 //unit = words with 4 bytes per word ; 2004 --> 40 micro-second trace + 4 word header
+
+#define EVENT_BUILD_TIME 109 // 100 = 1 micro-second ; should be < L + G ~ 5.04 us (note 0.08 us scale factor in set file)
+
+#define RAWE_REBIN_FACTOR 2.0 // Rebin 32k pixie16 spectra to something smaller to fit better into 8k.
+
+/////////////////////
+// RAW EVENT TYPES //
+/////////////////////
+struct subevent
+{
+    int chn; 
+    int sln;
+    int crn;
+    int id;
+    int hlen;
+    int elen;
+    int trlen;          //number of samples
+    int trwlen;         //number of words (two samples per word)
+    int fcode;          //pileup flag
+    long long int time;
+    int ctime;
+    int ctimef;
+    int energy;
+    int extra;
+    short tr[4096];
+    int esum[4];
+    int qsum[8];        
+}; 
+struct subevent subevt[MAX_ID]={0};
+int sevtmult=0;
+unsigned long long int sevtcount=0;
+unsigned long long int pileupcount=0;
+unsigned long long int evtcount=0;
+
+
+
+//////////////////////////////////////////
+// INPUT CALIBRATION AND MAP PARAMETERS //
+//////////////////////////////////////////
+
+float ecal[MAX_ID][2]={0};
+float tcal[MAX_ID][2]={0};
+float new_gain=1.0;
+
+char map2type[MAX_ID]={0};
+int map2det[MAX_ID]={0};
+int map2deti[MAX_ID]={0};
+float mapangles[MAX_ID][2]={0}; //theta and phi
+float mapanglesi[MAX_ID][2]={0}; //theta and phi
+float mapangles1[MAX_ID][2]={0}; //theta and phi
+float mapangles2[MAX_ID][2]={0};//theta and phi
+
+////////////////////
+// DETECTOR TYPES //
+////////////////////
+
+//G = Ge
+#define MAX_GE 16 // max number of Ge detectors
+#define MAX_GE_XTL 4 // max number of crystals per Ge detector
+#define MAX_GE_SEG 3 // max number of segments per Ge detector
+#define MAX_GE_BGO 1 // max number of BGO PMTs per Ge detector 
+#define GE_BGO_SUPPRESSION TRUE
+struct gdetector
+{    
+    int xmult;
+    int xid[MAX_GE_XTL]; 
+    int xe[MAX_GE_XTL];
+    long long int xt[MAX_GE_XTL]; 
+    int xct[MAX_GE_XTL];   
+    float xtheta[MAX_GE_XTL][4]; 
+    float xphi[MAX_GE_XTL][4];
+    bool xpileup[MAX_GE_XTL];
+    bool xsuppress[MAX_GE_XTL];
+    int xsubevtid[MAX_GE_XTL];
+
+    int smult;
+    int sid[MAX_GE_SEG]; 
+    int se[MAX_GE_SEG];
+    long long int st[MAX_GE_SEG];
+    int sct[MAX_GE_SEG];
+    float stheta[MAX_GE_SEG][4]; 
+    float sphi[MAX_GE_SEG][4];
+    bool spileup[MAX_GE_SEG];
+    bool ssuppress[MAX_GE_SEG];
+    int ssubevtid[MAX_GE_SEG];
+
+    int bgomult;
+    int bgoid[MAX_GE_BGO];
+    int bgoe[MAX_GE_BGO];
+    long long int bgot[MAX_GE_BGO];
+    int bgoct[MAX_GE_BGO];
+    float bgotheta[MAX_GE_XTL][4]; 
+    float bgophi[MAX_GE_XTL][4];
+    bool bgopileup[MAX_GE_BGO];
+    int bgosubevtid[MAX_GE_BGO];
+
+    int id;
+    int energy;
+    int edop;
+    long long int time;   
+    int ctime;  
+    float theta[3]; //det, xtl, or seg angle
+    float phi[3]; //det, xtl, or seg angle  
+    bool suppress; //at least one xtl was suppressed by bgo
+    bool pileup; //two or more unspressed xtls but at least one had pileup
+    bool nonprompt; //two or more unspressed xtls but at least one was non-prompt with first xtl
+    bool clean;
+}; 
+struct gdetector ge[MAX_GE]={0};
+int gmult=0;
+unsigned long long int gcount=0;
+
+//S = Si
+// .......
+
+///////////////////////////////////////
+// SPECTRA and FILE NAME DEFINITIONS // 
+///////////////////////////////////////
+//All spectra are considered two-dimensional arrays
+//Must add "write spectra" at end of file 
+
+//[Y-dim][X-dim]
+
+////////////////
+//Event Spectra
+////////////////
+
+#define HIT "hit.spn"
+int hit[2][4096]={0}; //first for all hits, second for pilup hits
+
+#define MULT "mult.spn" //total detector multiplicity for one event
+int mult[1][4096]={0};
+
+#define TDIFID "tdif.sec" //time diference between sequential events of a single channel ; 1 micro-second bins
+int tdifid[MAX_ID][8192]={0};
+
+#define E_RAW "e_raw.sec"
+int e_raw[MAX_ID][8192]={0};
+
+#define E_CAL "e_cal.sec"
+int e_cal[MAX_ID][8192]={0};
+
+#define TEVT_RAW "tevt_raw.sec" // 10 second bins
+int tevt_raw[MAX_ID][8192]={0};
+
+#define TEVT_CAL "tevt_cal.sec"
+int tevt_cal[MAX_ID][8192]={0}; // 10 second bins
+
+#define TCFD_RAW "tcfd_raw.sec"
+int tcfd_raw[MAX_ID][8192]={0};
+
+#define TDIF_RAW "tdif_raw.spn" // 10 nano-second bins
+int tdif_raw[MAX_ID][4096]={0};
+
+#define TDIF_CAL "tdif_cal.spn" // 10 nano-second bins
+int tdif_cal[MAX_ID][4096]={0};
+
+#define TDIF_CAL0_ETHRESH "tdif_cal0_ethresh.spn" //time difference relative to channel 0 ; 10 nano-second bins
+int tdif_cal0_ethresh[MAX_ID][4096]={0};
+
+#define IDID "idid.spn" //id vs id correlation hit pattern 
+int idid[MAX_ID][MAX_ID]={0};
+
+////////////////////////////
+//Detector Processed Spectra
+////////////////////////////
+
+//Ge
+#define GE_BGO_TDIF "ge_bgo_tdif.spn"
+int ge_bgo_tdif[MAX_GE][4096]={0};
+
+#define GE_XTL_TDIF "ge_xtl_tdif.spn"
+int ge_xtl_tdif[MAX_GE][4096]={0};
+
+#define GE_XTL_TDIF_ETHRESH "ge_xtl_tdif_ethresh.spn"
+int ge_xtl_tdif_ethresh[MAX_GE][4096]={0};
+
+#define GE_SPE_XTL "ge_spe_xtl.sec"
+int ge_spe_xtl[MAX_GE*MAX_GE_XTL][8192]={0};
+
+#define GE_SPE "ge_spe.sec"
+int ge_spe[MAX_GE][8192]={0};
+
+#define GE_SPE_CLEAN "ge_spe_clean.sec"
+int ge_spe_clean[MAX_GE][8192]={0};
+
+//trinity
+#define PID "pid.spn"
+int pid[4096][4096]={{0}};
+#define PID_EVSP "pid_evsp.spn"
+int pid_evsp[4096][4096]={{0}};
+#define PID_EVST "pid_evst.spn"
+int pid_evst[4096][4096]={{0}};
+#define PID_EVSTAU "pid_evstau.spn"
+int pid_evstau[4096][4096]={{0}};
+#define PID_EVSR "pid_evsr.spn"
+int pid_evsr[4096][4096]={{0}};
+#define PID_TAUVSR "pid_tauvsr.spn"
+int pid_tauvsr[4096][4096]={{0}};
+#define GAGGDT "gaggdt.spn"
+int gaggdt[4096][4096]={{0}};
+#define SPTRACE "sptrace.spn"
+int sptrace[1][4096]={{0}};
+
+//////////////////////
+//Final User Spectra
+//////////////////////
+
+//gamma-gamma
+#define GG_TDIF "gg_tdif.spn"
+int gg_tdif[4096][4096]={0};
+
+#define GG_PROMPT "gg_prompt.spn"
+int gg_prompt[4096][4096]={0};
+
+#define GG_NPROMPT "gg_nprompt.spn"
+int gg_nprompt[4096][4096]={0};
+
+
+
+
+///////////////////////
+// Write 2-byte data //
+///////////////////////
+void write_data2(char *filename, short *data, int xdim, int ydim, int overwrite) { //2byte per channel Write / Add to previous
+    
+    FILE *FP;
+    int i;
+    short *previous;
+    
+    if(!overwrite) {
+        //allocate memory for 1d-array for reading in rows of 2d Radware matrix
+        if ( ( previous = (short *)malloc(xdim * ydim * sizeof(short)) ) == NULL ) {
+            printf("\nError, memory not allocated.\n");
+            exit(1);
+        }
+      
+        //open previous spectra file  
+        if( (FP=fopen(filename, "r")) != NULL ){
+            fread(previous, sizeof(short)*xdim*ydim, 1, FP);        
+            fclose(FP);
+            //update spectra
+            for (i=0; i<xdim*ydim; i++) {
+                if(previous[i] < (powf(2,sizeof(short)*8.0)-2)) 
+                    data[i] = data[i] + previous[i];
+                }   
+            }
+        else{
+            printf("%s did not previously exist, creating ...\n", filename);       
+        }     
+   
+        //Deallocate previous data
+        free(previous);
+    }
+  
+    FP=fopen(filename, "w");
+    fwrite(data, sizeof(short)*xdim, ydim, FP);
+    fclose(FP);            
+}    
+
+
+///////////////////////
+// Write 4-byte data //
+///////////////////////
+void write_data4(char *filename, int *data, int xdim, int ydim, int overwrite) { //4byte per channel Write / Add to previous
+    
+    FILE *FP;
+    int i;
+    int *previous;
+   
+    if(!overwrite) {
+        //allocate memory for 1d-array for reading in rows of 2d Radware matrix
+        if ( ( previous = (int *)malloc(xdim * ydim * sizeof(int)) ) == NULL ) {
+            printf("\nError, memory not allocated.\n");
+            exit(1);
+        }
+      
+        //open previous spectra file  
+        if( (FP=fopen(filename, "r")) != NULL ){
+            fread(previous, sizeof(int)*xdim*ydim, 1, FP);        
+            fclose(FP);
+            //update spectra
+            for (i=0; i<xdim*ydim; i++) {
+                if(previous[i] < (powf(2,sizeof(int)*8.0)-2)) 
+                    data[i] = data[i] + previous[i];
+            }   
+        } 
+        else{
+            printf("%s did not previously exist, creating ...\n", filename);       
+        }   
+       
+        //Deallocate previous data
+        free(previous);
+    }
+  
+    FP=fopen(filename, "w");
+    fwrite(data, sizeof(int)*xdim, ydim, FP);
+    fclose(FP);            
+}    
+
+
+
+
+///////////////////////////////////
+// START OF MAIN FUNCTION        //
+///////////////////////////////////
+int main(int argc, char **argv) {
+  
+  int i=0, j=0, k=0;
+  float tempf=0;
+  int max1=0, min1=0;
+  int max2=0, min2=0;
+  int maxid1=-1, minid1=-1;
+  int maxid2=-1, minid2=-1;
+  div_t e_div;
+  lldiv_t lle_div;
+
+  int overwrite = 1;    
+
+  double etrace0,etrace1,btrace0,btrace1;
+  double ptrace0,ptrace1,ttrace0,ttrace1,tautrace0,tautrace1;
+  int dbcount = 0;
+  long long int strace[500];
+  memset(strace, 0, sizeof(strace));
+
+  //temp buffer for each sub event
+  unsigned int sub[MAX_SUB_LENGTH];
+  memset(sub, 0, sizeof(sub));
+  
+  //Reference time and difference for event building
+  long long int etime, tdif, idtime[MAX_ID]={0}, temptime;
+  
+   
+  // Check that the corrent number of arguments were provided.
+  if (argc<3)    {
+    printf("Incorrect number of arguments:\n");
+    printf("%s -op datafile calibrationfile mapfile \n", argv[0]);
+    printf("\n .... calibration file is optional\n");
+    printf(" .... map file is optional\n");
+    printf(" .... o for overwrite spectra\n");
+    printf(" .... u for update spectra\n");
+    printf(" .... p for print realtime stats\n");
+    return 1;
+  }
+    
+  if(strstr(argv[1], "u") != NULL) {
+    overwrite = 0;
+    printf("Updating Spectra\n");
+  }
+  else {
+    printf("Overwriting Spectra\n");
+  }  
+  
+  //open list-mode data file from PXI digitizer  
+  FILE *fpr;
+  long int fprsize,fprpos;
+  if ((fpr = fopen(argv[2], "r")) == NULL) {
+    fprintf(stderr, "Error, cannot open input file %s\n", argv[2]);
+    return 1;
+  }
+ 
+  //get file size
+  fseek(fpr, 0L, SEEK_END);
+  fprsize = ftell(fpr);
+  rewind(fpr);
+  
+   
+  //open debug file for streaming an 1d array
+  FILE *debugfile;
+  if ((debugfile = fopen(DEBUGFN, "w")) == NULL) {
+    fprintf(stderr, "Error, cannot open %s\n", DEBUGFN);
+    return 1;
+  }   
+   
+
+  //buffer for reading in text files for calibrations and maps below
+  char line[LINE_LENGTH];
+
+
+  //open energy and time calibration file (e.g., *.ca3 file)  
+  int calid=0; 
+  float caloffset=0.0, calgain=0.0;
+  int firstcal=0, necal=0, ntcal=0;
+  
+  FILE *fprcal;
+  
+  for (i=0; i<MAX_ID; i++) {
+    ecal[i][0] = caloffset;
+    ecal[i][1] = calgain;
+    tcal[i][0] = caloffset;
+    tcal[i][1] = calgain;     
+  }  
+  
+  if (argc >= 4) {
+
+    if ((fprcal = fopen(argv[3], "r")) == NULL) {
+        fprintf(stderr, "Error, cannot open input file %s\n", argv[3]);
+        return 1;
+    } 
+    
+    printf("%s loaded!\n", argv[3]);
+
+	while(fgets(line, LINE_LENGTH, fprcal) != NULL){
+	    calid=0; caloffset=0; calgain=0;
+        for(i=0; i<LINE_LENGTH; i++){
+            if(line[i]=='#'){
+                if(PRINT_CAL)printf("%s", line);
+                break;
+            }
+            else if(line[i]>=0){
+                if (firstcal==0) {
+                    sscanf(line,"%f\n", &new_gain);    
+	                if(PRINT_CAL) printf("%f\n", new_gain);                                       
+                    firstcal=1;
+                    break;
+                }    
+	            sscanf(line,"%d\t%f\t%f\n", &calid, &caloffset, &calgain);
+	            if(PRINT_CAL) printf("%d\t%.4f\t%.4f\n", calid, caloffset, calgain);
+		        if(calid >=0 && calid < MAX_ID) {
+                    ecal[calid][0] = caloffset;
+                    ecal[calid][1] = calgain;
+                    necal++;
+          		    break;
+			    }
+		        if(calid >=1000 && calid < 1000+MAX_ID) {
+                    tcal[calid-1000][0] = caloffset;
+                    tcal[calid-1000][1] = calgain;
+                    ntcal++;
+          		    break;
+			    }			    
+			    else {
+				    printf("Error in reading %s : bad id or format\n", argv[3]);
+				    return -1;
+			    }
+            }
+	        else if(line[i]=='\n'){
+                if(PRINT_CAL) printf("\n");
+                break;
+            }
+            else {
+                continue;
+            }
+
+        }
+        	memset(line, 0, LINE_LENGTH);
+ 	  }
+ 	  
+ 	fclose(fprcal);
+    printf("read %d energy calibrations\n", necal);
+    printf("read %d time calibrations\n", ntcal);
+  }
+  
+
+
+  //open ID->Detector map file (e.g., *.map file)  
+  int mapid=0, detid=0, detidi; 
+  char dettype=0;
+  float theta=0, phi=0, thetai=0, phii=0, theta1=0, phi1=0, theta2=0, phi2=0; 
+  FILE *fprmap;
+  int nmmap=0;
+  if (argc >= 5) {
+
+    if ((fprmap = fopen(argv[4], "r")) == NULL) {
+        fprintf(stderr, "Error, cannot open input file %s\n", argv[4]);
+        return 1;
+    } 
+
+    printf("%s loaded!\n", argv[4]);
+
+	while(fgets(line, LINE_LENGTH, fprmap) != NULL){
+	    mapid=0; dettype=0; thetai=0; phii=0; theta=0; phi=0; theta1=0; phi1=0; theta2=0; phi2=0; 
+        for(i=0; i<LINE_LENGTH; i++){
+            if(line[i]=='#'){
+                if(PRINT_MAP)printf("%s", line);
+                break;
+            }
+            else if(line[i]>=0){   
+	            sscanf(line,"%d\t%c\t%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", &mapid, &dettype, &detid, &detidi, &theta, &phi, &thetai, &phii, &theta1, &phi1, &theta2, &phi2);
+	            if(PRINT_MAP) printf("%d\t%c\t%d\t%d\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n", mapid, dettype, detid, detidi, theta, phi, thetai, phii, theta1, phi1, theta2, phi2);
+		        if(mapid >=0 && mapid < MAX_ID) {
+                    map2type[mapid] = dettype;
+                    map2det[mapid] = detid;
+                    map2deti[mapid] = detidi;
+                    mapangles[mapid][0] = theta;
+                    mapangles[mapid][1] = phi;
+                    mapanglesi[mapid][0]=thetai; 
+                    mapanglesi[mapid][1]=phii; 
+                    mapangles1[mapid][0]= theta1; 
+                    mapangles1[mapid][1]= phi1; 
+                    mapangles2[mapid][0]= theta2; 
+                    mapangles2[mapid][1]= phi2; 
+                    nmmap++;
+          		    break;
+			    }			    
+			    else {
+				    printf("Error in reading %s : bad id or format\n", argv[4]);
+				    return -1;
+			    }
+            }
+	        else if(line[i]=='\n'){
+                if(PRINT_MAP) printf("\n");
+                break;
+            }
+            else {
+                continue;
+            }
+
+        }
+        	memset(line, 0, LINE_LENGTH);
+ 	  }
+ 	  
+ 	fclose(fprmap);
+    printf("read %d id maps\n", nmmap);
+  }
+
+
+
+
+  /////////////////////
+  // MAIN WHILE LOOP //
+  /////////////////////
+  while (1) { //main while loop 
+
+
+      /////////////////////////////////
+      // UNPACK DATA AND EVENT BUILD //
+      /////////////////////////////////
+      etime=-1; tdif=-1; sevtmult=0;  
+      //memset(&subevt, 0, sizeof(subevt)); //not needed since everything is redefined (except maybe trace on pileup evts)
+      while (1) { //get subevents and event build for one "event" 
+        
+       // memset(&subevt[sevtmult], 0, sizeof(subevt[sevtmult])); //not needed since everything is redefined (except maybe trace on pileup evts)
+        
+        //read 4-byte header
+        if (fread(sub, sizeof(int)*HEADER_LENGTH, 1, fpr) != 1) break;
+        subevt[sevtmult].chn = sub[0] & 0xF;
+        subevt[sevtmult].sln = (sub[0] & 0xF0) >> 4;
+        subevt[sevtmult].crn = (sub[0] & 0xF00) >> 8;
+        subevt[sevtmult].id = subevt[sevtmult].crn*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (subevt[sevtmult].sln-BOARD_START)*MAX_CHANNELS_PER_BOARD + subevt[sevtmult].chn;   
+        subevt[sevtmult].hlen = (sub[0] & 0x1F000) >> 12;
+        subevt[sevtmult].elen = (sub[0] & 0x7FFE0000) >> 17;
+        subevt[sevtmult].fcode = (sub[0] & 0x80000000) >> 31;
+        subevt[sevtmult].time = ( (long long int)(sub[2] & 0xFFFF) << 32) + sub[1];
+        subevt[sevtmult].ctime = (sub[2] & 0x7FFF0000) >> 16;
+        subevt[sevtmult].ctimef = (sub[2] & 0x80000000) >> 31;
+        subevt[sevtmult].energy = (sub[3] & 0xFFFF);
+        subevt[sevtmult].trlen = (sub[3] & 0x7FFF0000) >> 16;
+        subevt[sevtmult].trwlen = subevt[sevtmult].trlen / 2;
+        subevt[sevtmult].extra = (sub[3] & 0x80000000) >> 31;       
+ 
+        //rebin raw trap energy from 32k to ....            
+        tempf = (float)subevt[sevtmult].energy/RAWE_REBIN_FACTOR;// + RAND;       
+        subevt[sevtmult].energy = (int)tempf;                        
+ 
+        //check lengths (sometimes all of the bits for trace length are turned on ...)
+       /* if (subevt[sevtmult].elen - subevt[sevtmult].hlen != subevt[sevtmult].trwlen) {
+            printf("SEVERE ERROR: event, header, and trace length inconsistencies found\n");
+            printf("event length = %d\n", subevt[sevtmult].elen);
+            printf("header length = %d\n", subevt[sevtmult].hlen);
+            printf("trace length = %d\n", subevt[sevtmult].trwlen);  
+            printf("Extra = %d\n", subevt[sevtmult].extra); 
+            printf("fcode = %d\n", subevt[sevtmult].fcode);              
+            //sleep(1);          
+            //return 0;
+        } */ 
+        
+       
+        //Set reference time for event building
+        if (etime == -1) {
+            etime = subevt[sevtmult].time;
+            tdif = 0;
+        }
+        else {
+            tdif = subevt[sevtmult].time - etime;
+            if (tdif < 0) {
+                printf("SEVERE ERROR: tdiff < 0, file must be time sorted\n");
+                printf("etime = %lld, time = %lld, and tdif = %lld\n", etime, subevt[sevtmult].time, tdif);                
+                return 0;   
+            }    
+        }    
+      
+        //Check for end of event, rewind, and break out of while loop
+        if (tdif > EVENT_BUILD_TIME) {
+            fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR); //fwrite/fread is buffered by system ; storing this in local buffer is no faster!
+            break;           
+        }    
+               
+        
+        //time between sequential events for a single channel ; useful for determining optimal event building time
+        temptime = (subevt[sevtmult].time - idtime[subevt[sevtmult].id])/100; //rebin to 1 micro-second
+        if ( temptime >= 0 && temptime < 8192) {
+            tdifid[subevt[sevtmult].id][temptime]++;
+        }    
+        idtime[subevt[sevtmult].id]=subevt[sevtmult].time; //store time for next subevent of channel    
+    
+        // total pileups
+        if (subevt[sevtmult].fcode==1) {
+            pileupcount++;
+        }
+        
+        //Histogram raw spectra
+        hit[0][subevt[sevtmult].id]++;
+        if (subevt[sevtmult].fcode==1) 
+        hit[1][subevt[sevtmult].id]++;     
+  
+        if (subevt[sevtmult].energy >= 0 && subevt[sevtmult].energy < 8192)
+        e_raw[subevt[sevtmult].id][subevt[sevtmult].energy]++;
+ 
+        if (subevt[sevtmult].time/1000000000 >= 0 && subevt[sevtmult].time/1000000000 < 8192) // rebin to 10 seconds
+        tevt_raw[subevt[sevtmult].id][subevt[sevtmult].time/1000000000]++; // rebin to 10 seconds
+
+        if (subevt[sevtmult].ctime >= 0 && subevt[sevtmult].ctime < 8192)
+        tcfd_raw[subevt[sevtmult].id][subevt[sevtmult].ctime]++;
+
+        if (tdif >= 0 && tdif < 4096 && sevtmult!=0)
+        tdif_raw[subevt[sevtmult].id][tdif]++;
+
+
+        //if CFD is enabled, ctime will be non-zero
+        //tempf = (float)subevt[sevtmult].ctime*10.0/32768.0;
+        //subevt[sevtmult].time = subevt[sevtmult].time + (long long int)tempf;    
+
+        //Calibrate energy and time
+        tempf = ((float)subevt[sevtmult].energy*ecal[subevt[sevtmult].id][1] + ecal[subevt[sevtmult].id][0])/new_gain;// + RAND;       
+        subevt[sevtmult].energy = (int)tempf; 
+        //subevt[sevtmult].time += (long long int)tcal[subevt[sevtmult].id][0];       
+	
+        //Histogram calibrated spectra
+        if (subevt[sevtmult].energy >= 0 && subevt[sevtmult].energy < 8192)
+        e_cal[subevt[sevtmult].id][subevt[sevtmult].energy]++;        
+        
+        if (subevt[sevtmult].time/1000000000 >= 0 && subevt[sevtmult].time/1000000000 < 8192)        
+        tevt_cal[subevt[sevtmult].id][subevt[sevtmult].time/1000000000]++;
+         
+        //continue on if no trace, esum, or qsum
+        if (subevt[sevtmult].hlen==HEADER_LENGTH && subevt[sevtmult].trwlen==0 ) {
+            sevtmult++;
+            continue;
+        }
+        
+        //more data than just the header; read entire sub event
+        fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR);
+        if (fread(sub, sizeof(int)*subevt[sevtmult].elen, 1, fpr) != 1) break;
+                              
+        //trace
+        k=0;
+        for (i = subevt[sevtmult].hlen; i < subevt[sevtmult].elen; i++) {      
+            subevt[sevtmult].tr[i - subevt[sevtmult].hlen + k] = sub[i] & 0x3FFF;
+            subevt[sevtmult].tr[i - subevt[sevtmult].hlen + k + 1] = (sub[i]>>16) & 0x3FFF;
+            k=k+1;
+        } 
+        
+     // if (subevt[sevtmult].id == 4 && subevt[sevtmult].fcode == 1) DB(subevt[sevtmult].tr);
+            
+        //continue if no esum or qsum   
+        if (subevt[sevtmult].hlen==HEADER_LENGTH) {
+            sevtmult++;        
+            continue;
+        }
+        
+        //esum
+        if (subevt[sevtmult].hlen==8 || subevt[sevtmult].hlen==16) { 
+            for (i=4; i < 8; i++) {
+                subevt[sevtmult].esum[i-4] = sub[i];
+            }
+        }
+    
+        //qsum
+        if (subevt[sevtmult].hlen==12) { 
+            for (i=4; i < 12; i++) {
+                subevt[sevtmult].qsum[i-4] = sub[i];
+            }
+        }
+    
+        //qsum
+        if (subevt[sevtmult].hlen==16) { 
+            for (i=8; i < 16; i++) {
+                subevt[sevtmult].qsum[i-8] = sub[i];
+            }
+        }    
+        
+        sevtmult++;
+     
+      } //end while loop for unpacking sub events and event building for one "event"
+      if (sevtmult==0) break; //end main WHILE LOOP when out of events 
+      mult[0][sevtmult]++; //Histogram raw sub event multiplicity 
+      sevtcount += sevtmult;
+      evtcount++; //event-built number
+      /////////////////////////////////////
+      // END UNPACK DATA AND EVENT BUILD //
+      /////////////////////////////////////
+
+  	//int GAGG1=58;
+	//int GAGG2=59;
+
+  	int GAGG1=60;
+	int GAGG2=70;
+
+      if (sevtmult == 2 && (subevt[0].id == GAGG1 || subevt[0].id == GAGG2) && (subevt[1].id == GAGG1 || subevt[1].id == GAGG2) && ((subevt[0].time - subevt[1].time) + 2000 > 1992) && ((subevt[0].time - subevt[1].time) + 2000 < 2005) ) {
+
+	if (evtcount < 200) DB(subevt[0].tr);
+
+        etrace0 = 0;
+        btrace0 = 0;
+        etrace1 = 0;
+        btrace1 = 0;
+	ptrace0 = 0;
+	ptrace1 = 0;
+	ttrace0 = 0;
+	ttrace1 = 0;
+	tautrace0 = 0;
+	tautrace1 = 0;
+
+	for (i=0; i<60; i++) {
+            btrace0 = btrace0 + subevt[0].tr[i];
+            btrace1 = btrace1 + subevt[1].tr[i];
+	}
+	btrace0 = btrace0/60.;
+        btrace1 = btrace1/60.;
+
+	for (i=60; i<173; i++) { //180-230
+            etrace0 = etrace0 + subevt[0].tr[i] - btrace0;
+            etrace1 = etrace1 + subevt[1].tr[i] - btrace1;
+	    tautrace0 = tautrace0 + (subevt[0].tr[i] - btrace0)*i;
+            tautrace1 = tautrace1 + (subevt[1].tr[i] - btrace1)*i;
+       	}
+	tautrace0 = tautrace0 / etrace0;
+        tautrace1 = tautrace1 / etrace1;
+        etrace0 = etrace0 / 2.;
+        etrace1 = etrace1 / 2.;
+
+
+        //peak sum
+        for (i=77; i<93; i++) { //180-223 // 197-213
+            ptrace0 = ptrace0 + subevt[0].tr[i] - btrace0;
+            ptrace1 = ptrace1 + subevt[1].tr[i] - btrace1;
+        }
+//        for (i=193; i<199; i++) { //180-223 // 197-213
+//            ptrace0 = ptrace0 - (subevt[0].tr[i] - btrace0);
+//            ptrace1 = ptrace1 - (subevt[1].tr[i] - btrace1);
+//        }
+        ptrace0 = ptrace0 / 2.;
+        ptrace1 = ptrace1 / 2.;
+	ptrace0 = (ptrace0 + ptrace1)/10;
+
+        //tail sum
+        for (i=101; i<160; i++) { //223-293 // 221-280
+            ttrace0 = ttrace0 + subevt[0].tr[i] - btrace0;
+            ttrace1 = ttrace1 + subevt[1].tr[i] - btrace1;
+	}
+        ttrace0 = ttrace0 / 2.;
+        ttrace1 = ttrace1 / 2.;
+	ttrace0 = (ttrace0 + ttrace1)/10;
+
+        if ((int)ptrace0 > 0 && (int)ptrace0 < 4096 && (int)ttrace0 > 0 && (int)ttrace0 < 4096) pid[(int)ptrace0][(int)ttrace0]++;
+        if ((int)etrace0 > 0 && (int)etrace0 < 4096 && (int)ptrace0 > 0 && (int)ptrace0 < 4096) pid_evsp[(int)etrace0][(int)ptrace0]++;
+        if ((int)etrace0 > 0 && (int)etrace0 < 4096 && (int)ttrace0 > 0 && (int)ttrace0 < 4096) pid_evst[(int)etrace0][(int)ttrace0]++;
+        if ((int)etrace0 > 0 && (int)etrace0 < 4096 && (int)tautrace0 > 0 && (int)tautrace0 < 4096) pid_evstau[(int)etrace0][(int)tautrace0]++;
+        if ((int)(etrace0) > 0 && (int)(etrace0) < 4096 && (int)((100.*ttrace0)/ptrace0) > 0 && (int)((100.*ttrace0)/ptrace0) < 4096) pid_evsr[(int)(etrace0)][(int)((100.*ttrace0)/ptrace0)]++;
+        if ((int)(tautrace0) > 0 && (int)(tautrace0) < 4096 && (int)((100.*ttrace0)/ptrace0) > 0 && (int)((100.*ttrace0)/ptrace0) < 4096) pid_tauvsr[(int)(tautrace0)][(int)((100.*ttrace0)/ptrace0)]++;
+
+
+
+/*
+	for (i=1; i<6; i++) {
+	    etrace0 = etrace0 + subevt[0].qsum[i];
+            etrace1 = etrace1 + subevt[1].qsum[i];
+	    //printf("subevt[0].qsum[%d]=%d\n", i, subevt[0].qsum[i]);
+	}
+	etrace0 = etrace0 - subevt[0].qsum[0]*0.5/1.8;
+        etrace1 = etrace1 - subevt[1].qsum[0]*0.5/1.8;
+
+	etrace0 = etrace0 / 1;
+        etrace1 = etrace1 / 1;
+*/
+
+	//return 0;
+
+        //printf("etrace = %f\n", etrace);
+        if ((int)etrace0 >0 && (int)etrace0 < 8192) e_raw[200-GAGG1+subevt[0].id][(int)etrace0]++;
+        if ((int)etrace1 >0 && (int)etrace1 < 8192) e_raw[200-GAGG1+subevt[1].id][(int)etrace1]++;
+        if ((int)etrace0 >0 && (int)etrace0 < 8192 && (int)ttrace0>500/10 && (int)ptrace0 > 978/10 && (int)ptrace0 < 1270/10 ) {
+		e_raw[203][(int)etrace0]++;
+        	//if (evtcount < 200) DB(subevt[0].tr);
+	}
+      //	if ((int)ttrace0<40 && (int)ptrace0 > 131) {
+                if ((int)etrace0>1350 && (int)etrace0<1400 && (int)((100.*ttrace0)/ptrace0) < 200 ) {
+//			DB(subevt[0].tr);
+			for (i=0; i<500; i++) {
+				strace[i]=strace[i]+subevt[0].tr[i];
+			}
+			dbcount++;
+		}
+      //  }
+
+        etrace1 = etrace0 + etrace1;
+        etrace1 = etrace1/2.0;
+        if ((int)etrace1 >0 && (int)etrace1 < 8192) {
+          e_raw[202][(int)etrace1]++;
+        }
+
+	if ((subevt[0].time - subevt[1].time) + 2000 > 0 && (subevt[0].time - subevt[1].time) + 2000 < 4096 && (int)etrace1 >0 && (int)etrace1 < 4096)
+	gaggdt[(subevt[0].time - subevt[1].time) + 2000][(int)etrace1]++;
+
+      }
+
+      //skip detector building below if no map file
+      if (argc >= 5) {        
+
+
+      ////////////////////////////////////////
+      // MAP SUB EVENTS INTO DETECTOR TYPES //
+      ////////////////////////////////////////
+      memset(&ge, 0, sizeof(ge)); //This is needed but could be replaced by setting suppress, pileup, nonprompt, clean, x/s/bmult to zero at start of loop!
+      
+      for (i=0; i<sevtmult; i++) { 
+        
+        for (j=0; j<sevtmult; j++) {
+            if (i!=j)
+            idid[subevt[i].id][subevt[j].id]++;
+        }
+        
+        //printf("i=%d, sevtmult=%d, subevt[i].id=%d, map2type[subevt[i].id]=%c, subevt[i].energy=%d, subevt[i].time=%lld\n", i, sevtmult, subevt[i].id, map2type[subevt[i].id], subevt[i].energy, subevt[i].time);
+        //fflush(stdout);
+          
+        //if CFD is enabled, ctime will be non-zero
+        //tempf = (float)subevt[i].ctime*10.0/32768.0;
+        //tempf = (float)subevt[i].ctime*1.0/32768.0; //should make no difference ; need decimals or convert time to 1 ns bins
+        //printf("%lld e=%d ", subevt[i].time, subevt[i].energy);
+        //subevt[i].time = subevt[i].time + (long long int)tempf;  
+        //printf("%lld and %f\n", subevt[i].time, tempf); 
+  
+               
+        //Histogram calibrated tdif spectra (Do here to keep subevents within an event time ordered during event build above)
+        if (i==0) etime = subevt[i].time;
+        tdif = abs(subevt[i].time - etime); 
+        if (tdif >= 0 && tdif < 4096 && i!=0) tdif_cal[subevt[i].id][tdif]++;
+        //if (tdif > 20 && tdif < 50) DB(subevt[i].tr);
+        
+        //tdif with respect to channel id 0
+        if (i!=0 && subevt[0].id==0 && subevt[0].energy >= 10 && subevt[0].energy <= 8000 && subevt[i].energy >= 10  && subevt[i].energy <= 8000) {
+            if (tdif >= 0 && tdif < 4096 ) tdif_cal0_ethresh[subevt[i].id][tdif]++;   
+        }            
+        
+
+        /////////////////////
+        // G = Ge Detector //
+        /////////////////////              
+        if ( map2type[subevt[i].id] == 'G' ) { //Keep G and ge or switch to C and clover/cl?
+        
+            if (map2deti[subevt[i].id] > 0 && map2deti[subevt[i].id] <= MAX_GE_XTL) { //Ge crystal 
+               if ( ge[map2det[subevt[i].id]].xmult >= MAX_GE_XTL) {               
+                    printf("SEVERE ERROR: Same Ge(xtl) twice within event build; Make event build time smaller!\n");    
+                    printf("ge[map2det[subevt[i].id]].xmult=%d, ge[map2det[subevt[i].id]].smult=%d, ge[map2det[subevt[i].id]].bgomult=%d\n", ge[map2det[subevt[i].id]].xmult, ge[map2det[subevt[i].id]].smult,ge[map2det[subevt[i].id]].bgomult);
+                    continue;
+                    //return -1;
+                }   
+                ge[map2det[subevt[i].id]].xid[ge[map2det[subevt[i].id]].xmult] = map2deti[subevt[i].id];               
+                ge[map2det[subevt[i].id]].xe[ge[map2det[subevt[i].id]].xmult] = subevt[i].energy;               
+                ge[map2det[subevt[i].id]].xt[ge[map2det[subevt[i].id]].xmult] = subevt[i].time;                           
+                ge[map2det[subevt[i].id]].xct[ge[map2det[subevt[i].id]].xmult] = subevt[i].ctime;  
+                ge[map2det[subevt[i].id]].xpileup[ge[map2det[subevt[i].id]].xmult] = subevt[i].fcode;  
+                ge[map2det[subevt[i].id]].xsubevtid[ge[map2det[subevt[i].id]].xmult] = i;  
+                        
+                ge[map2det[subevt[i].id]].xtheta[ge[map2det[subevt[i].id]].xmult][0] = mapangles[subevt[i].id][0];   
+                ge[map2det[subevt[i].id]].xtheta[ge[map2det[subevt[i].id]].xmult][1] = mapanglesi[subevt[i].id][0];               
+                ge[map2det[subevt[i].id]].xtheta[ge[map2det[subevt[i].id]].xmult][2] = mapangles1[subevt[i].id][0];               
+                ge[map2det[subevt[i].id]].xtheta[ge[map2det[subevt[i].id]].xmult][3] = mapangles2[subevt[i].id][0];               
+                ge[map2det[subevt[i].id]].xphi[ge[map2det[subevt[i].id]].xmult][0] = mapangles[subevt[i].id][1];               
+                ge[map2det[subevt[i].id]].xphi[ge[map2det[subevt[i].id]].xmult][1] = mapanglesi[subevt[i].id][1];               
+                ge[map2det[subevt[i].id]].xphi[ge[map2det[subevt[i].id]].xmult][2] = mapangles1[subevt[i].id][1];               
+                ge[map2det[subevt[i].id]].xphi[ge[map2det[subevt[i].id]].xmult][3] = mapangles2[subevt[i].id][1]; 
+
+                ge[map2det[subevt[i].id]].xmult++;   
+            }                                                                
+            if (map2deti[subevt[i].id] > MAX_GE_XTL && map2deti[subevt[i].id] <= MAX_GE_XTL + MAX_GE_SEG) { //Ge segment
+               if ( ge[map2det[subevt[i].id]].smult >= MAX_GE_SEG ) {               
+                    printf("SEVERE ERROR: Same Ge(seg) twice within event build; Make event build time smaller!\n");    
+                    printf("ge[map2det[subevt[i].id]].xmult=%d, ge[map2det[subevt[i].id]].smult=%d, ge[map2det[subevt[i].id]].bgomult=%d\n", ge[map2det[subevt[i].id]].xmult, ge[map2det[subevt[i].id]].smult,ge[map2det[subevt[i].id]].bgomult);
+                    continue;
+                    //return -1;
+                }   
+                ge[map2det[subevt[i].id]].sid[ge[map2det[subevt[i].id]].smult] = map2deti[subevt[i].id];               
+                ge[map2det[subevt[i].id]].se[ge[map2det[subevt[i].id]].smult] = subevt[i].energy;               
+                ge[map2det[subevt[i].id]].st[ge[map2det[subevt[i].id]].smult] = subevt[i].time;                           
+                ge[map2det[subevt[i].id]].sct[ge[map2det[subevt[i].id]].smult] = subevt[i].ctime;  
+                ge[map2det[subevt[i].id]].spileup[ge[map2det[subevt[i].id]].smult] = subevt[i].fcode;  
+                ge[map2det[subevt[i].id]].ssubevtid[ge[map2det[subevt[i].id]].smult] = i;  
+                            
+                ge[map2det[subevt[i].id]].smult++;               
+            }            
+            if (map2deti[subevt[i].id] > MAX_GE_XTL + MAX_GE_SEG) { //BGO 
+               if ( ge[map2det[subevt[i].id]].bgomult >= MAX_GE_BGO ) {               
+                    printf("SEVERE ERROR: Same Ge(bgo) twice within event build; Make event build time smaller!\n");    
+                    printf("ge[map2det[subevt[i].id]].xmult=%d, ge[map2det[subevt[i].id]].smult=%d, ge[map2det[subevt[i].id]].bgomult=%d\n", ge[map2det[subevt[i].id]].xmult, ge[map2det[subevt[i].id]].smult,ge[map2det[subevt[i].id]].bgomult);
+                    continue;
+                    //return -1;
+                }   
+                ge[map2det[subevt[i].id]].bgoid[ge[map2det[subevt[i].id]].bgomult] = map2deti[subevt[i].id];               
+                ge[map2det[subevt[i].id]].bgoe[ge[map2det[subevt[i].id]].bgomult] = subevt[i].energy;               
+                ge[map2det[subevt[i].id]].bgot[ge[map2det[subevt[i].id]].bgomult] = subevt[i].time;                           
+                ge[map2det[subevt[i].id]].bgoct[ge[map2det[subevt[i].id]].bgomult] = subevt[i].ctime;              
+                ge[map2det[subevt[i].id]].bgopileup[ge[map2det[subevt[i].id]].bgomult] = subevt[i].fcode;
+                ge[map2det[subevt[i].id]].bgosubevtid[ge[map2det[subevt[i].id]].bgomult] = i;  
+
+                ge[map2det[subevt[i].id]].bgomult++;               
+            } 
+                 
+                
+        } //end G
+
+        ///////////////////
+        // S=Si Detector //
+        ///////////////////  
+        
+      } // end i loop over sevtmult     
+      ////////////////////////////////////////////
+      // END MAP SUB EVENTS INTO DETECTOR TYPES //
+      ////////////////////////////////////////////
+
+
+      ////////////////////////////
+      // PROCESS DETECTOR TYPES //
+      ////////////////////////////
+   
+
+      /////////////////////
+      // G = Ge Detector //
+      /////////////////////      
+      gmult = 0;
+      for (i=0; i<MAX_GE; i++) {
+        
+        max1 = -1; max2 = -1;
+        maxid1 = -1; maxid2 = -1;
+        
+        ge[i].id = i;
+        //Addback and Compton Suppression 
+	for (j=0; j<ge[i].xmult; j++) {
+            
+            //compton suppression per crystal 
+            if (GE_BGO_SUPPRESSION == TRUE) {
+                for (k=0; k<ge[i].bgomult; k++) {
+                    tdif = abs( ge[i].xt[j] - ge[i].bgot[k] );
+                    if (tdif >= 0 && tdif <= 4096) ge_bgo_tdif[i][tdif]++;
+                    if ( (tdif < 50 && ge[i].bgoe[k] > 10) || ge[i].bgopileup[k]==TRUE) { //need to fix bgo pileup with trace analysis 
+                        ge[i].xsuppress[j] = TRUE; 
+                        ge[i].suppress = TRUE;                        
+                    }  
+                }
+            }      
+            
+            //addback     
+            if (ge[i].xsuppress[j] == FALSE) {       
+                if (ge[i].xpileup[j] == TRUE) {
+                    ge[i].pileup = TRUE;     
+                    continue;
+                }              
+                //xtl spectra
+                if (ge[i].xe[j] > 0 && ge[i].xe[j] < 8192 && ge[i].id >= 1) ge_spe_xtl[(ge[i].id-1)*MAX_GE_XTL + ge[i].xid[j]][ge[i].xe[j]]++;                        
+                tdif = abs( ge[i].xt[j] - ge[i].time );
+                if (tdif >= 0 && tdif <= 4096 && ge[i].time != 0 ) ge_xtl_tdif[i][tdif]++;           
+                                       
+                if (ge[i].xe[j] > 50 && ge[i].xe[j] < 4000) {                  
+                    if (tdif >= 0 && tdif <= 4096 && ge[i].time != 0 ) ge_xtl_tdif_ethresh[i][tdif]++;           
+                    if (tdif < 20 || ge[i].time == 0) {  
+                        ge[i].energy = ge[i].energy + ge[i].xe[j];
+                        if (max1 < ge[i].xe[j]) {
+                            max1 = ge[i].xe[j];
+                            maxid1 = j;
+                            ge[i].time = ge[i].xt[j];
+                            ge[i].ctime = ge[i].xct[j];
+                        } 
+                    }
+                    else {
+                        ge[i].nonprompt = TRUE; // the first time will become the adopted value / event   
+                    }        
+                } 
+            }  
+           
+        }
+        
+        if (max1 == -1) continue;
+        
+        //Segmentation Position and Compton Suppression
+        for (j=0; j<ge[i].smult; j++) {
+
+            //compton suppression per segment 
+            if (GE_BGO_SUPPRESSION == TRUE) {            
+                for (k=0; k<ge[i].bgomult; k++) {
+                    tdif = abs( ge[i].st[j] - ge[i].bgot[k] );
+                    if (tdif < 50 && ge[i].bgoe[k] > 10) {
+                        ge[i].ssuppress[j] = TRUE;    
+                    }
+                }
+            }      
+
+            //segment
+            if (ge[i].ssuppress[j] == FALSE && ge[i].se[j] > 0 && ge[i].se[j] < 10000) {                    
+                if (max2 < ge[i].se[j]) {
+                    max2 = ge[i].se[j];
+                    maxid2 = j;
+                }                                           
+            }  
+
+        }    
+        
+        //Angle assignments
+        ge[i].theta[0] = ge[i].xtheta[maxid1][0];                       //detector center
+        ge[i].phi[0] = ge[i].xphi[maxid1][0];
+        ge[i].theta[1] = ge[i].xtheta[maxid1][1];                       //crystal center
+        ge[i].phi[1] = ge[i].xphi[maxid1][1];       
+        if (ge[i].sid[maxid2] == 6) {                                   // side channel C  
+            ge[i].theta[2] = ge[i].xtheta[maxid1][2];
+            ge[i].phi[2] = ge[i].xphi[maxid1][2];  
+        }
+        else if (ge[i].sid[maxid2] == 5 || ge[i].sid[maxid2] == 7) {    // side channel L/R
+            ge[i].theta[2] = ge[i].xtheta[maxid1][3];
+            ge[i].phi[2] = ge[i].xphi[maxid1][3];             
+        }  
+        else {      
+            ge[i].theta[2] = ge[i].xtheta[maxid1][1];                   // side channel failure --> crystal center
+            ge[i].phi[2] = ge[i].xphi[maxid1][1];  
+        }        
+        
+        
+        //clean addback
+        if (ge[i].suppress == FALSE && ge[i].pileup == FALSE && ge[i].nonprompt == FALSE) {
+            ge[i].clean=TRUE; //maybe clean should not include suppress == FALSE?
+        }    
+        
+        //Ge spectra
+        if (ge[i].energy > 0 && ge[i].energy < 8192) ge_spe[ge[i].id][ge[i].energy]++;         
+        if (ge[i].energy > 0 && ge[i].energy < 8192 && ge[i].clean == TRUE) ge_spe_clean[ge[i].id][ge[i].energy]++;
+        
+        //copy data and increment counters
+        ge[gmult]=ge[i];
+        gmult++;        
+        gcount++;
+        
+      } //end G
+
+     
+
+      ///////////////////////////
+      // END PROCESS DETECTORS //
+      ///////////////////////////
+
+
+
+      ////////////////////////
+      // FINAL USER SPECTRA //
+      ////////////////////////
+        
+        
+      //gamma-gamma time and energy 
+      for (i=0; i<gmult; i++) {
+        for (j=0; j<gmult; j++) {
+            if (i!=j) {
+                tdif = ge[i].time - ge[j].time + 2000;
+                //time difference matrix    
+                if (tdif >= 0 && tdif < 4096 && ge[i].energy >= 0 && ge[i].energy < 4096 && ge[j].energy >= 0 && ge[j].energy < 4096) {
+                    if (ge[i].energy < ge[j].energy)
+                    gg_tdif[ge[i].energy][tdif]++;
+                    else 
+                    gg_tdif[ge[j].energy][tdif]++;
+                }   
+                //prompt gamma-gamma
+                if (tdif >= 1920 && tdif <= 2080) {
+                    if (ge[i].energy >= 0 && ge[i].energy < 4096 && ge[j].energy >= 0 && ge[j].energy < 4096) 
+                    gg_prompt[ge[i].energy][ge[j].energy]++;
+                }                   
+                //non-prompt gamma-gamma (mult by 0.2555)
+                if ( (tdif >= 1532 && tdif <= 1846) || (tdif >= 2138 && tdif <= 2452 ) ) {
+                    if (ge[i].energy >= 0 && ge[i].energy < 4096 && ge[j].energy >= 0 && ge[j].energy < 4096) 
+                    gg_nprompt[ge[i].energy][ge[j].energy]++;                
+                }    
+            }
+        }
+      }  
+
+
+      ////////////////////////////
+      // END FINAL USER SPECTRA //
+      ////////////////////////////
+
+      } //end argc >= 5 condition 
+
+      //event stats, print status every 10000 events
+      lle_div=lldiv(evtcount,10000);
+      if ( lle_div.rem == 0 && strstr(argv[1], "p") != NULL) {
+        fprpos = ftell(fpr);
+        tempf = (float)fprsize/(1024.*1024.*1024.);
+        printf("Total SubEvents: \x1B[32m%llu \x1B[31m(%d%% pileup)\x1B[0m\nTotal Events: \x1B[32m%llu (%.1f <mult>)\x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\n\033[3A\r", sevtcount, (int)((100*pileupcount)/sevtcount), evtcount, (float)sevtcount/(float)evtcount, (100*fprpos/fprsize), tempf);
+      }      
+
+          
+  } // end main while loop 
+  /////////////////////////
+  // END MAIN WHILE LOOP //
+  /////////////////////////
+  fprpos = ftell(fpr);
+  tempf = (float)fprsize/(1024.*1024.*1024.);
+  printf("Total SubEvents: \x1B[32m%llu \x1B[31m(%d%% pileup)\x1B[0m\nTotal Events: \x1B[32m%llu (%.1f <mult>)\x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\n\033[3A\r", sevtcount, (int)((100*pileupcount)/sevtcount), evtcount, (float)sevtcount/(float)evtcount, (100*fprpos/fprsize), tempf);
+           
+
+  
+  
+  
+  
+  
+  ////////////////////
+  // WRITE SPECTRA //
+  ///////////////////
+  printf("\n\n\n\nWriting Spectra to Disk ...\n");
+     
+  //Event Spectra   
+  write_data4(HIT, *hit, 4096, 2, overwrite);
+  write_data4(MULT, *mult, 4096, 1, overwrite);
+  write_data4(TDIFID, *tdifid, MAX_ID, 8192, overwrite);
+  write_data4(E_RAW, *e_raw, MAX_ID, 8192, overwrite);
+  write_data4(E_CAL, *e_cal, MAX_ID, 8192, overwrite);
+  write_data4(TEVT_RAW, *tevt_raw, MAX_ID, 8192, overwrite);
+  write_data4(TEVT_CAL, *tevt_cal, MAX_ID, 8192, overwrite);
+  write_data4(TCFD_RAW, *tcfd_raw, MAX_ID, 8192, overwrite);
+  write_data4(TDIF_RAW, *tdif_raw, MAX_ID, 4096, overwrite);
+  write_data4(TDIF_CAL, *tdif_cal, MAX_ID, 4096, overwrite);
+ // write_data4(TDIF_CAL0_ETHRESH, *tdif_cal0_ethresh, MAX_ID, 4096, overwrite);
+//  write_data4(IDID, *idid, MAX_ID, MAX_ID, overwrite);
+
+  //Detector Processed Spectra
+  //Ge
+//  write_data4(GE_BGO_TDIF, *ge_bgo_tdif, MAX_GE, 4096, overwrite);
+//  write_data4(GE_XTL_TDIF, *ge_xtl_tdif, MAX_GE, 4096, overwrite);
+//  write_data4(GE_XTL_TDIF_ETHRESH, *ge_xtl_tdif_ethresh, MAX_GE, 4096, overwrite);
+//  write_data4(GE_SPE_XTL, *ge_spe_xtl, MAX_GE*MAX_GE_XTL, 8192, overwrite);
+//  write_data4(GE_SPE, *ge_spe, MAX_GE, 8192, overwrite);
+//  write_data4(GE_SPE_CLEAN, *ge_spe_clean, MAX_GE, 8192, overwrite);
+  //trinity
+  write_data4(PID, *pid, 4096, 4096, overwrite);
+  write_data4(PID_EVSP, *pid_evsp, 4096, 4096, overwrite);
+  write_data4(PID_EVST, *pid_evst, 4096, 4096, overwrite);
+  write_data4(PID_EVSTAU, *pid_evstau, 4096, 4096, overwrite);
+  write_data4(PID_EVSR, *pid_evsr, 4096, 4096, overwrite);
+  write_data4(PID_TAUVSR, *pid_tauvsr, 4096, 4096, overwrite);
+  write_data4(GAGGDT, *gaggdt, 4096, 4096, overwrite);
+
+  for (i=0; i<500; i++) {
+  	if (strace[i]>0) sptrace[0][i] = strace[i]/dbcount;
+  }
+  write_data4(SPTRACE, *sptrace, 1, 4096, overwrite);
+
+  //Final User Spectra
+  //gamma-gamma
+//  write_data4(GG_TDIF, *gg_tdif, 4096, 4096, overwrite);
+ // write_data4(GG_PROMPT, *gg_prompt, 4096, 4096, overwrite);
+ // write_data4(GG_NPROMPT, *gg_nprompt, 4096, 4096, overwrite);
+
+  
+  fclose(fpr);
+  fclose(debugfile);
+  
+  return 0;
+}
+
+