From b922b4474abf9d5572b1c9f527ec5184d4042760 Mon Sep 17 00:00:00 2001 From: "Ryan@Ubuntu" Date: Tue, 14 Dec 2021 10:35:05 -0500 Subject: [PATCH] basic of the pixie2root.cpp --- .gitignore | 8 + Analyzer.C | 189 ++++ Analyzer.h | 115 +++ AutoFit.C | 2701 ++++++++++++++++++++++++++++++++++++++++++++++++ makefile | 11 +- mapping.h | 28 + pixie2root.cpp | 474 +++++++++ scan.c | 1208 ++++++++++++++++++++++ 8 files changed, 4731 insertions(+), 3 deletions(-) create mode 100644 Analyzer.C create mode 100644 Analyzer.h create mode 100644 AutoFit.C create mode 100644 mapping.h create mode 100644 pixie2root.cpp create mode 100644 scan.c 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 +#include +#include +#include +#include +#include +#include + +//############################################ 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 +#include +#include +#include + +#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 +#include +#include +#include +#include +#include + +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 SplitStrAF(std::string tempLine, std::string splitter, int shift = 0){ + + std::vector 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;pGetNpar(); + 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 energy, height, sigma, lowE, highE ; +vector 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 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 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 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 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 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 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 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 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 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 xPeakList; +vector yPeakList; +vector xBGList; +vector 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 +#include +#include +#include +#include + +#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\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 )\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 )\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 +#include +#include +#include +#include + +#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= 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=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=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= 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= 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 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= 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 )\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 )\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; +} + +