diff --git a/armory/AnalysisLib.h b/armory/AnalysisLib.h deleted file mode 100644 index 4ff335c..0000000 --- a/armory/AnalysisLib.h +++ /dev/null @@ -1,296 +0,0 @@ -#ifndef ANALYSIS_LIB_H -#define ANALYSIS_LIB_H - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -namespace AnalysisLib { - -std::vector SplitStr(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 - while( secStr.substr(0, 1) == " ") secStr = secStr.substr(1); - - ///check if secStr is end with space - while( secStr.back() == ' ') secStr = secStr.substr(0, secStr.size()-1); - - output.push_back(secStr); - ///printf(" |%s---\n", secStr.c_str()); - - }while(pos != std::string::npos ); - - return output; -}; - - -//************************************** TCutG - -TObjArray * LoadListOfTCut(TString fileName, TString cutName = "cutList"){ - - if( fileName == "" ) return nullptr; - - TObjArray * cutList = nullptr; - - TFile * fCut = new TFile(fileName); - bool isCutFileOpen = fCut->IsOpen(); - if(!isCutFileOpen) { - printf( "Failed to open rdt-cutfile 1 : %s\n" , fileName.Data()); - }else{ - cutList = (TObjArray *) fCut->FindObjectAny(cutName); - - if( cutList ){ - int numCut = cutList->GetEntries(); - printf("=========== found %d cutG in %s \n", numCut, fCut->GetName()); - - for(int i = 0; i < numCut ; i++){ - printf("cut name : %s , VarX: %s, VarY: %s, numPoints: %d \n", - cutList->At(i)->GetName(), - ((TCutG*)cutList->At(i))->GetVarX(), - ((TCutG*)cutList->At(i))->GetVarY(), - ((TCutG*)cutList->At(i))->GetN() - ); - } - } - } - - return cutList; -} - -TCutG * LoadSingleTCut( TString fileName, TString cutName = "cutEZ"){ - - if( fileName == "" ) return nullptr; - TCutG * cut = nullptr; - - TFile * fCut = new TFile(fileName); - bool isCutFileOpen = fCut->IsOpen(); - if( !isCutFileOpen) { - printf( "Failed to open E-Z cutfile : %s\n" , fileName.Data()); - }else{ - cut = (TCutG *) fCut->FindObjectAny(cutName); - if( cut != NULL ) { - printf("Found EZ cut| name : %s, VarX: %s, VarY: %s, numPoints: %d \n", - cut->GetName(), - cut->GetVarX(), - cut->GetVarY(), - cut->GetN() - ); - } - } - - return cut; -} - - -//************************************** Others -std::vector> combination(std::vector arr, int r){ - - std::vector> output; - - int n = arr.size(); - std::vector v(n); - std::fill(v.begin(), v.begin()+r, 1); - do { - //for( int i = 0; i < n; i++) { printf("%d ", v[i]); }; printf("\n"); - - std::vector temp; - for (int i = 0; i < n; ++i) { - if (v[i]) { - //printf("%.1f, ", arr[i]); - temp.push_back(arr[i]); - } - } - //printf("\n"); - - output.push_back(temp); - - } while (std::prev_permutation(v.begin(), v.end())); - - return output; -} - -double* sumMeanVar(std::vector data){ - - int n = data.size(); - double sum = 0; - for( int k = 0; k < n; k++) sum += data[k]; - double mean = sum/n; - double var = 0; - for( int k = 0; k < n; k++) var += pow(data[k] - mean,2); - - static double output[3]; - output[0] = sum; - output[1] = mean; - output[2] = var; - - return output; -} - -double* fitSlopeIntercept(std::vector dataX, std::vector dataY){ - - double * smvY = sumMeanVar(dataY); - double sumY = smvY[0]; - double meanY = smvY[1]; - - double * smvX = sumMeanVar(dataX); - double sumX = smvX[0]; - double meanX = smvX[1]; - double varX = smvX[2]; - - int n = dataX.size(); - double sumXY = 0; - for( int j = 0; j < n; j++) sumXY += dataX[j] * dataY[j]; - - double slope = ( sumXY - sumX * sumY/n ) / varX; - double intercept = meanY - slope * meanX; - - static double output[2]; - output[0] = slope; - output[1] = intercept; - - return output; - -} - -std::vector> FindMatchingPair(std::vector enX, std::vector enY){ - - //output[0] = fitEnergy; - //output[1] = refEnergy; - - int nX = enX.size(); - int nY = enY.size(); - - std::vector fitEnergy; - std::vector refEnergy; - - if( nX > nY ){ - - std::vector> output = combination(enX, nY); - - double * smvY = sumMeanVar(enY); - double sumY = smvY[0]; - double meanY = smvY[1]; - double varY = smvY[2]; - - double optRSquared = 0; - double absRSqMinusOne = 1; - int maxID = 0; - - for( int k = 0; k < (int) output.size(); k++){ - - double * smvX = sumMeanVar(output[k]); - double sumX = smvX[0]; - double meanX = smvX[1]; - double varX = smvX[2]; - - double sumXY = 0; - for( int j = 0; j < nY; j++) sumXY += output[k][j] * enY[j]; - - double rSq = abs(sumXY - sumX*sumY/nY)/sqrt(varX*varY); - - //for( int j = 0; j < nY ; j++){ printf("%.1f, ", output[k][j]); }; printf("| %.10f\n", rSq); - - if( abs(rSq-1) < absRSqMinusOne ) { - absRSqMinusOne = abs(rSq-1); - optRSquared = rSq; - maxID = k; - } - } - - fitEnergy = output[maxID]; - refEnergy = enY; - - printf(" R^2 : %.20f\n", optRSquared); - - //calculation fitting coefficient - //double * si = fitSlopeIntercept(fitEnergy, refEnergy); - //printf( " y = %.4f x + %.4f\n", si[0], si[1]); - - }else if( nX < nY ){ - - std::vector> output = combination(enY, nX); - - - double * smvX = sumMeanVar(enX); - double sumX = smvX[0]; - double meanX = smvX[1]; - double varX = smvX[2]; - - double optRSquared = 0; - double absRSqMinusOne = 1; - int maxID = 0; - - for( int k = 0; k < (int) output.size(); k++){ - - double * smvY = sumMeanVar(output[k]); - double sumY = smvY[0]; - double meanY = smvY[1]; - double varY = smvY[2]; - - double sumXY = 0; - for( int j = 0; j < nX; j++) sumXY += output[k][j] * enX[j]; - - double rSq = abs(sumXY - sumX*sumY/nX)/sqrt(varX*varY); - - //for( int j = 0; j < nX ; j++){ printf("%.1f, ", output[k][j]); }; printf("| %.10f\n", rSq); - - if( abs(rSq-1) < absRSqMinusOne ) { - absRSqMinusOne = abs(rSq-1); - optRSquared = rSq; - maxID = k; - } - } - - fitEnergy = enX; - refEnergy = output[maxID]; - printf(" R^2 : %.20f\n", optRSquared); - - }else{ - fitEnergy = enX; - refEnergy = enY; - - //if nX == nY, ther could be cases that only partial enX and enY are matched. - - } - - printf("fitEnergy = ");for( int k = 0; k < (int) fitEnergy.size() ; k++){ printf("%7.2f, ", fitEnergy[k]); }; printf("\n"); - printf("refEnergy = ");for( int k = 0; k < (int) refEnergy.size() ; k++){ printf("%7.2f, ", refEnergy[k]); }; printf("\n"); - - std::vector> haha; - haha.push_back(fitEnergy); - haha.push_back(refEnergy); - - return haha; - -} - - -} - -#endif \ No newline at end of file diff --git a/armory/AutoFit.C b/armory/AutoFit.C deleted file mode 100644 index 4b94781..0000000 --- a/armory/AutoFit.C +++ /dev/null @@ -1,2861 +0,0 @@ -/*************************************************** - * This is a root macro for Auto fitting - * - * Created by Tsz Leung (Ryan) TANG, around 2019. - * updated 01-04-2022 - * - * contact goluckyryan@gmail.com - ***************************************************/ - -#ifndef AutoFit_C -#define AutoFit_C - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -namespace AutoFit{ - -//Global fit paramaters -std::vector BestFitMean; -std::vector BestFitCount; -std::vector BestFitSigma; - -TString recentFitMethod; - -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("\n"); - printf("------- Utility : \n"); - printf(" SaveFitPara() - Save the inital/Best fit parameters.\n"); - printf(" ShowFitMerhod() - Show this menual.\n"); - printf("---------------------------------------------------------\n"); -} - -void AutoFit(){ - ShowFitMethod(); -} - -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 = std::max(0., (1+2*cos(ang))/3.); - double g = std::max(0., (1 - cos(ang) + sqrt(3)* sin(ang))/3.); - double b = std::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"); -} - -std::vector energy, height, sigma, lowE, highE ; -std::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()); - std::ifstream file; - file.open(fitParaFile.Data()); - - if( !file){ - printf("\ncannot read file : %s \n", fitParaFile.Data()); - return 0; - } - - while( file.good()) { - std::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()); - - std::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"); - - recentFitMethod = "fitGaussPol"; - - 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); - - BestFitCount.clear(); - BestFitMean.clear(); - BestFitSigma.clear(); - - BestFitCount.push_back(paraA[0]); - BestFitMean.push_back(paraA[1]); - BestFitSigma.push_back(paraA[2]); - -} - -//######################################## -//#### fit 2 gauss + pol-1 // not updated -//######################################## -std::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"); - - recentFitMethod = "fit2GaussP1"; - - std::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); - - - BestFitCount.clear(); - BestFitMean.clear(); - BestFitSigma.clear(); - - for( int i = 0; i < 2; i++){ - BestFitCount.push_back(paraA[3*i]/ bw); - BestFitMean.push_back(paraA[3*i+1]); - BestFitSigma.push_back(paraA[3*i+2]); - } - - 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"); - - recentFitMethod = "fitGF3Pol"; - - 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 - std::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 -//############################################## -std::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"); - - recentFitMethod = "fitAuto"; - - 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 ); - std::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); - - std::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")); - - BestFitCount.clear(); - BestFitMean.clear(); - BestFitSigma.clear(); - - 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])); - - BestFitCount.push_back(paraA[3*i]/ bw); - BestFitMean.push_back(paraA[3*i+1]); - BestFitSigma.push_back(paraA[3*i+2]); - } - - return exPos; - -} - - -//######################################## -//###### NOT tested -//######################################## -std::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(); - } - - recentFitMethod = "fitNGF3"; - - 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 ); - std::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); - - std::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"); - - recentFitMethod = "fitNGauss"; - - 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")); - - BestFitCount.clear(); - BestFitMean.clear(); - BestFitSigma.clear(); - - 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])); - - BestFitCount.push_back(paraA[3*i]/ bw); - BestFitMean.push_back(paraA[3*i+1]); - BestFitSigma.push_back(paraA[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"); - - recentFitMethod = "fitNGaussSub"; - - 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")); - - BestFitCount.clear(); - BestFitMean.clear(); - BestFitSigma.clear(); - - 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])); - - BestFitCount.push_back(paraA[3*i]/ bw); - BestFitMean.push_back(paraA[3*i+1]); - BestFitSigma.push_back(paraA[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"); - - recentFitMethod = "fitNGaussPol"; - - 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")); - - BestFitCount.clear(); - BestFitMean.clear(); - BestFitSigma.clear(); - - 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])); - - BestFitCount.push_back(paraA[3*i]/ bw); - BestFitMean.push_back(paraA[3*i+1]); - BestFitSigma.push_back(paraA[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"); - - recentFitMethod = "fitNGaussPolSub"; - - 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")); - - BestFitCount.clear(); - BestFitMean.clear(); - BestFitSigma.clear(); - - 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])); - - BestFitCount.push_back(paraA[3*i]/ bw); - BestFitMean.push_back(paraA[3*i+1]); - BestFitSigma.push_back(paraA[3*i+2]); - } - - -} - -//################################################# -//#### fit N-Gauss with pol-n BG using mouse click -//################################################# - -int nClick = 0; -bool peakFlag = 1; -std::vector xPeakList; -std::vector yPeakList; -std::vector xBGList; -std::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(bool isBestFit = true, TString fileName = "AutoFit_para.txt"){ - - if( xPeakList.size() == 0 && BestFitMean.size() == 0 ){ - printf(" no fit paramaters availible. \n"); - return; - } - - if( recentFitMethod == "fitGF3Pol" || recentFitMethod == "fitGF3" ){ - printf(" Not support for GF3 fitting. \n"); - return; - } - - 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"); - - - if ( xPeakList.size() == 0 || isBestFit ){ - for( int i = 0 ; i < BestFitMean.size() ; i++){ - fprintf(file_out, "%.3f %.3f %.3f 0 %.3f 0 %.0f\n", - BestFitMean[i], - BestFitMean[i] - 5*BestFitSigma[i], - BestFitMean[i] + 5*BestFitSigma[i], - BestFitSigma[i], - BestFitCount[i]); - } - }else{ - 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, double meanRange = 0){ - - printf("=========================================================\n"); - printf("======== fit n-Gauss + Pol-%d BG using mouse click =====\n", degPol ); - printf("==========================================================\n"); - - recentFitMethod = "clickFitNGaussPol"; - - 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); - - if( meanRange <= 0 ) { - 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 ); - }else{ - fit->SetParLimits(3*i+1, xPeakList[i] - meanRange/2. , xPeakList[i] + meanRange/2. ); - } - 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")); - - BestFitCount.clear(); - BestFitMean.clear(); - BestFitSigma.clear(); - - 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])); - - BestFitCount.push_back(paraA[3*i]/ bw); - BestFitMean.push_back(paraA[3*i+1]); - BestFitSigma.push_back(paraA[3*i+2]); - } - -} - - - -void clickFitNGaussPolSub(TH1F * hist, int degPol, double sigmaMax = 0, double meanRange = 0){ - - printf("=========================================================\n"); - printf("= fit n-Gauss + Pol-%d BG using mouse click (method-2) =\n", degPol ); - printf("==========================================================\n"); - - recentFitMethod = "clickFitNGaussPolSub"; - - 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); - - if( meanRange <= 0 ) { - 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 ); - }else{ - fit->SetParLimits(3*i+1, xPeakList[i] - meanRange/2. , xPeakList[i] + meanRange/2. ); - } - 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")); - - BestFitCount.clear(); - BestFitMean.clear(); - BestFitSigma.clear(); - - 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])); - - BestFitCount.push_back(paraA[3*i]/ bw); - BestFitMean.push_back(paraA[3*i+1]); - BestFitSigma.push_back(paraA[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/armory/EventBuilder.cpp b/armory/EventBuilder.cpp deleted file mode 100644 index 13a0f0f..0000000 --- a/armory/EventBuilder.cpp +++ /dev/null @@ -1,351 +0,0 @@ -#include "SolReader.h" -#include -#include - -#include "TFile.h" -#include "TTree.h" -#include "TMath.h" -#include "TString.h" -#include "TMacro.h" -//#include "TClonesArray.h" // plan to save trace as TVector with TClonesArray -//#include "TVector.h" - -#define MAX_MULTI 64 -#define MAX_TRACE_LEN 2500 - -#define tick2ns 8 // 1 tick = 8 ns - -SolReader ** reader; -Hit ** hit; - -std::vector> idList; - -unsigned long totFileSize = 0; -unsigned long processedFileSize = 0; - -std::vector activeFileID; -std::vector groupIndex; -std::vector> group; // group[i][j], i = group ID, j = group member) - -void findEarliestTime(int &fileID, int &groupID){ - - unsigned long firstTime = 0; - for( int i = 0; i < (int) activeFileID.size(); i++){ - int id = activeFileID[i]; - if( i == 0 ) { - firstTime = hit[id]->timestamp; - fileID = id; - groupID = i; - //printf("%d | %d %lu %d | %d \n", id, reader[id]->GetBlockID(), hit[id]->timestamp, hit[id]->channel, (int) activeFileID.size()); - continue; - } - if( hit[id]->timestamp <= firstTime) { - firstTime = hit[id]->timestamp; - fileID = id; - groupID = i; - //printf("%d | %d %lu %d | %d \n", id, reader[id]->GetBlockID(), hit[id]->timestamp, hit[id]->channel, (int) activeFileID.size()); - } - } - -} - -unsigned long long evID = 0; -unsigned int multi = 0; -unsigned short bd[MAX_MULTI] = {0}; -unsigned short sn[MAX_MULTI] = {0}; -unsigned short ch[MAX_MULTI] = {0}; -unsigned short e[MAX_MULTI] = {0}; -unsigned short e2[MAX_MULTI] = {0}; //for PSD energy short -unsigned long long e_t[MAX_MULTI] = {0}; -unsigned short e_f[MAX_MULTI] = {0}; -unsigned short lowFlag[MAX_MULTI] = {0}; -unsigned short highFlag[MAX_MULTI] = {0}; -int traceLen[MAX_MULTI] = {0}; -int trace[MAX_MULTI][MAX_TRACE_LEN] = {0}; - -void fillData(int &fileID, const bool &saveTrace){ - bd[multi] = idList[fileID][1]; - sn[multi] = idList[fileID][3]; - ch[multi] = hit[fileID]->channel; - e[multi] = hit[fileID]->energy; - e2[multi] = hit[fileID]->energy_short; - e_t[multi] = hit[fileID]->timestamp; - e_f[multi] = hit[fileID]->fine_timestamp; - lowFlag[multi] = hit[fileID]->flags_low_priority; - highFlag[multi] = hit[fileID]->flags_high_priority; - - if( saveTrace ){ - traceLen[multi] = hit[fileID]->traceLenght; - for( int i = 0; i < TMath::Min(traceLen[multi], MAX_TRACE_LEN); i++){ - trace[multi][i] = hit[fileID]->analog_probes[0][i]; - } - } - - multi++; - reader[fileID]->ReadNextBlock(); -} - -void printEvent(){ - printf("==================== evID : %llu\n", evID); - for( int i = 0; i < multi; i++){ - printf(" %2d | %d %d | %llu %d \n", i, bd[i], ch[i], e_t[i], e[i] ); - } - printf("==========================================\n"); -} - -//^################################################################################## -int main(int argc, char ** argv){ - - printf("=======================================================\n"); - printf("=== SOLARIS Event Builder sol --> root ===\n"); - printf("=======================================================\n"); - - if( argc <= 3){ - printf("%s [outfile] [timeWindow] [saveTrace] [sol-1] [sol-2] ... \n", argv[0]); - printf(" outfile : output root file name\n"); - printf(" timeWindow : number of tick, 1 tick = %d ns.\n", tick2ns); - printf(" saveTrace : 1 = save trace, 0 = no trace\n"); - printf(" sol-X : the sol file(s)\n"); - return -1; - } - - // for( int i = 0; i < argc; i++){ - // printf("%d | %s\n", i, argv[i]); - // } - - TString outFileName = argv[1]; - int timeWindow = abs(atoi(argv[2])); - const bool saveTrace = atoi(argv[3]); - - const int nFile = argc - 4; - TString inFileName[nFile]; - for( int i = 0 ; i < nFile ; i++){ - inFileName[i] = argv[i+4]; - } - - //*======================================== setup reader - reader = new SolReader*[nFile]; - hit = new Hit *[nFile]; - - for( int i = 0 ; i < nFile ; i++){ - reader[i] = new SolReader(inFileName[i].Data()); - hit[i] = reader[i]->hit; //TODO check is file open propertly - reader[i]->ReadNextBlock(); // read the first block - } - - //*======================================== group files - idList.clear(); - for( int i = 0; i < nFile; i++){ - TString fn = inFileName[i]; - - int pos = fn.Last('/'); // path - fn.Remove(0, pos+1); - - pos = fn.First('_'); // expName; - fn.Remove(0, pos+1); - - pos = fn.First('_'); // runNum; - fn.Remove(0, pos+1); - - pos = fn.First('_'); // digiID - TString f1 = fn; - int digiID = f1.Remove(pos).Atoi(); - fn.Remove(0, pos+1); - - pos = fn.Last('_'); // digi serial num - f1 = fn; - int digisn = f1.Remove(pos).Atoi(); - fn.Remove(0, pos+1); - - pos = fn.First('.'); // get the file id; - int indexID = fn.Remove(pos).Atoi(); - - int fileID = i; - std::vector haha = {fileID, digiID, indexID, digisn}; - idList.push_back(haha); - } - - // sort by digiID - std::sort(idList.begin(), idList.end(), [](const std::vector& a, const std::vector& b){ - if (a[1] == b[1]) { - return a[2] < b[2]; - } - return a[1] < b[1]; - }); - - group.clear(); // group[i][j], i is the group Index = digiID - int last_id = 0; - std::vector kaka; - for( int i = 0; i < (int) idList.size() ; i++){ - if( i == 0 ) { - kaka.clear(); - last_id = idList[i][1]; - kaka.push_back(idList[i][0]); - continue; - } - - if( idList[i][1] != last_id ) { - last_id = idList[i][1]; - group.push_back(kaka); - kaka.clear(); - kaka.push_back(idList[i][0]); - }else{ - kaka.push_back(idList[i][0]); - } - } - group.push_back(kaka); - - printf(" out file : \033[1;33m%s\033[m\n", outFileName.Data()); - printf(" Event building time window : %d tics = %d nsec \n", timeWindow, timeWindow*tick2ns); - printf(" Save Trace ? %s \n", saveTrace ? "Yes" : "No"); - printf(" Number of input file : %d \n", nFile); - for( int i = 0; i < nFile; i++){ - printf(" %2d| %5.1f MB| %s \n", i, reader[i]->GetFileSize()/1024./1024., inFileName[i].Data()); - totFileSize += reader[i]->GetFileSize(); - } - printf("------------------------------------\n"); - for( int i = 0; i < (int) group.size(); i++){ - printf("Group %d :", i); - for( int j = 0; j < (int) group[i].size(); j ++){ - printf("%d, ", group[i][j]); - } - printf("\n"); - } - printf("------------------------------------\n"); - - //*======================================== setup tree - TFile * outRootFile = new TFile(outFileName, "recreate"); - outRootFile->cd(); - - TTree * tree = new TTree("tree", outFileName); - - tree->Branch("evID", &evID, "event_ID/l"); - tree->Branch("multi", &multi, "multi/i"); - tree->Branch("bd", bd, "board[multi]/s"); - tree->Branch("sn", sn, "sn[multi]/s"); - tree->Branch("ch", ch, "channel[multi]/s"); - tree->Branch("e", e, "energy[multi]/s"); - tree->Branch("e2", e2, "energy_short[multi]/s"); - tree->Branch("e_t", e_t, "timestamp[multi]/l"); - tree->Branch("e_f", e_f, "fine_timestamp[multi]/s"); - tree->Branch("lowFlag", lowFlag, "lowFlag[multi]/s"); - tree->Branch("highFlag", highFlag, "highFlag[multi]/s"); - - if( saveTrace){ - tree->Branch("tl", traceLen, "traceLen[multi]/I"); - tree->Branch("trace", trace, Form("trace[multi][%d]/I", MAX_TRACE_LEN)); - } - - //*=========================================== build event - - //@---- using file from group[i][0] first - - //--- find earlist time among the files - activeFileID.clear(); - groupIndex.clear(); //the index of each group - - for(int i = 0; i < (int) group.size(); i++) { - groupIndex.push_back(0); - activeFileID.push_back(group[i][0]); - } - - int fileID = 0; - int groupID = 0; - findEarliestTime(fileID, groupID); - fillData(fileID, saveTrace); - - unsigned long firstTimeStamp = hit[fileID]->timestamp; - unsigned long lastTimeStamp = 0; - - int last_precentage = 0; - while((activeFileID.size() > 0)){ - - findEarliestTime(fileID, groupID); - if( reader[fileID]->IsEndOfFile() ){ - groupIndex[groupID] ++; - if( groupIndex[groupID] < (int) group[groupID].size() ){ - activeFileID[groupID] = group[groupID][groupIndex[groupID]]; - fileID = activeFileID[groupID]; - }else{ - activeFileID.erase(activeFileID.begin() + groupID); - } - } - - if( hit[fileID]->timestamp - e_t[0] < timeWindow ){ - fillData(fileID, saveTrace); - }else{ - outRootFile->cd(); - tree->Fill(); - evID ++; - - multi = 0; - fillData(fileID, saveTrace); - } - - ///========= calculate progress - processedFileSize = 0; - for( int p = 0; p < (int) group.size(); p ++){ - for( int q = 0; q <= groupIndex[p]; q++){ - if( groupIndex[p] < (int) group[p].size() ){ - int id = group[p][q]; - processedFileSize += reader[id]->GetFilePos(); - } - } - } - double percentage = processedFileSize * 100/ totFileSize; - if( percentage >= last_precentage ) { - printf("Processed : %llu, %.0f%% | %lu/%lu | ", evID, percentage, processedFileSize, totFileSize); - for( int i = 0; i < (int) activeFileID.size(); i++) printf("%d, ", activeFileID[i]); - printf(" \n\033[A\r"); - last_precentage = percentage + 1.0; - } - }; ///====== end of event building loop - - processedFileSize = 0; - for( int p = 0; p < (int) group.size(); p ++){ - for( int q = 0; q < (int) group[p].size(); q++){ - int id = group[p][q]; - processedFileSize += reader[id]->GetFilePos(); - } - } - double percentage = processedFileSize * 100/ totFileSize; - printf("Processed : %llu, %.0f%% | %lu/%lu \n", evID, percentage, processedFileSize, totFileSize); - - lastTimeStamp = hit[fileID]->timestamp; - //*=========================================== save file - outRootFile->cd(); - tree->Fill(); - evID ++; - tree->Write(); - - //*=========================================== Save timestamp as TMacro - TMacro timeStamp; - TString str; - str.Form("%lu", firstTimeStamp); timeStamp.AddLine( str.Data() ); - str.Form("%lu", lastTimeStamp); timeStamp.AddLine( str.Data() ); - timeStamp.Write("timeStamp"); - - unsigned int numBlock = 0; - for( int i = 0; i < nFile; i++){ - //printf("%d | %8ld | %10u/%10u\n", i, reader[i]->GetBlockID() + 1, reader[i]->GetFilePos(), reader[i]->GetFileSize()); - numBlock += reader[i]->GetBlockID() + 1; - } - - printf("===================================== done. \n"); - printf("Number of Block Scanned : %u\n", numBlock); - printf(" Number of Event Built : %lld\n", evID); - printf(" Output Root File Size : %.2f MB\n", outRootFile->GetSize()/1024./1024.); - printf(" first timestamp : %lu \n", firstTimeStamp); - printf(" last timestamp : %lu \n", lastTimeStamp); - unsigned long duration = lastTimeStamp - firstTimeStamp; - printf(" total duration : %lu = %.2f sec \n", duration, duration * tick2ns * 1.0 / 1e9 ); - printf("===================================== end of summary. \n"); - - - //^############## delete new - for( int i = 0; i < nFile; i++) delete reader[i]; - delete [] reader; - outRootFile->Close(); - - return 0; -} \ No newline at end of file diff --git a/armory/GeneralSort.C b/armory/GeneralSort.C deleted file mode 100644 index b8729a9..0000000 --- a/armory/GeneralSort.C +++ /dev/null @@ -1,205 +0,0 @@ -#define GeneralSort_cxx - -#include -#include -#include - -#include "GeneralSort.h" - -#include -#include -#include -#include -#include - -Long64_t processedEntry = 0; -float lastPercentage = 0; - -//^############################################################## -Bool_t GeneralSort::Process(Long64_t entry){ - - if( entry < 1 ) printf("============================== start processing data\n"); - - ///initialization - for( int i = 0; i < mapping::nDetType; i++){ - for( int j = 0; j < mapping::detNum[i]; j++){ - eE[i][j] = TMath::QuietNaN(); - eT[i][j] = 0; - - if( isTraceExist && traceMethod > 0){ - teE[i][j] = TMath::QuietNaN(); - teT[i][j] = TMath::QuietNaN(); - teR[i][j] = TMath::QuietNaN(); - } - } - } - - multi = 0; - b_event_ID->GetEntry(entry); - b_multi->GetEntry(entry); - b_bd->GetEntry(entry); - b_ch->GetEntry(entry); - b_e->GetEntry(entry); - b_e_t->GetEntry(entry); - - for( int i = 0 ; i < multi; i++){ - int detID = mapping::map[bd[i]][ch[i]]; - int detType = mapping::FindDetTypeIndex(detID); - int low = (i == 0 ? 0 : mapping::detMaxID[detType-1]); - int reducedDetID = detID - low; - - eE[detType][reducedDetID] = e[i] * mapping::detParity[detType]; - eT[detType][reducedDetID] = e_t[i]; - - } - - //@===================================== Trace - if( isTraceExist && traceMethod >= 0 ){ - - b_tl->GetEntry(entry); - b_trace->GetEntry(entry); - - int countTrace = 0; - - arr->Clear("C"); - - for( int i = 0; i < multi; i++){ - int detID = mapping::map[bd[i]][ch[i]]; - - - int traceLength = tl[i]; - gTrace = (TGraph*) arr->ConstructedAt(countTrace, "C"); - gTrace->Clear(); - gTrace->Set(traceLength); - - gTrace->SetTitle(Form("ev:%llu,nHit:%d,id:%d,len:%d", evID, i, detID, traceLength)); - countTrace ++; - - for( int k = 0 ; k < traceLength; k++){ - gTrace->SetPoint(k, k, trace[i][k]); - } - - //***=================== fit - if( traceMethod == 1){ - - int detType = mapping::FindDetTypeIndex(detID); - //TODO use a blackList - //if( mapping::detTypeName[detType] != "rdt") continue; - - //TODO try custom build fiting algorithm. May be faster? - gFit = new TF1("gFit", fitFunc, 0, traceLength, numPara); - gFit->SetLineColor(6); - gFit->SetRange(0, traceLength); - - gFit->SetParameter(0, e[i]); - gFit->SetParameter(1, 100); //triggerTime //TODO how to not hardcode? - gFit->SetParameter(2, 10); //rise time //TODO how to not hardcode? - gFit->SetParameter(3, trace[i][0]); //base line - gFit->SetParameter(4, 100); // decay //TODO how to not hardcode? - gFit->SetParameter(5, -0.01); // pre-rise slope //TODO how to not hardcode? - - gFit->SetParLimits(1, 85, 125); //raneg for the trigger time - gFit->SetParLimits(5, -2, 0); - - gTrace->Fit("gFit", "QR", "", 0, traceLength); - - int low = (i == 0 ? 0 : mapping::detMaxID[detType-1]); - int reducedDetID = detID - low; - - teE[detType][reducedDetID] = gFit->GetParameter(0); - teT[detType][reducedDetID] = gFit->GetParameter(1); - teR[detType][reducedDetID] = gFit->GetParameter(2); - - delete gFit; - gFit = nullptr; - } - - //***=================== Trapezoid filter - if( traceMethod == 2){ - //TODO - } - - } - - } - - if( !isParallel){ - processedEntry ++; - float percentage = processedEntry*100/NumEntries; - if( percentage >= lastPercentage ) { - TString msg; msg.Form("%lu", NumEntries); - int len = msg.Sizeof(); - printf("Processed : %*lld, %3.0f%% | Elapsed %6.1f sec | expect %6.1f sec\n\033[A\r", len, entry, percentage, stpWatch.RealTime(), stpWatch.RealTime()/percentage*100); - stpWatch.Start(kFALSE); - lastPercentage = percentage + 1.0; - if( lastPercentage >= 100) printf("\n"); - - } - } - - newSaveTree->Fill(); - - return kTRUE; -} - -//^############################################################## -void GeneralSort::Terminate(){ - - printf("=============================== %s\n", __func__); - - DecodeOption(); - - if( !isParallel){ - stpWatch.Start(kFALSE); - saveFile->cd(); - newSaveTree->Print("toponly"); - newSaveTree->Write(); - saveFile->Close(); - } - - //get entries - saveFile = TFile::Open(saveFileName); - if( saveFile->IsOpen() ){ - TTree * tree = (TTree*) saveFile->FindObjectAny("gen_tree"); - int validCount = tree->GetEntries(); - - saveFile->Close(); - - printf("=========================================================================\n"); - PrintTraceMethod(); - printf("----- saved as \033[1;33m%s\033[0m. valid event: %d\n", saveFileName.Data() , validCount); - printf("=========================================================================\n"); - } -} - -//^############################################################## -void GeneralSort::Begin(TTree * tree){ - - printf( "=================================================================\n"); - printf( "===================== SOLARIS GeneralSort.C =================\n"); - printf( "=================================================================\n"); - - mapping::PrintMapping(); - - DecodeOption(); - if(!isParallel) { - tree->GetEntriesFast(); - stpWatch.Start(); - } - -} - -void GeneralSort::SlaveBegin(TTree * /*tree*/){ - -} - -void GeneralSort::SlaveTerminate(){ - if( isParallel){ - printf("%s::SaveTree\n", __func__); - saveFile->cd(); - newSaveTree->Write(); - fOutput->Add(proofFile); - saveFile->Close(); - } - -} \ No newline at end of file diff --git a/armory/GeneralSort.h b/armory/GeneralSort.h deleted file mode 100644 index 632ed7f..0000000 --- a/armory/GeneralSort.h +++ /dev/null @@ -1,362 +0,0 @@ -#ifndef GeneralSort_h -#define GeneralSort_h - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -//^######################################### Skip list for trace fitting -//TODO - -/*********************************======= - -the sequence of each method - -1) Begin (master, stdout) - 2) SlaveBegin (slave) - 3) Init - 4) Notify - 5) Process (looping each event) - 6) SlaveTerminate -7) Terminate - -// variable in the Master will not pass to Slave - -******************************************/ - -/// in Process_Sort, copy the Mapping.h to ~/.proof/working/ -#include "../working/Mapping.h" - -//^######################################### FIT FUNCTION -const int numPara = 6; -double fitFunc(double * x, double * par){ - /// par[0] = A - /// par[1] = t0 - /// par[2] = rise time - /// par[3] = baseline - /// par[4] = decay - /// par[5] = pre-rise slope - - if( x[0] < par[1] ) return par[3] + par[5] * (x[0]-par[1]); - - return par[3] + par[0] * (1 - TMath::Exp(- (x[0] - par[1]) / par[2]) ) * TMath::Exp(- (x[0] - par[1]) / par[4]); -} - -//^######################################### TRAPEZOID -TGraph * TrapezoidFilter(TGraph * trace){ - ///Trapezoid filter https://doi.org/10.1016/0168-9002(94)91652-7 - - //TODO how to not hard code? - int baseLineEnd = 80; - int riseTime = 10; //ch - int flatTop = 20; - float decayTime = 2000; - - TGraph * trapezoid = new TGraph(); - trapezoid->Clear(); - - ///find baseline; - double baseline = 0; - for( int i = 0; i < baseLineEnd; i++){ - baseline += trace->Eval(i); - } - baseline = baseline*1./baseLineEnd; - - int length = trace->GetN(); - - double pn = 0.; - double sn = 0.; - for( int i = 0; i < length ; i++){ - - double dlk = trace->Eval(i) - baseline; - if( i - riseTime >= 0 ) dlk -= trace->Eval(i - riseTime) - baseline; - if( i - flatTop - riseTime >= 0 ) dlk -= trace->Eval(i - flatTop - riseTime) - baseline; - if( i - flatTop - 2*riseTime >= 0) dlk += trace->Eval(i - flatTop - 2*riseTime) - baseline; - - if( i == 0 ){ - pn = dlk; - sn = pn + dlk*decayTime; - }else{ - pn = pn + dlk; - sn = sn + pn + dlk*decayTime; - } - trapezoid->SetPoint(i, i, sn / decayTime / riseTime); - } - return trapezoid; -} - -TStopwatch stpWatch; - -//^######################################### Class definition -// Header file for the classes stored in the TTree if any. -class GeneralSort : 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; - Int_t multi; - Int_t bd[100]; //[multi] - Int_t ch[100]; //[multi] - Int_t e[100]; //[multi] - ULong64_t e_t[100]; //[multi] - UShort_t lowFlag[100]; //[multi] - UShort_t highFlag[100]; //[multi] - Int_t tl[100]; //[multi] - Int_t trace[100][2500]; //[multi] - - // List of branches - TBranch *b_event_ID; //! - TBranch *b_multi; //! - TBranch *b_bd; //! - TBranch *b_ch; //! - TBranch *b_e; //! - TBranch *b_e_t; //! - TBranch *b_lowFlag; //! - TBranch *b_highFlag; //! - TBranch *b_tl; //! - TBranch *b_trace; //! - - GeneralSort(TTree * /*tree*/ =0) : fChain(0) { - printf("constructor :: %s\n", __func__); - - isTraceExist = false; - traceMethod = 0; // -1 = ignore trace, 0 = no trace fit, 1 = fit, 2 = trapezoid - - isParallel = false; - - saveFile = nullptr; - newSaveTree = nullptr; - - eE = nullptr; - eT = nullptr; - - arr = nullptr; - gTrace = nullptr; - gFit = nullptr; - arrTrapezoid = nullptr; - gTrapezoid = nullptr; - - teE = nullptr; - teT = nullptr; - teR = nullptr; - - } - virtual ~GeneralSort() { } - 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(GeneralSort,0); - - bool isTraceExist; - int traceMethod; - - void SetTraceMethod(int methodID) {traceMethod = methodID;} - void PrintTraceMethod(); - - void SetUpTree(); - void DecodeOption(); - bool isParallel; - - unsigned long NumEntries; - - TString saveFileName; - TFile * saveFile; //! - TProofOutputFile * proofFile; //! - TTree * newSaveTree; //! - - Float_t ** eE; //! - ULong64_t ** eT; //! - - //trace - TClonesArray * arr ;//! - TGraph * gTrace; //! - TF1 * gFit; //! - TClonesArray * arrTrapezoid ;//! - TGraph * gTrapezoid; //! - - //trace energy, trigger time, rise time - Float_t **teE; //! - Float_t **teT; //! - Float_t **teR; //! - -}; - -#endif - -#ifdef GeneralSort_cxx - -//^############################################################## -void GeneralSort::SetUpTree(){ - - printf("%s\n", __func__); - - if( isParallel){ - proofFile = new TProofOutputFile(saveFileName, "M"); - saveFile = proofFile->OpenFile("RECREATE"); - }else{ - saveFile = new TFile(saveFileName,"RECREATE"); - } - - newSaveTree = new TTree("gen_tree", "Tree After GeneralSort"); - newSaveTree->SetDirectory(saveFile); - newSaveTree->AutoSave(); - - newSaveTree->Branch( "evID", &evID, "EventID/l"); // simply copy - - eE = new Float_t * [mapping::nDetType]; - eT = new ULong64_t * [mapping::nDetType]; - - for( int i = 0 ; i < mapping::nDetType; i++){ - eE[i] = new Float_t[mapping::detNum[i]]; - eT[i] = new ULong64_t[mapping::detNum[i]]; - - for( int j = 0; j < mapping::detNum[i]; j++){ - eE[i][j] = TMath::QuietNaN(); - eT[i][j] = 0; - } - - newSaveTree->Branch( mapping::detTypeName[i].c_str(), eE[i], Form("%s[%d]/F", mapping::detTypeName[i].c_str(), mapping::detNum[i])); - newSaveTree->Branch( (mapping::detTypeName[i]+"_t").c_str(), eT[i], Form("%s_Timestamp[%d]/l", mapping::detTypeName[i].c_str(), mapping::detNum[i])); - } - - - if( isTraceExist && traceMethod >= 0){ - - arr = new TClonesArray("TGraph"); - - newSaveTree->Branch("trace", arr, 256000); - arr->BypassStreamer(); - - if( traceMethod > 0 ){ - - teE = new Float_t * [mapping::nDetType]; - teT = new Float_t * [mapping::nDetType]; - teR = new Float_t * [mapping::nDetType]; - - for( int i = 0 ; i < mapping::nDetType; i++){ - teE[i] = new Float_t[mapping::detNum[i]]; - teT[i] = new Float_t[mapping::detNum[i]]; - teR[i] = new Float_t[mapping::detNum[i]]; - - for( int j = 0; j Branch( ("w" + mapping::detTypeName[i]).c_str(), teE[i], Form("trace_%s[%d]/F", mapping::detTypeName[i].c_str(), mapping::detNum[i])); - newSaveTree->Branch( ("w" + mapping::detTypeName[i]+"T").c_str(), teT[i], Form("trace_%s_time[%d]/l", mapping::detTypeName[i].c_str(), mapping::detNum[i])); - newSaveTree->Branch( ("w" + mapping::detTypeName[i]+"R").c_str(), teR[i], Form("trace_%s_rise[%d]/l", mapping::detTypeName[i].c_str(), mapping::detNum[i])); - } - - } - - } - newSaveTree->Print("toponly"); //very important, otherwise the mac will blow up. -} - -//^############################################################## -void GeneralSort::DecodeOption(){ - TString option = GetOption(); - if( option != ""){ - TObjArray * tokens = option.Tokenize(","); - traceMethod = ((TObjString*) tokens->At(0))->String().Atoi(); - saveFileName = ((TObjString*) tokens->At(1))->String(); - isParallel = ((TObjString*) tokens->At(2))->String().Atoi(); - }else{ - traceMethod = -1; - saveFileName = "temp.root"; - isParallel = false; - } - - printf("|%s| %d %s %d \n", option.Data(), traceMethod, saveFileName.Data(), isParallel); - -} - -//^############################################################## -void GeneralSort::Init(TTree *tree){ - - // Set branch addresses and branch pointers - if (!tree) return; - fChain = tree; - fChain->SetMakeClass(1); - - fChain->SetBranchAddress("evID", &evID, &b_event_ID); - fChain->SetBranchAddress("multi", &multi, &b_multi); - fChain->SetBranchAddress("bd", bd, &b_bd); - fChain->SetBranchAddress("ch", ch, &b_ch); - fChain->SetBranchAddress("e", e, &b_e); - fChain->SetBranchAddress("e_t", e_t, &b_e_t); - fChain->SetBranchAddress("lowFlag", lowFlag, &b_lowFlag); - fChain->SetBranchAddress("highFlag", highFlag, &b_highFlag); - - TBranch * br = (TBranch *) fChain->GetListOfBranches()->FindObject("tl"); - if( br == NULL ){ - printf(" ++++++++ no Trace.\n"); - isTraceExist = false; - }else{ - printf(" ++++++++ Found Trace.\n"); - isTraceExist = true; - fChain->SetBranchAddress("tl", tl, &b_tl); - fChain->SetBranchAddress("trace", trace, &b_trace); - } - - NumEntries = fChain->GetEntries(); - printf( " ========== total Entry : %ld\n", NumEntries); - - //########################### Get Option - DecodeOption(); - - if( isTraceExist ){ - PrintTraceMethod(); - }else{ - printf("++++++++ no Trace found\n"); - } - - SetUpTree(); - - printf("---- end of Init %s\n ", __func__); - -} - -Bool_t GeneralSort::Notify(){ - return kTRUE; -} - -void GeneralSort::PrintTraceMethod(){ - const char* traceMethodStr; - switch(traceMethod) { - case -1 : traceMethodStr = "Ignore Trace"; break; - case 0 : traceMethodStr = "Copy"; break; - case 1 : traceMethodStr = "Fit"; break; - case 2 : traceMethodStr = "Trapezoid"; break; - default: traceMethodStr = "Unknown"; break; - } - printf("\033[1;33m ===== Trace method ? %s \033[m\n", traceMethodStr); -} - -#endif // #ifdef GeneralSort_cxx diff --git a/armory/GeneralSortAgent.C b/armory/GeneralSortAgent.C deleted file mode 100644 index 63db067..0000000 --- a/armory/GeneralSortAgent.C +++ /dev/null @@ -1,56 +0,0 @@ - -#include "TTree.h" -#include "TProof.h" -#include "TChain.h" -#include "TMacro.h" -#include "TFile.h" - -void GeneralSortAgent(Int_t runNum, int nWorker = 1, int traceMethod = -1){ - - TString name; - name.Form("../root_data/run%03d.root", runNum); - - TChain * chain = new TChain("tree"); - chain->Add(name); - - chain->GetListOfFiles()->Print(); - - printf("\033[1;33m---------------------total number of Events %llu\033[0m\n", chain->GetEntries()); - - if( chain->GetEntries() == 0 ) return; - - //this is the option for TSelector, the first one is traceMethod, 2nd is save fileName; - TString option; - - if( abs(nWorker) == 1){ - - option.Form("%d,../root_data/gen_run%03d.root,%d", traceMethod, runNum, 0); - chain->Process("../armory/GeneralSort.C+", option); - - }else{ - - TProof * p = TProof::Open("", Form("workers=%d", abs(nWorker))); - p->ShowCache(); - printf("----------------------------\n"); - - chain->SetProof(); - option.Form("%d,../root_data/gen_run%03d.root,%d", traceMethod, runNum, 1); - chain->Process("../armory/GeneralSort.C+", option); - } - - //========== open the output root and copy teh timestamp Marco - - TFile * f1 = new TFile(Form("../root_data/run%03d.root", runNum), "READ"); - TMacro * m = (TMacro* ) f1->Get("timeStamp"); - m->AddLine(Form("%d", runNum)); - - TFile * f2 = new TFile(Form("../root_data/gen_run%03d.root", runNum), "UPDATE"); - f2->cd(); - m->Write("timeStamp"); - - f1->Close(); - f2->Close(); - - delete chain; - -} \ No newline at end of file diff --git a/armory/Hit.h b/armory/Hit.h deleted file mode 100644 index 1a945b1..0000000 --- a/armory/Hit.h +++ /dev/null @@ -1,287 +0,0 @@ -#ifndef HIT_H -#define HIT_H - -#include -#include -#include -#include - -#define MaxTraceLenght 8100 - -enum DataFormat{ - - ALL = 0x00, - OneTrace = 0x01, - NoTrace = 0x02, - Minimum = 0x03, - MiniWithFineTime = 0x04, - Raw = 0x0A, - -}; - -namespace DPPType{ - - const std::string PHA = "DPP_PHA"; - const std::string PSD = "DPP_PSD"; - -}; - -class Hit { - public: - - unsigned short dataType; - std::string DPPType; - - ///============= for dpp-pha - uint8_t channel; // 6 bit - uint16_t energy; // 16 bit - uint16_t energy_short; // 16 bit, only for PSD - uint64_t timestamp; // 48 bit - uint16_t fine_timestamp; // 10 bit - uint16_t flags_low_priority; // 12 bit - uint16_t flags_high_priority; // 8 bit - size_t traceLenght; // 64 bit - uint8_t downSampling; // 8 bit - bool board_fail; - bool flush; - uint8_t analog_probes_type[2]; // 3 bit for PHA, 4 bit for PSD - uint8_t digital_probes_type[4]; // 4 bit for PHA, 5 bit for PSD - int32_t * analog_probes[2]; // 18 bit - uint8_t * digital_probes[4]; // 1 bit - uint16_t trigger_threashold; // 16 bit - size_t event_size; // 64 bit - uint32_t aggCounter; // 32 bit - - ///============= for raw - uint8_t * data; - size_t dataSize; /// number of byte of the data, size/8 = word [64 bits] - uint32_t n_events; - - bool isTraceAllZero; - - Hit(){ - Init(); - } - - ~Hit(){ - ClearMemory(); - } - - void Init(){ - DPPType = DPPType::PHA; - dataType = DataFormat::ALL; - - channel = 0; - energy = 0; - energy_short = 0; - timestamp = 0; - fine_timestamp = 0; - downSampling = 0; - board_fail = false; - flush = false; - flags_low_priority = 0; - flags_high_priority = 0; - trigger_threashold = 0; - event_size = 0; - aggCounter = 0; - analog_probes[0] = NULL; - analog_probes[1] = NULL; - digital_probes[0] = NULL; - digital_probes[1] = NULL; - digital_probes[2] = NULL; - digital_probes[3] = NULL; - - analog_probes_type[0] = 0xFF; - analog_probes_type[1] = 0xFF; - digital_probes_type[0] = 0xFF; - digital_probes_type[1] = 0xFF; - digital_probes_type[2] = 0xFF; - digital_probes_type[3] = 0xFF; - data = NULL; - - isTraceAllZero = true; // indicate trace are all zero - } - - void ClearMemory(){ - if( data != NULL ) delete data; - - if( analog_probes[0] != NULL) delete analog_probes[0]; - if( analog_probes[1] != NULL) delete analog_probes[1]; - - if( digital_probes[0] != NULL) delete digital_probes[0]; - if( digital_probes[1] != NULL) delete digital_probes[1]; - if( digital_probes[2] != NULL) delete digital_probes[2]; - if( digital_probes[3] != NULL) delete digital_probes[3]; - - isTraceAllZero = true; - } - - void SetDataType(unsigned int type, std::string dppType){ - dataType = type; - DPPType = dppType; - ClearMemory(); - - if( dataType == DataFormat::Raw){ - data = new uint8_t[20*1024*1024]; - }else{ - analog_probes[0] = new int32_t[MaxTraceLenght]; - analog_probes[1] = new int32_t[MaxTraceLenght]; - - digital_probes[0] = new uint8_t[MaxTraceLenght]; - digital_probes[1] = new uint8_t[MaxTraceLenght]; - digital_probes[2] = new uint8_t[MaxTraceLenght]; - digital_probes[3] = new uint8_t[MaxTraceLenght]; - - isTraceAllZero = true; - - } - } - - void ClearTrace(){ - if( isTraceAllZero ) return; // no need to clear again - - for( int i = 0; i < MaxTraceLenght; i++){ - analog_probes[0][i] = 0; - analog_probes[1][i] = 0; - - digital_probes[0][i] = 0; - digital_probes[1][i] = 0; - digital_probes[2][i] = 0; - digital_probes[3][i] = 0; - } - isTraceAllZero = true; - } - - void PrintEnergyTimeStamp(){ - printf("ch: %2d, energy: %u, timestamp: %lu ch, traceLenght: %lu\n", channel, energy, timestamp, traceLenght); - } - - std::string AnaProbeType(uint8_t probeType){ - - if( DPPType == DPPType::PHA){ - switch(probeType){ - case 0: return "ADC"; - case 1: return "Time filter"; - case 2: return "Energy filter"; - default : return "none"; - } - }else if (DPPType == DPPType::PSD){ - switch(probeType){ - case 0: return "ADC"; - case 9: return "Baseline"; - case 10: return "CFD"; - default : return "none"; - } - }else{ - return "none"; - } - } - - std::string DigiProbeType(uint8_t probeType){ - - if( DPPType == DPPType::PHA){ - switch(probeType){ - case 0: return "Trigger"; - case 1: return "Time filter armed"; - case 2: return "Re-trigger guard"; - case 3: return "Energy filter baseline freeze"; - case 4: return "Energy filter peaking"; - case 5: return "Energy filter peaking ready"; - case 6: return "Energy filter pile-up guard"; - case 7: return "Event pile-up"; - case 8: return "ADC saturation"; - case 9: return "ADC saturation protection"; - case 10: return "Post-saturation event"; - case 11: return "Energy filter saturation"; - case 12: return "Signal inhibit"; - default : return "none"; - } - }else if (DPPType == DPPType::PSD){ - switch(probeType){ - case 0: return "Trigger"; - case 1: return "CFD Filter Armed"; - case 2: return "Re-trigger guard"; - case 3: return "ADC Input Baseline freeze"; - case 20: return "ADC Input OverThreshold"; - case 21: return "Charge Ready"; - case 22: return "Long Gate"; - case 7: return "Pile-Up Trig."; - case 24: return "Short Gate"; - case 25: return "Energy Saturation"; - case 26: return "Charge over-range"; - case 27: return "ADC Input Neg. OverThreshold"; - default : return "none"; - } - - }else{ - return "none"; - } - } - - std::string HighPriority(uint16_t prio){ - std::string output; - - bool pileup = prio & 0x1; - //bool pileupGuard = (prio >> 1) & 0x1; - //bool eventSaturated = (prio >> 2) & 0x1; - //bool postSatEvent = (prio >> 3) & 0x1; - //bool trapSatEvent = (prio >> 4) & 0x1; - //bool SCA_Event = (prio >> 5) & 0x1; - - output = std::string("Pile-up: ") + (pileup ? "Yes" : "No"); - - return output; - } - - //TODO LowPriority - - void PrintAll(){ - - switch(dataType){ - case DataFormat::ALL : printf("============= Type : ALL\n"); break; - case DataFormat::OneTrace : printf("============= Type : OneTrace\n"); break; - case DataFormat::NoTrace : printf("============= Type : NoTrace\n"); break; - case DataFormat::MiniWithFineTime : printf("============= Type : Min with FineTimestamp\n"); break; - case DataFormat::Minimum : printf("============= Type : Minimum\n"); break; - case DataFormat::Raw : printf("============= Type : Raw\n"); return; break; - default : return; - } - - printf("ch : %2d (0x%02X), fail: %d, flush: %d\n", channel, channel, board_fail, flush); - if( DPPType == DPPType::PHA ) printf("energy: %u, timestamp: %lu, fine_timestamp: %u \n", energy, timestamp, fine_timestamp); - if( DPPType == DPPType::PSD ) printf("energy: %u, energy_S : %u, timestamp: %lu, fine_timestamp: %u \n", energy, energy_short, timestamp, fine_timestamp); - printf("flag (high): 0x%02X, (low): 0x%03X, traceLength: %lu\n", flags_high_priority, flags_low_priority, traceLenght); - printf("Agg counter : %u, trigger Thr.: %u, downSampling: %u \n", aggCounter, trigger_threashold, downSampling); - printf("AnaProbe Type: %s(%u), %s(%u)\n", AnaProbeType(analog_probes_type[0]).c_str(), analog_probes_type[0], - AnaProbeType(analog_probes_type[1]).c_str(), analog_probes_type[1]); - printf("DigProbe Type: %s(%u), %s(%u), %s(%u), %s(%u)\n", DigiProbeType(digital_probes_type[0]).c_str(), digital_probes_type[0], - DigiProbeType(digital_probes_type[1]).c_str(), digital_probes_type[1], - DigiProbeType(digital_probes_type[2]).c_str(), digital_probes_type[2], - DigiProbeType(digital_probes_type[3]).c_str(), digital_probes_type[3]); - } - - void PrintTrace(unsigned short ID){ - for(unsigned short i = 0; i < (unsigned short)traceLenght; i++){ - if( ID == 0 ) printf("%4d| %6d\n", i, analog_probes[0][i]); - if( ID == 1 ) printf("%4d| %6d\n", i, analog_probes[1][i]); - if( ID == 2 ) printf("%4d| %u\n", i, digital_probes[0][i]); - if( ID == 3 ) printf("%4d| %u\n", i, digital_probes[1][i]); - if( ID == 4 ) printf("%4d| %u\n", i, digital_probes[2][i]); - if( ID == 5 ) printf("%4d| %u\n", i, digital_probes[3][i]); - } - } - - void PrintAllTrace(){ - for(unsigned short i = 0; i < (unsigned short)traceLenght; i++){ - printf("%4d| %6d %6d %1d %1d %1d %1d\n", i, analog_probes[0][i], - analog_probes[1][i], - digital_probes[0][i], - digital_probes[1][i], - digital_probes[2][i], - digital_probes[3][i]); - } - } - -}; - -#endif \ No newline at end of file diff --git a/armory/Makefile b/armory/Makefile deleted file mode 100644 index 4de3a66..0000000 --- a/armory/Makefile +++ /dev/null @@ -1,12 +0,0 @@ -CC=g++ -CFLAG= -g -ROOTFLAG=`root-config --cflags --glibs` - -all: EventBuilder - -EventBuilder: EventBuilder.cpp SolReader.h Hit.h - $(CC) $(CFLAG) EventBuilder.cpp -o EventBuilder ${ROOTFLAG} - - -clean: - -rm EventBuilder \ No newline at end of file diff --git a/armory/Monitor_Util.C b/armory/Monitor_Util.C deleted file mode 100644 index fa817a7..0000000 --- a/armory/Monitor_Util.C +++ /dev/null @@ -1,359 +0,0 @@ -#ifndef Utilities -#define Utilities - -#include -#include - -//This file runs after on Monitor.C -//This file is parasite on Monitor.C - -int canvasSize[2] = {2000, 1200}; - -void listDraws(void) { - printf("------------------- List of Plots -------------------\n"); - printf(" newCanvas() - Create a new Canvas\n"); - printf("-----------------------------------------------------\n"); - printf(" rawID() - Raw \033[0;31me\033[0m, \033[0;31mring\033[0m, \033[0;31mxf\033[0m, \033[0;31mxn\033[0m vs detID\n"); - printf(" rawe() - Raw \033[0;31me\033[0m for all %d detectors\n", numDet); - printf(" rawxf() - Raw \033[0;31mxf\033[0m for all %d detectors\n", numDet); - printf(" rawxn() - Raw \033[0;31mxn\033[0m for all %d detectors\n", numDet); - printf(" rawxfVxn() - Raw \033[0;31mxf\033[0m vs. \033[0;31mxn\033[0m for all %d detectors\n", numDet); - printf(" raweVxs() - Raw \033[0;31me\033[0m vs. Raw \033[0;31mxs = xf + xn\033[0m for all %d detectors\n", numDet); - printf(" raweVx() - Raw \033[0;31me\033[0m vs. RAW \033[0;31mx\033[0m for all %d detectors\n", numDet); - printf("-----------------------------------------------------\n"); - printf(" eVxsCal() - Raw \033[0;31me\033[0m vs. Corrected \033[0;31mxs\033[0m for all %d detectors\n", numDet); - printf(" ecal() - Calibrated \033[0;31me\033[0m for all %d detectors\n", numDet); - printf("xfCalVxnCal() - Calibrated \033[0;31mxf\033[0m vs. \033[0;31mxn\033[0m for all %d detectors\n", numDet); - printf(" eCalVxCal() - Cal \033[0;31me\033[0m vs. \033[0;31mx\033[0m for all %d detectors\n", numDet); - printf("-----------------------------------------------------\n"); - printf(" recoils() - Raw DE vs. E Recoil spectra\n"); - //printf(" elum() - Luminosity Energy Spectra\n"); - //printf(" ic() - Ionization Chamber Spectra\n"); - printf("-----------------------------------------------------\n"); - printf(" eCalVz() - Energy vs. Z\n"); - printf(" eCalVzRow() - Energy vs. Z for each row\n"); - printf(" excite() - Excitation Energy\n"); - printf(" ExThetaCM() - Ex vs ThetaCM\n"); - printf(" ExVxCal() - Ex vs X for all %d detectors\n", numDet); - printf("-----------------------------------------------------\n"); - printf(" ShowFitMethod() - Shows various fitting methods \n"); - printf(" RDTCutCreator() - Create RDT Cuts [May need to edit]\n"); - printf(" Check_rdtGate() - Check RDT Cuts. \n"); - printf(" readTrace() - read trace from gen_runXXX.root \n"); - printf(" readRawTrace() - read trace from runXXX.root \n"); - printf(" Check1D() - Count Integral within a range\n"); - printf("-----------------------------------------------------\n"); - printf(" %s\n", canvasTitle.Data()); - printf("-----------------------------------------------------\n"); -} - -int xD, yD; -void FindBesCanvasDivision(int nPad){ - for( int i = TMath::Sqrt(nPad); i >= 2 ; i--){ - if( nPad % i == 0 ) { - yD = i; - xD = nPad/i; - break; - } - } -} - -int nCanvas=0; -void newCanvas(int sizeX = 800, int sizeY = 600, int posX = 0, int posY = 0){ - TString name; name.Form("cNewCanvas%d | %s", nCanvas, canvasTitle.Data()); - TCanvas * cNewCanvas = new TCanvas(name, name, posX, posY, sizeX, sizeY); - nCanvas++; - cNewCanvas->cd(); -} - -void rawID(){ - TCanvas * cRawID = (TCanvas *) gROOT->FindObjectAny("cRawID"); - if( cRawID == NULL ) cRawID = new TCanvas("cRawID", Form("Raw e, Ring, xf, xn vs ID | %s", canvasTitle.Data()), canvasSize[0], canvasSize[1]); - cRawID->Clear();cRawID->Divide(2,2); - cRawID->cd(1); cRawID->cd(1)->SetGrid(); heVID->Draw("colz"); - cRawID->cd(2); cRawID->cd(2)->SetGrid(); hMultiHit->Draw(); - cRawID->cd(3); cRawID->cd(3)->SetGrid(); hxfVID->Draw("colz"); - cRawID->cd(4); cRawID->cd(4)->SetGrid(); hxnVID->Draw("colz"); -} - -void rawe(Bool_t isLogy = false) { - TCanvas *cRawE = (TCanvas *) gROOT->FindObjectAny("cRawE"); - if( cRawE == NULL ) cRawE = new TCanvas("cRawE",Form("E raw | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - cRawE->Clear();cRawE->Divide(numCol,numRow); - for (Int_t i=0; i < numDet; i++) { - cRawE->cd(i+1); - cRawE->cd(i+1)->SetGrid(); - if( isLogy ) cRawE->cd(i+1)->SetLogy(); - he[i]->Draw(""); - } -} - -void rawxf(Bool_t isLogy = false) { - TCanvas *cRawXf = (TCanvas *) gROOT->FindObjectAny("cRawXf"); - if( cRawXf == NULL ) cRawXf = new TCanvas("cRawXf",Form("Xf raw | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - cRawXf->Clear();cRawXf->Divide(numCol,numRow); - for (Int_t i=0; icd(i+1); - cRawXf->cd(i+1)->SetGrid(); - if( isLogy ) cRawXf->cd(i+1)->SetLogy(); - hxf[i]->Draw(""); - } -} - -void rawxn(Bool_t isLogy = false) { - TCanvas *cRawXn = (TCanvas *) gROOT->FindObjectAny("cRawXn"); - if( cRawXn == NULL ) cRawXn = new TCanvas("cRawXn",Form("Xn raw | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - cRawXn->Clear();cRawXn->Divide(numCol,numRow); - for (Int_t i=0; icd(i+1); - cRawXn->cd(i+1)->SetGrid(); - if( isLogy ) cRawXn->cd(i+1)->SetLogy(); - hxn[i]->Draw(""); - } -} - -void rawxfVxn(void) { - TCanvas *cxfxn = (TCanvas *) gROOT->FindObjectAny("cxfxn"); - if( cxfxn == NULL ) cxfxn = new TCanvas("cxfxn",Form("XF vs. XN | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - cxfxn->Clear(); cxfxn->Divide(numCol,numRow); - for (Int_t i=0;icd(i+1); - cxfxn->cd(i+1)->SetGrid(); - hxfVxn[i]->Draw("col"); - } -} - -void raweVxs(void) { - TCanvas *cxfxne = (TCanvas *) gROOT->FindObjectAny("cxfxne"); - if( cxfxne == NULL ) cxfxne = new TCanvas("cxfxne",Form("E - XF+XN | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - cxfxne->Clear(); cxfxne->Divide(numCol,numRow); - TLine line(0,0, 4000, 4000); line.SetLineColor(2); - for (Int_t i=0;icd(i+1); - cxfxne->cd(i+1)->SetGrid(); - heVxs[i]->Draw("col"); - line.Draw("same"); - } -} - -void raweVx(void) { - TCanvas *ceVx = (TCanvas *) gROOT->FindObjectAny("ceVx"); - if(ceVx == NULL) ceVx = new TCanvas("ceVx",Form("E vs. X = (xf-xn)/e | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - ceVx->Clear(); ceVx->Divide(numCol,numRow); - for (Int_t i=0;icd(i+1); heVx[i]->Draw("col"); - } -} - -void eVxsCal(void) { - TCanvas *cxfxneC = (TCanvas *) gROOT->FindObjectAny("cxfxneC"); - if(cxfxneC == NULL)cxfxneC = new TCanvas("cxfxneC",Form("Raw E - Corrected XF+XN | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - cxfxneC->Clear(); cxfxneC->Divide(numCol,numRow); - TLine line(0,0, 4000, 4000); line.SetLineColor(2); - for (Int_t i=0;icd(i+1); - cxfxneC->cd(i+1)->SetGrid(); - heVxsCal[i]->Draw("col"); - line.Draw("same"); - } -} - -void ecal(void) { - TCanvas *cEC = (TCanvas *) gROOT->FindObjectAny("cEC"); - if(cEC == NULL) cEC = new TCanvas("cEC",Form("E corrected | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - cEC->Clear();cEC->Divide(numCol,numRow); - for (Int_t i=0; icd(i+1); - cEC->cd(i+1)->SetGrid(); - heCal[i]->Draw(""); - } - - TCanvas *cEC2 = (TCanvas *) gROOT->FindObjectAny("cEC2"); - if(cEC2 == NULL) cEC2 = new TCanvas("cEC2",Form("E corrected | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - cEC2->Clear(); - heCalID->Draw("colz"); -} - -void xfCalVxnCal(void) { - TCanvas *cxfxnC = (TCanvas *) gROOT->FindObjectAny("cxfxnC"); - if(cxfxnC == NULL) cxfxnC = new TCanvas("cxfxnC",Form("XF vs XN corrected | %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - cxfxnC->Clear(); cxfxnC->Divide(numCol,numRow); - for (Int_t i=0;icd(i+1); - cxfxnC->cd(i+1)->SetGrid(); - hxfCalVxnCal[i]->Draw("col"); - } -} - -void eCalVxCal(void) { - TCanvas *cecalVxcal = (TCanvas *) gROOT->FindObjectAny("cecalVxcal"); - if( cecalVxcal == NULL ) cecalVxcal = new TCanvas("cecalVxcal",Form("ECALVXCAL | %s",canvasTitle.Data()),canvasSize[0], canvasSize[1]); - cecalVxcal->Clear(); cecalVxcal->Divide(numCol,numRow); - for (Int_t i=0;icd(i+1); - heCalVxCal[i]->SetMarkerStyle(7); - heCalVxCal[i]->Draw(""); - } -} - -void recoils(bool isLogz = false) { - TCanvas *crdt = (TCanvas *) gROOT->FindObjectAny("crdt"); - if( crdt == NULL ) crdt = new TCanvas("crdt",Form("raw RDT | %s", canvasTitle.Data()),1700, 0, 1000,1000); - crdt->Clear();crdt->Divide(2,2); - - if( isLogz ) crdt->cd(1)->SetLogz(); crdt->cd(1); hrdt2D[0]->Draw("col"); - if( isLogz ) crdt->cd(2)->SetLogz(); crdt->cd(2); hrdt2D[1]->Draw("col"); - if( isLogz ) crdt->cd(3)->SetLogz(); crdt->cd(3); hrdt2D[3]->Draw("col"); - if( isLogz ) crdt->cd(4)->SetLogz(); crdt->cd(4); hrdt2D[2]->Draw("col"); - - TCanvas *crdtID = (TCanvas *) gROOT->FindObjectAny("crdtID"); - if( crdtID == NULL ) crdtID = new TCanvas("crdtID",Form("raw RDT ID | %s", canvasTitle.Data()),0,0, 500, 500); - crdtID->Clear(); - if( isLogz ) crdtID->SetLogz(); - hrdtID->Draw("colz"); - - TCanvas *crdtS = (TCanvas *) gROOT->FindObjectAny("crdtS"); - if( crdtS == NULL ) crdtS = new TCanvas("crdtS",Form("raw RDT | %s", canvasTitle.Data()),600, 0, 1000, 1000); - crdtS->Clear(); crdtS->Divide(2,4); - for( int i = 0; i < 8; i ++){ - crdtS->cd(i+1); - if( isLogz ) crdtS->cd(i+1)->SetLogy(); - hrdt[i]->Draw(""); - } -} - -// void elum(void) { -// TCanvas *celum = (TCanvas *) gROOT->FindObjectAny("celum"); -// if( celum == NULL ) celum = new TCanvas("celum",Form("ELUM | %s", canvasTitle.Data()),1000,1000); -// celum->Clear(); celum->Divide(4,4); -// for( int i = 0 ; i < 16 ; i++){ -// celum->cd(i+1); -// helum[i]->Draw(""); -// } - -// TCanvas *celumID = (TCanvas *) gROOT->FindObjectAny("celumID"); -// if( celumID == NULL ) celumID = new TCanvas("celumID",Form("ELUM-ID | %s", canvasTitle.Data()),1100, 0, 500,500); -// celumID->Clear(); -// helumID->Draw("colz"); - -// } - - -// void ic(){ -// TCanvas *cic = (TCanvas *) gROOT->FindObjectAny("cic"); -// if( cic == NULL ) cic = new TCanvas("cic",Form("Ionization Chamber | %s", canvasTitle.Data() ),1200,800); - -// cic->Clear(); cic->SetGrid(0); cic->Divide(3,2); -// gStyle->SetOptStat("neiou"); - -// cic->cd(1); hic0->Draw(); -// cic->cd(2); hic1->Draw(); -// cic->cd(3); hic2->Draw(); -// cic->cd(4); hic01->Draw("colz"); -// cic->cd(5); hic02->Draw("colz"); -// cic->cd(6); hic12->Draw("colz"); -// } - -void eCalVz(void) { - TCanvas *cecalVz = (TCanvas *) gROOT->FindObjectAny("cecalVz"); - if( cecalVz == NULL ) cecalVz = new TCanvas("cevalVz",Form("ECALVZ : %s", canvasTitle.Data()),1000,650); - cecalVz->Clear(); cecalVz->Divide(2,1); - gStyle->SetOptStat("neiou"); - cecalVz->cd(1);heCalVz->Draw("col"); - cecalVz->cd(2);heCalVzGC->Draw("col"); -} - -void eCalVzRow() { - TCanvas *cecalVzRow = (TCanvas *) gROOT->FindObjectAny("cecalVzRow"); - if( cecalVzRow == NULL ) cecalVzRow = new TCanvas("cevalVzRow",Form("eCal - Z : %s", canvasTitle.Data()),canvasSize[0], canvasSize[1]); - FindBesCanvasDivision(numRow); - cecalVzRow->Clear(); cecalVzRow->Divide(xD,yD); - gStyle->SetOptStat("neiou"); - - for(int row = 0; row < numRow; row ++){ - cecalVzRow->cd(row+1); - cecalVzRow->cd(row+1)->SetGrid(); - hecalVzRow[row]->Draw("colz"); - } -} - -void excite(void) { - TCanvas *cex = (TCanvas *) gROOT->FindObjectAny("cex"); - if( cex == NULL ) cex = new TCanvas("cex",Form("EX : %s", canvasTitle.Data()),0, 0, 1000,650); - cex->Clear(); - gStyle->SetOptStat("neiou"); - hEx->Draw(""); - - - TCanvas *cexI = (TCanvas *) gROOT->FindObjectAny("cexI"); - if( cexI == NULL ) cexI = new TCanvas("cexI",Form("EX : %s", canvasTitle.Data()),500, 0, 1600,1000); - cexI->Clear();cexI->Divide(numCol,numRow); - gStyle->SetOptStat("neiou"); - for( int i = 0; i < numDet; i++){ - cexI->cd(i+1); - hExi[i]->Draw(""); - } -} - -void ExThetaCM(void) { - TCanvas *cExThetaCM = (TCanvas *) gROOT->FindObjectAny("cExThetaCM"); - if( cExThetaCM == NULL ) cExThetaCM = new TCanvas("cExThetaCM",Form("EX - ThetaCM | %s", canvasTitle.Data()),650,650); - cExThetaCM->Clear(); - gStyle->SetOptStat("neiou"); - hExThetaCM->Draw("colz"); -} - -void ExVxCal(TString drawOpt = "") { - TCanvas *cExVxCal = (TCanvas *) gROOT->FindObjectAny("cExVxCal"); - if( cExVxCal == NULL ) cExVxCal = new TCanvas("cExVxCal",Form("EX | %s", canvasTitle.Data()),1600,1000); - cExVxCal->Clear(); - gStyle->SetOptStat("neiou"); - - cExVxCal->Divide(numCol,numRow); - for( int i = 0; i < numDet; i++){ - cExVxCal->cd(i+1); - if( drawOpt == "" )hExVxCal[i]->SetMarkerStyle(7); - hExVxCal[i]->Draw(drawOpt); - } -} - -void Count1DH(TString name, TH1F * hist, TCanvas * canvas, int padID, double x1, double x2, Color_t color){ - - int k1 = hist->FindBin(x1); - int k2 = hist->FindBin(x2); - - int hight = 0 ; - for( int i = k1; i < k2 ; i ++){ - int temp = hist->GetBinContent(i); - if( temp > hight ) hight = temp; - } - hight = hight * 1.2; - int max = hist->GetMaximum(); - - canvas->cd(padID); - - if( color != 0 ){ - TBox box; - box.SetFillColorAlpha(color, 0.1); - box.DrawBox(x1, 0, x2, hight); - } - - int count = hist->Integral(k1, k2); - - TLatex text; - text.SetTextFont(82); - text.SetTextSize(0.06); - if( color != 0 ){ - text.SetTextColor(color); - text.DrawLatex(x1, hight, Form("%d", count)); - }else{ - text.DrawLatex((x1+x2)/2., max, Form("%d", count)); - } - - printf(" %s : %d \n", name.Data(), count); - -} - - - -#endif diff --git a/armory/Process_BasicConfig b/armory/Process_BasicConfig deleted file mode 100644 index ccc16de..0000000 --- a/armory/Process_BasicConfig +++ /dev/null @@ -1,73 +0,0 @@ -#!/bin/bash -l - -############################################## -# -# This script define color, PCID, and dataPath -# -############################################## - -if [ ! -z $RED ]; then - echo "Process_BasicConfig already loaded." - return -fi - -RED='\033[1;31m' -YELLOW='\033[1;33m' -ORANGE='\033[0;33m' -GREEN='\033[1;32m' -BLUE='\033[0;34m' -CYAN='\033[0;36m' -NC='\033[0m' #no color -LRED='\033[1;91m' - -############## need to distingish mac and daq -Arch="$(uname -s)" -PCName="$(hostname)" -PCID=-1 #if PCID == 1 (DAQ), 2 (MAC), -1(OTHER) - -#------ Set up data folder, check disk space -echo -e "${YELLOW} ##################### Check computer name and arch. ${NC}" -echo "PC name : ${PCName}" -echo "Archetech: ${Arch}" - -if [ ${Arch} == "Linux" ] && [ ${PCName} == "solaris-daq" ]; then - - PCID=1 - - pathsSetting=${HOME}/SOLARIS_QT6_DAQ/programSettings.txt - if [ -e ${pathsSetting} ]; then - #echo "Found DAQ programSettings.txt for paths settings" - - analysisPath=$(cat ${pathsSetting} | head -n 2 | tail -n 1) - - if [ ! "${analysisPath}" = "$SOLARISANADIR" ]; then - echo "The analysisPath from ${analysisPath} is different from present folder $SOLARISANADIR. Abort." - exit - fi - - rawDataPathParent=$(cat ${pathsSetting} | head -n 3 | tail -n 1) - rootDataPathParent=$(cat ${pathsSetting} | head -n 4 | tail -n 1) - - databaseIP=$(cat ${pathsSetting} | head -n 6 | tail -n 1) - databaseName=$(cat ${pathsSetting} | head -n 7 | tail -n 1) - - #echo ${rawDataPathParent} - #echo ${rootDataPathParent} - #echo ${databaseIP} - #echo ${databaseName} - - else - - echo "${RED} Cannot found DAQ programSettings.txt for path settings ${NC}" - echo "Seek Ryan for help" - exit - - fi - -fi - -if [ ${Arch} == "Darwin" ] && [ ${PCName} == "SOLARISs-Mac-Studio.local" ]; then - PCID=2 - rawDataPathParent=${HOME}/experimentalData/ - rootDataPathParent=${HOME}/experimentalData/ -fi \ No newline at end of file diff --git a/armory/Process_Download b/armory/Process_Download deleted file mode 100755 index 042be0c..0000000 --- a/armory/Process_Download +++ /dev/null @@ -1,105 +0,0 @@ -#!/bin/bash - -if [ $# -eq 0 ] || [ $1 == "-help" ]; then - echo "$./process_Download [RunNum]" - echo " RunNum = run number" - echo " * if RunNum = all, sync all" - exit 1 -fi; -RUN=$1 - -runNum=${RUN#0} #remove zero -RUN=$(printf '%03d' $runNum) ##add back the zero - -source $SOLARISANADIR/armory/Process_BasicConfig -source $SOLARISANADIR/working/expName.sh - -MacRawDataPath=$rawDataPathParent/$expName/ -IP=solarisdaq # defined at ~/.ssh/config -USR=solaris - -if [ ${RUN} == "all" ]; then - - if [ ${PCID} -eq 2 ]; then - - #=========== Ping to check the connectivity - echo "Checking $IP connetivity, max wait for 3 sec...." - ping -c 3 -W 3 $IP #> /dev/null - if [ $? -ne 0 ]; then - echo -e "$RED !!!!!!! $IP is not alive $NC" - else - #============ Get the raw data - rsync -avuht --progress $USR@$IP:$rawDataPath/$expName_* $MacRawDataPath/data/. - - echo -e "$YELLOW======== rsync RunTimeStamp.dat $NC" - rsync -avuht --progress $USR@$IP:$rawDataPath/$RunTimeStamp* $MacRawDataPath/data/. - - echo -e "$YELLOW======== rsync expName.sh $NC" - rsync -avuht --progress $USR@$IP:Analysis/working/expName.sh $SOLARISANADIR/working/. - fi - else - echo -e "$RED############### Only in SOLARIS MAC can donwload data. skip.$NC" - fi - - echo -e "$YELLOW=============================================$NC" - tail -10 $MacRawDataPath/data/RunTimeStamp.dat - echo -e "$YELLOW=============================================$NC" - - exit 1 -fi - -#just in case people put %03d as RUN -if [ "${RUN:0:2}" == "00" ]; then - RUN=${RUN:2:1} -elif [ "${RUN:0:1}" == "0" ]; then - RUN=${RUN:1:2} -else - RUN=$(printf '%d' $RUN) -fi -RUN=$(printf '%03d' ${RUN}) - -####################################### -#################### Download raw data -echo -e "${RED}######################### Download raw data: run ${RUN}${NC}" -if [ ${PCID} -eq 2 ]; then - - #=========== Ping to check the connectivity - echo "Checking $IP connetivity, max wait for 3 sec...." - ping -c 3 $IP -W 3 #> /dev/null - if [ $? -ne 0 ]; then - echo -e "$RED !!!!!!! $IP is not alive $NC" - else - #============ Get the raw data - echo -e "================= RUN $RUN: Get the raw data `date`" - - rsync -avuht --progress $USR@$IP:$rawDataPath/$expName_$RUN* $MacRawDataPath/data/. - - echo -e "$YELLOW======== rsync RunTimeStamp.dat $NC" - rsync -avuht --progress $USR@$IP:$rawDataPath/$RunTimeStamp* $MacRawDataPath/data/. - - echo -e "$YELLOW======== rsync expName.sh $NC" - rsync -avuht --progress $USR@$IP:Analysis/working/expName.sh $SOLARISANADIR/working/. - fi -else - echo -e "$RED############### Only in SOLARIS MAC can donwload data. skip.$NC" -fi - -echo -e "$YELLOW=============================================$NC" -tail -10 $MacRawDataPath/data/RunTimeStamp.dat -echo -e "$YELLOW=============================================$NC" - -count=`ls -1 $SOLARISANADIR/data_raw/${expName}_${RUN}_*.sol 2>/dev/null | wc -l` -echo -e "========== Number of Files : ${count}${NC}" -if [ ${count} -eq 0 ]; then - echo "============================================" - echo "==== RAW Files of RUN-${RUN} not found! " - echo "============================================" - - isRunDataExist=false - exit 1 -else - echo -e "${YELLOW}" - du -hc $SOLARISANADIR/data_raw/${expName}_${RUN}_*.sol - echo -e "$NC=============================================" - isRunDataExist=true -fi diff --git a/armory/Process_EventBuilder b/armory/Process_EventBuilder deleted file mode 100755 index 942808c..0000000 --- a/armory/Process_EventBuilder +++ /dev/null @@ -1,108 +0,0 @@ -#!/bin/bash - -if [ -z $SOLARISANADIR ]; then - echo "###### env variable SOLARISANADIR not defined. Abort. Please run the SOLARIS.sh." - echo "better add \"source \" into .bashrc" - exit -fi - -if [ $# -ne 3 ] || [ $1 == "-help" ]; then - echo "$ Process_EventBuilder [RunNum] [EventBuild] [timeWin]" - echo " RunNum = run number" - echo " EventBld = 2/1/0/-1/-2 || 2 = with Trace" - echo " timeWin = number of tick for an event " - echo "" - exit 1 -fi; - -RUN=$1 -EventBld=$2 -timeWin=$3 - -source ${SOLARISANADIR}/armory/Process_BasicConfig -source ${SOLARISANADIR}/working/expName.sh - -runNum=${RUN#0} #remove zero -RUN=$(printf '%03d' $runNum) ##add back the zero - -rawDataPath=$SOLARISANADIR/data_raw -rootDataPath=$SOLARISANADIR/root_data - -rawDataPattern="$rawDataPath/${expName}_${RUN}_*.sol" -rootDataName="$rootDataPath/run$RUN.root" - -dir=$(pwd) -cd ${SOLARISANADIR}/armory -make -cd ${dir} - -#==== check raw data exist -isRawDataExist=`ls -1 ${rawDataPattern}* 2>/dev/null | wc -l` - -if [ ! $isRawDataExist -gt 0 ]; then - echo -e "${LRED}################# Run Data $rawDataPattern not exist. Abort. ${NC}" - exit -fi - -echo -e "${CYAN} ============== list of files ${NC}" -\du -h ${rawDataPattern}* - -totSize=$(\du -hc ${rawDataPattern}* | tail -n 1 | awk '{print $1}') -echo -e "${CYAN} ============== total file size : ${totSize}${NC}" - - -#==== check raw data timeStamp -if [ ${Arch} == "Linux" ]; then - rawDataTime=`stat -c "%Z" ${rawDataPattern}* | sort -rn | head -1` -else - rawDataTime=`stat -f "%Sm" -t "%Y%m%d%H%M%S" $rawDataPattern | sort -rn | head -1` -fi - -#==== check if root data exist -isRootDataExist=`ls -1 $rootDataName 2>/dev/null | wc -l` - -#==== if root data exist, check timeStamp -if [ ${isRootDataExist} -gt 0 ]; then - if [ ${Arch} == "Linux" ]; then - rootDataTime=`stat -c "%Z" $rootDataName | sort -rn | head -1` - else - rootDataTime=`stat -f "%Sm" -t "%Y%m%d%H%M%S" $rootDataName | sort -rn | head -1` - fi -else - rootDataTime=0 -fi - - -if [ ${EventBld} -eq 0 ]; then - - echo -e "${LRED}>>>>>>>>>>>>>>>>>>>>> Event Building Skipped by user. ${NC}" - -elif [ ${EventBld} -ge 1 ]; then - - if [ ${rawDataTime} -ge ${rootDataTime} ]; then - - echo -e "${LRED}>>>>>>>>>>>>>>>>>>>>> Event Building $(date) ${NC}" - if [ ${EventBld} -eq 1 ]; then - EventBuilder $rootDataName ${timeWin} 0 $rawDataPattern - elif [ ${EventBld} -eq 2 ]; then - EventBuilder $rootDataName ${timeWin} 1 $rawDataPattern - fi - echo -e "${LRED}<<<<<<<<<<<<<<<< Done Event Building $(date) ${NC}" - - else - echo -e "${GREEN} root data are newer than raw data. No need to merged again.${NC}" - echo -e "${GREEN} You can Force merging using option -${EventBld}, ${ORANGE} see ./process_run.sh -help${NC}" - echo -e "${LRED}>>>>>>>>>>>>>>>>>>>>> Event Building Skipped. ${NC}" - fi - -else - - echo -e "${LRED}>>>>>>>>>>>>>>> Force Event Building $(date) ${NC}" - if [ ${EventBld} -eq -1 ]; then - EventBuilder $rootDataName ${timeWin} 0 $rawDataPattern - elif [ ${EventBld} -eq -2 ]; then - EventBuilder $rootDataName ${timeWin} 1 $rawDataPattern - fi - echo -e "${LRED}<<<<<<<<<<<<<<<< Done Event Building $(date) ${NC}" - -fi \ No newline at end of file diff --git a/armory/Process_MultiRuns b/armory/Process_MultiRuns deleted file mode 100755 index fd0257f..0000000 --- a/armory/Process_MultiRuns +++ /dev/null @@ -1,96 +0,0 @@ -#!/bin/bash -l - -if [ -z $SOLARISANADIR ]; then - echo "###### env variable SOLARISANADIR not defined. Abort. Please run the SOLARIS.sh." - echo "better add \"source \" into .bashrc" - exit -fi - -if [ $# -le 2 ] || [ $1 == "-help" ]; then - echo "$./process_MultiRuns [RunNum1] [RunNum2] [EventBuild] [GeneralSort] [numMacTerminal]" - echo " RunNum1 = start run number" - echo " RunNum2 = stop run number" - echo " EventBld = 2/1/0/-1/-2 || 2 = with Trace" - echo " GeneralSort = n/0/-n || n = number of worker" - echo " TraceMethod = -1/0/1/2 || -1 no trace, 0 save trace, 1 fit, 2 trapezoid" - echo " numMacTerminal = n ^ || in order to speed up in Mac " - echo " * negative option = force " - echo " ^ only for mac, and it will override GeneralSort to be 1 worker. " - exit 1 -fi; - - -runID1=$1 -runID2=$2 -buidEvt=$3 -nWorker=0 -traceMethod=0 -nTerminal=0 - -if [ $# -ge 4 ]; then nWorker=$4; fi -if [ $# -ge 5 ]; then traceMethod=$5; fi -if [ $# -ge 6 ]; then nTerminal=$6; fi - -source ${SOLARISANADIR}/armory/Process_BasicConfig -source ${SOLARISANADIR}/working/expName.sh - - -if [ $PCID -eq 2 ]; then - if [ $nTerminal -ge 2 ]; then - if [[ $nWorker -ge 0 ]]; then - nWorker=1; - else - nWorker=-1; - fi - fi -else - $nTerminal=0; -fi - - -if [ $nTerminal -eq 0 ]; then - - for i in $(seq $runID1 $runID2); do - Process_Run $i $buidEvt $nWorker $traceMethod 0 - done - -else - - if [ $PCID -ne 2 ]; then - exit 1 - fi - - # Calculate the size of each segment - segment_size=$(( ($runID2 - $runID1 + 1) / $nTerminal )) - - # Iterate through the segments - for i in $(seq 0 $(( $nTerminal - 1 ))); do - start=$(( $runID1 + ($i * $segment_size) )) - end=$(( $start + $segment_size - 1 )) - - # Handle the last segment, which may be larger - if [[ $i -eq $(( $nTerminal - 1 )) ]]; then - end=$runID2 - fi - - echo "Segment $(( $i + 1 )): [$start, $end]" - - profile_name="Homebrew" - width=700 - height=500 - x_pos=200 - y_pos=200 - -osascript <\" into .bashrc" - exit -fi - - -if [ $# -eq 0 ] || [ $1 == "-help" ]; then - echo "$ Process_Run [RunNum] [EventBuild] [GeneralSort] [TraceMethod] [Monitor]" - echo " RunNum = run number / \"lastRun\" " - echo " EventBld = 2/1/0/-1/-2 || 2 = with Trace" - echo " GeneralSort = n/0/-n || n = number of worker" - echo " TraceMethod = -1/0/1/2 || -1 no trace, 0 save trace, 1 fit, 2 trapezoid(not implemented)" - echo " Monitor = 2/1/0 || 1 = single run, 2 = using the list in ChainMonitors.C" - echo "" - echo " * negative option = force (except for TraceMethod and Monitor)." - echo " * Defult timeWindow for Event builder is 100 tick = 800 ns." - echo "" - exit 1 -fi; - -RUN=$1 -runNum=$1 -EventBld=2 -nWorker=1 -TraceMethod=-1 -isMonitor=1 - -if [ $# -ge 2 ]; then EventBld=$2; fi -if [ $# -ge 3 ]; then nWorker=$3; fi -if [ $# -ge 4 ]; then TraceMethod=$4; fi -if [ $# -ge 5 ]; then isMonitor=$5; fi - -if [ "$RUN" == "lastRun" ]; then - RUN=$runID -fi - -RUN=${RUN##*(0)} #remove zero -RUN=$(printf '%03d' $RUN) ##add back the zero - -################################### Setting display -echo "#################################################" -echo "### Process_Run #####" -echo "#################################################" -echo "### RunID : ${RUN}" -echo "### Event Builder : $EventBld" -echo "### General Sort : $nWorker worker(s)" -echo "### Trace Method : $TraceMethod" -echo "### Monitor : $isMonitor" -echo "#################################################" - -source ${SOLARISANADIR}/armory/Process_BasicConfig -source ${SOLARISANADIR}/working/expName.sh - -if [ "$PWD" != "${SOLARISANADIR}/working" ]; then - echo "============= go to the Working directory" - cd "${SOLARISANADIR}/working" -fi - - -#################################### CHECK IS RUN DATA EXIST -isRunDataExist=true - -#################################### EVENT BUILDER -source Process_Download $RUN - -if [ $isRunDataExist ]; then - source Process_EventBuilder $RUN $EventBld $timeWin -fi - -#################################### GeneralSort - -if [ $isRunDataExist ]; then - source Process_Sort $RUN $nWorker $TraceMethod -fi -#################################### Monitor - -if [ $isMonitor -eq 0 ]; then - echo -e "${LRED}>>>>>>>>>>>>>>>>>>>>> Monitor Skipped by user. ${NC}" -elif [ $isMonitor -eq 1 ]; then - root -l "ChainMonitors.C($RUN)" -elif [ $isMonitor -eq 2 ]; then - root -l "ChainMonitors.C" -fi - - - - diff --git a/armory/Process_Sort b/armory/Process_Sort deleted file mode 100755 index a397393..0000000 --- a/armory/Process_Sort +++ /dev/null @@ -1,91 +0,0 @@ -#!/bin/bash -l - -if [ -z $SOLARISANADIR ]; then - echo "###### env variable SOLARISANADIR not defined. Abort. Please run the SOLARIS.sh." - echo "better add \"source \" into .bashrc" - exit -fi - -if [ $# -eq 0 ] || [ $1 == "-help" ]; then - echo "$ Process_Sort [RunNum] [GeneralSort] [TraceMethod]" - echo " RunNum = run number" - echo " GeneralSort = n/0/-n || n = number of worker" - echo " TraceMethod = -1/0/1/2 || -1 no trace, 0 save trace, 1 fit, 2 trapezoid" - echo "" - exit 1 -fi; - -RUN=$1 -nWorker=$2 -TraceMethod=$3 - -source $SOLARISANADIR/armory/Process_BasicConfig -source $SOLARISANADIR/working/expName.sh - -runNum=${RUN#0} #remove zero -RUN=$(printf '%03d' $runNum) ##add back the zero - -if [ ${nWorker} -eq 0 ]; then - echo -e "${LRED}>>>>>>>>>>>>>>>>>>>>> GeneralSort Skipped by user. ${NC}" -else - source $ROOTSYS/bin/thisroot.sh - - #--------- Check is runXXX.root exist - rootDataPath=$SOLARISANADIR/root_data - rootDataName="$rootDataPath/run$RUN.root" - isRootDataExist=`ls -1 $rootDataName 2>/dev/null | wc -l` - - #==== if root data exist, check timeStamp - if [ $isRootDataExist -gt 0 ]; then - if [ ${Arch} == "Linux" ]; then - rootDataTime=`stat -c "%Z" $rootDataName | sort -rn | head -1` - else - rootDataTime=`stat -f "%Sm" -t "%Y%m%d%H%M%S" $rootDataName | sort -rn | head -1` - fi - else - rootDataTime=0 - echo -e "$LRED ################# run$RUN.root does not exist. Abort. $NC" - exit 1 - fi - - #-------- check gen_root timestamp - genRootDataName="$rootDataPath/gen_run$RUN.root" - isGenRootDataExist=`ls -1 $genRootDataName 2>/dev/null | wc -l` - - #----- if gen_runXXX.data exist, check timeStamp - if [ $isGenRootDataExist -gt 0 ]; then - if [ ${Arch} == "Linux" ]; then - genRootDataTime=`stat -c "%Z" $genRootDataName | sort -rn | head -1` - else - genRootDataTime=`stat -f "%Sm" -t "%Y%m%d%H%M%S" $genRootDataName | sort -rn | head -1` - fi - else - genRootDataTime=0 - fi - - mkdir -p ~/.proof/working - cp ${SOLARISANADIR}/working/Mapping.h ~/.proof/working/. - - if [ $nWorker -le -1 ]; then - echo -e "${LRED}>>>>>>>>>>>>>>> Force GeneralSort $(date) ${NC}" - root -l -q -b "$SOLARISANADIR/armory/GeneralSortAgent.C($runNum, $nWorker, $TraceMethod)" - echo -e "${LRED}<<<<<<<<<<<<<<<< Done GeneralSort $(date) ${NC}" - fi - - if [ $nWorker -ge 1 ]; then - - if [ $rootDataTime -ge $genRootDataTime ]; then - - echo -e "${LRED}>>>>>>>>>>>>>>>>>>>>> GeneralSort $(date) ${NC}" - root -l -q -b "$SOLARISANADIR/armory/GeneralSortAgent.C($runNum, $nWorker, $TraceMethod)" - echo -e "${LRED}<<<<<<<<<<<<<<<< Done GeneralSort $(date) ${NC}" - - else - echo -e "${GREEN} gen_run$RUN.root is newer than run$RUN.root. No need to GeneralSort again.${NC}" - echo -e "${GREEN} You can Force GeneralSort using option -${nWorker}, ${ORANGE} see Process_Run -help${NC}" - echo -e "${LRED}>>>>>>>>>>>>>>>>>>>>> GeneralSort Skipped. ${NC}" - fi - - fi - -fi \ No newline at end of file diff --git a/armory/SolReader.h b/armory/SolReader.h deleted file mode 100644 index 8761666..0000000 --- a/armory/SolReader.h +++ /dev/null @@ -1,258 +0,0 @@ -#ifndef SOLREADER_H -#define SOLREADER_H - -#include /// for FILE -#include -#include -#include -#include -#include // time in nano-sec - -#include "Hit.h" - -class SolReader { - private: - FILE * inFile; - unsigned int inFileSize; - unsigned int filePos; - unsigned int totNumBlock; - - unsigned short blockStartIdentifier; - unsigned int numBlock; - bool isScanned; - - void init(); - - std::vector blockPos; - - public: - SolReader(); - SolReader(std::string fileName, unsigned short dataType); - ~SolReader(); - - void OpenFile(std::string fileName); - int ReadNextBlock(int isSkip = 0); // opt = 0, noraml, 1, fast - int ReadBlock(unsigned int index, bool verbose = false); - - void ScanNumBlock(); - - bool IsEndOfFile() const {return (filePos >= inFileSize ? true : false);} - unsigned int GetBlockID() const {return numBlock - 1;} - unsigned int GetNumBlock() const {return numBlock;} - unsigned int GetTotalNumBlock() const {return totNumBlock;} - unsigned int GetFilePos() const {return filePos;} - unsigned int GetFileSize() const {return inFileSize;} - - void RewindFile(); - - Hit * hit; - -}; - -void SolReader::init(){ - inFileSize = 0; - numBlock = 0; - filePos = 0; - totNumBlock = 0; - hit = new Hit(); - - isScanned = false; - - blockPos.clear(); - -} - -SolReader::SolReader(){ - init(); -} - -SolReader::SolReader(std::string fileName, unsigned short dataType = 0){ - init(); - OpenFile(fileName); - hit->SetDataType(dataType, DPPType::PHA); -} - -SolReader::~SolReader(){ - if( !inFile ) fclose(inFile); - delete hit; -} - -inline void SolReader::OpenFile(std::string fileName){ - inFile = fopen(fileName.c_str(), "r"); - if( inFile == NULL ){ - printf("Cannot open file : %s \n", fileName.c_str()); - }else{ - fseek(inFile, 0L, SEEK_END); - inFileSize = ftell(inFile); - rewind(inFile); - } -} - -inline int SolReader::ReadBlock(unsigned int index, bool verbose){ - if( isScanned == false) return -1; - if( index >= totNumBlock )return -1; - fseek(inFile, 0L, SEEK_SET); - - if( verbose ) printf("Block index: %u, File Pos: %u byte\n", index, blockPos[index]); - - fseek(inFile, blockPos[index], SEEK_CUR); - - filePos = blockPos[index]; - - numBlock = index; - - return ReadNextBlock(); -} - -inline int SolReader::ReadNextBlock(int isSkip){ - if( inFile == NULL ) return -1; - if( feof(inFile) ) return -1; - if( filePos >= inFileSize) return -1; - - fread(&blockStartIdentifier, 2, 1, inFile); - - if( (blockStartIdentifier & 0xAA00) != 0xAA00 ) { - printf("header fail.\n"); - return -2 ; - } - - if( ( blockStartIdentifier & 0xF ) == DataFormat::Raw ){ - hit->SetDataType(DataFormat::Raw, ((blockStartIdentifier >> 1) & 0xF) == 0 ? DPPType::PHA : DPPType::PSD); - } - hit->dataType = blockStartIdentifier & 0xF; - hit->DPPType = ((blockStartIdentifier >> 4) & 0xF) == 0 ? DPPType::PHA : DPPType::PSD; - - if( hit->dataType == DataFormat::ALL){ - if( isSkip == 0 ){ - fread(&hit->channel, 1, 1, inFile); - fread(&hit->energy, 2, 1, inFile); - if( hit->DPPType == DPPType::PSD ) fread(&hit->energy_short, 2, 1, inFile); - fread(&hit->timestamp, 6, 1, inFile); - fread(&hit->fine_timestamp, 2, 1, inFile); - fread(&hit->flags_high_priority, 1, 1, inFile); - fread(&hit->flags_low_priority, 2, 1, inFile); - fread(&hit->downSampling, 1, 1, inFile); - fread(&hit->board_fail, 1, 1, inFile); - fread(&hit->flush, 1, 1, inFile); - fread(&hit->trigger_threashold, 2, 1, inFile); - fread(&hit->event_size, 8, 1, inFile); - fread(&hit->aggCounter, 4, 1, inFile); - }else{ - fseek(inFile, hit->DPPType == DPPType::PHA ? 31 : 33, SEEK_CUR); - } - fread(&hit->traceLenght, 8, 1, inFile); - if( isSkip == 0){ - fread(hit->analog_probes_type, 2, 1, inFile); - fread(hit->digital_probes_type, 4, 1, inFile); - fread(hit->analog_probes[0], hit->traceLenght*4, 1, inFile); - fread(hit->analog_probes[1], hit->traceLenght*4, 1, inFile); - fread(hit->digital_probes[0], hit->traceLenght, 1, inFile); - fread(hit->digital_probes[1], hit->traceLenght, 1, inFile); - fread(hit->digital_probes[2], hit->traceLenght, 1, inFile); - fread(hit->digital_probes[3], hit->traceLenght, 1, inFile); - }else{ - fseek(inFile, 6 + hit->traceLenght*(12), SEEK_CUR); - } - - }else if( hit->dataType == DataFormat::OneTrace){ - if( isSkip == 0 ){ - fread(&hit->channel, 1, 1, inFile); - fread(&hit->energy, 2, 1, inFile); - if( hit->DPPType == DPPType::PSD ) fread(&hit->energy_short, 2, 1, inFile); - fread(&hit->timestamp, 6, 1, inFile); - fread(&hit->fine_timestamp, 2, 1, inFile); - fread(&hit->flags_high_priority, 1, 1, inFile); - fread(&hit->flags_low_priority, 2, 1, inFile); - }else{ - fseek(inFile, hit->DPPType == DPPType::PHA ? 14 : 16, SEEK_CUR); - } - fread(&hit->traceLenght, 8, 1, inFile); - if( isSkip == 0){ - fread(&hit->analog_probes_type[0], 1, 1, inFile); - fread(hit->analog_probes[0], hit->traceLenght*4, 1, inFile); - }else{ - fseek(inFile, 1 + hit->traceLenght*4, SEEK_CUR); - } - - }else if( hit->dataType == DataFormat::NoTrace){ - if( isSkip == 0 ){ - fread(&hit->channel, 1, 1, inFile); - fread(&hit->energy, 2, 1, inFile); - if( hit->DPPType == DPPType::PSD ) fread(&hit->energy_short, 2, 1, inFile); - fread(&hit->timestamp, 6, 1, inFile); - fread(&hit->fine_timestamp, 2, 1, inFile); - fread(&hit->flags_high_priority, 1, 1, inFile); - fread(&hit->flags_low_priority, 2, 1, inFile); - }else{ - fseek(inFile, hit->DPPType == DPPType::PHA ? 14 : 16, SEEK_CUR); - } - - }else if( hit->dataType == DataFormat::MiniWithFineTime){ - if( isSkip == 0 ){ - fread(&hit->channel, 1, 1, inFile); - fread(&hit->energy, 2, 1, inFile); - if( hit->DPPType == DPPType::PSD ) fread(&hit->energy_short, 2, 1, inFile); - fread(&hit->timestamp, 6, 1, inFile); - fread(&hit->fine_timestamp, 2, 1, inFile); - }else{ - fseek(inFile, hit->DPPType == DPPType::PHA ? 11 : 13, SEEK_CUR); - } - - }else if( hit->dataType == DataFormat::Minimum){ - if( isSkip == 0 ){ - fread(&hit->channel, 1, 1, inFile); - fread(&hit->energy, 2, 1, inFile); - if( hit->DPPType == DPPType::PSD ) fread(&hit->energy_short, 2, 1, inFile); - fread(&hit->timestamp, 6, 1, inFile); - }else{ - fseek(inFile, hit->DPPType == DPPType::PHA ? 9 : 11, SEEK_CUR); - } - - }else if( hit->dataType == DataFormat::Raw){ - fread(&hit->dataSize, 8, 1, inFile); - if( isSkip == 0){ - fread(hit->data, hit->dataSize, 1, inFile); - }else{ - fseek(inFile, hit->dataSize, SEEK_CUR); - } - } - - numBlock ++; - filePos = ftell(inFile); - return 0; -} - -void SolReader::RewindFile(){ - rewind(inFile); - filePos = 0; - numBlock = 0; -} - -void SolReader::ScanNumBlock(){ - if( inFile == NULL ) return; - if( feof(inFile) ) return; - - numBlock = 0; - blockPos.clear(); - - blockPos.push_back(0); - - while( ReadNextBlock(1) == 0){ - blockPos.push_back(filePos); - printf("%u, %.2f%% %u/%u\n\033[A\r", numBlock, filePos*100./inFileSize, filePos, inFileSize); - } - - totNumBlock = numBlock; - numBlock = 0; - isScanned = true; - printf("\nScan complete: number of data Block : %u\n", totNumBlock); - rewind(inFile); - filePos = 0; - - //for( int i = 0; i < totNumBlock; i++){ - // printf("%7d | %u \n", i, blockPos[i]); - //} - -} - -#endif \ No newline at end of file diff --git a/armory/readRawTrace.C b/armory/readRawTrace.C deleted file mode 100644 index f2630f9..0000000 --- a/armory/readRawTrace.C +++ /dev/null @@ -1,144 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "../working/Mapping.h" - -void readRawTrace(TString fileName, int minDetID = 0, int maxDetID = 1000){ - -/**///============================================================== - - TFile * f1 = new TFile (fileName, "read"); - TTree * tree = (TTree *) f1->Get("tree"); - - if( tree == NULL ) { - printf("===== Are you using gen_runXXX.root ? please use runXXX.root\n"); - return; - } - - int totnumEntry = tree->GetEntries(); - printf( "========== total Entry : %d \n", totnumEntry); - - TCanvas * cRead = new TCanvas("cRead", "Read Raw Trace", 0, 1500, 800, 300); - cRead->Divide(1,1); - for( int i = 1; i <= 2 ; i++){ - cRead->cd(i)->SetGrid(); - } - cRead->SetGrid(); - - gStyle->SetOptFit(0); - -/**///============================================================== - Int_t bd[200]; - Int_t ch[200]; - Int_t numHit; - Int_t trace[200][2500]; - Int_t traceLength[200]; - tree->SetBranchAddress("bd", bd); - tree->SetBranchAddress("ch", ch); - tree->SetBranchAddress("multi", &numHit); - tree->SetBranchAddress("trace", trace); - tree->SetBranchAddress("tl", traceLength); - - TLatex text ; - text.SetNDC(); - text.SetTextFont(62); - text.SetTextSize(0.06); - text.SetTextColor(2); - - bool breakFlag = false; - bool lastEvFlag = false; - std::vector oldEv; - int evPointer = 0; - - TGraph * graph = new TGraph(); - - for( int ev = 0; ev < totnumEntry; ev++){ - - if( lastEvFlag ) { - ev = oldEv[evPointer]; - lastEvFlag = false; - } - tree->GetEntry(ev); - - int countJ = 0; - - for( int j = 0; j < numHit ; j++){ - - int detID = mapping[bd[j]][ch[j]]; - - if( !(minDetID <= detID && detID <= maxDetID ) ) continue; - - if( countJ == 0 ) printf("-------------------------------- ev : %d, evPointer : %d| num Trace : %d\n", ev, evPointer, numHit); - - printf("nHit: %d, id : %d, trace Length : %u ( enter = next , q = stop, w = last)\n", j, detID, traceLength[j]); - - graph->Clear(); - graph->Set(traceLength[j]); - for( int k = 0; k < traceLength[j] ; k++){ - graph->SetPoint(k, k, trace[j][k]); - //printf("%4d, %5d |", k, trace[j][k]); - //if( k % 10 ==0 ) printf("\n"); - } - - graph->SetTitle(Form("ev: %d, nHit : %d, id : %d, trace Length : %u\n", ev, j, detID, traceLength[j])); - - cRead->cd(1); - cRead->Clear(); - graph->Draw("Al"); - //g->GetXaxis()->SetRangeUser(0, g->GetN()); - //g->GetYaxis()->SetRangeUser(7500, 35000); - - cRead->Update(); - gSystem->ProcessEvents(); - - - char s[80]; - fgets(s, sizeof s, stdin); - - if( s[0] == 113 ) { // 'q' = 113 - breakFlag = true; - break; - }else if( s[0] == 119 ) { // 'w' = 119 - - if( j > 0 || countJ > 0 ) { - j = j - 2; - } - - if( (j == 0 && (int)oldEv.size() >= 0 && evPointer > 0 ) || countJ == 0 ) { - if( evPointer > 0 ) evPointer --; - if( evPointer == 0 ) printf(" the first event!!! \n"); - lastEvFlag = true; - //printf(" ev : %d, evPt : %d \n", oldEv[evPointer], evPointer); - break; - } - } - - countJ ++; - - } - if( breakFlag ) break; - - if( lastEvFlag == false && countJ > 0 ){ - if( oldEv.size() == evPointer ) oldEv.push_back(ev); - evPointer ++; - } - - - } - - //gROOT->ProcessLine(".q"); - -} diff --git a/armory/readTrace.C b/armory/readTrace.C deleted file mode 100644 index 980aeb7..0000000 --- a/armory/readTrace.C +++ /dev/null @@ -1,171 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -void readTrace(TString fileName, int minDetID = 0, int maxDetID = 1000, bool isGoodOnly = false){ - -/**///============================================================== - - TFile * f1 = new TFile (fileName, "read"); - TTree * tree = (TTree *) f1->Get("gen_tree"); - - int totnumEntry = tree->GetEntries(); - printf( "========== total Entry : %d \n", totnumEntry); - - TCanvas * cRead = new TCanvas("cRead", "Read Trace", 0, 1500, 800, 300); - cRead->Divide(1,1); - for( int i = 1; i <= 2 ; i++){ - cRead->cd(i)->SetGrid(); - } - cRead->SetGrid(); - - gStyle->SetOptFit(0); - -/**///============================================================== - TClonesArray * arr = new TClonesArray("TGraph"); - tree->SetBranchAddress("trace", &arr); - - //TODO UP-Down arrow for pervious-next control - TLine timeVLine; - - TLatex text ; - text.SetNDC(); - text.SetTextFont(62); - text.SetTextSize(0.06); - text.SetTextColor(2); - - bool breakFlag = false; - bool lastEvFlag = false; - std::vector oldEv; - int evPointer = 0; - - for( int ev = 0; ev < totnumEntry; ev++){ - - if( lastEvFlag ) { - ev = oldEv[evPointer]; - lastEvFlag = false; - } - tree->GetEntry(ev); - - int nTrace = arr->GetEntriesFast(); - - int countJ = 0; - - for( int j = 0; j < nTrace ; j++){ - - TGraph * g = (TGraph*) arr->At(j); - - TString gTitle; - gTitle = g->GetTitle(); - ///printf("TGraph Title : %s\n", gTitle.Data()); - - int detID = 0; - int pos = gTitle.Index("id:"); - gTitle.Remove(0, pos+3); - gTitle.Remove(3); - detID = gTitle.Atoi(); - - if( !(minDetID <= detID && detID <= maxDetID ) ) continue; - - if( countJ == 0 ) printf("-------------------------------- ev : %d, evPointer : %d| num Trace : %d\n", ev, evPointer, nTrace); - - - TF1 * gFit = (TF1 *) g->FindObject("gFit"); - - TString kTitle; - - if( gFit != NULL ){ - double base, time = 0, riseTime, energy, chiSq; - energy = gFit->GetParameter(0); - time = gFit->GetParameter(1); - riseTime = gFit->GetParameter(2); - base = gFit->GetParameter(3); - chiSq = gFit->GetChisquare()/gFit->GetNDF(); - int kind = gFit->GetLineColor(); - int det = gFit->GetLineStyle(); - - ///if( !(minDetID <= det && det <= maxDetID ) ) continue; - - if( isGoodOnly){ - if( abs(energy) < 1.5* g->GetRMS(2) ) continue; - if( time > gFit->GetXmax() || time < gFit->GetXmin()) continue; - if( time > 200 || time < 20) continue; - if( riseTime > gFit->GetXmax()/7. || riseTime < 0 ) continue; - } - - //if( energy < 400 ) continue; - - kTitle = Form("(%3d,%d), base: %8.1f, rise: %6.2f, time: %6.1f, energy: %8.1f | chi2: %8.1f, RMS: %6.1f", - det, kind, base, riseTime, time, energy, chiSq, g->GetRMS(2)); - - printf("%s |(q = break, w = last one)", kTitle.Data()); - - - timeVLine.SetX1(time); - timeVLine.SetX2(time); - timeVLine.SetY1(-1000); - timeVLine.SetY2(20000); - timeVLine.SetLineColor(4); - } - - cRead->cd(1); - //cRead->Clear(); - g->Draw("AC"); - //g->GetXaxis()->SetRangeUser(0, g->GetN()); - //g->GetYaxis()->SetRangeUser(7500, 35000); - timeVLine.Draw("same"); - if( gFit != NULL ) text.DrawText(0.11, 0.85, kTitle.Data()); - - cRead->Update(); - gSystem->ProcessEvents(); - - - char s[80]; - fgets(s, sizeof s, stdin); - - if( s[0] == 113 ) { // 'q' = 113 - breakFlag = true; - break; - }else if( s[0] == 119 ) { // 'w' = 119 - - if( j > 0 || countJ > 0 ) { - j = j - 2; - } - - if( (j == 0 && (int)oldEv.size() >= 0 && evPointer > 0 ) || countJ == 0 ) { - if( evPointer > 0 ) evPointer --; - if( evPointer == 0 ) printf(" the first event!!! \n"); - lastEvFlag = true; - //printf(" ev : %d, evPt : %d \n", oldEv[evPointer], evPointer); - break; - } - } - - countJ ++; - - } - if( breakFlag ) break; - - if( lastEvFlag == false && countJ > 0 ){ - if( oldEv.size() == evPointer ) oldEv.push_back(ev); - evPointer ++; - } - - - } - - //gROOT->ProcessLine(".q"); - -}