commit e3aa858ff0f0017cbac0cb5ed106f3e17c6456b9 Author: ryan@pauli Date: Tue May 24 18:21:12 2022 -0400 initial commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..8f7a6fb --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +*.pcm +*.d +*.so +*.root + +legacy_code +nscl2pixie +nscl2pixie_haha diff --git a/Cali_gamma_histogram.C b/Cali_gamma_histogram.C new file mode 100644 index 0000000..534b933 --- /dev/null +++ b/Cali_gamma_histogram.C @@ -0,0 +1,267 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "armory/AnalysisLibrary.h" + +vector Cali_gamma_histogram(TH1F * hist, int opt = 0, float threshold = 0.1, int peakDensity = 10){ + /**///======================================================== load tree + + printf("============================================================= \n"); + printf("====================== Cali_gamma.C ========================= \n"); + printf("============================================================= \n"); + + /**///======================================================== Browser or Canvas + + vector sol; + + if( hist->GetEntries() < 100 ){ + sol.push_back(0); + sol.push_back(1); + return sol; + } + + Int_t Div[2] = {2,2}; //x,y + Int_t size[2] = {600,300}; //x,y + + TCanvas * cAlpha = (TCanvas *) gROOT->FindObjectAny("cAlpha"); + if( cAlpha == NULL ){ + cAlpha = new TCanvas("cAlpha", "cAlpha", 0, 0, size[0]*Div[0], size[1]*Div[1]); + } + cAlpha->Clear(); + cAlpha->Divide(Div[0],Div[1]); + + for( int i = 1; i <= Div[0]*Div[1] ; i++){ + cAlpha->cd(i)->SetGrid(); + } + + gStyle->SetOptStat(0); + gStyle->SetStatY(1.0); + gStyle->SetStatX(0.99); + gStyle->SetStatW(0.2); + gStyle->SetStatH(0.1); + + if(cAlpha->GetShowEditor() )cAlpha->ToggleEditor(); + if(cAlpha->GetShowToolBar() )cAlpha->ToggleToolBar(); + + /**///========================================================= Analysis + + cAlpha->cd(1); + hist->Draw(); + cAlpha->Update(); + gSystem->ProcessEvents(); + + int dummy = 0; + int temp = 0; + + int nPeaks = 10; + vector energy; + vector refEnergy; + vector fitEnergy; + + + printf("---- finding peak using TSepctrum Class...\n"); + + TSpectrum * spec = new TSpectrum(20); + nPeaks = spec->Search(hist, peakDensity, "", threshold); + printf("---- found %2d peaks | ", nPeaks); + + double * xpos = spec->GetPositionX(); + double * ypos = spec->GetPositionY(); + + cAlpha->Update(); + gSystem->ProcessEvents(); + + vector height; + + int * inX = new int[nPeaks]; + TMath::Sort(nPeaks, xpos, inX, 0 ); + 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("%7.2f, ", energy[j]); + printf("\n"); + + + //------------ 3, correction + int option = 0; + if ( opt == 0 ) { + printf("========== which detector to be the reference?\n"); + printf("-1 = manual reference\n"); + printf("-2 = use 228Th, first 5 strongest peaks \n"); + printf("-3 = use 207Bi, 3 peaks \n"); + printf("-4 = use 152Eu, 11 peaks \n"); + printf("-9 = stop \n"); + printf("your choice = "); + temp = scanf("%d", &option); + }else{ + option = opt; + } + + if( option == -9 ) { + printf("------ stopped by user.\n"); + return sol; + } + + + + //======== fill reference energy + if(option == -1){ + for( int i = 0; i < nPeaks; i++){ + float eng; + printf("%2d-th peak energy %f ch (-1 to skip, -2 to skip all):", i, energy[i]); + temp = scanf("%f", &eng); + if( eng >= 0 ) { + refEnergy.push_back(eng); + fitEnergy.push_back(energy[i]); + printf(" input: %f \n", eng); + }else if ( eng == -1 ){ + printf(" input: skipped \n"); + }else if ( eng == -2 ){ + break; + } + }; + } + + if( option == -2 ){ + fitEnergy = energy; + + refEnergy.clear(); + refEnergy.push_back(5.423); + refEnergy.push_back(5.685); + refEnergy.push_back(6.288); + refEnergy.push_back(6.778); + refEnergy.push_back(8.785); + } + + if( option == -3 ){ + + fitEnergy = energy; + + refEnergy.clear(); + refEnergy.push_back(0.569702); + refEnergy.push_back(1.063662); + refEnergy.push_back(1.770237); + } + + if( option == -4 ){ + fitEnergy = energy; + + refEnergy.clear(); + refEnergy.push_back( 121.783); + refEnergy.push_back( 244.699); + refEnergy.push_back( 344.281); + refEnergy.push_back( 411.115); + refEnergy.push_back( 443.965); + refEnergy.push_back( 778.903); + refEnergy.push_back( 867.390); + refEnergy.push_back( 964.055); + refEnergy.push_back(1085.842); + refEnergy.push_back(1112.087); + refEnergy.push_back(1408.022); + } + + //==================== adjusting energy + int n = refEnergy.size(); + for( int k = 0; k < n; k++) printf("%2d-th peak : %f \n", k, refEnergy[k]); + + nPeaks = energy.size(); + + vector k1 = fitEnergy; + vector k2 = refEnergy; + //===== when nPeaks != refEnergy.size(), need to matching the two vector size by checking the r-squared. + if( option != -1 ){ + vector> haha = FindMatchingPair(fitEnergy, refEnergy); + k1 = haha[0]; + k2 = haha[1]; + } + + TGraph * graph = new TGraph(min(n, nPeaks), &k1[0], &k2[0] ); + cAlpha->cd(2); + graph->Draw("A*"); + gSystem->ProcessEvents(); + + TF1 * fit = new TF1("fit", "pol1" ); + graph->Fit("fit", "q"); + + double a0 = fit->GetParameter(0); + double a1 = fit->GetParameter(1); + double a2 = 0.0 ; //fit->GetParameter(2); + + printf("----- a0: %9.6f, a1: %9.6f (%14.8f) \n", a0, a1, 1./a1); + printf(" %13.10f\t %13.10f \n", a0, a1); + + sol.clear(); + sol.push_back(a0); + sol.push_back(a1); + + //====== Plot residue + vector resi; + for( int i = 0 ; i < min(n, nPeaks); i++){ + resi.push_back(k2[i] - fit->Eval(k1[i])); + } + TGraph * graphResi = new TGraph(min(n, nPeaks), &k2[0], &resi[0] ); + + cAlpha->cd(4); + graphResi->Draw("A*"); + graphResi->GetYaxis()->SetTitle("residue [keV]"); + graphResi->GetXaxis()->SetTitle("cali energy [keV]"); + gSystem->ProcessEvents(); + + //====== Plot adjusted spectrum + + cAlpha->cd(3); + + int bin = hist->GetNbinsX(); + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + + double xMinC, xMaxC; + if( a1 > 0 ) { + xMinC = xMin*xMin*a2 + xMin*a1 + a0; + xMaxC = xMax*xMax*a2 + xMax*a1 + a0; + }else{ + xMaxC = xMin*xMin*a2 + xMin*a1 + a0; + xMinC = xMax*xMax*a2 + xMax*a1 + a0; + } + + + TH1F * calH = new TH1F("calH", "calibrated energy", bin, xMinC, xMaxC); + + /* + FILE * paraOut; + TString filename; + filename.Form("hist%s.dat", hist->GetName()); + paraOut = fopen (filename.Data(), "w+"); + */ + for( int i = 1; i <= bin; i ++){ + int y = hist->GetBinContent(i); + int x = hist->GetBinCenter(i); + + calH->Fill(x*x*a2 + x*a1+a0, y); + + //fprintf(paraOut, "%9.6f, %9d\n", x*x*a2 + x*a1+a0, y); + + } + + // fflush(paraOut); + // fclose(paraOut); + + calH->Draw(); + gSystem->ProcessEvents(); + + return sol; + +} diff --git a/DecayFinder.C b/DecayFinder.C new file mode 100644 index 0000000..139ba6e --- /dev/null +++ b/DecayFinder.C @@ -0,0 +1,246 @@ +#define DecayFinder_cxx + +#include "DecayFinder.h" +#include +#include +#include +#include +#include +#include +#include +#include + +//===================== Settings +Bool_t debug = false; + +TString pidCutFileName = "PIDCuts.root"; + +double distThreshold = 0.2; +double forwardTime = 250; //ms +double backwardTime = 100; //ms + +//===================== histograms + +TH2F * hPID; +TH1F ** hDecay ; +TH2F * hDist; +TH1F ** hvetoF; + +TH2F * hIonsPos; +TH2F * hBetaPos; + +//===================== Parameters + +int tick2ns = 8; +double tick2ms = tick2ns / 1e6; + +TFile * cutFile; +TCutG * cut = NULL; +TObjArray * cutList; +int numCut = 0; + +int numMatch = 0; + +//===================== Begins +void DecayFinder::Begin(TTree * /*tree*/){ + + TString option = GetOption(); + + hPID = new TH2F("hPID", "PID; TOF; dE", 400, -260, -220, 400, 0, 6500); + + hDist = new TH2F("hDist", "dist; dx; dy", 100, -distThreshold, distThreshold, 100, -distThreshold, distThreshold); + + hIonsPos = new TH2F("hIonsPos", "Ions Pos; x; y", 200, 0, 1, 200, 0, 1); + hBetaPos = new TH2F("hBetaPos", "Beta Pos; x; y", 200, 0, 1, 200, 0, 1); + + //-------------- GetCut; + cutFile = new TFile(pidCutFileName, "UPDATE"); + bool listExist = cutFile->GetListOfKeys()->Contains("cutList"); + if( !listExist ) { + cutList = new TObjArray(); + }else{ + cutList = (TObjArray*) cutFile->FindObjectAny("cutList"); + numCut = cutList->GetLast()+1; + printf("------------ found %d cuts in %s \n", numCut, pidCutFileName.Data()); + for( int k = 0; k < numCut; k++){ + if( cutList->At(k) != NULL ){ + printf("found a cut at %2d \n", k); + }else{ + printf(" No cut at %2d \n", k); + } + } + } + + hDecay = new TH1F *[numCut]; + hvetoF = new TH1F *[numCut]; + + for( int i = 0; i < numCut; i++){ + hDecay[i] = new TH1F(Form("hDecay%02d",i), Form("Decay cut-%02d ; [ms]; count", i), 50, -backwardTime, forwardTime); + hvetoF[i] = new TH1F(Form("hvetoF%02d",i), Form("veto-F cut-%02d ; [ms]; count", i), 100, 0, 6000); + } + + printf("=====================================\n"); + +} + +double dist(double x1, double x2, double y1, double y2){ + return sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)); +} + +//===================== Process +Bool_t DecayFinder::Process(Long64_t entry){ + + if( entry > (Long64_t)1e6 ) return kTRUE; + + b_TOF->GetEntry(entry); + b_energy->GetEntry(entry); + b_crossTime->GetEntry(entry); + b_crossEnergy->GetEntry(entry); + b_flag->GetEntry(entry); + b_vetoFlag->GetEntry(entry); + b_xIons->GetEntry(entry); + b_yIons->GetEntry(entry); + b_dyIonsTime->GetEntry(entry); + b_veto_r->GetEntry(entry); + b_veto_f->GetEntry(entry); + + ULong64_t ImplantEntry = 0; + ULong64_t ImplantTime = 0; + double ImplantX = TMath::QuietNaN(); + double ImplantY = TMath::QuietNaN(); + + if( debug && entry % 10 == 0 ) printf("------------- %llu\n", entry); + + if( flag & 1 ) hPID->Fill(TOF, energy); + + int cutID = -1; + for(int i = 0; i < numCut; i++){ + cut = (TCutG*) cutList->At(i); + if( cut->IsInside(TOF, energy) ) cutID = i; + } + + if( cutID == -1 ) return kTRUE; + if( (flag & 1) == 0 ) return kTRUE; /// no beam + //if( (flag & 2) == 0 ) return kTRUE; /// no Ions + if( (flag & 4) == 0 ) return kTRUE; /// no beta + + //if( veto_f > 0 ) return kTRUE; + if( veto_r > 0 ) return kTRUE; + if( dyIonsTime[0] == 0 ) return kTRUE; + + ImplantEntry = entry; + ImplantTime = dyIonsTime[0]; + ImplantX = xIons; + ImplantY = yIons; + + hIonsPos->Fill( xIons, yIons); + + if( debug ) printf("===========%8lld, %13llu, %f, %f\n", ImplantEntry, ImplantTime, ImplantX, ImplantY); + + for( ULong64_t k = ImplantEntry - 10000; k < totNumEntry ; k++){ + b_flag->GetEntry(k); + b_vetoFlag->GetEntry(k); + b_dyBetaTime->GetEntry(k); + b_crossTime->GetEntry(k); + b_crossEnergy->GetEntry(k); + b_veto_f->GetEntry(k); + b_veto_r->GetEntry(k); + b_xBeta->GetEntry(k); + b_yBeta->GetEntry(k); + + if ( k == ImplantEntry ) continue; + + //if( crossTime > 0 ) continue; + if( crossEnergy > 0 ) continue; + //if( (flag & 2) == 2 ) continue; /// has Ions + if( (flag & 4) == 0 ) continue; /// no Beta + if( !TMath::IsNaN(veto_f) ) continue; + if( !TMath::IsNaN(veto_r) ) continue; + if( dyBetaTime[0] == 0) continue; + if( dist(ImplantX, xBeta, ImplantY, yBeta) > distThreshold ) continue; + + if( k < ImplantEntry && ( ImplantTime - dyBetaTime[0] ) * tick2ms > backwardTime) continue; + if( k > ImplantEntry && ( dyBetaTime[0] - ImplantTime ) * tick2ms > forwardTime ) break; + + clock.Stop("timer"); + Double_t time = clock.GetRealTime("timer"); + clock.Start("timer"); + printf(" %llu[%.2f%%], searched next %lld events | match case %d | expect : %5.2f min\r", + ImplantEntry, ImplantEntry*100./totNumEntry, k - ImplantEntry, + numMatch, totNumEntry*time/(ImplantEntry+1.)/60.); + + if( debug ) printf(" %8lld, %13llu, %f, %f\n", k, dyBetaTime[0], xBeta, yBeta); + if( debug ) printf("+++++++++++%8s, %13llu | %.f ms\n", "+", dyBetaTime[0], (dyBetaTime[0] - ImplantTime)*tick2ms); + + hDist->Fill(ImplantX - xBeta, ImplantY - yBeta); + + numMatch ++; + double dT = 0; + if( dyBetaTime[0] > ImplantTime ){ + dT = (dyBetaTime[0] - ImplantTime)*tick2ms; + }else{ + dT = (ImplantTime - dyBetaTime[0])*tick2ms; + dT = -dT; + } + + hDecay[cutID]->Fill(dT); + hvetoF[cutID]->Fill(crossEnergy); + hBetaPos->Fill(xBeta, yBeta); + + } + + return kTRUE; +} + + +void DecayFinder::Terminate(){ + + //============== Save the results + TFile * haha = new TFile("results.root", "recreate"); + haha->cd(); + for( int i = 0; i < numCut; i++){ + hDecay[i]->Write(); + } + hPID->Write(); + haha->Close(); + + //============== Plot result + gStyle->SetOptStat("neiou"); + TCanvas * cDecay = new TCanvas("cDecay", "Decay", 1000, 1000); + cDecay->Divide(2,2); + + int padID = 1; cDecay->cd(padID); + hPID->Draw("colz"); + for( int i = 0; i < numCut; i++){ + cut = (TCutG*) cutList->At(i); + cut->SetLineColor(i+2); + cut->Draw("same"); + } + + padID++; cDecay->cd(padID); //cDecay->cd(padID)->SetLogy(); + double yMax = 0; + for( int i = 0; i < numCut; i++){ + if( hDecay[i]->GetMaximum() > yMax ) yMax = hDecay[i]->GetMaximum(); + } + for( int i = 0; i < numCut; i++){ + hDecay[i]->SetLineColor(i+2); + hDecay[i]->SetMaximum(yMax *1.2); + hDecay[i]->Draw(i == 0 ? "" : "same"); + } + + //padID++; cDecay->cd(padID); + //hDist->Draw("colz"); + // + //padID++; cDecay->cd(padID); + //for( int i = 0; i < numCut; i++){ + // hvetoF[i]->SetLineColor(i+2); + // hvetoF[i]->Draw(i == 0 ? "" : "same"); + //} + + padID++; cDecay->cd(padID); + hIonsPos->Draw("colz"); + + padID++; cDecay->cd(padID); + hBetaPos->Draw("colz"); + +} diff --git a/DecayFinder.h b/DecayFinder.h new file mode 100644 index 0000000..93aee38 --- /dev/null +++ b/DecayFinder.h @@ -0,0 +1,136 @@ +////////////////////////////////////////////////////////// +// This class has been automatically generated on +// Mon May 16 17:09:53 2022 by ROOT version 6.24/06 +// from TTree tree/tree +// found on file: run-0238_00_02.root +////////////////////////////////////////////////////////// + +#ifndef DecayFinder_h +#define DecayFinder_h + +#include +#include +#include +#include +#include + +// Header file for the classes stored in the TTree if any. + +class DecayFinder : 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 eventID; + //Int_t runID; + Double_t TOF; + Double_t energy; + UInt_t crossEnergy; + ULong64_t crossTime; + Short_t flag; + Short_t vetoFlag; + Double_t xIons; + Double_t yIons; + Double_t xBeta; + Double_t yBeta; + ULong64_t dyIonsTime[4]; + ULong64_t dyBetaTime[4]; + Double_t veto_f; + Double_t veto_r; + + // List of branches + TBranch *b_eventID; //! + //TBranch *b_runID; //! + TBranch *b_TOF; //! + TBranch *b_energy; //! + TBranch *b_crossTime; //! + TBranch *b_crossEnergy; //! + TBranch *b_flag; //! + TBranch *b_vetoFlag; //! + TBranch *b_xIons; //! + TBranch *b_yIons; //! + TBranch *b_xBeta; //! + TBranch *b_yBeta; //! + TBranch *b_dyIonsTime; //! + TBranch *b_dyBetaTime; //! + TBranch *b_veto_f; //! + TBranch *b_veto_r; //! + + DecayFinder(TTree * /*tree*/ =0) : fChain(0) { } + virtual ~DecayFinder() { } + 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(DecayFinder,0); + + ULong64_t totNumEntry; + + //clock + TBenchmark clock; + +}; + +#endif + +#ifdef DecayFinder_cxx +void DecayFinder::Init(TTree *tree){ + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("eventID", &eventID, &b_eventID); + //fChain->SetBranchAddress("runID", &runID, &b_runID); + fChain->SetBranchAddress("TOF", &TOF, &b_TOF); + fChain->SetBranchAddress("energy", &energy, &b_energy); + fChain->SetBranchAddress("crossTime", &crossTime, &b_crossTime); + fChain->SetBranchAddress("crossEnergy", &crossEnergy, &b_crossEnergy); + fChain->SetBranchAddress("flag", &flag, &b_flag); + fChain->SetBranchAddress("vetoFlag", &vetoFlag, &b_vetoFlag); + fChain->SetBranchAddress("xIons", &xIons, &b_xIons); + fChain->SetBranchAddress("yIons", &yIons, &b_yIons); + fChain->SetBranchAddress("xBeta", &xBeta, &b_xBeta); + fChain->SetBranchAddress("yBeta", &yBeta, &b_yBeta); + fChain->SetBranchAddress("dyIonsTime", dyIonsTime, &b_dyIonsTime); + fChain->SetBranchAddress("dyBetaTime", dyBetaTime, &b_dyBetaTime); + fChain->SetBranchAddress("veto_f", &veto_f, &b_veto_f); + fChain->SetBranchAddress("veto_r", &veto_r, &b_veto_r); + + totNumEntry = fChain->GetEntries(); + + //==== clock + clock.Reset(); + clock.Start("timer"); + +} + +Bool_t DecayFinder::Notify(){ + + return kTRUE; +} + +void DecayFinder::SlaveTerminate(){ + +} + + +void DecayFinder::SlaveBegin(TTree * /*tree*/){ + TString option = GetOption(); +} + + +#endif // #ifdef DecayFinder_cxx diff --git a/PIDCutChecker.C b/PIDCutChecker.C new file mode 100644 index 0000000..137754b --- /dev/null +++ b/PIDCutChecker.C @@ -0,0 +1,54 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +void PIDCutChecker(TH2F* hist, + TString saveFileName = "PIDCuts.root", + bool isLogz = false + ){ + + printf("================ Graphic Cut Checker for PID ============== \n"); + + gStyle->SetOptStat("neiou"); + + TCanvas * cCutCreator = new TCanvas("cCutCreator", "PID Cut Creator", 100, 100, 800, 800); + if( !cCutCreator->GetShowToolBar() ) cCutCreator->ToggleToolBar(); + + cCutCreator->Update(); + if( isLogz ) cCutCreator->cd()->SetLogz(); + hist->Draw("colz"); + + TFile * cutFile; + TCutG * cut = NULL; + TObjArray * cutList; + + ///===================== load the cutFile and load the cutList + + int numCut = 0; + cutFile = new TFile(saveFileName, "UPDATE"); + bool listExist = cutFile->GetListOfKeys()->Contains("cutList"); + if( !listExist ) { + cutList = new TObjArray(); + }else{ + cutList = (TObjArray*) cutFile->FindObjectAny("cutList"); + numCut = cutList->GetLast()+1; + printf("----- found %d cuts \n", numCut); + for( int k = 0; k < numCut; k++){ + if( cutList->At(k) != NULL ){ + printf("found a cut at %2d \n", k); + cut = (TCutG*) cutList->At(k); + cut->Draw("same"); + }else{ + printf(" No cut at %2d \n", k); + } + } + } + cCutCreator->Modified(); cCutCreator->Update(); + +} diff --git a/PIDCutCreator.C b/PIDCutCreator.C new file mode 100644 index 0000000..bc29cdf --- /dev/null +++ b/PIDCutCreator.C @@ -0,0 +1,90 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +void PIDCutCreator(TH2F* hist, + TString saveFileName = "PIDCuts.root", + bool isLogz = false + ){ + + printf("================ Graphic Cut Creator for PID ============== \n"); + + gStyle->SetOptStat("neiou"); + + TCanvas * cCutCreator = new TCanvas("cCutCreator", "PID Cut Creator", 100, 100, 800, 800); + if( !cCutCreator->GetShowToolBar() ) cCutCreator->ToggleToolBar(); + + cCutCreator->Update(); + if( isLogz ) cCutCreator->cd()->SetLogz(); + hist->Draw("colz"); + + TFile * cutFile; + TCutG * cut = NULL; + TObjArray * cutList; + + ///===================== load the cutFile and load the cutList + + int numCut = 0; + cutFile = new TFile(saveFileName, "UPDATE"); + if( cutFile == NULL ){ + cutFile = new TFile(saveFileName, "recreate"); + } + bool listExist = cutFile->GetListOfKeys()->Contains("cutList"); + if( !listExist ) { + cutList = new TObjArray(); + }else{ + cutList = (TObjArray*) cutFile->FindObjectAny("cutList"); + numCut = cutList->GetLast()+1; + printf("----- found %d cuts \n", numCut); + for( int k = 0; k < numCut; k++){ + if( cutList->At(k) != NULL ){ + printf("found a cut at %2d \n", k); + cut = (TCutG*) cutList->At(k); + cut->Draw("same"); + }else{ + printf(" No cut at %2d \n", k); + } + } + } + cCutCreator->Modified(); cCutCreator->Update(); + + int countCut = 0; + while (true){ + + printf("======== make a graphic cut on the plot (double click to stop), %2d-th cut: ", countCut ); + + gPad->WaitPrimitive(); + + cut = (TCutG*) gROOT->FindObject("CUTG"); + + if( cut == NULL ){ + printf(" stopped by user.\n"); + break; + } + + TString name; name.Form("cut%dd", countCut); + cut->SetName(name); + cut->SetVarX("TOF"); + cut->SetVarY("dE"); + cut->SetTitle(Form("cut%2d", countCut)); + cut->SetLineColor(countCut+1); + printf(" cut-%d \n", countCut); + + cutList->AddAtAndExpand(cut, countCut); + countCut ++; + + } + + ///================= Summary + + cutList->Write("cutList", TObject::kSingleKey); + printf("====> saved %d cuts into %s\n", countCut, saveFileName.Data()); + gROOT->ProcessLine(".q"); + +} diff --git a/armory/AnalysisLibrary.h b/armory/AnalysisLibrary.h new file mode 100644 index 0000000..b4f7eb8 --- /dev/null +++ b/armory/AnalysisLibrary.h @@ -0,0 +1,296 @@ +#ifndef Analysis_Library_h +#define Analysis_Library_h + +#include +#include +#include +#include +#include +#include +#include +#include + +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; +} + +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; + +} + + +std::vector> LoadCorrectionParameters(TString corrFile, bool show=false){ + + printf(" load correction parameters : %s", corrFile.Data()); + std::ifstream file; + file.open(corrFile.Data()); + + std::vector> corr; + corr.clear(); + + std::vector detCorr; + detCorr.clear(); + + if( file.is_open() ){ + + while( file.good() ){ + + std::string line; + getline(file, line); + + if( line.substr(0,1) == "#" ) continue; + if( line.substr(0,2) == "//" ) continue; + if( line.size() == 0 ) continue; + + //printf("%s\n", line.c_str()); + std::vector temp = SplitStr(line, " "); + + detCorr.clear(); + for( int i = 0; i < (int) temp.size() ; i++){ + detCorr.push_back(std::stod(temp[i])); + } + corr.push_back(detCorr); + } + + file.close(); + + printf(".... done\n"); + if( show ){ + printf("===== correction parameters \n"); + for( int i = 0; i < (int) corr.size(); i++){ + printf("det : %2d | ", i ); + int len = (int) corr[i].size(); + for( int j = 0; j < len - 1 ; j++){ + printf("%14.6f, ", corr[i][j]); + } + printf("%14.6f\n", corr[i][len-1]); + } + } + + }else{ + printf(".... fail\n"); + std::vector temp = {0, 1}; + for( int i = 0; i < 36; i++){ + corr.push_back(temp); + } + } + + return corr; +} + +#endif + diff --git a/armory/Analyzer_Utili.c b/armory/Analyzer_Utili.c new file mode 100644 index 0000000..4549f22 --- /dev/null +++ b/armory/Analyzer_Utili.c @@ -0,0 +1,195 @@ +void listDraws(void) { + printf("------------------- List of Plots -------------------\n"); + printf(" newCanvas() - Create a new Canvas\n"); + printf("-----------------------------------------------------\n"); + printf(" eVID() - e vs ID\n"); + printf(" drawE() - e for all %d detectors\n", NCRYSTAL); + //printf(" drawGG() - Gamma - Gamma Coincident for all %d detectors\n", NCRYSTAL); + printf("-----------------------------------------------------\n"); + printf(" energyCalibration() - Calibrate energy \n"); + printf("-----------------------------------------------------\n"); +} + +int nCanvas=0; +void newCanvas(int sizeX = 800, int sizeY = 600, int posX = 0, int posY = 0){ + TString name; name.Form("cNewCanvas%d", nCanvas); + TCanvas * cNewCanvas = new TCanvas(name, name, posX, posY, sizeX, sizeY); + nCanvas++; + cNewCanvas->cd(); +} + +void eVID(bool cal = false, bool logz = false){ + TCanvas * cRawID = (TCanvas *) gROOT->FindObjectAny("cRawID"); + if( cRawID == NULL ) cRawID = new TCanvas("cRawID", "raw ID", 1000, 800); + + if( cRawID->GetShowEventStatus() == 0 ) cRawID->ToggleEventStatus(); + + cRawID->cd(1)->SetGrid(); + if( logz ) cRawID->cd(1)->SetLogz(); + cal ? heCalVID->Draw("colz") : heVID->Draw("colz"); +} + +void drawE(int CloverID = -1, bool cali = false, bool isLogy = false, double xMin = 0, double xMax = 0){ + + int nCrystalPerClover = 4; + int nClover = NCRYSTAL / nCrystalPerClover; + + if( CloverID >= nClover ) { + printf("Clover-ID > nClover = %d. \n", nClover); + return; + } + + int size = 300; + + TCanvas *cRawE = (TCanvas *) gROOT->FindObjectAny("cRawE"); + if( cRawE == NULL ) cRawE = new TCanvas("cRawE", cali ? "Cal e" : "Raw e", size * nClover, size * nCrystalPerClover); + + if( cRawE->GetShowEventStatus() == 0 ) cRawE->ToggleEventStatus(); + + cRawE->Clear(); + if( CloverID >= 0 ) { + nClover = 1; + cRawE->Divide(nClover, 1); + }else{ + cRawE->Divide(nClover, nCrystalPerClover, 0); + } + + ///find max y + double maxY = 0; + int nDet = nClover*nCrystalPerClover; + for( int i = (CloverID < 0 ? 0 : nCrystalPerClover*CloverID) ; i < (CloverID < 0 ? nDet : nCrystalPerClover*CloverID + nDet) ; i++){ + int mBin = cali ? heCal[i]->GetMaximumBin() : he[i]->GetMaximumBin(); + double max = cali ? heCal[i]->GetBinContent(mBin) : he[i]->GetBinContent(mBin); + if( max > maxY ) maxY = max; + } + maxY = maxY * 1.1; + ///printf("max Y : %f \n", maxY); + + for (Int_t i = 0; i < nClover; i++) { + + int hID = nCrystalPerClover * CloverID + i ; + if( cali ) { + heCal[hID]->SetMaximum(maxY); + heCal[hID]->Draw("same"); + }else{ + he[hID]->SetMaximum(maxY); + he[hID]->Draw("same"); + } + + + + ///for( Int_t j = 0; j < nCrystalPerClover; j++){ + /// int canvasID = CloverID < 0 ? nClover*j+ i + 1 : j + 1; + /// cRawE->cd(canvasID); + /// cRawE->cd(canvasID)->SetGrid(); + /// cRawE->cd(canvasID)->SetTickx(2); + /// cRawE->cd(canvasID)->SetTicky(2); + /// cRawE->cd(canvasID)->SetBottomMargin(0.06); + /// if ( i == nClover -1 ) cRawE->cd(canvasID)->SetRightMargin(0.002); + /// if( isLogy ) cRawE->cd(canvasID)->SetLogy(); + /// int hID = CloverID < 0 ? nCrystalPerClover*i+ j : nCrystalPerClover * CloverID + j ; + /// if( cali ) { + /// if ( xMin != 0 || xMax != 0 ) heCal[hID]->GetXaxis()->SetRangeUser(xMin, xMax); + /// heCal[hID]->SetMaximum(maxY); + /// heCal[hID]->Draw(""); + /// }else{ + /// if ( xMin != 0 || xMax != 0 ) he[hID]->GetXaxis()->SetRangeUser(xMin, xMax); + /// he[hID]->SetMaximum(maxY); + /// he[hID]->Draw(""); + /// } + ///} + } + + cRawE->SetCrosshair(1); + +} + +/** +void drawGG(){ + + int nCrystal = 4; + int numCol = NCRYSTAL / nCrystal; + + int size = 300; + + TCanvas *cGG = (TCanvas *) gROOT->FindObjectAny("cGG"); + if( cGG == NULL ) cGG = new TCanvas("cGG", "Gamma - Gamma Coin.", size * NCRYSTAL, size * NCRYSTAL); + cGG->Clear();cGG->Divide(NCRYSTAL, NCRYSTAL); + + for( int i = 0; i < NCRYSTAL; i ++){ + for( int j = i+1; j < NCRYSTAL; j ++){ + cGG->cd( NCRYSTAL * i + j +1 ); + hgg[i][j]->Draw("colz"); + } + } + +} +*/ + +void energyCalibration(int detID = -1, int BG = 10, double threshold = 0.1, double sigmaMax = 5, int peakDensity = 10){ + + TCanvas *cCal = (TCanvas *) gROOT->FindObjectAny("cCal"); + if( cCal == NULL ) cCal = new TCanvas("cCal", "Energy Calibration", 1000, 0, 1000, 600); + cCal->Clear(); + + cCal->Divide(2,1); + cCal->SetGrid(); + + vector refEnergy = {121.738, + 244.699, + 344.281, + 411.115, + 443.965, + 778.903, + 867.390, + 964.055, + 1085.842, + ///1089.700, + 1112.087, + 1408.022}; + + double a0[NCRYSTAL]; + double a1[NCRYSTAL]; + + for( int i = 0 ; i < NCRYSTAL; i++){ + if( detID >= 0 && i != detID ) continue; + + cCal->cd(1); + he[i]->Draw(); + vector peaks = fitAuto(he[i], BG, threshold, sigmaMax, peakDensity); + vector> output = FindMatchingPair(peaks, refEnergy); + + vector haha1 = output[0]; + vector haha2 = output[1]; + + TGraph * graph = new TGraph(haha1.size(), &haha1[0], &haha2[0] ); + cCal->cd(2); + graph->Draw("A*"); + + TF1 * fit = new TF1("fit", "pol1" ); + graph->Fit("fit", ""); + + a0[i] = fit->GetParameter(0); + a1[i] = fit->GetParameter(1); + + if( detID < 0 ) { + printf("%2d | a0: %14.10f, a1: %14.10f\n", i, a0[i], a1[i]); + }else{ + printf("%2d | a0, a1 = %14.10f\t%14.10f\n", i, a0[i], a1[i]); + } + } + + if( detID < 0 ){ + FILE * paraOut; + TString filename; + filename.Form("correction_e_auto.dat"); + paraOut = fopen (filename.Data(), "w+"); + printf("=========== save e-correction parameters to %s \n", filename.Data()); + for( int i = 0; i < NCRYSTAL; i++){ + fprintf(paraOut, "%14.10f\t%14.10f\n", a0[i], a1[i]); + } + fflush(paraOut); + fclose(paraOut); + } + +} diff --git a/armory/AutoFit.C b/armory/AutoFit.C new file mode 100644 index 0000000..70cb7fe --- /dev/null +++ b/armory/AutoFit.C @@ -0,0 +1,2787 @@ +/************************************************************ +# +# Created by Tsz Leung (Ryan) Tang +# goluckyryan@gmail.com +# rtang@fsu.edu +# +*************************************************************/ + +#ifndef AutoFit_C +#define AutoFit_C + +#include +#include +#include +#include +#include +#include + +void showFitMethod(){ + printf("----------------------- Method of Fitting ---------------\n"); + printf("---------------------------------------------------------\n"); + printf(" fitAuto() - estimate BG, find peak, and fit n-Gauss \n"); + printf(" fitGaussPol() - fit 1 Gauss + pol-n BG\n"); + printf(" fit2GaussP1() - fit 2 Gauss + pol-1 BG \n"); + printf(" fitGF3Pol() - fit GF3 + pol-n BG \n"); + //printf(" fitNGF3() - fit n-GF3, estimated BG \n"); + printf(" fitNGauss() - fit n-Gauss, estimated BG, need input\n"); + printf(" fitNGaussSub() - fit estimated BG with Pol, subtract, fit n-Guass\n"); + printf(" fitNGaussPol() - fit n-Gauss + pol-n BG \n"); + printf(" fitNGaussPolSub() - subtract Pol-n BG, fit n-Gauss \n"); + printf("\n"); + printf(" fitDecay() - fit A * Exp(-t/T) + B\n"); + printf("\n"); + printf("------- Mouse click Fit : \n"); + printf(" clickFitNGaussPol() - fit n-Gauss + pol-n BG \n"); + printf(" clickFitNGaussPolSub() - Fit Pol-n BG, subtract, fit n-Gauss\n"); + printf(" saveFitPara() - Save the initial guess parameters.\n"); + printf("---------------------------------------------------------\n"); +} + +std::vector SplitStrAF(std::string tempLine, std::string splitter, int shift = 0){ + + std::vector output; + + size_t pos; + do{ + pos = tempLine.find(splitter); // fine splitter + if( pos == 0 ){ //check if it is splitter again + tempLine = tempLine.substr(pos+1); + continue; + } + + std::string secStr; + if( pos == std::string::npos ){ + secStr = tempLine; + }else{ + secStr = tempLine.substr(0, pos+shift); + tempLine = tempLine.substr(pos+shift); + } + + //check if secStr is begin with space + if( secStr.substr(0, 1) == " "){ + secStr = secStr.substr(1); + } + + output.push_back(secStr); + //printf(" |%s---\n", secStr.c_str()); + + }while(pos != std::string::npos ); + + return output; +} + +TColor RGBWheel(double ang){ + + ang = ang * TMath::DegToRad(); + + double r = max(0., (1+2*cos(ang))/3.); + double g = max(0., (1 - cos(ang) + sqrt(3)* sin(ang))/3.); + double b = max(0., (1 - cos(ang) - sqrt(3)* sin(ang))/3.); + + TColor col(r,g,b); + + return col; + +} + +int nPeaks = 16; +Double_t nGauss(Double_t *x, Double_t *par) { + Double_t result = 0; + for (Int_t p=0;pGetNpar(); + int count = totPar/numParPerPeak; + printf("ID "); + for( int i = 0; i < numParPerPeak; i++){ + printf("par%d ", i); + } + printf("\n"); + + for( int i = 0; i < count ; i++){ + printf("%3d ", i); + for( int j = 0; j < numParPerPeak; j++){ + int parID = numParPerPeak * i + j; + printf("%.3f(%.3f) ", fit->GetParameter(parID), fit->GetParError(parID)); + } + printf("\n"); + } +} + +void GoodnessofFit(TH1F * hist, TF1 * fit){ + + int nBin = hist->GetNbinsX(); + int effBin = 0; + double mean = 0; + double ysq = 0; + double SSR = 0; + double chisq = 0; //with estimated error be sqrt(y) + double Xsq = 0; // for Pearson's chi-sq test + for( int i = 1; i <= nBin; i++){ + + double e = hist->GetBinError(i); + if( e > 0 ) { + effBin ++; + double y = hist->GetBinContent(i); + double x = hist->GetBinCenter(i); + double ybar = fit->Eval(x); + ysq += y*y; + mean + y; + SSR += (y - ybar)*(y-ybar); + chisq += (y - ybar)*(y-ybar)/e/e; + + + if( ybar > e ) { + Xsq += (y - ybar)*(y-ybar)/ybar; + }else{ + Xsq += (y - ybar)*(y-ybar)/e; + } + //printf(" %d | x : %f, y : %f, ybar : %f , X-sq : %f\n", i, x, y, ybar, Xsq); + } + } + mean = mean / nBin; + double SSTotal = ysq + mean*mean; + + int npar = fit->GetNpar(); + int ndf = effBin - npar; + printf("#################################################\n"); + printf("## Goodness of Fit. ##\n"); + printf("#################################################\n"); + printf(" eff. Bin(%d) - numPar(%d) = ndf = %d \n", effBin, npar, ndf); + + printf("============== Regression analysis\n"); + printf("----------------- R-sq \n"); + printf(" SSTotal = %f \n", SSTotal); + printf(" SSR = %f \n", SSR); + printf(" MSR = %f <-- is it similar to sample variance?\n", SSR/ndf); + double Rsq = 1 - SSR/SSTotal; + printf(" R-sq = %f \n", Rsq ); + + printf("----------------- Chi-sq \n"); + printf(" Chi-sq = %f \n", chisq); + printf(" rd. Chi-sq = %f \n", chisq/ndf); + printf("ROOT Chi-Sq = %f , NDF = %d \n", fit->GetChisquare(), fit->GetNDF()); + //================ chi-sq test + printf("============== Hypothesis testing\n"); + printf(" Null Hypothesis : the fitting model is truth. \n"); + printf(" * p-value = prob. that Null Hypo. is truth. \n"); + printf(" * the Pearson's test in here only for background free data \n"); + printf(" Pearson's X-sq = %.2f \n", Xsq); + double p = 1 - TMath::Prob(Xsq, ndf); + printf(" Pearson's p-value = %.2f %s 0.05 | %s\n", p, p < 0.05 ? "<": ">", p < 0.05 ? "reject" : "cannot reject"); + double pchi = 1 - TMath::Prob(chisq, ndf); + printf(" Chi-sq p-value = %.2f %s 0.05 | %s\n", pchi, pchi < 0.05 ? "<": ">", pchi < 0.05 ? "reject" : "cannot reject"); + double pRoot = 1- fit->GetProb(); + printf("ROOT Chi-sq p-value = %.2f %s 0.05 | %s\n", pRoot, pRoot < 0.05 ? "<": ">", pRoot < 0.05 ? "reject" : "cannot reject"); + printf("################################################\n"); +} + +vector energy, height, sigma, lowE, highE ; +vector energyFlag, sigmaFlag; + +bool loadFitParameters(TString fitParaFile){ + + energy.clear(); energyFlag.clear(); + sigma.clear(); sigmaFlag.clear(); + + lowE.clear(); highE.clear(); + + height.clear(); + + bool paraRead = false; + + printf("====================================================================== \n"); + printf("----- loading fit parameters from : %s", fitParaFile.Data()); + ifstream file; + file.open(fitParaFile.Data()); + + if( !file){ + printf("\ncannot read file : %s \n", fitParaFile.Data()); + return 0; + } + + while( file.good()) { + string tempLine; + getline(file, tempLine); + + if( tempLine.substr(0, 1) == "#" ) continue; + if( tempLine.substr(0, 2) == "//" ) continue; + if( tempLine.size() < 5 ) continue; + + ///printf("%s\n", tempLine.c_str()); + + vector temp = SplitStrAF(tempLine, " "); + + if( temp.size() < 7 ) continue; + + energy.push_back( atof(temp[0].c_str())); + lowE.push_back( atof(temp[1].c_str())); + highE.push_back( atof(temp[2].c_str())); + energyFlag.push_back(atoi(temp[3].c_str())); + sigma.push_back( atof(temp[4].c_str())); + sigmaFlag.push_back( atoi(temp[5].c_str())); + height.push_back( atof(temp[6].c_str())); + + } + + printf("... done.\n"); + + int n = energy.size(); + TString limStr = "(fixed)"; + printf("%2s| %34s | %10s \n", "ID", "Peak [MeV]", "Sigma [MeV]"); + for( int j = 0; j < n; j ++){ + if( energyFlag[j] == 0 ) limStr.Form("(limited, %6.3f - %6.3f)", lowE[j], highE[j]); + printf("%2d| %7.4f %-26s | %7.4f (%5s) \n", j, energy[j], limStr.Data(), sigma[j], sigmaFlag[j] == 1 ? "fixed" : "free"); + } + + paraRead = true; + + printf("====================================================================== \n"); + + return paraRead; + +} + +TCanvas * NewCanvas(TString name, TString title, int divX, int divY, int padSizeX, int padSizeY){ + TCanvas * output = NULL; + if( gROOT->FindObjectAny(name) == NULL ){ + output = new TCanvas(name, title, divX * padSizeX, divY * padSizeY); + }else{ + output = (TCanvas *) gROOT->FindObjectAny(name) ; + output->Clear(); + } + output->Divide(divX, divY); + return output; +} + +void ScaleAndDrawHist(TH1F * hist, double xMin, double xMax){ + + if ( xMin != xMax ) hist->GetXaxis()->SetRangeUser(xMin, xMax); + int maxBin = hist->GetMaximumBin(); + double ymax = hist->GetBinContent(maxBin); + hist->GetYaxis()->SetRangeUser(0, 1.1 * ymax); + hist->Draw(); + +} + +void PlotResidual(TH1F * hist, TF1 * fit){ + + TH1F * hRes = (TH1F*) hist->Clone(); + hRes->GetListOfFunctions()->Clear(); + hRes->SetTitle("Residual"); + hRes->SetName("hRes"); + hRes->SetYTitle("Hist - fit"); + hRes->Sumw2(0); + hRes->Sumw2(1); + double xMin = hRes->GetXaxis()->GetXmin(); + double xMax = hRes->GetXaxis()->GetXmax(); + fit->SetRange(xMin, xMax); + hRes->Add(fit, -1); + hRes->Draw(); + +} + + +//######################################## +//### Fit a Gauss + Pol-n +//######################################## +void fitGaussPol(TH1F * hist, double mean, double sigmaMax, int degPol, double xMin, double xMax, TString optStat = ""){ + + printf("=========================================================\n"); + printf("================ fit 1-Gauss + Pol-%d BG ================\n", degPol); + printf(" * mean Range +- 5 sigmaMax \n"); + printf(" * inital parameters of the polynomial is random/pow(10, i) \n"); + printf("==========================================================\n"); + + gStyle->SetOptStat(optStat); + TCanvas * cFitGaussPol = NewCanvas("cFitGaussPol", Form("fit Gauss & Pol-%d | fitGaussPol", degPol), 1, 2, 800, 350); + cFitGaussPol->cd(1); + + ScaleAndDrawHist(hist, xMin, xMax); + + const int nPar = 3 + degPol + 1; + + TString funcExp = "[0] * TMath::Gaus(x, [1], [2], 1)"; + for( int i = 0; i < degPol+1 ; i++){ + funcExp += Form(" + [%d]*TMath::Power(x,%d)", i+3 , i); + } + TF1 * fit = new TF1("fit", funcExp.Data(), xMin, xMax); + + double * para = new double[nPar]; + para[0] = 100 * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[1] = mean; + para[2] = sigmaMax/2.; + for( int i = 0 ; i < degPol+1; i++){ + para[3 + i ] = gRandom->Rndm()/TMath::Power(10, i); + } + + fit->SetLineWidth(2); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + fit->SetParLimits(0, 0, 1e9); + fit->SetParLimits(1, mean - 5*sigmaMax, mean + 5 * sigmaMax); + fit->SetParLimits(2, 0, sigmaMax); + + hist->Fit("fit", "Rq"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + double bw = hist->GetBinWidth(1); + + printf("histogram name : %s \n====== Gaussian:\ncount: %8.0f(%3.0f)\nmean : %8.4f(%8.4f)\nsigma: %8.4f(%8.4f) \n", + hist->GetName(), + paraA[0] / bw, paraE[0] /bw, + paraA[1], paraE[1], + paraA[2], paraE[2]); + + printf("------- the polnomail BG:\n"); + for(int i = 0 ; i < degPol+1; i++){ + printf("%2d | %8.4f(%8.4f) \n", i, paraA[3+i], paraE[3+i]); + } + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.DrawLatex(0.12, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + text.DrawLatex(0.12, 0.75,Form("count: %4.0f(%3.0f)", paraA[0] / bw, paraE[0] /bw)); + text.DrawLatex(0.12, 0.70,Form("E_{x}: %6.3f(%5.3f) MeV", paraA[1], paraE[1])); + text.DrawLatex(0.12, 0.65,Form("#sigma: %3.0f(%3.0f) keV", paraA[2] * 1000., paraE[2] * 1000.)); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.60, 0.85 - 0.05*i ,Form("%3s: %8.3f(%8.3f)", Form("p%d", i), paraA[3+i], paraE[3+i])); + } + + TString expression = "[0] "; + for( int j = 1; j < degPol + 1; j++){ + expression += Form(" + [%d]*TMath::Power(x, %d)", j, j); + } + TF1 * g0 = new TF1("g0", expression.Data(), xMin, xMax); + for( int j = 0; j < degPol + 1 ; j++){ + g0->SetParameter(j, paraA[3+j]); + } + g0->SetLineColor(4); + g0->Draw("same"); + + TF1* f0 = new TF1("f0", "[0] * TMath::Gaus(x, [1], [2], 1)", xMin, xMax); + f0->SetParameter(0, paraA[0]); + f0->SetParameter(1, paraA[1]); + f0->SetParameter(2, paraA[2]); + f0->SetLineColor(2); + f0->SetNpx(1000); + f0->Draw("same"); + +// GoodnessofFit(hist, fit); + + cFitGaussPol->cd(2); + PlotResidual(hist, fit); + +} + +//######################################## +//#### fit 2 gauss + pol-1 // not updated +//######################################## +vector fit2GaussP1(TH1F * hist, double mean1, double sigma1, + double mean2, double sigma2, + double xMin, double xMax, TString optStat = "", bool newCanvas = false){ + + + printf("=========================================================\n"); + printf("================ fit 2-Gauss + Pol-1 BG ================\n" ); + printf(" NOT updated. It works, but the code is old \n"); + printf("==========================================================\n"); + + + vector output; + output.clear(); + + gStyle->SetOptStat(optStat); + TCanvas * cFit2GaussP1 = NewCanvas("cFit2GaussP1", "fit Gauss & P1 | fit2GaussP1", 1, 1, 800, 350); + cFit2GaussP1->cd(1); + + ScaleAndDrawHist(hist, xMin, xMax); + + TF1 * fit = new TF1("fit", "[0] * TMath::Gaus(x, [1], [2], 1) + [3] * TMath::Gaus(x, [4], [5], 1) + [6] + [7]*x", xMin, xMax); + + double * para = new double[8]; + para[0] = 20 * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[1] = mean1; + para[2] = sigma1; + para[3] = 100 * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[4] = mean2; + para[5] = sigma2; + para[6] = 1; + para[7] = 0; + + fit->SetLineWidth(2); + fit->SetLineColor(2); + fit->SetNpx(1000); + fit->SetParameters(para); + + hist->Fit("fit", "Rq"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + double bw = hist->GetBinWidth(1); + + printf("%7s ====== count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + hist->GetName(), + paraA[0] / bw, paraE[0] /bw, + paraA[1], paraE[1], + paraA[2], paraE[2]); + printf("%7s ====== count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + "", + paraA[3] / bw, paraE[3] /bw, + paraA[4], paraE[4], + paraA[5], paraE[5]); + + output.push_back( paraA[0]/bw); + output.push_back( paraE[0]/bw); + output.push_back( paraA[1]); + output.push_back( paraE[1]); + output.push_back( paraA[2]); + output.push_back( paraE[2]); + + output.push_back( paraA[3]/bw); + output.push_back( paraE[3]/bw); + output.push_back( paraA[4]); + output.push_back( paraE[4]); + output.push_back( paraA[5]); + output.push_back( paraE[5]); + + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + text.DrawLatex(0.15, 0.75,Form("count: %4.0f(%3.0f), E_{x}: %6.3f(%5.3f) MeV, #sigma: %3.0f(%3.0f) keV ", + paraA[0] / bw, paraE[0] /bw, + paraA[1], paraE[1], + paraA[2] * 1000., paraE[2] * 1000.)); + text.DrawLatex(0.15, 0.7, Form("count: %4.0f(%3.0f), E_{x}: %6.3f(%5.3f) MeV, #sigma: %3.0f(%3.0f) keV ", + paraA[3] / bw, paraE[3] /bw, + paraA[4], paraE[4], + paraA[5] * 1000., paraE[5] * 1000.)); + + text.DrawLatex(0.15, 0.6, Form("Line : %6.3f(%5.3f) + %6.3f(%5.3f)x ", + paraA[6], paraE[6], + paraA[7], paraE[7])); + + GoodnessofFit(hist, fit); + return output; +} + + +//######################################## +//#### fit for gamma + pol-n BG +//######################################## +void fitGF3Pol(TH1F * hist, double mean, double sigmaMax, double ratio, double beta, double step, int degPol, double xMin, double xMax, TString optStat = ""){ + + printf("=========================================================\n"); + printf("================ fit GF1 + Pol-%d BG ================\n", degPol); + printf(" * mean Range = xMin, xMax \n"); + printf(" * inital parameters of the polynomial is random/pow(10, i) \n"); + printf("==========================================================\n"); + + + gStyle->SetOptStat(optStat); + + gStyle->SetOptStat(optStat); + TCanvas * cFitGF3Pol = NewCanvas("cFitGF3Pol", Form("fit GF3 + pol-%d | fitGF3Pol", degPol), 1, 2, 800, 350); + cFitGF3Pol->cd(1); + + ScaleAndDrawHist(hist, xMin, xMax); + + nPeaks = 1; + nPols = degPol; + int nPar = 6*nPeaks + degPol + 1; + + TF1 * fit = new TF1("fit", nGF3Pol, xMin, xMax, nPar); + + fit->Print(); + + double * para = new double[nPar]; + para[0] = hist->GetMaximum() *4; + para[1] = mean; + para[2] = sigmaMax/2.; + para[3] = ratio ; + para[4] = beta ; + para[5] = step ; + for( int i = 0 ; i < degPol + 1; i++){ + para[6+i] = gRandom->Rndm()/TMath::Power(10, i); + } + + fit->SetLineWidth(2); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + fit->SetParLimits(0, 0, 1e9); + fit->SetParLimits(1, xMin, xMax); + fit->SetParLimits(2, 0.00000001, sigmaMax); + fit->SetParLimits(3, 0, 0.5); + fit->SetParLimits(4, 1, 400); + fit->SetParLimits(5, 0, 0.5); + + hist->Fit("fit", "Rq"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + double chisquare = fit->GetChisquare(); + int ndf = fit->GetNDF(); + double bw = hist->GetBinWidth(1); + + printf("histogram : %s \n", hist->GetName()); + printf("========= The Gaussian \n"); + printf("count: %8.0f(%3.0f)\n", paraA[0] / bw, paraE[0] /bw); + printf("mean : %8.4f(%8.4f)\n", paraA[1], paraE[1]); + printf("sigma: %8.4f(%8.4f)\n", paraA[2], paraE[2]); + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(1); + + text.DrawLatex(0.12, 0.65, Form("count : %5.0f(%5.0f)", paraA[0]/bw, paraE[0]/bw)); + text.DrawLatex(0.12, 0.60, Form(" mean : %5.3f(%5.3f) keV", paraA[1], paraE[1])); + text.DrawLatex(0.12, 0.55, Form("sigma : %5.3f(%5.3f) keV", paraA[2], paraE[2])); + text.DrawLatex(0.12, 0.50, Form(" FWHM : %5.3f(%5.3f) keV", paraA[2] *2.355, paraE[2]*2.355)); + + text.DrawLatex(0.12, 0.40, Form("#chi^2/ndf : %5.3f", chisquare/ndf)); + + //GoodnessofFit(hist, fit); + + /// 0 1 2 3 4 5 + string label[8] = {"Area", "mean", "sigma", "ratio", "beta", "step"}; + printf("---------- The detail\n"); + for(int i = 0 ; i < 6 ; i++){ + printf("%d | %8s | %f (%f) \n", i, label[i].c_str(), paraA[i], paraE[i]); + text.DrawLatex(0.65, 0.85-0.05*i, Form("%6s: %5.3f(%5.3f)", label[i].c_str(), paraA[i], paraE[i])); + } + for(int i = 6 ; i < nPar; i++){ + printf("%d | %8s | %f (%f) \n", i, Form("p%d", i-6), paraA[i], paraE[i]); + text.DrawLatex(0.65, 0.85-0.05*i, Form("%6s: %5.3f (%5.3f) \n", Form("p%d", i-6), paraA[i], paraE[i])); + } + + /// norm * (1-ratio)* TMath::Gaus(x[0], mean, sigma, 1) + TF1 * g0 = new TF1("g0", "[0] * (1.0-[3]) * TMath::Gaus(x, [1], [2], 1)", xMin, xMax); + g0->SetParameter(0, paraA[0]); + g0->SetParameter(1, paraA[1]); + g0->SetParameter(2, paraA[2]); + g0->SetParameter(3, paraA[3]); + g0->SetNpx(1000); + g0->SetLineColor(kRed); + + /// norm * ratio * exp( sigma * sigma/2/beta/beta)* exp((x[0]-mean)/beta) * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) + sigma/beta/sqrt(2)) ; + TF1 * g1 = new TF1("g1", "[0] * [3] * exp( [2] * [2]/2/[4]/[4]) / (2* [4])* exp((x-[1])/[4]) * TMath::Erfc( (x-[1])/([2] * sqrt(2)) + [2]/[4]/sqrt(2)) ", xMin, xMax); + g1->SetParameter(0, paraA[0]); + g1->SetParameter(1, paraA[1]); + g1->SetParameter(2, paraA[2]); + g1->SetParameter(3, paraA[3]); + g1->SetParameter(4, paraA[4]); + g1->SetNpx(1000); + g1->SetLineColor(kGreen +3); + + /// norm * step * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) ); + TF1 * g2 = new TF1("g2", "[0] * [3] * TMath::Erfc( (x-[1])/([2] * sqrt(2)) );", xMin, xMax); + g2->SetParameter(0, paraA[0]); + g2->SetParameter(1, paraA[1]); + g2->SetParameter(2, paraA[2]); + g2->SetParameter(3, paraA[5]); + g2->SetNpx(1000); + g2->SetLineColor(kViolet); + + + TString expression = "[0] "; + for( int j = 1; j < degPol + 1; j++){ + expression += Form(" + [%d]*TMath::Power(x, %d)", j, j); + } + TF1 * g3 = new TF1("g3", expression.Data(), xMin, xMax); + for( int j = 0; j < degPol + 1 ; j++){ + g3->SetParameter(j, paraA[6+j]); + } + g3->SetLineColor(kBlue); + g3->Draw("same"); + + + g0->Draw("same"); + g1->Draw("same"); + g2->Draw("same"); + g3->Draw("same"); + + cFitGF3Pol->cd(2); + PlotResidual(hist, fit); + +} + +//############################################## +//##### Auto Fit n-Gauss with estimated BG +//############################################## +vector fitAuto(TH1F * hist, int bgEst = 10, + double peakThreshold = 0.05, + double sigmaMax = 0, + int peakDensity = 4, + TString optStat = ""){ + + printf("================================================================\n"); + printf("========== Auto Fit n-Gauss with estimated BG ==================\n"); + printf(" * bgEst = parameter of BG estimation, larger BG, more linear \n"); + printf(" * peakThreshold = precentage of the highest peak that will count \n"); + printf(" * sigmaMax = maximum sigma, if -1, program try to find the sigma \n"); + printf(" * peakDensity = peak will closer when the number is larger "); + printf(" \n"); + printf(" after peaks found, the i-th peaks will be limited by the mid-point\n"); + printf(" by the (i-1)-th peak and the i-th peak, and the mid-point of the\n"); + printf(" i-th peak and (i+1)-th peak \n"); + printf(" i.e. [peak(i-1)+peak(i)]/2 < limit of peak(i) < [peak(i)+peak(i+1)]/2 \n"); + printf("================================================================\n"); + + gStyle->SetOptStat(optStat); + TCanvas *cFitAuto = NewCanvas("cFitAuto","Auto Fitting | fitAuto", 1, 4, 800, 300); + cFitAuto->cd(1); + + ScaleAndDrawHist(hist, 0, 0); + + TH1F * specS = (TH1F*) hist->Clone(); + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + int xBin = hist->GetXaxis()->GetNbins(); + + TString titleH; + titleH.Form("fitted spectrum (BG=%d); Ex [MeV]; Count / %4.0f keV", bgEst, (xMax-xMin)*1000./xBin ); + specS->SetTitle(titleH); + specS->SetName("specS"); + ///specS->GetXaxis()->SetTitleSize(0.06); + ///specS->GetYaxis()->SetTitleSize(0.06); + ///specS->GetXaxis()->SetTitleOffset(0.7); + ///specS->GetYaxis()->SetTitleOffset(0.6); + + //=================== find peak and fit + gStyle->SetOptFit(0); + TSpectrum * peak = new TSpectrum(50); + nPeaks = peak->Search(hist, peakDensity, "", peakThreshold); + + if( bgEst > 0 ) { + printf("============= estimating background...\n"); + TH1 * h1 = peak->Background(hist, bgEst); + h1->Draw("same"); + printf("============= substracting the linear background...\n"); + specS->Sumw2(); + specS->Add(h1, -1.); + } + + cFitAuto->cd(2)->SetGrid(); + cFitAuto->cd(2); + specS->Draw("hist"); + + //========== Fitting + printf("============= Fitting....."); + printf(" found %d peaks \n", nPeaks); + + double * xpos = peak->GetPositionX(); + double * ypos = peak->GetPositionY(); + + int * inX = new int[nPeaks]; + TMath::Sort(nPeaks, xpos, inX, 0 ); + vector energy, height; + for( int j = 0; j < nPeaks; j++){ + energy.push_back(xpos[inX[j]]); + height.push_back(ypos[inX[j]]); + } + for( int j = 0; j < nPeaks; j++){ + printf(" energy : %f , %f \n", energy[j], height[j]); + } + + + if( sigmaMax == 0 ){ + printf("------------- Estimate sigma for each peak \n"); + sigma.clear(); + int binMin = hist->FindBin(xMin); + int binMax = hist->FindBin(xMax); + for( int i = 0; i < nPeaks ; i++){ + int b0 = hist->FindBin(energy[i]); + double sMin = (xMax-xMin)/5., sMax = (xMax-xMin)/5.; + //---- backward search, stop when + for( int b = b0-1 ; b > binMin ; b-- ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y < (height[i])/2. ) { + sMin = energy[i] - hist->GetBinCenter(b); + break; + } + } + //---- forward search, stop when + for( int b = b0+1 ; b < binMax ; b++ ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y < (height[i])/2. ) { + sMax = hist->GetBinCenter(b) - energy[i]; + break; + } + } + + double temp = TMath::Min(sMin, sMax); + /// When there are multiple peaks closely packed : + if( i > 0 && temp > 2.5 * sigma.back() ) temp = sigma.back(); + sigma.push_back(temp); + + printf("%2d | x : %8.2f | sigma(est) %f \n", i, energy[i], sigma[i]); + } + }else if( sigmaMax < 0 ){ + printf("========== use user input sigma : %f (fixed)\n", abs(sigmaMax)); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(abs(sigmaMax)); + }else if( sigmaMax > 0 ){ + printf("========== use user input sigma : %f/2. \n", sigmaMax/2.); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(sigmaMax/2.); + } + + + + int numParPerPeak = 3; + const int n = numParPerPeak * nPeaks; + double * para = new double[n]; + for(int i = 0; i < nPeaks ; i++){ + para[numParPerPeak*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[numParPerPeak*i+1] = energy[i]; + if( sigmaMax == 0 ){ + para[numParPerPeak*i+2] = sigma[i]; + }else if(sigmaMax < 0 ){ + para[numParPerPeak*i+2] = abs(sigmaMax); + }else if(sigmaMax > 0 ){ + para[numParPerPeak*i+2] = sigmaMax/2.; + } + } + + TF1 * fit = new TF1("fit", nGauss, xMin, xMax, 3 * nPeaks ); + fit->SetLineWidth(2); + fit->SetLineColor(2); + fit->SetNpx(1000); + fit->SetParameters(para); + + if( nPeaks > 1 ){ + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(numParPerPeak*i+0, 0, 1e+9); + double de1 = 1, de2 = 1; + if( i == 0 ){ + de2 = (energy[i+1] - energy[i])/2.; + de1 = de2; + }else if( i < nPeaks -1 ){ + de1 = (energy[i] - energy[i-1])/2.; + de2 = (energy[i+1] - energy[i])/2.; + }else{ + de1 = (energy[i] - energy[i-1])/2.; + de2 = de1; + } + + fit->SetParLimits(numParPerPeak*i+1, energy[i] - de1 , energy[i] + de2); + if( sigmaMax== 0 ) fit->SetParLimits(numParPerPeak*i+2, 0, 1.5*sigma[i]); // add 50% margin of sigma + if( sigmaMax < 0 ) fit->FixParameter(numParPerPeak*i+2, abs(sigmaMax)); + if( sigmaMax > 0 ) fit->SetParLimits(numParPerPeak*i+2, 0, sigmaMax); + } + }else{ + fit->SetParLimits(0, 0, 1e+9); + fit->SetParLimits(2, 0, sigmaMax); + } + + specS->Fit("fit", "q"); + + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + //======== calculate reduce chi-squared + //GoodnessofFit(specS, fit); + + double bw = specS->GetBinWidth(1); + + vector exPos; + for(int i = 0; i < nPeaks ; i++){ + exPos.push_back(paraA[numParPerPeak*i+1]); + printf("%2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[numParPerPeak*i] / bw, paraE[numParPerPeak*i] /bw, + paraA[numParPerPeak*i+1], paraE[numParPerPeak*i+1], + paraA[numParPerPeak*i+2], paraE[numParPerPeak*i+2]); + + } + cFitAuto->Update(); + + //draw the indivual fit + fit->Draw("same"); + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[numParPerPeak*i]); + gFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]); + gFit[i]->SetParameter(2, paraA[numParPerPeak*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + specS->Draw("hist same"); + + //======== print text on plot + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + cFitAuto->cd(3); + PlotResidual(specS, fit); + + cFitAuto->cd(4); + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + } + + return exPos; + +} + + +//######################################## +//###### NOT tested +//######################################## +vector fitNGF3(TH1 * hist, int bgEst = 10, + double peakThreshold = 0.1, + double sigmaMax = 20, + int peakDensity = 4, + TString optStat = "", bool newPlot = true){ + + TCanvas *cFitAuto = NULL; + if( newPlot ){ + cFitAuto = new TCanvas("cFitAuto","Auto Fitting", 100, 100, 800,800); + cFitAuto->Divide(1,2); + + gStyle->SetOptStat(optStat); + cFitAuto->cd(1); + hist->Draw(); + } + + TH1F * specS = (TH1F*) hist->Clone(); + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + int xBin = hist->GetXaxis()->GetNbins(); + + TString titleH; + titleH.Form("fitted spectrum (BG=%d); Ex [MeV]; Count / %4.0f keV", bgEst, (xMax-xMin)*1000./xBin ); + specS->SetTitle(titleH); + specS->SetName("specS"); + ///specS->GetXaxis()->SetTitleSize(0.06); + ///specS->GetYaxis()->SetTitleSize(0.06); + ///specS->GetXaxis()->SetTitleOffset(0.7); + ///specS->GetYaxis()->SetTitleOffset(0.6); + + //=================== find peak and fit + gStyle->SetOptFit(0); + TSpectrum * peak = new TSpectrum(50); + nPeaks = peak->Search(hist, peakDensity, "", peakThreshold); + + if( bgEst > 0 ) { + printf("============= estimating background...\n"); + TH1 * h1 = peak->Background(hist, bgEst); + h1->Draw("same"); + printf("============= substracting the linear background...\n"); + specS->Sumw2(); + specS->Add(h1, -1.); + } + + if( newPlot ){ + cFitAuto->cd(2)->SetGrid(); + cFitAuto->cd(2); + } + specS->Draw("hist"); + + + //========== Fitting + if( newPlot ){ + printf("============= Fitting....."); + printf(" found %d peaks \n", nPeaks); + } + double * xpos = peak->GetPositionX(); + double * ypos = peak->GetPositionY(); + + int * inX = new int[nPeaks]; + TMath::Sort(nPeaks, xpos, inX, 0 ); + vector energy, height; + for( int j = 0; j < nPeaks; j++){ + energy.push_back(xpos[inX[j]]); + height.push_back(ypos[inX[j]]); + } + if( newPlot ){ + for( int j = 0; j < nPeaks; j++){ + printf(" energy : %f , %f \n", energy[j], height[j]); + } + } + + int numParPerPeak = 6; + const int n = numParPerPeak * nPeaks; + double * para = new double[n]; + for(int i = 0; i < nPeaks ; i++){ + para[numParPerPeak*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[numParPerPeak*i+1] = energy[i]; + para[numParPerPeak*i+2] = sigmaMax; + para[numParPerPeak*i+3] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()) * 0.1; + para[numParPerPeak*i+4] = sigmaMax; + para[numParPerPeak*i+5] = -4; + } + + TF1 * fit = new TF1("fit", nGF3, xMin, xMax, numParPerPeak * nPeaks ); + fit->SetLineWidth(2); + fit->SetLineColor(2); + fit->SetNpx(1000); + fit->SetParameters(para); + + if( nPeaks > 1 ){ + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(numParPerPeak*i+0, 0, 1e+9); + double de1 = 1, de2 = 1; + if( i == 0 ){ + de2 = (energy[i+1] - energy[i])/2.; + de1 = de2; + }else if( i < nPeaks -1 ){ + de1 = (energy[i] - energy[i-1])/2.; + de2 = (energy[i+1] - energy[i])/2.; + }else{ + de1 = (energy[i] - energy[i-1])/2.; + de2 = de1; + } + + fit->SetParLimits(numParPerPeak*i+1, energy[i] - de1 , energy[i] + de2); + fit->SetParLimits(numParPerPeak*i+2, 0, sigmaMax * 5); + + fit->SetParLimits(numParPerPeak*i+3, 0, 1e+9); + fit->SetParLimits(numParPerPeak*i+4, 0, sigmaMax); + fit->SetParLimits(numParPerPeak*i+5, -10, -2); + + } + }else{ + fit->SetParLimits(0, 0, 1e+9); + fit->SetParLimits(2, 0, sigmaMax); + fit->SetParLimits(3, 0, 1e+9); + fit->SetParLimits(4, 0, sigmaMax); + fit->SetParLimits(5, -10, -2); + } + specS->Fit("fit", "q"); + + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + //======== calculate reduce chi-squared + if( newPlot ) GoodnessofFit(specS, fit); + + double bw = specS->GetBinWidth(1); + + vector exPos; + for(int i = 0; i < nPeaks ; i++){ + exPos.push_back(paraA[numParPerPeak*i+1]); + double totCount = paraA[numParPerPeak*i] + paraA[numParPerPeak*i+3]; + double totCountErr = sqrt(paraE[numParPerPeak*i]*paraE[numParPerPeak*i] + paraE[numParPerPeak*i+3]*paraE[numParPerPeak*i+3]); + printf("%2d , count: %8.0f(%3.0f)+%8.0f(%3.0f)=%8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f), skewneww: %4.1f(%4.1f) \n", + i, + paraA[numParPerPeak*i] / bw, paraE[numParPerPeak*i] /bw, + paraA[numParPerPeak*i+3] / bw, paraE[numParPerPeak*i+3] /bw, + totCount / bw, totCountErr /bw, + paraA[numParPerPeak*i+1], paraE[numParPerPeak*i+1], + paraA[numParPerPeak*i+2], paraE[numParPerPeak*i+2], + paraA[numParPerPeak*i+5], paraE[numParPerPeak*i+5]); + + //PrintPar(fit, numParPerPeak); + } + if( newPlot ) cFitAuto->Update(); + + //draw the indivual fit + fit->Draw("same"); + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + TF1 ** kFit = new TF1 *[nn]; + TF1 ** zFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1) + [3] * TMath::Gaus(x,[1],[4], 1) * ( 1 + TMath::Erf( [5]*(x-[1])/sqrt(2)/[4] ))", xMin, xMax); + gFit[i]->SetParameter(0, paraA[numParPerPeak*i]); + gFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]); + gFit[i]->SetParameter(2, paraA[numParPerPeak*i+2]); + gFit[i]->SetParameter(3, paraA[numParPerPeak*i+3]); + gFit[i]->SetParameter(4, paraA[numParPerPeak*i+4]); + gFit[i]->SetParameter(5, paraA[numParPerPeak*i+5]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + + kFit[i] = new TF1(Form("kFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1) * ( 1 + TMath::Erf( [3]*(x-[1])/sqrt(2)/[2] ))", xMin, xMax); + kFit[i]->SetParameter(0, paraA[numParPerPeak*i+3]); + kFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]); + kFit[i]->SetParameter(2, paraA[numParPerPeak*i+4]); + kFit[i]->SetParameter(3, paraA[numParPerPeak*i+5]); + kFit[i]->SetLineColor(i+1); + kFit[i]->SetNpx(1000); + kFit[i]->SetLineWidth(1); + kFit[i]->Draw("same"); + + zFit[i] = new TF1(Form("zFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + zFit[i]->SetParameter(0, paraA[numParPerPeak*i]); + zFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]); + zFit[i]->SetParameter(2, paraA[numParPerPeak*i+2]); + zFit[i]->SetLineColor(i+1); + zFit[i]->SetNpx(1000); + zFit[i]->SetLineWidth(1); + zFit[i]->Draw("same"); + + } + + specS->Draw("hist same"); + + return exPos; + +} + +//######################################## +//###### fir N Gauss with estimated BG +//######################################## +void fitNGauss(TH1F * hist, int bgEst = 10, TString fitFile = "AutoFit_para.txt", TString optStat = ""){ + + printf("================================================================\n"); + printf("================ fit N-Gauss with estimated BG ================\n"); + printf(" * bgEst = larger of bgEst, more linear the estimated BG \n"); + printf(" * need the file input \n"); + printf(" \n"); + printf(" 1) The histogram will be subtracted by the estimated BG. \n"); + printf(" 2) n-Gauss will then be fitted the BG subtracted histogram \n"); + printf("================================================================\n"); + + + bool isParaRead = loadFitParameters(fitFile); + if( !isParaRead ) { + printf("Please provide a valid input file\n"); + return; + } + + nPeaks = energy.size(); + + gStyle->SetOptStat(optStat); + TCanvas *cFitNGauss = NewCanvas("cFitNGauss","Fit n-Gauss | fitNGauss", 1,4, 800, 300);; + cFitNGauss->cd(1); + + ScaleAndDrawHist(hist, 0, 0); + + TH1F * specS = (TH1F*) hist->Clone(); + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + int xBin = hist->GetXaxis()->GetNbins(); + + TString titleH; + titleH.Form("fitNGauss (BG = %2d); Ex [MeV]; Count / %4.0f keV", bgEst, (xMax-xMin)*1000./xBin ); + specS->SetTitle(titleH); + specS->SetName("specS"); + + //=================== find peak and fi + + gStyle->SetOptFit(0); + //cFitNGauss->cd(2)->SetGrid(); + cFitNGauss->cd(2); + + if( bgEst > 0 ) { + printf("============= estimating background...\n"); + TSpectrum * peak = new TSpectrum(50); + TH1 * h1 = peak->Background(hist, bgEst); + cFitNGauss->cd(1); + h1->Draw("same"); + cFitNGauss->cd(2); + printf("============= substracting the estimated background...\n"); + specS->Sumw2(); + specS->Add(h1, -1.); + } + + specS->Draw("hist"); + + //========== Fitting + printf("============= Fitting %d-Gauss..... \n", nPeaks); + + const int n = 3 * nPeaks; + double * para = new double[n]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = energy[i]; + para[3*i+2] = sigma[i]/2.; + } + + TF1 * fit = new TF1("fit", nGauss, xMin, xMax, 3* nPeaks ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + + if( energyFlag[i] == 1 ) { + fit->FixParameter(3*i+1, energy[i]); + }else{ + fit->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fit->FixParameter(3*i+2, sigma[i]); + }else{ + fit->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + specS->Fit("fit", "q"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + //======== calculate reduce chi-squared + //GoodnessofFit(specS, fit); + + double bw = specS->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + //draw the indivual fit + specS->Draw("hist"); + fit->Draw("same"); + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + specS->Draw("hist same"); + //specS->Draw("E same"); + + + //======== print text on plot + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.06); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + cFitNGauss->cd(3); + PlotResidual(specS, fit); + + cFitNGauss->cd(4); + + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + } + + cFitNGauss->Update(); +} + + +//######################################## +//#### not updated +//######################################## +void fitNGaussSub(TH1F * hist, int bgEst = 10, int degPol = 1, TString fitFile = "AutoFit_para.txt", TString optStat = ""){ + + + printf("==================================================================\n"); + printf("======== fit N-Gauss with estimated BG (method-2) ================\n"); + printf(" * bgEst = larger of bgEst, more linear the estimated BG \n"); + printf(" * degPol = degree of polynomial \n"); + printf(" * need the file input \n"); + printf(" \n"); + printf(" 1) A BG is estimated, and then the BG is fitted by a polynomial. \n"); + printf(" 2) The histogram will be subtracted by the polynomial. \n"); + printf(" 3) n-Gauss will then be fitted the subtracted histogram \n"); + printf("================================================================\n"); + + + bool isParaRead = loadFitParameters(fitFile); + if( !isParaRead ) { + printf("Please provide a valid input file\n"); + return; + } + + nPeaks = energy.size(); + nPols = degPol; + + gStyle->SetOptStat(optStat); + TCanvas *cFitNGaussSub = NewCanvas("cFitNGaussSub","Fit n-Gauss, replace estimated BG with Pol-n | fitNGauss2", 1, 4, 800, 300 ); + //if(! cFitNGaussSub->GetShowEventStatus() ) cFitNGaussSub->ToggleEventStatus(); + + cFitNGaussSub->cd(1); + ScaleAndDrawHist(hist, 0, 0); + + TH1F * specS = (TH1F*) hist->Clone(); + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + int xBin = hist->GetXaxis()->GetNbins(); + + TString titleH; + titleH.Form("fitNGauss2 (replace Est. BG with Pol-%d) (BG = %2d); Ex [MeV]; Count / %4.0f keV", degPol, bgEst, (xMax-xMin)*1000./xBin ); + specS->SetTitle(titleH); + specS->SetName("specS"); + + printf("============= estimating background...\n"); + TSpectrum * peak = new TSpectrum(50); + TH1 * h1 = peak->Background(hist, bgEst); + + printf("============= fit the est-background with a polynomial function...\n"); + + TString polExp = "[0]"; + for( int i = 1; i < degPol + 1; i++){ + polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + } + TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax); + + bg->SetParameter(0, 50); + bg->SetParameter(0, 0); + bg->SetLineColor(2); + bg->SetNpx(1000); + h1->Fit("bg", "R", ""); + + hist->Draw(); + bg->Draw("same"); + + + //======== print text on plot + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + + const Double_t * paraAt = bg->GetParameters(); + const Double_t * paraEt = bg->GetParErrors(); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraAt[i], paraEt[i])); + } + + gStyle->SetOptFit(0); +// cFitNGaussSub->cd(2)->SetGrid(); + cFitNGaussSub->cd(2); + + printf("============= substracting the polynomial background...\n"); + specS->Sumw2(); + specS->Add(bg, -1.); + specS->Draw("hist"); + + //========== Fitting + printf("============= Fitting..... \n"); + + const int n = 3 * nPeaks; + double * para = new double[n]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = energy[i]; + para[3*i+2] = sigma[i]/2.; + } + + TF1 * fit = new TF1("fit", nGauss, xMin, xMax, 3* nPeaks ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + + if( energyFlag[i] == 1 ) { + fit->FixParameter(3*i+1, energy[i]); + }else{ + fit->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fit->FixParameter(3*i+2, sigma[i]); + }else{ + fit->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + specS->Fit("fit", "q"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + GoodnessofFit(specS, fit); + + double bw = specS->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + + //draw the indivual fit + specS->Draw("hist"); + fit->Draw("same"); + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + specS->Draw("hist same"); + //specS->Draw("E same"); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + + cFitNGaussSub->cd(3); + PlotResidual(specS, fit); + + cFitNGaussSub->cd(4); + + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + } + +} + + +//######################################## +//#### fit N-Gauss with pol-n BG +//######################################## +void fitNGaussPol(TH1F * hist, int degPol, TString fitFile = "AutoFit_para.txt",double xMin = 0, double xMax = 0, TString optStat = ""){ + + printf("================================================================\n"); + printf("================ fit N-Gauss with Pol-%1d BG ==================\n", degPol); + printf(" * degPol = degree of polynomial \n"); + printf(" * need the file input \n"); + printf(" * xMin, xMax = if left empty, full range will be used\n"); + printf(" \n"); + printf(" 1) The histogram will be fitted by n-Gauss + Pol \n"); + printf("================================================================\n"); + + + bool isParaRead = loadFitParameters(fitFile); + if( !isParaRead ) { + printf("Please provide a valid input file\n"); + return; + } + + gStyle->SetOptStat(optStat); + nPeaks = energy.size(); + nPols = degPol; + + TCanvas * cFitNGaussPol = NewCanvas("cFitNGaussPol", Form("Fitting with n-Gauss + pol-%d | fitNGaussPol", degPol), 1, 3, 800, 300); + //if(! cFitNGaussPol->GetShowEventStatus() ) cFitNGaussPol->ToggleEventStatus(); + cFitNGaussPol->cd(1); + + ScaleAndDrawHist(hist, xMin, xMax); + + if( xMin == xMax){ + xMin = hist->GetXaxis()->GetXmin(); + xMax = hist->GetXaxis()->GetXmax(); + } + int xBin = hist->GetXaxis()->GetNbins(); + + + printf("============= find the polynomial background ..... \n"); + + int nPar = 3 * nPeaks + nPols + 1; + double * para = new double[nPar]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = energy[i]; + para[3*i+2] = sigma[i]/2.; + } + + for(int i = 0; i < nPols + 1; i++){ + //if( ggg == NULL ){ + para[3*nPeaks+i] = gRandom->Rndm()/TMath::Power(10, 3*i); + //}else{ + // para[3*nPeaks+i] = gPara[i]; + //} + } + + TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + if( energyFlag[i] == 1 ) { + fit->FixParameter(3*i+1, energy[i]); + }else{ + fit->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fit->FixParameter(3*i+2, sigma[i]); + }else{ + fit->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + hist->Fit("fit", "Rq"); + + //=========== get the polynomial part + const Double_t* paraA = fit->GetParameters(); + const Double_t* paraE = fit->GetParErrors(); + + TString polExp = "[0]"; + for( int i = 1; i < degPol + 1; i++){ + polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + } + + TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax); + for( int i = 0; i < degPol + 1; i++){ + bg->SetParameter(i, paraA[3*nPeaks+i]); + } + bg->SetNpx(1000); + + for( int i = 0; i < degPol + 1; i++){ + printf("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i]); + } + printf("====================================\n"); + + cFitNGaussPol->cd(1); + bg->Draw("same"); + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(1); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i])); + } + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + + //GoodnessofFit(specS, fit); + + double bw = hist->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + cFitNGaussPol->Update(); + + cFitNGaussPol->cd(2); + PlotResidual(hist, fit); + + cFitNGaussPol->cd(3); + + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + } + +} + + + + +//######################################## +//#### fit N-Gauss with pol-n BG +//######################################## +void fitNGaussPolSub(TH1F * hist, int degPol, TString fitFile = "AutoFit_para.txt",double xMin = 0, double xMax = 0, TString optStat = ""){ + + printf("================================================================\n"); + printf("========= fit N-Gauss with Pol-%1d BG / 2nd method ============\n", degPol); + printf(" * degPol = degree of polynomial \n"); + printf(" * need the file input \n"); + printf(" * xMin, xMax = if left empty, full range will be used\n"); + printf(" \n"); + printf(" 1) The histogram will be fitted by n-Gauss + Pol -> to get the estimated BG \n"); + printf(" 2) The histogram will be subtracted by the Pol BG. \n"); + printf(" 3) n-Gauss will then be fitted the BG subtracted histogram \n"); + printf("================================================================\n"); + + + bool isParaRead = loadFitParameters(fitFile); + if( !isParaRead ) { + printf("Please provide a valid input file\n"); + return; + } + + gStyle->SetOptStat(optStat); + nPeaks = energy.size(); + nPols = degPol; + + TCanvas * cFitNGaussPolSub = NewCanvas("cFitNGaussPolSub", Form("Fitting with n-Gauss + pol-%d (method-2) | fitGaussPol2", degPol), 1, 4, 800, 300); + //if(! cFitNGaussPol->GetShowEventStatus() ) cFitNGaussPol->ToggleEventStatus(); + cFitNGaussPolSub->cd(1); + + ScaleAndDrawHist(hist, xMin, xMax); + + if( xMin == xMax){ + xMin = hist->GetXaxis()->GetXmin(); + xMax = hist->GetXaxis()->GetXmax(); + } + int xBin = hist->GetXaxis()->GetNbins(); + + printf("============= find the polynomial background ..... \n"); + + int nPar = 3 * nPeaks + nPols + 1; + double * para = new double[nPar]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = energy[i]; + para[3*i+2] = sigma[i]/2.; + } + + for(int i = 0; i < nPols + 1; i++){ + para[3*nPeaks+i] = gRandom->Rndm()/TMath::Power(10, 3*i); + } + + TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + if( energyFlag[i] == 1 ) { + fit->FixParameter(3*i+1, energy[i]); + }else{ + fit->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fit->FixParameter(3*i+2, sigma[i]); + }else{ + fit->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + hist->Fit("fit", "nq"); + + //=========== get the polynomial part and substract + const Double_t* paraAt = fit->GetParameters(); + const Double_t* paraEt = fit->GetParErrors(); + + TString polExp = "[0]"; + for( int i = 1; i < degPol + 1; i++){ + polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + } + + TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax); + for( int i = 0; i < degPol + 1; i++){ + bg->SetParameter(i, paraAt[3*nPeaks+i]); + } + bg->SetNpx(1000); + + for( int i = 0; i < degPol + 1; i++){ + printf("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraAt[3*nPeaks+i], paraEt[3*nPeaks+i]); + } + printf("====================================\n"); + + cFitNGaussPolSub->cd(1); + bg->Draw("same"); + + TH1F * specS = (TH1F*) hist->Clone(); + TString titleH; + titleH.Form("pol-%d BG Subtracted spectrum (fitNGaussPol-2); Ex [MeV]; Count / %4.0f keV", degPol, (xMax-xMin)*1000./xBin ); + specS->SetTitle(titleH); + specS->SetName("specS"); + + //=================== find peak and fit + + gStyle->SetOptFit(0); + ///cFitNGaussPol->cd(2)->SetGrid(); + cFitNGaussPolSub->cd(2); + + printf("============= substracting the polynomial background...\n"); + specS->Sumw2(); + specS->Add(bg, -1.); + specS->Draw("hist"); + + //======= fit again + printf("============= fitting the subtracted spectrum.... \n"); + + nPar = 3* nPeaks; + TF1 * fita = new TF1("fita", nGauss, xMin, xMax, nPar ); + fita->SetLineWidth(3); + fita->SetLineColor(1); + fita->SetNpx(1000); + fita->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fita->SetParLimits(3*i , 0, 1e9); + + if( energyFlag[i] == 1 ) { + fita->FixParameter(3*i+1, energy[i]); + }else{ + fita->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fita->FixParameter(3*i+2, sigma[i]); + }else{ + fita->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + specS->Fit("fita", "q"); + + const Double_t* paraE = fita->GetParErrors(); + const Double_t* paraA = fita->GetParameters(); + + //GoodnessofFit(specS, fit); + + double bw = specS->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + + //draw the indivual fit + specS->Draw("hist"); + fita->Draw("same"); + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + specS->Draw("hist same"); + //specS->Draw("E same"); + + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(1); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i])); + } + + double chi2 = fita->GetChisquare(); + int ndf = fita->GetNDF(); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + + cFitNGaussPolSub->cd(3); + PlotResidual(specS, fita); + + cFitNGaussPolSub->cd(4); + + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + } + +} + +//######################################## +//###### fit A * Exp(-t/T) + B +//######################################## +void fitDecay(TH1F * hist){ + + gStyle->SetOptStat(""); + TCanvas * cFitDecay = NULL; + if( gROOT->FindObjectAny("cFitDecay") == NULL ){ + cFitDecay = new TCanvas("cFitDecay", Form("fit A * Exp(-t/T) + B | fitDecay"), 800, 1000); + }else{ + delete gROOT->FindObjectAny("cFitDecay") ; + cFitDecay = new TCanvas("cFitDecay", Form("fit A * Exp(-t/T) + B | fitDecay"), 800, 1000); + } + + cFitDecay->Divide(1, 2); + + double yMax = hist->GetBinContent(hist->GetMaximumBin()); + hist->SetMaximum(yMax*1.2); + cFitDecay->cd(1); hist->Draw(); + double estHalf = hist->GetBinCenter(hist->FindLastBinAbove(yMax/2.)) *1.3; + + //printf("estimator : A = %.2f, T = %.2f \n", yMax, estHalf); + + TF1 * fit = new TF1("fit", "[0]*TMath::Exp(-x/[1]/TMath::Log(2))+[2]"); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + + fit->SetParameter(0, yMax ); + fit->SetParameter(1, estHalf ); + fit->SetParameter(2, 1); + + fit->SetParLimits(0, 0, 1e9); + fit->SetParLimits(1, 0, 1e9); + fit->SetParLimits(2, 0, 1e9); + + gStyle->SetOptFit(000000); + hist->Fit("fit", "q"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(1); + + text.DrawLatex(0.6, 0.8, "A e^{#left(-#frac{t}{T Ln(2)}#right)}+B"); + text.DrawLatex(0.6, 0.75, Form("A : %.f(%.f)", paraA[0], paraE[0])); + text.DrawLatex(0.6, 0.7, Form("T : %.f(%.f)", paraA[1], paraE[1])); + text.DrawLatex(0.6, 0.65, Form("B : %.f(%.f)", paraA[2], paraE[2])); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.DrawLatex(0.6, 0.6, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + + TF1 * fite = new TF1("fite", "[0]*TMath::Exp(-x/[1])" ); + fite->SetParameter(0, paraA[0]); + fite->SetParameter(1, paraA[1]); + fite->Draw("same"); + + TF1 * fitbg = new TF1("fitbg", "pol0" ); + fitbg->SetParameter(0, paraA[2]); + fitbg->Draw("same"); + + cFitDecay->cd(2); + PlotResidual(hist, fit); + + +} + +//################################################# +//#### fit N-Gauss with pol-n BG using mouse click +//################################################# + +int nClick = 0; +bool peakFlag = 1; +vector xPeakList; +vector yPeakList; +vector xBGList; +vector yBGList; + +TH1F * tempHist; +int markerStyle = 23; +int markerColor = 2; +int markerBGColor = 4; + +void Clicked() { + int event = gPad->GetEvent(); + if (event != 11) return; + TObject *select = gPad->GetSelected(); + if (!select) return; + + TH1F *h = (TH1F*)select; + int px = gPad->GetEventX(); + double xd = gPad->AbsPixeltoX(px); + float x = gPad->PadtoX(xd); + + if( peakFlag ) { + xPeakList.push_back(x); + }else{ + xBGList.push_back(x); + } + int b = tempHist->FindBin(x); + double y = tempHist->GetBinContent(b); + if( peakFlag ){ + yPeakList.push_back(y); + }else{ + yBGList.push_back(y); + } + // add marker in the histogram + TMarker * mark = new TMarker(x,y, markerStyle); + mark->SetMarkerColor(markerColor); + tempHist->GetListOfFunctions()->Add(mark); + + printf("%2d | x : %8.2f , y : %.0f \n", nClick, x, y); + + nClick ++; + +} + + +void saveFitPara(TString fileName = "AutoFit_para.txt"){ + printf("Save to : %s \n", fileName.Data()); + FILE * file_out; + file_out = fopen (fileName.Data(), "w+"); + + fprintf(file_out, "# for n-Gauss fit, can use \"#\", or \"//\" to comment out whole line\n"); + fprintf(file_out, "# peak low high fixed? sigma_Max fixed? hight\n"); + + for( int i = 0 ; i < xPeakList.size() ; i++){ + fprintf(file_out, "%.3f %.3f %.3f 0 %.3f 0 %.0f\n", + xPeakList[i], + xPeakList[i] - 5*sigma[i], + xPeakList[i] + 5*sigma[i], + sigma[i], + yPeakList[i]); + } + fclose(file_out); +} + +void clickFitNGaussPol(TH1F * hist, int degPol, double sigmaMax = 0){ + + printf("=========================================================\n"); + printf("======== fit n-Gauss + Pol-%d BG using mouse click =====\n", degPol ); + printf("==========================================================\n"); + + gStyle->SetOptStat(""); + gStyle->SetOptFit(0); + + TCanvas * cClickFitNGaussPol = NULL; + if( gROOT->FindObjectAny("cClickFitGaussPol") == NULL ){ + cClickFitNGaussPol = new TCanvas("cClickFitNGaussPol", Form("fit Gauss & Pol-%d by mouse click | clickFitGaussPol", degPol), 1400, 1200); + }else{ + delete gROOT->FindObjectAny("cClickFitNGaussPol") ; + cClickFitNGaussPol = new TCanvas("cClickFitNGaussPol", Form("fit Gauss & Pol-%d by mouse click | clickFitGaussPol", degPol), 1400, 1200); + } + cClickFitNGaussPol->Divide(1, 4); + for(int i = 1; i <= 4 ; i++){ + cClickFitNGaussPol->cd(i)->SetGrid(0,0); + } + + if(! cClickFitNGaussPol->GetShowEventStatus() ) cClickFitNGaussPol->ToggleEventStatus(); + cClickFitNGaussPol->cd(1); + + + hist->GetListOfFunctions()->Clear(); + ScaleAndDrawHist(hist, 0, 0); + TH1F* hspec = (TH1F*) hist->Clone(); + hspec->Sumw2(); + cClickFitNGaussPol->Update(); + cClickFitNGaussPol->Draw(); + + TLatex helpMsg; + helpMsg.SetNDC(); + helpMsg.SetTextFont(82); + helpMsg.SetTextSize(0.06); + helpMsg.SetTextColor(kRed); + helpMsg.DrawLatex(0.15, 0.8, "Click for peaks: (Double-click / x / Ctrl to end) "); + + printf("--------- Click the histogram for peaks: (Double-click / x / Ctrl to end) \n"); + nClick = 0; + xPeakList.clear(); + yPeakList.clear(); + markerColor = 2; + markerStyle = 23; + peakFlag = 1; + cClickFitNGaussPol->cd(1)->SetCrosshair(1); + cClickFitNGaussPol->cd(1)->AddExec("ex", "Clicked()"); + tempHist = hist; + TObject * obj ; + do{ + obj = gPad->WaitPrimitive(); + if( obj == NULL ) break; + }while( obj != NULL); + + if( degPol >= 0 ){ + printf("--------- Click the histogram for Background: (Double-click / x / Ctrl to end) \n"); + printf(" * when no input, program will estimate \n"); + + cClickFitNGaussPol->cd(1)->Clear(); + hist->Draw(); + helpMsg.SetTextColor(markerBGColor); + helpMsg.DrawLatex(0.15, 0.8, "Click for BG: (Double-click / x / Ctrl to end) "); + helpMsg.DrawLatex(0.15, 0.75, "* when no input, program will estimate"); + cClickFitNGaussPol->Update(); + cClickFitNGaussPol->Draw(); + + nClick = 0; + xBGList.clear(); + yBGList.clear(); + markerColor = markerBGColor; + markerStyle = 33; + peakFlag = 0; + cClickFitNGaussPol->cd(1)->AddExec("ex", "Clicked()"); + do{ + obj = gPad->WaitPrimitive(); + if( obj == NULL ) break; + }while( obj != NULL); + + } + + cClickFitNGaussPol->cd(1)->DeleteExec("ex"); + cClickFitNGaussPol->cd(1)->SetCrosshair(0); + cClickFitNGaussPol->cd(1)->Clear(); + hist->Draw(); + + tempHist = NULL; + + cClickFitNGaussPol->Update(); + cClickFitNGaussPol->Draw(); + + nPols = degPol; + double xMin = hspec->GetXaxis()->GetXmin(); + double xMax = hspec->GetXaxis()->GetXmax(); + + TString polExp = "[0]"; + for( int i = 1; i < degPol + 1; i++) polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + TF1 *bg = new TF1("bg", polExp.Data(), xMin, xMax); + bg->SetNpx(1000); + bg->SetLineColor(markerBGColor); + bg->SetLineWidth(1); + if( xBGList.size() > 0 ) { + printf("---------------- fit the BG with Pol-%d \n", nPols); + TGraph * gBG = new TGraph((int) xBGList.size(), &xBGList[0], &yBGList[0]); + for( int i = 0; i < degPol + 1; i++) bg->SetParameter(i, gRandom->Rndm()/TMath::Power(10, 3*i)); + gBG->Fit("bg", "Rq"); + bg->Draw("same"); + //printf("--------------- Subtracting the BG \n"); + //hspec->Add(bg, -1); + }else{ + for( int i = 0; i < nPols+1; i++) bg->SetParameter(i, 0.); + } + + nPeaks = (int) xPeakList.size(); + if( sigmaMax == 0 ){ + printf("------------- Estimate sigma for each peak \n"); + sigma.clear(); + int binMin = hist->FindBin(xMin); + int binMax = hist->FindBin(xMax); + for( int i = 0; i < nPeaks ; i++){ + int b0 = hist->FindBin(xPeakList[i]); + double estBG = bg->Eval(xPeakList[i]); + double sMin = (xMax-xMin)/5., sMax = (xMax-xMin)/5.; + //---- backward search, stop when + for( int b = b0-1 ; b > binMin ; b-- ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) { + sMin = xPeakList[i] - hist->GetBinCenter(b); + break; + } + } + //---- forward search, stop when + for( int b = b0+1 ; b < binMax ; b++ ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) { + sMax = hist->GetBinCenter(b) - xPeakList[i]; + break; + } + } + + double temp = TMath::Min(sMin, sMax); + /// When there are multiple peaks closely packed : + if( i > 0 && temp > 2.5 * sigma.back() ) temp = sigma.back(); + sigma.push_back(temp); + + printf("%2d | x : %8.2f | sigma(est) %f \n", i, xPeakList[i], sigma[i]); + + } + + //---- use the mean of the sigma + double sigma0 = 0; + for( int i = 0; i < nPeaks ; i++) sigma0 += sigma[i]; + sigma0 = sigma0/(nPeaks+1); + for( int i = 0; i < nPeaks ; i++) sigma[i] = sigma0; + printf("========== use the mean sigma : %f \n", sigma0); + }else if( sigmaMax < 0 ){ + printf("========== use user input sigma : %f (fixed)\n", abs(sigmaMax)); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(abs(sigmaMax)); + }else if( sigmaMax > 0 ){ + printf("========== use user input sigma : %f/2. \n", sigmaMax/2.); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(sigmaMax/2.); + } + + printf("-------------- Fit the spectrum with %d-Gauss + Pol-%d\n", nPeaks, nPols ); + cClickFitNGaussPol->cd(2); + hspec->Draw("hist"); + + int nPar = 3 * nPeaks + nPols + 1; + double * para = new double[nPar]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = yPeakList[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = xPeakList[i]; + if( sigmaMax == 0){ + para[3*i+2] = sigma[i]; + }else if(sigmaMax < 0 ){ + para[3*i+2] = abs(sigmaMax); + }else if(sigmaMax > 0 ){ + para[3*i+2] = sigmaMax/2.; + } + } + for(int i = 0; i < nPols+1 ; i++){ + para[3*nPeaks+i] = bg->GetParameter(i); + } + + TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //limit parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + + double dE1, dE2; + if( i == 0 ){ + dE2 = xPeakList[i+1] - xPeakList[i]; + dE1 = dE2; + }else if ( i == nPeaks-1 ){ + dE1 = xPeakList[i] - xPeakList[i-1]; + dE2 = dE1; + }else{ + dE1 = xPeakList[i] - xPeakList[i-1]; + dE2 = xPeakList[i+1] - xPeakList[i]; + } + + fit->SetParLimits(3*i+1, xPeakList[i] - dE1 , xPeakList[i] + dE2 ); + if( sigmaMax== 0 ) fit->SetParLimits(3*i+2, 0, 1.5*sigma[i]); // add 50% margin of sigma + if( sigmaMax < 0 ) fit->FixParameter(3*i+2, abs(sigmaMax)); + if( sigmaMax > 0 ) fit->SetParLimits(3*i+2, 0, sigmaMax); + } + + hspec->Fit("fit", "Rq"); + fit->Draw("same"); + + //=========== get the polynomial part + const Double_t* paraA = fit->GetParameters(); + const Double_t* paraE = fit->GetParErrors(); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + double bw = hspec->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i])); + } + + TF1 ** gFit = new TF1 *[nPeaks]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + if( degPol >= 0 ){ + TF1 *bg2 = new TF1("bg", polExp.Data(), xMin, xMax); + bg2->SetNpx(1000); + bg2->SetLineColor(markerBGColor); + bg2->SetLineWidth(1); + for( int i = 0; i < nPols + 1; i++){ + bg2->SetParameter(i, paraA[3*nPeaks+i]); + } + bg2->Draw("same"); + } + + cClickFitNGaussPol->cd(3); + PlotResidual(hspec, fit); + + cClickFitNGaussPol->cd(4); + + text.SetTextSize(0.05); + text.SetTextColor(2); + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + } +} + + + +void clickFitNGaussPolSub(TH1F * hist, int degPol, double sigmaMax = 0){ + + printf("=========================================================\n"); + printf("= fit n-Gauss + Pol-%d BG using mouse click (method-2) =\n", degPol ); + printf("==========================================================\n"); + + gStyle->SetOptStat(""); + gStyle->SetOptFit(0); + + TCanvas * cClickFitNGaussPolsub = NULL; + if( gROOT->FindObjectAny("cClickFitGaussPolsub") == NULL ){ + cClickFitNGaussPolsub = new TCanvas("cClickFitNGaussPol", Form("fit Gauss & Pol-%d by mouse click (sub) | clickFitGaussPolsub", degPol), 1400, 1200); + }else{ + delete gROOT->FindObjectAny("cClickFitNGaussPolsub") ; + cClickFitNGaussPolsub = new TCanvas("cClickFitNGaussPolsub", Form("fit Gauss & Pol-%d by mouse click (sub) | clickFitGaussPolsub", degPol), 1400, 1200); + } + cClickFitNGaussPolsub->Divide(1, 4); + for(int i = 1; i <= 4 ; i++){ + cClickFitNGaussPolsub->cd(i)->SetGrid(0,0); + } + + if(! cClickFitNGaussPolsub->GetShowEventStatus() ) cClickFitNGaussPolsub->ToggleEventStatus(); + cClickFitNGaussPolsub->cd(1); + + + hist->GetListOfFunctions()->Clear(); + ScaleAndDrawHist(hist, 0, 0); + TH1F* hspec = (TH1F*) hist->Clone(); + hspec->Sumw2(); + cClickFitNGaussPolsub->Update(); + cClickFitNGaussPolsub->Draw(); + + TLatex helpMsg; + helpMsg.SetNDC(); + helpMsg.SetTextFont(82); + helpMsg.SetTextSize(0.06); + helpMsg.SetTextColor(kRed); + helpMsg.DrawLatex(0.15, 0.8, "Click for peaks: (Double-click / x / Ctrl to end) "); + + + printf("--------- Click the histogram for peaks: (Double-click / x / Ctrl to end) \n"); + nClick = 0; + xPeakList.clear(); + yPeakList.clear(); + markerColor = 2; + markerStyle = 23; + peakFlag = 1; + cClickFitNGaussPolsub->cd(1)->SetCrosshair(1); + cClickFitNGaussPolsub->cd(1)->AddExec("ex", "Clicked()"); + tempHist = hist; + TObject * obj ; + do{ + obj = gPad->WaitPrimitive(); + if( obj == NULL ) break; + }while( obj != NULL); + + if( degPol >= 0 ){ + printf("--------- Click the histogram for Background: (Double-click / x / Ctrl to end) \n"); + + cClickFitNGaussPolsub->cd(1)->Clear(); + hist->Draw(); + helpMsg.SetTextColor(markerBGColor); + helpMsg.DrawLatex(0.15, 0.8, "Click for BG: (Double-click / x / Ctrl to end) "); + cClickFitNGaussPolsub->Update(); + cClickFitNGaussPolsub->Draw(); + + nClick = 0; + xBGList.clear(); + yBGList.clear(); + markerColor = markerBGColor; + markerStyle = 33; + peakFlag = 0; + do{ + obj = gPad->WaitPrimitive(); + if( obj == NULL ) break; + }while( obj != NULL); + } + cClickFitNGaussPolsub->cd(1)->DeleteExec("ex"); + cClickFitNGaussPolsub->cd(1)->SetCrosshair(0); + cClickFitNGaussPolsub->cd(1)->Clear(); + hist->Draw(); + + tempHist = NULL; + + if( xBGList.size() == 0 ) helpMsg.DrawLatex(0.15, 0.75, " No BG defined ! fitting could be problematics. "); + + cClickFitNGaussPolsub->Update(); + cClickFitNGaussPolsub->Draw(); + + nPols = degPol; + double xMin = hspec->GetXaxis()->GetXmin(); + double xMax = hspec->GetXaxis()->GetXmax(); + + TString polExp = "[0]"; + for( int i = 1; i < degPol + 1; i++) polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + TF1 *bg = new TF1("bg", polExp.Data(), xMin, xMax); + bg->SetNpx(1000); + bg->SetLineColor(markerBGColor); + bg->SetLineWidth(1); + if( xBGList.size() > 0 ) { + printf("---------------- fit the BG with Pol-%d \n", nPols); + + TGraph * gBG = new TGraph((int) xBGList.size(), &xBGList[0], &yBGList[0]); + + for( int i = 0; i < degPol + 1; i++) bg->SetParameter(i, gRandom->Rndm()/TMath::Power(10, 3*i)); + + gBG->Fit("bg", "Rq"); + bg->Draw("same"); + + printf("--------------- Subtracting the BG \n"); + hspec->Add(bg, -1); + }else{ + for( int i = 0; i < nPols+1; i++) bg->SetParameter(i, 0.); + } + + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.06); + + double chi2BG = bg->GetChisquare(); + int ndfBG = bg->GetNDF(); + + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2BG/ndfBG)); + + //=========== get the polynomial BG + const Double_t* paraAt = bg->GetParameters(); + const Double_t* paraEt = bg->GetParErrors(); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.2f(%8.2f)\n", Form("p%d", i), paraAt[i], paraEt[i])); + } + + nPeaks = (int) xPeakList.size(); + if( sigmaMax == 0 ){ + printf("------------- Estimate sigma for each peak \n"); + sigma.clear(); + int binMin = hist->FindBin(xMin); + int binMax = hist->FindBin(xMax); + + for( int i = 0; i < nPeaks ; i++){ + int b0 = hist->FindBin(xPeakList[i]); + double estBG = bg->Eval(xPeakList[i]); + double sMin = (xMax-xMin)/5., sMax = (xMax-xMin)/5.; + //---- backward search, stop when + for( int b = b0-1 ; b > binMin ; b-- ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) { + sMin = xPeakList[i] - hist->GetBinCenter(b); + break; + } + } + //---- forward search, stop when + for( int b = b0+1 ; b < binMax ; b++ ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) { + sMax = hist->GetBinCenter(b) - xPeakList[i]; + break; + } + } + + double temp = TMath::Min(sMin, sMax); + /// When there are multiple peaks closely packed : + if( i > 0 && temp > 2.5 * sigma.back() ) temp = sigma.back(); + sigma.push_back(temp); + + printf("%2d | x : %8.2f | sigma(est) %f \n", i, xPeakList[i], sigma[i]); + + } + }else if( sigmaMax < 0 ){ + printf("========== use user input sigma : %f (fixed)\n", abs(sigmaMax)); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(abs(sigmaMax)); + }else if( sigmaMax > 0 ){ + printf("========== use user input sigma : %f/2. \n", sigmaMax/2.); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(sigmaMax/2.); + } + + printf("-------------- Fit the subtracted spectrum with %d-Gauss\n", nPeaks ); + cClickFitNGaussPolsub->cd(2); + hspec->Draw("hist"); + + int nPar = 3 * nPeaks; + double * para = new double[nPar]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = yPeakList[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = xPeakList[i]; + if( sigmaMax == 0){ + para[3*i+2] = sigma[i]; + }else if(sigmaMax < 0 ){ + para[3*i+2] = abs(sigmaMax); + }else if(sigmaMax > 0 ){ + para[3*i+2] = sigmaMax/2.; + } + } + + TF1 * fit = new TF1("fit", nGauss, xMin, xMax, nPar ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //limit parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + + double dE1, dE2; + if( i == 0 ){ + dE2 = xPeakList[i+1] - xPeakList[i]; + dE1 = dE2; + }else if ( i == nPeaks-1 ){ + dE1 = xPeakList[i] - xPeakList[i-1]; + dE2 = dE1; + }else{ + dE1 = xPeakList[i] - xPeakList[i-1]; + dE2 = xPeakList[i+1] - xPeakList[i]; + } + + fit->SetParLimits(3*i+1, xPeakList[i] - dE1 , xPeakList[i] + dE2 ); + if( sigmaMax== 0 ) fit->SetParLimits(3*i+2, 0, 1.5*sigma[i]); // add 50% margin of sigma + if( sigmaMax < 0 ) fit->FixParameter(3*i+2, abs(sigmaMax)); + if( sigmaMax > 0 ) fit->SetParLimits(3*i+2, 0, sigmaMax); + } + + hspec->Fit("fit", "Rq"); + fit->Draw("same"); + + //=========== get the fit parameters + const Double_t* paraA = fit->GetParameters(); + const Double_t* paraE = fit->GetParErrors(); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + double bw = hspec->GetBinWidth(1); + + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + TF1 ** gFit = new TF1 *[nPeaks]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + cClickFitNGaussPolsub->cd(3); + PlotResidual(hspec, fit); + + cClickFitNGaussPolsub->cd(4); + + text.SetTextSize(0.05); + text.SetTextColor(2); + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + } +} + + +//######################################## +//#### fit N-Gauss with pol-1 BG at fixed for < 0 +//######################################## +void fitSpecial(TH1F * hist, TString fitFile = "AutoFit_para.txt"){ + + printf("================================================================\n"); + printf("================ Special fit for h074_82Kr =====================\n"); + printf(" * need the file input \n"); + printf("================================================================\n"); + + + bool isParaRead = loadFitParameters(fitFile); + if( !isParaRead ) { + printf("Please provide a valid input file\n"); + return; + } + + gStyle->SetOptStat(""); + gStyle->SetOptFit(0); + nPeaks = energy.size(); + nPols = 1; + + TCanvas * cFitSpecial = NewCanvas("cFitSpecial", "Fitting for h074_82Kr", 1, 3, 800, 300); + //if(! cFitSpecial->GetShowEventStatus() ) cFitSpecial->ToggleEventStatus(); + cFitSpecial->cd(1); + + ScaleAndDrawHist(hist, 0, 0); + + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + int xBin = hist->GetXaxis()->GetNbins(); + + printf("============= find the pol-1 background ..... \n"); + + TF1 * fitpol1 = new TF1("fitpol1", "pol1", xMin, -0.3); + fitpol1->SetParameter(0, gRandom->Rndm()); + fitpol1->SetParameter(1, gRandom->Rndm()); + + hist->Fit("fitpol1", "Rq"); + + double x0 = fitpol1->GetParameter(0); + double x1 = fitpol1->GetParameter(1); + + + int nPar = 3 * nPeaks + nPols + 1; + double * para = new double[nPar]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = energy[i]; + para[3*i+2] = sigma[i]/2.; + } + + para[3*nPeaks+0] = x0; + para[3*nPeaks+1] = x1; + + + TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + if( energyFlag[i] == 1 ) { + fit->FixParameter(3*i+1, energy[i]); + }else{ + fit->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fit->FixParameter(3*i+2, sigma[i]); + }else{ + fit->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + fit->FixParameter(3*nPeaks + 0, x0); + fit->FixParameter(3*nPeaks + 1, x1); + + hist->Fit("fit", "Rq"); + + //=========== get the polynomial part + const Double_t* paraA = fit->GetParameters(); + const Double_t* paraE = fit->GetParErrors(); + + TString polExp = "[0]"; + for( int i = 1; i < nPols + 1; i++){ + polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + } + + TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax); + for( int i = 0; i < nPols + 1; i++){ + bg->SetParameter(i, paraA[3*nPeaks+i]); + } + bg->SetNpx(1000); + + for( int i = 0; i < nPols + 1; i++){ + printf("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i]); + } + printf("====================================\n"); + + cFitSpecial->cd(1); + bg->Draw("same"); + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(1); + + for( int i = 0; i < nPols + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i])); + } + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + + //GoodnessofFit(specS, fit); + + double bw = hist->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + cFitSpecial->Update(); + + cFitSpecial->cd(2); + PlotResidual(hist, fit); + + cFitSpecial->cd(3); + + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + } + +} + + + + +#endif diff --git a/armory/DataBlock.h b/armory/DataBlock.h new file mode 100644 index 0000000..ecfe1b2 --- /dev/null +++ b/armory/DataBlock.h @@ -0,0 +1,110 @@ + +#ifndef DATABLOCK_H +#define DATABLOCK_H + +#include +#include +#include +#include +#include + +#include "TSystem.h" +#include "TObject.h" +#include "TFile.h" +#include "TTree.h" +#include "TString.h" + +#define MAX_TRACE_LENGHT 16000 + +class DataBlock{ + +public: + UShort_t ch; + UShort_t slot; + UShort_t crate; + UShort_t headerLength; /// headerLength > 4, more data except tarce. + UShort_t eventLength; /// eventLength = headerLength + trace + Bool_t pileup; + ULong64_t time; + Bool_t cfd_forced; + Bool_t cfd_source; + UShort_t cfd; + UShort_t energy; + UShort_t trace_length; + Bool_t trace_out_of_range; + + Int_t trailing; + Int_t leading; + Int_t gap; + Int_t baseline; + Int_t QDCsum[8]; + + ULong64_t eventID; + + UShort_t trace[MAX_TRACE_LENGHT]; + + DataBlock(){ + Clear(); + }; + + ~DataBlock(){}; + + void Clear(){ + ch = 0; + slot = 0; + crate = 0; + headerLength = 0; + eventLength = 0; + pileup = false; + time = 0; + cfd_forced = false; + cfd_source = false; + cfd = 0; + energy = 0; + trace_length = 0; + trace_out_of_range = 0; + eventID = 0; + ClearQDC(); + ClearTrace(); + } + + void ClearQDC(){ + trailing = 0; + leading = 0; + gap = 0; + baseline = 0; + for( int i = 0; i < 8; i++) QDCsum[i] = -1; + } + + void ClearTrace(){ + for( int i = 0 ; i < MAX_TRACE_LENGHT; i++) trace[i] = 0; + } + + + void Print(int opt = 0){ + printf("============== eventID : %llu\n", eventID); + printf("Crate: %d, Slot: %d, Ch: %d \n", crate, slot, ch); + printf("HeaderLength: %d, Event Length: %d, energy: %d, timeStamp: %llu\n", headerLength, eventLength, energy, time); + printf("trace_length: %d, pile-up:%d\n", trace_length, pileup); + printf("CFD : %d , Forced : %d , Source : %d\n", cfd, cfd_forced, cfd_source); + if( headerLength > 4 ){ + if( headerLength > 12 ){ + printf(" trailing : %d\n", trailing); + printf(" leading : %d\n", leading); + printf(" gap : %d\n", gap); + printf(" baseLine : %d\n", baseline); + } + if( opt >= 1 ){ + printf(" QDCsum : \n"); + for( int i = 0; i < 8; i++) printf(" %-10d\n", QDCsum[i]); + } + } + if( opt >= 2 && eventLength > headerLength ){ + printf(" trace:\n"); + for( int i = 0 ; i < trace_length ; i++)printf("%3d| %-10d\n",i, trace[i]); + } + } + +}; + +#endif diff --git a/armory/EventBuilder_evt b/armory/EventBuilder_evt new file mode 100755 index 0000000..5c51149 Binary files /dev/null and b/armory/EventBuilder_evt differ diff --git a/armory/EventBuilder_evt.cpp b/armory/EventBuilder_evt.cpp new file mode 100644 index 0000000..18310bf --- /dev/null +++ b/armory/EventBuilder_evt.cpp @@ -0,0 +1,270 @@ +/**********************************************************/ +/* */ +/* Modified by Ryan From */ +/* */ +/* PXI SCAN CODE -- J.M. Allmond (ORNL) -- July 2016 */ +/* */ +/**********************************************************/ + +#include +#include +#include +#include +#include + +#include "TFile.h" +#include "TTree.h" +#include "TMath.h" +#include "TBenchmark.h" + +#define RAND ((float) rand() / ((unsigned int) RAND_MAX + 1)) // random number in interval (0,1) + +#define MAX_CRATES 2 +#define MAX_BOARDS_PER_CRATE 13 +#define MAX_CHANNELS_PER_BOARD 16 +#define BOARD_START 2 + +#define MAX_ID MAX_CRATES*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + +#define HEADER_LENGTH 4 //unit = words with 4 bytes per word +#define MAX_SUB_LENGTH 2016 //unit = words with 4 bytes per word ; 2004 --> 40 micro-second trace + 4 word header + +#define RAWE_REBIN_FACTOR 2.0 // Rebin 32k pixie16 spectra to something smaller to fit better into 8k. + +#include "../mapping.h" + +#include "../armory/evtReader.h" + +unsigned long long int dataCount=0; +unsigned long long int pileUpCount=0; +unsigned long long int evtCount=0; +int tick2ns = 4; + +/////////////////////////////////// +// START OF MAIN FUNCTION // +/////////////////////////////////// +int main(int argc, char **argv) { + + printf("=====================================\n"); + printf("=== evt.to --> root ===\n"); + printf("=====================================\n"); + + // Check that the corrent number of arguments were provided. + if (argc <= 3) { + printf("Incorrect number of arguments:\n"); + printf("%s [outFile] [timeWindow] [to File1] [to File2] .... \n", argv[0]); + printf(" outFile : output root file name\n"); + printf(" timeWindow : number of tick, 1 tick = %d ns. default = 100 \n", tick2ns); + return 1; + } + + TString outFileName = argv[1]; + int timeWindow = atoi(argv[2]); + int nFile = argc - 3; + TString inFileName[nFile]; + for( int i = 0 ; i < nFile ; i++){ + inFileName[i] = argv[i+3]; + } + + printf("====================================\n"); + + evtReader * evt = new evtReader(); + DataBlock * data = evt->data; + + printf(" Number of input file : %d \n", nFile); + printf(" out file : \033[1;31m%s\033[m\n", outFileName.Data()); + printf(" Event building time window : %d tics = %d nsec \n", timeWindow, timeWindow*tick2ns); + + + TFile * outRootFile = new TFile(outFileName, "recreate"); + outRootFile->cd(); + TTree * tree = new TTree("tree", outFileName); + + unsigned long long evID = 0; + int multi = 0; + int id[MAX_ID] = {0}; + double e[MAX_ID] = {TMath::QuietNaN()}; + unsigned long long e_t[MAX_ID] = {0}; + int cfd[MAX_ID] = {0}; + bool pileup[MAX_ID] = {0}; + int qdc[MAX_ID][8] = {0}; + int multiClover = 0 ; /// this is total multiplicity for all crystal + int runID = 0; // date-run-fileID + int multiBeam = 0; + //unsigned short pileup[MAXMULTI]; + + tree->Branch("evID", &evID, "event_ID/l"); + tree->Branch("multi", &multi, "multi/I"); + tree->Branch("detID", id, "detID[multi]/I"); + tree->Branch("e", e, "e[multi]/D"); + tree->Branch("e_t", e_t, "e_timestamp[multi]/l"); + //tree->Branch("pileup", pileup, "pileup[multi]/O"); + tree->Branch("cfd", cfd, "cfd[multi]/I"); + tree->Branch("qdc", qdc, "qdc[multi][8]/I"); + tree->Branch("multiClover", &multiClover, "multiplicity_Clover/I"); + tree->Branch("multiBeam", &multiBeam, "multiplicity_Beam/I"); + tree->Branch("runID", &runID, "runID/I"); + + int countGP = 0; //gamma-particle coincident + double totalDataSize = 0; + int outFileCount = 0; + + for( int i = 0; i < nFile; i++){ + + evt->OpenFile(inFileName[i], false); + if( evt->IsOpen() == false ) continue; + + printf("==================================================== %d / %d\n", i+1, nFile); + printf("\033[1;31m%s \033[m\n", inFileName[i].Data()); + + int pos = inFileName[i].Last('/'); + TString temp = inFileName[i]; + temp.Remove(0, pos+1); + pos = temp.Last('-'); + temp.Remove(0, pos+1); + pos = temp.Last('.'); + temp.Remove(pos); + runID = atoi(temp.Data()); + + long long int etime = -1; + long long int tdif = -1; + + while ( evt->IsEndOfFile() == false ) { //main while loop + + if ( evt->ReadBlock() == -1) break; + + if ( data->crate == 0 ) continue; + + //Set reference time for event building + if (etime == -1) { + etime = data->time; + tdif = 0; + multi = 0; + multiClover = 0; + multiBeam = 0; + }else { + tdif = data->time - etime; + } + + //Check for end of event, rewind, and break out of while loop + if (tdif > timeWindow) { + + //Gate + //if( multiClover > 2 && multiBeam ==0 ) { + // + outRootFile->cd(); + tree->Fill(); + //printf("------------------\n"); + // + // countGP++; + //} + + evID ++; + + //clear data + etime = data->time; + tdif = 0; + multi = 0; + multiClover = 0; + multiBeam = 0; + + int haha = data->crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data->slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + data->ch; + id[multi] = mapping[haha]; + e[multi] = data->energy; + e_t[multi] = data->time; + cfd[multi] = data->cfd_forced == 0 ? (data->cfd - 16384*(data->cfd_source == true ? 1 : 0))/2 : 0; + pileup[multi] = data->pileup; + for( int i = 0; i < 8; i++) qdc[multi][i] = data->QDCsum[i]; + + if( 0 <= id[multi] && id[multi] < 100 ) multiClover += 1; + if( id[multi] == 300 || id[multi] == 301 || id[multi] == 307 || id[multi] == 308 ) multiBeam += 1; + //printf("id: %d, multiBeam: %d, %llu, tdif %llu\n", id[multi], multiBeam, e_t[multi], tdif); + multi++ ; + }else{ + //if within time window, fill array; + int haha = data->crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data->slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + data->ch; + id[multi] = mapping[haha]; + e[multi] = data->energy; + e_t[multi] = data->time; + cfd[multi] = data->cfd_forced == 0 ? (data->cfd - 16384*(data->cfd_source == true ? 1 : 0))/2 : 0; + //printf(" %d , %d , %d, %d \n", data->cfd_forced, data->cfd_source, data->cfd, cfd[multi]); + pileup[multi] = data->pileup; + for( int i = 0; i < 8; i++) qdc[multi][i] = data->QDCsum[i]; + if( 0 <= id[multi] && id[multi] < 100 ) multiClover += 1; + if( id[multi] == 300 || id[multi] == 301 || id[multi] == 307 || id[multi] == 308 ) multiBeam += 1; + //printf("id: %d, multiBeam: %d, %llu, tdif %llu\n", id[multi], multiBeam, e_t[multi], tdif); + multi++ ; + + } + + // total pileups + if (data->pileup == 1) { + pileUpCount++; + } + + evt->PrintStatus(10000); + + } // end main while loop + + evt->PrintStatus(1); + printf("\n\n\n"); + printf(" total number of event built : %llu\n", evID); + //printf(" total number of Gamma - GAGG coincdient : %d (%.3f %%)\n", countGP, countGP*1.0/evID); + + outRootFile->cd(); + tree->Write(); + + totalDataSize += (evt->GetFileSize())/1024./1024./1024.; + + double rootFileSize = outRootFile->GetSize()/1024./1024./1024. ; // in GB + printf(" ----------- root file size : %.3f GB\n", rootFileSize); + printf(" ---------- total read size : %.3f GB\n", totalDataSize); + printf(" ----------- reduction rate : %.3f %%\n", rootFileSize*100./totalDataSize); + + evt->CloseFile(); + + /* + if( rootFileSize > 3.0 ) { + break; + }*/ + + ///try to open a new root file when file size > 2 GB + /*if( rootFileSize > 2.0 ) { + + outRootFile->Close(); + delete outRootFile; + delete tree; + outFileCount += 1; + + if( outFileCount > 5 ) break; + + TString outFileName2 = outFileName; + outFileName2.Insert(outFileName.Sizeof() - 6, Form("_%03d",outFileCount)); + outRootFile = new TFile( outFileName2, "recreate"); + + tree = new TTree("tree", "tree"); + tree->Branch("evID", &evID, "event_ID/l"); + tree->Branch("multi", &multi, "multi/I"); + tree->Branch("detID", id, "detID[multi]/I"); + tree->Branch("e", e, "e[multi]/D"); + tree->Branch("e_t", e_t, "e_timestamp[multi]/l"); + tree->Branch("qdc", qdc, "qdc[multi][8]/I"); + tree->Branch("multiClover", &multiClover, "multiplicity_crystal/I"); + tree->Branch("multiBeam", &multiBeam, "multiplicity_GAGG/I"); + tree->Branch("runID", &runID, "runID/I"); + + }*/ + + } + + + outRootFile->Close(); + + printf("\n\n\n==================== finished.\r\n"); + + //printf(" number of Gamma - GAGG coincdient : %d\n", countGP); + + return 0; +} + + diff --git a/armory/EventBuilder_seperated.cpp b/armory/EventBuilder_seperated.cpp new file mode 100644 index 0000000..45fd3cf --- /dev/null +++ b/armory/EventBuilder_seperated.cpp @@ -0,0 +1,206 @@ +#include +#include +#include +#include +#include + +#include "TFile.h" +#include "TTree.h" +#include "TMath.h" +#include "TBenchmark.h" +#include "TStopwatch.h" +#include "TTreeIndex.h" + +#include "../mapping.h" + +Int_t eventID = 0 ; +double e[NCRYSTAL]; +ULong64_t e_t[NCRYSTAL]; +double bgo[NBGO]; +ULong64_t bgo_t[NBGO]; +Short_t other[NOTHER]; +Short_t multi; + +void ClearTreeData(){ + + for( int i = 0; i < NCRYSTAL; i++){ + e[i] = TMath::QuietNaN(); + e_t[i] = 0; + //pileup[i] = 0; + //hit[i] = 0; + } + for( int i = 0; i < NBGO; i++) { + bgo[i] = TMath::QuietNaN(); + bgo_t[i] = 0 ; + } + for( int i = 0; i < NOTHER; i++) { + other[i] = TMath::QuietNaN(); + } + multi = 0; +} + +int main(int argn, char **argv){ + printf("=====================================\n"); + printf("=== Event Builder ===\n"); + printf("=====================================\n"); + + if (argn != 2 && argn != 3 && argn != 4 ) { + printf("Usage :\n"); + printf("%s [_raw.root File] \n", argv[0]); + printf(" timeWindows : default = 100 \n"); + printf(" SaveFileName : default is *.root \n"); + return 1; + } + + TString inFileName = argv[1]; // need to check name + int timeWindow = 100; + if( argn >= 3 ) timeWindow = atoi(argv[2]); + + printf(">>> Opening input %s \n", inFileName.Data()); + TFile * inFile = new TFile(inFileName, "READ"); + if( inFile->IsOpen() == false ) { + printf("!!!! cannot open file %s \n", inFileName.Data()); + return 0; + } + + TTree * tree = (TTree *) inFile->Get("tree"); + + Long64_t evID; + UShort_t detID; + UShort_t energy; + ULong64_t energy_t; + + TBranch *b_data_ID; //! + TBranch *b_ID; //! + TBranch *b_energy; //! + TBranch *b_energy_timestamp; //! + + tree->SetBranchAddress("evID", &evID, &b_data_ID); + tree->SetBranchAddress("id", &detID, &b_ID); + tree->SetBranchAddress("e", &energy, &b_energy); + tree->SetBranchAddress("e_t", &energy_t, &b_energy_timestamp); + + Long64_t totnumEntry = tree->GetEntries(); + + printf( "total Entry : %lld \n", totnumEntry); + + printf(">>> Buidling Index using the timestamp\n"); + tree->BuildIndex("e_t"); + TTreeIndex *in = (TTreeIndex*) tree->GetTreeIndex(); + Long64_t * index = in->GetIndex(); + + ULong64_t time0; //time-0 for each event + int timeDiff; + + + TString outFileName = inFileName; + outFileName.Remove(inFileName.First("_raw")); + outFileName.Append(".root"); + if( argn >=4 ) outFileName = argv[3]; + + printf(">>> out File name : %s\n", outFileName.Data()); + + printf(">>> Create output tree\n"); + TFile * saveFile = new TFile(outFileName, "recreate"); + saveFile->cd(); + TTree * newtree = new TTree("tree", "tree"); + + newtree->Branch("evID", &eventID, "event_ID/l"); + newtree->Branch("e", e, Form("e[%d]/D", NCRYSTAL)); + newtree->Branch("e_t", e_t, Form("e_timestamp[%d]/l", NCRYSTAL)); + //newtree->Branch("p", pileup, Form("pile_up_flag[%d]/s", NCRYSTAL)); + //newtree->Branch("hit", hit, Form("hit[%d]/s", NCRYSTAL)); + + newtree->Branch("bgo", bgo, Form("BGO_e[%d]/D", NBGO)); + newtree->Branch("bgo_t", bgo_t, Form("BGO_timestamp[%d]/l", NBGO)); + + newtree->Branch("other", other, Form("other_e[%d]/D", NOTHER)); + + newtree->Branch("multi", &multi, "multiplicity_crystal/I"); + + ClearTreeData(); + + printf("================== Start processing....\n"); + Float_t Frac = 0.1; ///Progress bar + TStopwatch StpWatch; + StpWatch.Start(); + eventID = 0; + + for( Long64_t entry = 0; entry < totnumEntry; entry++){ + + + /*********** Progress Bar ******************************************/ + if (entry>totnumEntry*Frac-1) { + TString msg; msg.Form("%llu", totnumEntry/1000); + int len = msg.Sizeof(); + printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n", + Frac*100, len, entry/1000,totnumEntry/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac); + StpWatch.Start(kFALSE); + Frac+=0.1; + } + + entry = index[entry]; + + b_ID->GetEntry(entry); + b_energy->GetEntry(entry); + b_energy_timestamp->GetEntry(entry); + + if( time0 == 0) time0 = energy_t; + timeDiff = (int) (energy_t - time0); + + if( timeDiff < timeWindow ) { + + if ( detID < NCRYSTAL ){ + e[detID] = energy; + e_t[detID] = energy_t; + multi++; + } + if ( 100 <= detID && detID < 100 + NBGO ){ + bgo[detID-100] = energy; + bgo_t[detID-100] = energy_t; + } + if ( 200 <= detID && detID < 200 + NOTHER){ + other[detID-200] = energy; + } + + //printf("%d | %3d %6d %10llu, %3d\n", multi, detID, energy, energy_t, timeDiff); + + }else{ + //---- end of event + eventID ++; + + saveFile->cd(); + newtree->Fill(); + + ClearTreeData(); + + /// fill 1st data of an event + time0 = energy_t; + timeDiff = 0; + + if ( detID < NCRYSTAL ){ + e[detID] = energy; + e_t[detID] = energy_t; + multi = 1; + } + if ( 100 <= detID && detID < 100 + NBGO ){ + bgo[detID-100] = energy; + bgo_t[detID-100] = energy_t; + } + if ( 200 <= detID && detID < 200 + NOTHER){ + other[detID-200] = energy; + } + + } + + } + + printf("============================== finished.\n"); + + saveFile->cd(); + newtree->Write(); + saveFile->Close(); + + printf(" total number of event Built : %d \n", eventID); + +} diff --git a/armory/LoadCorrection.h b/armory/LoadCorrection.h new file mode 100644 index 0000000..a17b609 --- /dev/null +++ b/armory/LoadCorrection.h @@ -0,0 +1,111 @@ +#ifndef Analysis_Library_h +#define Analysis_Library_h + +#include +#include +#include +#include +#include +#include +#include +#include + +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; +} + +std::vector> LoadCorrectionParameters(TString corrFile, bool show=false){ + + printf(" load correction parameters : %s", corrFile.Data()); + std::ifstream file; + file.open(corrFile.Data()); + + std::vector> corr; + corr.clear(); + + std::vector detCorr; + detCorr.clear(); + + if( file.is_open() ){ + + while( file.good() ){ + + std::string line; + getline(file, line); + + if( line.substr(0,1) == "#" ) continue; + if( line.substr(0,2) == "//" ) continue; + if( line.size() == 0 ) continue; + + //printf("%s\n", line.c_str()); + std::vector temp = SplitStr(line, " "); + + detCorr.clear(); + for( int i = 0; i < (int) temp.size() ; i++){ + detCorr.push_back(std::stod(temp[i])); + } + corr.push_back(detCorr); + } + + file.close(); + + printf(".... done\n"); + if( show ){ + printf("===== correction parameters \n"); + for( int i = 0; i < (int) corr.size(); i++){ + printf("det : %2d | ", i ); + int len = (int) corr[i].size(); + for( int j = 0; j < len - 1 ; j++){ + printf("%14.6f, ", corr[i][j]); + } + printf("%14.6f\n", corr[i][len-1]); + } + } + + }else{ + printf(".... fail\n"); + std::vector temp = {0, 1}; + for( int i = 0; i < 36; i++){ + corr.push_back(temp); + } + } + + return corr; +} + +#endif + diff --git a/armory/MergeEVT b/armory/MergeEVT new file mode 100755 index 0000000..81c1e10 Binary files /dev/null and b/armory/MergeEVT differ diff --git a/armory/MergeEVT.cpp b/armory/MergeEVT.cpp new file mode 100644 index 0000000..10773b7 --- /dev/null +++ b/armory/MergeEVT.cpp @@ -0,0 +1,135 @@ +#include +#include +#include +#include +#include +#include "TFile.h" +#include "TTree.h" +#include "TString.h" +#include "TMath.h" +#include "TBenchmark.h" +#include + +#define MAX_CRATES 2 +#define MAX_BOARDS_PER_CRATE 13 +#define MAX_CHANNELS_PER_BOARD 16 +#define BOARD_START 2 + +#include "../mapping.h" +#include "../armory/DataBlock.h" +#include "../armory/evtReader.h" + +//############################################# +// main +//############################################# +int main(int argn, char **argv) { + + printf("=====================================\n"); + printf("=== evt --> _raw.root ===\n"); + printf("=====================================\n"); + + if (argn < 3 ) { + printf("Usage :\n"); + printf("%s [outFile] [evt1] [evt2] [evt3] ..... \n", argv[0]); + printf("e.g.: \n"); + printf("%s hahaha_raw.root haha-000.evt haha-001.evt haha-002.evt\n", argv[0]); + printf("%s hahaha_raw.root `ls haha-*.evt`\n", argv[0]); + return 1; + } + + TString outFileName = argv[1]; + int nFiles = argn-2; + TString inFileName[nFiles]; + for( int i = 0; i < nFiles ; i++){ + inFileName[i] = argv[i+2]; + printf(" in file - %2d: %s\n", i, inFileName[i].Data()); + } + + printf(" out file: %s\n", outFileName.Data()); + + + evtReader * evt = new evtReader(); + DataBlock * data = evt->data; + short detID; + printf("====================================\n"); + + //====== ROOT file + TFile * outFile = new TFile(outFileName, "recreate"); + TTree * tree = new TTree("tree", "tree"); + + tree->Branch("evID", &data->eventID, "data_ID/L"); + tree->Branch("detID", &detID, "detID/s"); + tree->Branch("e", &data->energy, "crystal_energy/s"); + tree->Branch("cfd", &data->cfd, "cfd/s"); + tree->Branch("e_t", &data->time, "crystal_timestamp/l"); + tree->Branch("p", &data->pileup, "pileup/O"); + // tree->Branch("trace_length", &data->trace_length, "trace_length/s"); + // tree->Branch("trace", data->trace, "trace[trace_length]/s"); + + TBenchmark gClock; + gClock.Reset(); + gClock.Start("timer"); + + //========================================= + //========================================= + //========================================= + //========================================= + for( int i = 0; i < nFiles; i++){ + + evt->OpenFile(inFileName[i]); + if( evt->IsOpen() == false ) continue; + + Long64_t measureCount = 0; + printf("\033[1;31mProcessing file: %s\033[0m\n", inFileName[i].Data()); + TBenchmark clock2; + clock2.Reset(); + clock2.Start("timer"); + + evt->ScanNumberOfBlock(); + Long64_t nBlock = evt->GetNumberOfBlock(); + + //=============== Read File + while( evt->IsEndOfFile() == false ){ + + evt->ReadBlock(); + evt->PrintStatus(nBlock/10); + + int id = data->crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data->slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + data->ch; + if( id > 352 || id < 0 ) continue; + detID = mapping[id]; + + //cern fill tree + outFile->cd(); + tree->Fill(); + } + + clock2.Stop("timer"); + double time = clock2.GetRealTime("timer"); + float tempf = (float)evt->GetFilePos()/(1024.*1024.*1024.); + printf(" measurements: \x1B[32m%lld \x1B[0m | %.3f GB\n", evt->GetBlockID(), tempf); + printf(" Time used:%3.0f min %5.2f sec\n", TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + printf(" Root file size so far: %.4f GB\n", outFile->GetSize()/1024./1024./1024.); + + } + + gClock.Stop("timer"); + double time = gClock.GetRealTime("timer"); + gClock.Start("timer"); + float tempf = (float)evt->GetFilePos()/(1024.*1024.*1024.); + printf("Total measurements: \x1B[32m%lld \x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[A\r", + evt->GetBlockID()+1, (100*evt->GetFilePos()/evt->GetFileSize()), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + + + //cern save root + outFile->cd(); + double totRootSize = outFile->GetSize()/1024./1024./1024.; + tree->Write(); + outFile->Close(); + + gClock.Stop("timer"); + time = gClock.GetRealTime("timer"); + printf("\n==================== finished.\r\n"); + printf("Total time spend : %3.0f min %5.2f sec\n", TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + printf(" File size of %s : %.3f GB \n", outFileName.Data(), totRootSize); + +} diff --git a/armory/ev22txt b/armory/ev22txt new file mode 100755 index 0000000..47d8916 Binary files /dev/null and b/armory/ev22txt differ diff --git a/armory/ev22txt.cpp b/armory/ev22txt.cpp new file mode 100644 index 0000000..7c44366 --- /dev/null +++ b/armory/ev22txt.cpp @@ -0,0 +1,81 @@ +#include +#include +#include +#include + +using namespace std; + +int main(int argc, char **argv){ + + printf("=====================================\n"); + printf("=== ev2 --> txt ===\n"); + printf("=====================================\n"); + + if ( argc != 2 && argc != 3 ){ + printf("Usage: \n"); + printf("%10s [ev2 file] [max number of event]\n", argv[0]); + } + + string inFileName = argv[1]; + string outFileName = argv[1]; + int found = outFileName.find('.'); + int len = outFileName.length(); + outFileName.erase(found, len-found); + outFileName.append(".txt"); + + long long maxNumEvent = 0; // if maxNumEvent <= 0, all event + if( argc >= 3 ) maxNumEvent = atoi(argv[2]); + + + //open list-mode data file from PXI digitizer + FILE *inFile = fopen(argv[1], "r"); + long int inFileSize,inFilepos; + if ( inFile == NULL) { + printf("Error, cannot open input file %s\n", argv[1]); + return 1; + } + + //get file size + fseek(inFile, 0L, SEEK_END); + inFileSize = ftell(inFile); + rewind(inFile); + long int inFilePos = 0; + + printf(" in file : %s | file size : %f MB\n", inFileName.c_str(), inFileSize/1024./1024.); + //printf(" out file : %s \n", outFileName.c_str()); + if( maxNumEvent <= 0 ){ + printf(" max number of event : ALL \n"); + }else{ + printf(" max number of event : %lld \n", maxNumEvent); + } + printf("-------------------------------------\n"); + + long long nEvent = 0; + short * multi = new short[1]; + short * id = new short[1]; + short * energy = new short[1] ; + + while ( inFilePos < inFileSize || feof(inFile) ) { + + printf("------------------- %lld\n", nEvent); + + fread(multi, 1, 1, inFile); + printf(" number of det : %d \n", multi[0]); + + for( int i = 0; i < multi[0] ; i++){ + + fread(id, 1, 1, inFile); + fread(energy, 2, 1, inFile); + + printf("%2d| %3d, %8d \n", i, id[0], energy[0]); + + } + fread(multi, 1, 1, inFile); + + nEvent++; + + if( maxNumEvent > 0 && nEvent > maxNumEvent ) break; + + } + +} diff --git a/armory/evt2hist b/armory/evt2hist new file mode 100755 index 0000000..51bd0ed Binary files /dev/null and b/armory/evt2hist differ diff --git a/armory/evt2hist.cpp b/armory/evt2hist.cpp new file mode 100644 index 0000000..380fbe3 --- /dev/null +++ b/armory/evt2hist.cpp @@ -0,0 +1,344 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "TSystem.h" +#include "TObject.h" +#include "TFile.h" +#include "TTree.h" +#include "TString.h" +#include "TMath.h" +#include "TGraph.h" +#include "TLatex.h" +#include "TBenchmark.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TApplication.h" +#include "TCanvas.h" +#include "TClonesArray.h" + +#include "../mapping.h" +#include "../armory/AnalysisLibrary.h" + +#include "../armory/DataBlock.h" +#include "../armory/evtReader.h" + +#define sleepTime 5000 ///sleep for 5 sec +//#############################TODO +// 1) multiple file +// 2) Change to GUI +// 4) eventBuilding +// 3) last 10% data + + +//################################## user setting +int eRange[2] = {50, 1000}; ///min, max + +bool PIDFlag = false; +int GAGGID = 209; +//############################################# +// main +//############################################### +int main(int argn, char **argv) { + + if (argn < 2 || argn > 7 ) { + printf("Usage :\n"); + printf("%s [evt File] [E corr] [raw E threshold] [Hist File] [Root File] [display data]\n", argv[0]); + printf(" [E corr] : correction file for gamma energy \n"); + printf(" [raw E threshold] : min raw E, default 10 ch \n"); + printf(" [Hist File] : if provided, saved histograms in root. if value = 1, auto fileName. \n"); + printf(" [Root File] : if provided, saved energy, timestamp, QDC, and trace in to root. if value = 1, auto fileName. \n"); + printf(" [display data] : number of event from the first one to display. default 0. \n"); + return 1; + } + + TString inFileName = argv[1]; + + TString corrFile = ""; + std::vector> eCorr; + if( argn >= 3 ){ + corrFile = argv[2]; + eCorr.clear(); + eCorr = LoadCorrectionParameters(corrFile); + } + + int rawEnergyThreshold = 10; + if( argn >= 4 ) rawEnergyThreshold = atoi(argv[3]); + + TString histFileName = ""; ///save gamma hist for calibration + if( argn >= 5 ) histFileName = argv[4]; + if( histFileName == "1" ){ + histFileName = inFileName; + histFileName.Remove(0, inFileName.Last('/')+1); + histFileName.Remove(histFileName.First('.')); + histFileName.Append("_hist.root"); + } + if( histFileName == "0" ) histFileName = ""; + + TString rootFileName = ""; ///save data into root + if( argn >= 6 ) rootFileName = argv[5]; + if( rootFileName == "1" ){ + rootFileName = inFileName; + rootFileName.Remove(0, inFileName.Last('/')+1); + rootFileName.Remove(rootFileName.First('.')); + rootFileName.Append("_raw.root"); + } + if( rootFileName == "0" ) rootFileName = ""; + + int maxDataDisplay = 0; + if( argn >= 7 ) maxDataDisplay = atoi(argv[6]); + + TBenchmark gClock; + gClock.Reset(); + gClock.Start("timer"); + + printf("====================================\n"); + + evtReader * evt = new evtReader(); + evt->OpenFile(inFileName); + if( evt->IsOpen() == false ) return -404; + DataBlock * data = evt->data; + + printf(" in file: \033[1;31m%s\033[m\n", inFileName.Data()); + printf(" Gamma energy correction file : %s\n", corrFile == "" ? "Not provided." : corrFile.Data()); + printf(" raw E threshold : %d ch\n", rawEnergyThreshold); + if( histFileName != "" ) printf(" Save histograms to %s\n", histFileName.Data()); + if( rootFileName != "" ) printf(" Save root to %s\n", rootFileName.Data()); + printf("--------------------------------\n"); + printf("Scanning the evt file... \n"); + evt->ScanNumberOfBlock(); + printf("go to 90%% of data \n"); + evt->JumptoPrecent(9); + + //================ ROOT tree + TFile * fFile = NULL; + TTree * tree = NULL; + + short detID; + + if( rootFileName != "" ){ + fFile = new TFile(rootFileName, "RECREATE"); + tree = new TTree("tree", "tree"); + + tree->Branch("headerLenght", &data->headerLength, "HeaderLength/s"); + tree->Branch("detID", &detID, "detID/s"); + tree->Branch("e", &data->energy, "energy/s"); + tree->Branch("e_t", &data->time, "timestamp/l"); + tree->Branch("p", &data->pileup, "pileup/O"); + tree->Branch("qdc", data->QDCsum, "QDC_sum[8]/I"); + tree->Branch("trace_length", &data->trace_length, "trace_length/s"); + tree->Branch("trace", data->trace, "trace[trace_length]/s"); + } + + //================ Historgrams + TH1F * he[NCRYSTAL]; + for( int i = 0 ; i < NCRYSTAL; i++){ + he[i] = new TH1F(Form("he%02d", i), Form("e-%2d", i), eRange[1]-eRange[0], eRange[0], eRange[1]); + switch (i % 4){ + case 0 : he[i]->SetLineColor(2); break; + case 1 : he[i]->SetLineColor(4); break; + case 2 : he[i]->SetLineColor(1); break; + case 3 : he[i]->SetLineColor(kGreen+3); break; + } + } + + + TH2F * hPID ; + if( PIDFlag ) hPID = new TH2F(Form("hPID%d", GAGGID), Form("GAGG - %d; tail; peak", GAGGID), 400, -10, 600, 400, -50, 1000); + + TGraph * gTrace = new TGraph(); + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(2); + + //================ Set Canvas + TApplication * app = new TApplication ("app", &argn, argv); + + TCanvas * canvas = new TCanvas("fCanvas", "Online Spectrum", 1800, 2000); + + canvas->Divide(3, TMath::Ceil(NCLOVER/3.), 0); + canvas->SetCrosshair(1); + for( int i = 0; i < 9 ; i++){ + canvas->cd(i+1)->SetBottomMargin(0.1); + canvas->cd(i+1)->SetRightMargin(0.002); + } + + ///TCanvas * cTrace = new TCanvas("cTrace", "Trace", 100, 100, 1000, 500); + TCanvas * cPID; + if( PIDFlag ) cPID = new TCanvas("cPID", "PID", 100, 100, 500, 500); + + //=============== Read File + int sleepCount = 0; + + while ( true ){ + + if( evt->ReadBlock()== -1 ) { + break; + //printf("\n\n\nReached the end of file, wait %d sec to see any update.\n", sleepTime); + //sleep( sleepTime ); + //evt->UpdateFileSize(); + //sleepCount ++; + //if( sleepCount > 1 ) { + // printf("waited for %d sec. exit.\n", 4* sleepTime); + // break; + //}else{ + // continue; + //} + } + sleepCount = 0; + + if( evt->GetBlockID() < maxDataDisplay ) { + printf("----------------------event Length: %u, fpos: %lu byte (%ld words)\n", data->eventLength, evt->GetFilePos(), evt->GetFilePos()/4); + data->Print(); + } + + //==== Fill Histogram + + int haha = data->crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data->slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + data->ch; + detID = mapping[haha]; + + if( 0 <= detID && detID < 100 && data->energy > rawEnergyThreshold ){ + if( corrFile != ""){ + ///========= apply correction + int order = (int) eCorr[detID].size(); + double eCal = 0; + for( int k = 0; k < order ; k++){ + eCal += eCorr[detID][k] * TMath::Power(data->energy, k); + } + he[detID]->Fill(eCal); + }else{ + he[detID]->Fill(data->energy); + } + } + + ///============ QDC + if( PIDFlag && detID == GAGGID && (data->headerLength < data->eventLength) ){ + double bg = (data->QDCsum[0] + data->QDCsum[1])/60.; + double peak = data->QDCsum[3]/20. - bg; + double tail = data->QDCsum[5]/55. - bg; + hPID->Fill( tail , peak); + } + + //===== Trace + ///if( data->trace_length > 0 ) { + /// gTrace->Clear(); + /// gTrace->Set(data->trace_length); + /// gTrace->SetTitle(Form("eventID : %llu, detID: %d", data->eventID, data->detID)); + /// + /// for( int i = 0; i < data->trace_length; i++) gTrace->SetPoint(i, i, data->trace[i]); + ///} + + if( rootFileName != "" ){ + fFile->cd(); + tree->Fill(); + } + + //==== event stats, print status every 10000 events + evt->PrintStatus(10000); + + //==== Plot Canvas + gClock.Stop("timer"); + int time = TMath::Floor(gClock.GetRealTime("timer")*1000); // in millisec + gClock.Start("timer"); + if( time % 1000 == 0 || time < 10){ + + ///==== for clover + for( int i = 0; i < NCLOVER; i++){ + double maxY = 0; + double y = 0; + for( int j = 0; j < 4; j++){ + int mBin = he[4*i+j]->GetMaximumBin(); + y = he[4*i+j]->GetBinContent(mBin); + if( maxY < y ) maxY = y; + } + for( int j = 0; j < 4; j++){ + canvas->cd(i+1); + he[4*i+j]->GetYaxis()->SetRangeUser(0, maxY*1.2); + if ( j == 0) { + he[4*i]->Draw(); + }else{ + he[4*i+j]->Draw("same"); + } + } + } + canvas->Modified(); + canvas->Update(); + + ///==== for trace + ///if( data->trace_length > 0 ){ + /// cTrace->cd(); + /// gTrace->Draw("AL"); + /// + /// for( int i = 0; i < 8; i++){ + /// text.DrawLatex(0.2, 0.8-0.05*i, Form("%d", data->QDCsum[i])); + /// } + /// cTrace->Modified(); + /// cTrace->Update(); + ///} + + ///=== for GAGG PID + if( PIDFlag ) { + cPID->cd(); + cPID->SetLogz(); + hPID->Draw("colz"); + cPID->Modified(); + cPID->Update(); + } + + gSystem->ProcessEvents(); + } + }//---- end of file loop + + + for( int i = 0; i < NCLOVER; i++){ + double maxY = 0; + double y = 0; + for( int j = 0; j < 4; j++){ + int mBin = he[4*i+j]->GetMaximumBin(); + y = he[4*i+j]->GetBinContent(mBin); + if( maxY < y ) maxY = y; + } + for( int j = 0; j < 4; j++){ + canvas->cd(i+1); + he[4*i+j]->GetYaxis()->SetRangeUser(0, maxY*1.2); + if ( j == 0) { + he[4*i]->Draw(); + }else{ + he[4*i+j]->Draw("same"); + } + } + } + canvas->Modified(); + canvas->Update(); + + gSystem->ProcessEvents(); + + evt->PrintStatus(1); + + printf("\n\n\n============= reached end of file\n"); + + if( histFileName != "" ) { + printf(" save gamma histograms : \033[1;3m%s\033[m\n", histFileName.Data()); + TFile * fHist = new TFile(histFileName, "RECREATE"); + for( int i = 0; i < NCRYSTAL; i++) he[i]->Write("", TObject::kOverwrite); + fHist->Close(); + } + + if( rootFileName != "" ){ + printf(" save into Root : \033[1;3m%s\033[m\n", rootFileName.Data()); + fFile->cd(); + tree->Write(); + fFile->Close(); + } + + printf("Crtl+C to end program.\n"); + + app->Run(); + +} diff --git a/armory/evtReader.h b/armory/evtReader.h new file mode 100644 index 0000000..8da3310 --- /dev/null +++ b/armory/evtReader.h @@ -0,0 +1,373 @@ +#ifndef EVTREADER_H +#define EVTREADER_H + +#include /// for FILE +#include +#include +#include +#include + +#include "TString.h" +#include "TBenchmark.h" + +#include "../armory/DataBlock.h" + +#define MAX_CRATES 2 +#define MAX_BOARDS_PER_CRATE 13 +#define MAX_CHANNELS_PER_BOARD 16 +#define BOARD_START 2 + +//TODO load the file into RAM, and read from the RAM + +class evtReader{ + + public: + DataBlock * data; + + private: + FILE * inFile; + + long int inFileSize; + long int inFilePos; + bool endOfFile; + bool isOpened; + Long64_t blockID; + long int nBlock; + + unsigned int extraHeader[14]; + unsigned int traceBlock[MAX_TRACE_LENGHT/2]; + + TBenchmark gClock; + + long int inFilePosPrecent[10]; + Long64_t blockIDPrecent[10]; + + bool fromRAM; + int * pxidata; + long nWords; + + ///============================================ Methods + public: + + evtReader(); + evtReader(TString inFileName, bool load2RAM); + ~evtReader(); + + void OpenFile(TString inFileName, bool load2RAM); + void CloseFile(); + + void UpdateFileSize(); + bool IsEndOfFile(); + + bool IsOpen() {return isOpened;} + long int GetFilePos() {return inFilePos;} + long int GetFileSize() {return inFileSize;} + Long64_t GetBlockID() {return blockID;} + Long64_t GetNumberOfBlock() {return nBlock;} + + + int ReadBlock(int opt = 0); /// 0 = default, fill data + /// 1 = no fill data + /// 2 = fill data and print + + void ScanNumberOfBlock(); + void JumptoPrecent(int precent); ///this is offset by 1 block + void PrintStatus(int mod); + +}; + + +//========================== implementation + +evtReader::evtReader(){ + inFile = 0; + data = new DataBlock(); + + inFileSize = 0; + inFilePos = 0; + + nBlock = 0; + blockID = -1; + endOfFile = false; + isOpened = false; + + fromRAM = false; + pxidata = NULL; + nWords = 0; +} + + +evtReader::~evtReader(){ + fclose(inFile); + delete inFile; + delete data; + delete pxidata; +} + + +evtReader::evtReader(TString inFileName, bool load2RAM = false){ + inFile = 0; + data = new DataBlock(); + + inFileSize = 0; + inFilePos = 0; + + nBlock = 0; + blockID = -1; + endOfFile = false; + isOpened = false; + fromRAM = false; // true until loaded to RAM + pxidata = NULL; + nWords = 0; + + OpenFile(inFileName, load2RAM); +} + +void evtReader::OpenFile(TString inFileName, bool load2RAM = false){ + inFile = fopen(inFileName, "r"); + if( inFile == NULL ){ + printf("Cannot read file : %s \n", inFileName.Data()); + }else{ + fseek(inFile, 0L, SEEK_END); + inFileSize = ftell(inFile); + rewind(inFile); ///back to the File begining + + //TODO load the evt file to RAM + if( load2RAM ) { + pxidata = (int *) malloc(inFileSize); + if( pxidata == NULL ){ + printf("\nError, memory not allocated.\n"); + }else{ + printf("Allocated %.3f GB of memory to buffer native PXI data\n", (float)inFileSize/(1024.0*1024.0*1024.0)); + printf("Now loading data to RAM\n"); + fread(pxidata, sizeof(pxidata), inFileSize/sizeof(pxidata), inFile); + fromRAM = true; + fclose(inFile); + inFilePos = 0; + } + } + + data->Clear(); + + gClock.Reset(); + gClock.Start("timer"); + + isOpened = true; + } +}; + +void evtReader::CloseFile(){ + fclose(inFile); + isOpened = false; + data->Clear(); + inFileSize = 0; + inFilePos = 0; + nBlock = 0; + blockID = -1; + endOfFile = false; +}; + +void evtReader::UpdateFileSize(){ + if( inFile == NULL ) return; + fseek(inFile, 0L, SEEK_END); + inFileSize = ftell(inFile); + fseek(inFile, inFilePos, SEEK_SET); +} + +bool evtReader::IsEndOfFile() { + int haha = feof(inFile); + return haha > 0 ? true: false; +} + + +int evtReader::ReadBlock(int opt){ + + if( feof(inFile) ) return -1; + if( endOfFile ) return -1; + + unsigned int header[4]; ///read 4 header, unsigned int = 4 byte = 32 bits. + + if( fromRAM ){ + for( int i = 0; i < 4 ; i++){ + header[i] = pxidata[nWords]; nWords += 1; + } + }else{ + if ( fread(header, sizeof(header), 1, inFile) != 1 ) { + endOfFile = true; + return -1; + } + } + blockID ++; + + + if( opt == 0 || opt == 2){ + /// see the Pixie-16 user manual, Table4-2 + data->eventID = blockID; + data->ch = header[0] & 0xF ; + data->slot = (header[0] >> 4) & 0xF; + data->crate = (header[0] >> 8) & 0xF; + data->headerLength = (header[0] >> 12) & 0x1F; + data->eventLength = (header[0] >> 17) & 0x3FFF; + data->pileup = header[0] >> 31 ; + data->time = ((ULong64_t)(header[2] & 0xFFFF) << 32) + header[1]; + data->cfd_forced = header[2] >> 16 & 0x8000; + data->cfd_source = header[2] >> 16 & 0x4000; + data->cfd = header[2] >> 16 & 0x3FFF; // 0x3FFF for 250MHz , 0x7FFF for 100MHz + data->energy = (header[3] & 0xFFFF ); + data->trace_length = (header[3] >> 16) & 0x7FFF; + data->trace_out_of_range = header[3] >> 31; + + data->ClearQDC(); + data->ClearTrace(); + + ///======== read QDCsum + if( data->headerLength >= 4 ){ + if( fromRAM ){ + for( int i = 0; i < data->headerLength - 4 ; i++){ + extraHeader[i] = pxidata[nWords]; nWords += 1; + } + }else{ + fread(extraHeader, sizeof(unsigned int) * (data->headerLength-4), 1, inFile); + } + if( data->headerLength == 8 || data->headerLength == 16){ + data->trailing = extraHeader[0]; + data->leading = extraHeader[1]; + data->gap = extraHeader[2]; + data->baseline = extraHeader[3]; + } + if( data->headerLength == 12 || data->headerLength == 16){ + for( int i = 0; i < 8; i++){ + int startID = 0; + if( data->headerLength > 12) startID = 4; ///the 1st 4 words + data->QDCsum[i] = extraHeader[i+startID]; + } + } + }else{ + for( int i = 0 ; i < 8; i++){ data->QDCsum[i] = 0;} + data->trailing = 0; + data->leading = 0; + data->gap = 0; + data->baseline = 0; + } + ///====== read trace + if( data->eventLength > data->headerLength ){ + if( fromRAM ){ + for( int i = 0; i < data->trace_length / 2 ; i++){ + traceBlock[i] = pxidata[nWords]; nWords += 1; + } + }else{ + fread(traceBlock, sizeof(unsigned int) * ( data->trace_length / 2 ), 1, inFile); + } + for( int i = 0; i < data->trace_length/2 ; i++){ + data->trace[2*i+0] = traceBlock[i] & 0xFFFF ; + data->trace[2*i+1] = (traceBlock[i] >> 16 ) & 0xFFFF ; + } + + ///make QDC by trace + /** + if( data->headerLength == 4 || data->headerLength == 8 ) { + for( int i = 0; i < 8; i++){ data->QDCsum[i] = 0;} + for( int i = 0; i < data->trace_length; i++){ + if( 0 <= i && i < 31 ) data->QDCsum[0] += data->trace[i]; + if( 31 <= i && i < 60 ) data->QDCsum[1] += data->trace[i]; + if( 60 <= i && i < 75 ) data->QDCsum[2] += data->trace[i]; + if( 75 <= i && i < 95 ) data->QDCsum[3] += data->trace[i]; + if( 95 <= i && i < 105 ) data->QDCsum[4] += data->trace[i]; + if( 105 <= i && i < 160 ) data->QDCsum[5] += data->trace[i]; + if( 160 <= i && i < 175 ) data->QDCsum[6] += data->trace[i]; + if( 175 <= i && i < 200 ) data->QDCsum[7] += data->trace[i]; + } + }*/ + } + } + + if( opt == 1 ){ + + data->headerLength = (header[0] >> 12) & 0x1F; + data->eventLength = (header[0] >> 17) & 0x3FFF; + data->trace_length = (header[3] >> 16) & 0x7FFF; + + if( data->headerLength >= 4 ){ + if( fromRAM ){ + nWords += (data->headerLength-4); + }else{ + fread(extraHeader, sizeof(unsigned int) * (data->headerLength-4), 1, inFile); + } + } + if( data->eventLength > data->headerLength ){ + if( fromRAM ){ + nWords += ( data->trace_length / 2 ); + }else{ + fread(traceBlock, sizeof(unsigned int) * ( data->trace_length / 2 ), 1, inFile); + } + } + } + + if ( fromRAM ){ + inFilePos = nWords * 32; + }else{ + inFilePos = ftell(inFile); + } + + if( opt == 2) data->Print(); + + return 1; +} + + +void evtReader::ScanNumberOfBlock(){ + + nBlock = 0; + int count = 0; + while( ReadBlock(1) != -1 ){ + nBlock ++; + int haha = (inFilePos*10/inFileSize)%10; + if( haha == count ) { + inFilePosPrecent[count] = inFilePos; + blockIDPrecent[count] = blockID; + count++; + } + + PrintStatus(10000); + } + + printf("\n\n\n"); + printf("scan complete: number of data Block : %ld\n", nBlock); + + inFilePos = 0; + blockID = -1; + + rewind(inFile); ///back to the File begining + endOfFile = false; + +} + +void evtReader::JumptoPrecent(int precent){ + + if( precent < 0 || precent > 10 ) { + printf("input precent should be 0 to 10\n"); + return; + } + + fseek(inFile, inFilePosPrecent[precent], SEEK_SET); + blockID = blockIDPrecent[precent]; + +} + +void evtReader::PrintStatus(int mod){ + + ///==== event stats, print status every 10000 events + + if ( blockID % mod == 0 ) { + + UpdateFileSize(); + gClock.Stop("timer"); + double time = gClock.GetRealTime("timer"); + gClock.Start("timer"); + printf("Total measurements: \x1B[32m%llu \x1B[0m\nReading Pos: \x1B[32m %.3f/%.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[A\033[A\r", + blockID, inFilePos/(1024.*1024.*1024.), inFileSize/1024./1024./1024, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + } + +} + +#endif diff --git a/armory/makefile b/armory/makefile new file mode 100644 index 0000000..29a7fe2 --- /dev/null +++ b/armory/makefile @@ -0,0 +1,25 @@ +CC=g++ + +all: EventBuilder_evt evt2hist MergeEVT ev22txt EventBuilder + +#this is for eventbuild +EventBuilder_evt: ../armory/EventBuilder_evt.cpp ../armory/DataBlock.h ../armory/evtReader.h ../mapping.h + $(CC) ../armory/EventBuilder_evt.cpp -o EventBuilder_evt `root-config --cflags --glibs` + +#this is for online root +MergeEVT: ../armory/MergeEVT.cpp ../armory/DataBlock.h ../armory/evtReader.h ../mapping.h + $(CC) ../armory/MergeEVT.cpp -o MergeEVT `root-config --cflags --glibs` + +#this is for online spectrums +evt2hist: ../armory/evt2hist.cpp ../armory/DataBlock.h ../armory/evtReader.h ../mapping.h + $(CC) ../armory/evt2hist.cpp -o evt2hist `root-config --cflags --glibs` + + +ev22txt: ../armory/ev22txt.cpp + $(CC) ../armory/ev22txt.cpp -o ev22txt + +EventBuilder: ../armory/EventBuilder.cpp + $(CC) ../armory/EventBuilder.cpp -o EventBuilder `root-config --cflags --glibs` + +clean: + -rm xia2root to2root MergeEVT evt2hist pxi-time-order ev22txt EventBuilder test diff --git a/convertPixie.sh b/convertPixie.sh new file mode 100755 index 0000000..b4cd768 --- /dev/null +++ b/convertPixie.sh @@ -0,0 +1,87 @@ + +nCore=10; + +if [ $# -eq 0 ]; then + echo "Usage: " + echo " ./convertPixie.sh [RUNNUM]" + exit; +fi + +runNum=$1 + +nFile=$(ls -1 rawdata/run${runNum}/run*evt | wc -l) +files=$(ls -1 rawdata/run${runNum}/run*evt) + +nLoop=$((${nFile}/nCore+1)) + +echo "===== number of files in run-${runNum} = ${nFile}" +echo " number of loops = ${nLoop}" +echo "=================================================" + +i=0 +while [[ $i -lt ${nLoop} ]]; do + + j=0 + while [[ $j -lt ${nCore} ]]; do + SegNum=$((${nCore}*i+j)) + + if [[ ${SegNum} -ge ${nFile} ]]; then break; fi + + #patching the prefix ZERO + if [ $runNum -lt 10 ]; then + runPrefix="000" + elif [ $runNum -lt 100 ]; then + runPrefix="00" + elif [ $runNum -lt 1000 ]; then + runPrefix="0" + fi + + #patching the prefix ZERO for segment + if [ $SegNum -lt 10 ]; then + segPrefix="0" + else + segPrefix="" + fi + + fileName=rawdata/run${runNum}/run-${runPrefix}${runNum}-${segPrefix}${SegNum}.evt + outFileName=pixie_evt/pixie-${runPrefix}${runNum}-${segPrefix}${SegNum}.evt + if [[ -f ${outFileName} ]]; then + echo "${outFileName} exist. Skip." + else + echo ${fileName} + ./nscl2pixie_haha ${fileName} ${outFileName} & + fi + ((j += 1)) + done + + while true; do + nthings=`ps aux | grep "[n]scl2pixie_haha"| wc -l` + echo $(date)" || "$nthings" nscl2pixie_haha are running." + + + if [[ $nthings > 0 ]]; then + sleep 10; + else + break + fi + done + + + + ((i += 1)) +done + + + +#Pixiefiles=$(ls -1 evt_data/pixie-0${run}*.evt) + +#for a in ${Pixiefiles} +#do +# echo ${a} +# echo ${a:15:7} +# ./armory/MergeEVT root_data/run-${a:15:7}_raw.root ${a} +#done + + + + diff --git a/correction_Clover.dat b/correction_Clover.dat new file mode 100644 index 0000000..51b4fe2 --- /dev/null +++ b/correction_Clover.dat @@ -0,0 +1,52 @@ +-0.2745420387 0.1903055400 + 0.1893425548 0.1818331614 +-0.0367071701 0.1824341324 +-0.1063887098 0.1910650773 +-0.1850638424 0.1946547185 +-0.2585036758 0.1865086505 + 0.2311507914 0.1930849561 +-0.0611156448 0.1880617739 + 0.2138792456 0.1857668487 +-0.1000184518 0.1830644733 + 2.0604294051 0.1385512827 + 0.3180233659 0.1386769235 + 0.0000000000 0.0000000000 + 0.0000000000 0.0000000000 + 0.0000000000 0.0000000000 + 0.0000000000 0.0000000000 +-0.2208351518 0.1364541549 +-0.0306010261 0.1783805884 +-0.0037076212 0.1960560535 + 0.0642160243 0.1839858595 +-0.0511685756 0.1859194405 +-0.1172917454 0.1856488571 +-0.1286314737 0.1872948275 +-0.1364423663 0.1855266236 +-0.0008503082 0.1902808189 + 0.1662583256 0.1907490728 + 0.0557733175 0.1920862347 + 0.1879598992 0.1887119382 + 0.0000000000 0.0000000000 + 0.0000000000 0.0000000000 + 0.0000000000 0.0000000000 + 0.0000000000 0.0000000000 +-0.9876919208 0.4327488664 + 0.2788794959 0.4412627192 +-0.9159334953 0.4314262803 + 0.1617243872 0.1540839667 +-0.0582723977 0.2020260902 +-0.2204402251 0.1920815477 + 0.0312224501 0.1779034597 +-0.0839891831 0.1812223578 + 0.0842978355 0.1560970523 +-0.1401545082 0.1992025914 + 0.2404464347 0.1853884121 + 0.0029329553 0.1837522238 + 1.0202556280 0.1941341886 + 1.0659732514 0.1919353184 +-0.2565005696 0.1973878067 + 0.0787616555 0.1961591821 +-0.0159272801 0.1863936847 + 0.4103358572 0.1853258451 + 0.3394901400 0.1936672410 +-0.2622309571 0.1876968520 diff --git a/correction_e.dat b/correction_e.dat new file mode 120000 index 0000000..c1b58f6 --- /dev/null +++ b/correction_e.dat @@ -0,0 +1 @@ +degai.cal \ No newline at end of file diff --git a/degai.cal b/degai.cal new file mode 100644 index 0000000..1502651 --- /dev/null +++ b/degai.cal @@ -0,0 +1,64 @@ + -0.106851 0.190142 0.000000 + 0.016704 0.181833 0.000000 + -0.145046 0.182478 0.000000 + -0.249700 0.191049 0.000000 + + -0.249730 0.194651 0.000000 + -0.267311 0.186496 0.000000 + 0.203408 0.193090 0.000000 + -0.317856 0.188170 0.000000 + + 0.305632 0.185745 0.000000 + -0.230991 0.183063 0.000000 + -0.043614 0.138764 0.000000 + 0.386491 0.138711 0.000000 + + 0.000000 1.000000 0.000000 + 0.000000 1.000000 0.000000 + 0.000000 1.000000 0.000000 + 0.000000 1.000000 0.000000 + + -0.353672 0.136451 0.000000 + -0.189959 0.178413 0.000000 + -0.063779 0.195860 0.000000 + -0.134865 0.183998 0.000000 + + -0.220694 0.185943 0.000000 + -0.165269 0.185650 0.000000 + -0.205702 0.187253 0.000000 + -0.139217 0.185499 0.000000 + + -0.031178 0.191581 0.000000 + -0.023608 0.190817 0.000000 + -0.100241 0.193328 0.000000 + 0.396967 0.188671 0.000000 + + 0.000000 1.000000 0.000000 + 0.000000 1.000000 0.000000 + 0.000000 1.000000 0.000000 + 0.000000 1.000000 0.000000 + + -0.974461 0.432854 0.000000 + -0.279092 0.441360 0.000000 + -0.698628 0.431157 0.000000 + 0.148492 0.154080 0.000000 + + -0.198020 0.202018 0.000000 + -0.281912 0.192079 0.000000 + -0.012972 0.177876 0.000000 + -0.337566 0.181243 0.000000 + + 0.056818 0.156063 0.000000 + -0.248622 0.199156 0.000000 + -0.046962 0.185395 0.000000 + 0.175072 0.183612 0.000000 + + -0.111547 0.195441 0.000000 + -0.328011 0.190841 0.000000 + -0.501598 0.197271 0.000000 + -0.121546 0.196111 0.000000 + + -0.147028 0.186377 0.000000 + 0.161169 0.185324 0.000000 + 0.157357 0.193662 0.000000 + -0.348731 0.187723 0.000000 diff --git a/mapping.h b/mapping.h new file mode 100644 index 0000000..7bd1f30 --- /dev/null +++ b/mapping.h @@ -0,0 +1,56 @@ +/************************************ +Clover : 0 - 99 +LaBr : 100 - 199 +Implat : 200 - 299, 200-209 low gain, 250-259 highgain, 290 = veto front, 291 = veto rear +beam line : 300 - 399, 300-309 Image Scint, 310-319 PPAC, 320-329 cross scint, 330-339 MSX +TAC : 400 - 499 + * *********************************/ +#ifndef MAPPING +#define MAPPING + +//==================== mapping + +#define NCLOVER 13 +#define NCRYSTAL NCLOVER*4 +#define NIMPLAT 10 +#define NVETO 2 +#define NLABR 15 +#define NBEAM 15 + +// 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 +int mapping[416] ={ -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-0 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-1 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-2 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-3 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-4 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-5 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-6 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-7 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-8 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-9 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-10 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-11 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-0, mod-12 + + 250, 200, -1, -1, 251, 252, 253, 254, 201, 202, 203, 204, -1, -1, -1, 500, //crate-1, mod-0 + -1, -1, -1, -1, -1, -1, -1, -1, 306, 310, 300, 301, 302, 303, 304, 305, //crate-1, mod-1 + 308, 309, -1, -1, 311, 312, 313, 314, 307, -1, -1, -1, -1, -1, -1, -1, //crate-1, mod-2 + 0, 1, 2, 3, 4, 5, 6, 7, 8, -1, -1, 11, -1, -1, -1, -1, //crate-1, mod-3 + 16, 17, 18, 19, 20, 21, 22, 23, -1, -1, -1, 27, -1, -1, -1, -1, //crate-1, mod-4 + 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, -1, -1, 46, 47, //crate-1, mod-5 + 48, 49, 50, 51, 9, 10, 24, 25, 26, -1, -1, -1, -1, -1, -1, -1, //crate-1, mod-6 + 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, -1, //crate-1, mod-7 + 290, 291, -1, -1, 44, 45, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-1, mod-8 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-1, mod-9 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-1, mod-10 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //crate-1, mod-11 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}; //crate-1, mod-12 + + +int mapDetID(int crate, int slot, int ch){ + int id = crate * 208 + 16*(slot-2) + ch; + //printf("id : %d, detID : %d \n", id, mapping[id]); + return mapping[id]; +} + +#endif diff --git a/nscl2pixie.c b/nscl2pixie.c new file mode 100644 index 0000000..e77800f --- /dev/null +++ b/nscl2pixie.c @@ -0,0 +1,288 @@ +#include +#include + +/* +Ring Item structure + +NSCLDAQ >= 11.0 +Ring Item Header + Size (4) + Type (4) (physics item, type=30) + +Ring Item Body + Ring Item Body Header + Size (4) + Timestamp (8) + Source ID (4) + Barrier Type (4) + + Body (for physics event) + Size (4) + Fragment #0 + Fragment Header (20) + Timestamp (8) + Source ID (4) + Payload size in bytes (4) + Barrier Type (4) + Fragment Payload + Ring Item Header (8) + Size (4) + Type (4) (physics item, type=30) + Ring Item Body Header (20) + Size (4) + Timestamp (8) + Source ID (4) + Barrier Type (4) + Ring Item Body (>8) + Body Size (4) - THIS IS IN 16-bit WORDS, NOT BYTES + Device Info (4) + Raw Data (XIA format) + Fragement #1 + Fragement Header (20) + Ring Item Header (8) + Ring Item Body Header (20) + Ring Item Body + . + . + . +*/ + +bool debug = false; + +int readRingItemHeader(FILE *file, unsigned int *ri_size) { + unsigned int ri_type; + if (fread(ri_size, (size_t)sizeof(int), 1, file) != 1) { + return -1; + } + if (fread(&ri_type, (size_t)sizeof(int), 1, file) != 1) { + return -1; + } + if (ri_type != 30 ) { + //non-physics item + if (ri_type == 1) { + if( debug) printf("\n============== BEGIN RUN ===============\n"); + } + else if (ri_type == 2) { + if( debug) printf("\n============== END RUN ===============\n"); + } + else if (ri_type == 3) { + if( debug) printf("\n============== PAUSE RUN ===============\n"); + } + else if (ri_type == 4) { + if( debug) printf("\n============== RESUME RUN ===============\n"); + } + else if (ri_type == 20) { + if( debug) printf("\nSCALERS FOUND..."); + } + else { + if( debug) printf("\nUntreated Ring item: %i %i\n", *ri_size, ri_type); + } + } + return ri_type; +} + +int readRingItemBodyHeader(FILE *file) { + unsigned int ribh_size; + unsigned long long int ribh_ts; + unsigned int ribh_sid; + unsigned int ribh_bt; + if (fread(&ribh_size, (size_t)sizeof(int), 1, file) != 1) { + return -1; + } + if (fread(&ribh_ts, (size_t)sizeof(unsigned long long int), 1, file) != 1) { + return -1; + } + if (fread(&ribh_sid, (size_t)sizeof(int), 1, file) != 1) { + return -1; + } + if (fread(&ribh_bt, (size_t)sizeof(int), 1, file) != 1) { + return -1; + } + + /* + std::cout << "Ring Item Body Header: " << std::endl; + std::cout << " size : " << ribh_size << std::endl; + std::cout << " timestamp : " << ribh_ts << std::endl; + std::cout << " source ID : " << ribh_sid << std::endl; + std::cout << " barrier type : " << ribh_bt << std::endl; + */ + return ribh_size; +} + +int readRingItemBody(FILE *file, unsigned int *rib_size) { + if (fread(rib_size, (size_t)sizeof(int), 1, file) != 1) { + return -1; + } + //std::cout << "Reading ring item body size = " << rib_size << std::endl; + *rib_size -= 4; + return *rib_size; +} + +int readNextFragment(FILE *file, unsigned int *rib_size) { + //this skips over the fragment header, the ring item header, and the ring item body header associated with the fragment + //reads the first 2 words of the ring item body (body size + device info) + //leaves the pointer right at the beginning of the actual PIXIE fragment + if (*rib_size > 48) { + //std::cout << "rib_size = " << rib_size << std::endl; + unsigned long long int frag_ts; + unsigned int frag_sid; + unsigned int frag_paysize; + unsigned int frag_bt; + + unsigned int frag_rih_size; + unsigned int frag_rih_type; + + unsigned int frag_ribh_size; + unsigned long long int frag_ribh_ts; + unsigned int frag_ribh_sid; + unsigned int frag_ribh_bt; + + if (fread(&frag_ts, (size_t)sizeof(unsigned long long int), 1, file) != 1) { return -1; } + if (fread(&frag_sid, (size_t)sizeof(int), 1, file) != 1) { return -1; } + if (fread(&frag_paysize, (size_t)sizeof(int), 1, file) != 1) { return -1; } + if (fread(&frag_bt, (size_t)sizeof(int), 1, file) != 1) { return -1; } + + if (fread(&frag_rih_size, (size_t)sizeof(int), 1, file) != 1) { return -1; } + if (fread(&frag_rih_type, (size_t)sizeof(int), 1, file) != 1) { return -1; } + + if (fread(&frag_ribh_size, (size_t)sizeof(int), 1, file) != 1) { return -1; } + if (fread(&frag_ribh_ts, (size_t)sizeof(unsigned long long int), 1, file) != 1) { return -1; } + if (fread(&frag_ribh_sid, (size_t)sizeof(int), 1, file) != 1) { return -1; } + if (fread(&frag_ribh_bt, (size_t)sizeof(int), 1, file) != 1) { return -1; } + + /* + std::cout << "Fragment Header: " << std::endl; + std::cout << " Timestamp: " << frag_ts << std::endl; + std::cout << " Source ID: " << frag_sid << std::endl; + std::cout << " Payload Size: " << frag_paysize << std::endl; + std::cout << " Barrier Type: " << frag_bt << std::endl; + + std::cout << "Fragment RIH: " << std::endl; + std::cout << " Size: " << frag_rih_size << std::endl; + std::cout << " Type: " << frag_rih_type << std::endl; + + + std::cout << "Fragment RIBH: " << std::endl; + std::cout << " Size: " << frag_ribh_size << std::endl; + std::cout << " Timestamp: " << frag_ribh_ts << std::endl; + std::cout << " Source ID: " << frag_ribh_sid << std::endl; + std::cout << " Barrier Type: " << frag_ribh_bt << std::endl; + */ + + unsigned int frag_size; //in 16 bit words not 32 bit words + unsigned int dev_info; + if (fread(&frag_size, (size_t)sizeof(int), 1, file) != 1) { + return -1; + } + if (fread(&dev_info, (size_t)sizeof(int), 1, file) != 1) { + return -1; + } + + /* + unsigned int ADC_freq = (dev_info & 0x0000ffff); + unsigned int ADC_res = (dev_info & 0xff0000) >> 16; + unsigned int mod_rev = (dev_info & 0xff000000) >> 24; + + std::cout << "dev_info = " << ADC_freq << " " << ADC_res << " " << mod_rev << std::endl; + + std::cout << "frag_size = " << frag_size << std::endl; + */ + + *rib_size -= (48 + frag_size*2); + if (*rib_size < 0 ) { + if( debug) printf("\nSEVERE ERROR, NAVIGATION THROUGH RING BUFFER BROKEN\n"); + return -1; + } + } + else if (*rib_size < 48) { + if( debug) printf("\nREMAINING DATA: \n"); + int rib_tmp = *rib_size; + int data; + while (rib_tmp > 0) { + fread(&data, (size_t)sizeof(int), 1, file); + if( debug) printf("%i\n", data); + rib_tmp -= 4; + } + } + else { + if( debug) printf("\nWarning! readNextFragment called when rib_size == 0\n"); + } + //std::cout << "readNextFragment returning " << rib_size << std::endl; + return *rib_size; +} + +int findNextFragment(FILE *file, unsigned int *rib_size) { + if (*rib_size > 0) { + return readNextFragment(file, rib_size); + } + + while (1) { + unsigned int ri_size; + int type = readRingItemHeader(file, &ri_size); + if (type == -1) { return -1; } + if (type != 30) { + fseek(file, ri_size-8, SEEK_CUR); //skip the rest of the non-physics event + continue; + } + int ribh_size = readRingItemBodyHeader(file); + if (ribh_size < 0) { return ribh_size; } + + readRingItemBody(file, rib_size); + if (*rib_size < 0) { return *rib_size; } + + return readNextFragment(file, rib_size); + } +} + +int copyPixieSubEvent(FILE *infile, FILE *outfile) { + unsigned int firstWords[4]; + if (fread(&firstWords, (size_t) sizeof(int)*4, (size_t) 1, infile) != 1) { + return -1; + } + fwrite(&firstWords, (size_t) sizeof(int)*4, (size_t) 1, outfile); + + int headerLength = (firstWords[0] & 0x1F000) >> 12; + int traceLength= (firstWords[3] & 0x7FFF0000) >> 16; + if (headerLength > 4) { + unsigned int otherWords[headerLength-4]; + if (fread(&otherWords, (size_t) 4, (size_t) headerLength-4, infile) != (size_t) headerLength-4) { + return -1; + } + fwrite(&otherWords, (size_t) 4, (size_t) headerLength-4, outfile); + } + if (traceLength > 0) { + unsigned short int trace[traceLength]; // = {0}; + if (fread(&trace[0], (size_t) 2, (size_t) traceLength, infile) != (size_t) traceLength) { + return -1; + } + fwrite(&trace[0], (size_t) 2, (size_t) traceLength, outfile); + } + return 1; +} + +int main(int argc, const char **argv) { + if (argc < 3) { + printf("usage: ./nscl2pixie [infile] [outfile] [debug=0]\n"); + return -1; + } + + FILE *infile = fopen(argv[1], "rb"); + FILE *outfile = fopen(argv[2], "wb"); + if( argc == 4 ) debug = atoi(argv[3]); + + unsigned int rib_size = 0; + long long unsigned int evts = 0; + int retval = 1; + while (1) { + retval = findNextFragment(infile, &rib_size); + if (retval < 0 ) { break; } + retval = copyPixieSubEvent(infile, outfile); + if (retval < 0 ) { break; } + ++evts; + } + + printf("%llu subevents copied\n", evts); + fclose(infile); + fclose(outfile); + +} diff --git a/peachCake.C b/peachCake.C new file mode 100644 index 0000000..5098ba0 --- /dev/null +++ b/peachCake.C @@ -0,0 +1,380 @@ +#define peachCake_cxx + +#include "peachCake.h" +#include +#include +#include +#include +#include +#include + +#include "./armory/AnalysisLibrary.h" +//============== histograms + +TH2F * hCloverID; +TH2F * hPID; +TH2F * hImplantIons; +TH2F * hImplantBeta; //high gain + +TH2F * hImplantX; +TH2F * hImplantY; + +TH2F * hImplantDyIons; +TH2F * hImplantDyBeta; + +TH1F * hDecay; +TH1F * hDecay_veto; + +TH1F * hFlag; +TH1F * hvetoFlag; + +//============ parameters +//TString cutFileName="PIDCuts.root"; +//TFile * cutFile; +//TCutG * cut = NULL; +//TObjArray * cutList; +//int numCut = 0; + +vector> eCorr; + + +void peachCake::Begin(TTree * /*tree*/){ + TString option = GetOption(); + + hCloverID = new TH2F("hCloverID", "Clover; ID; energy [keV]", 52, 0, 52, 400, 0, 2000); + hPID = new TH2F("hPID", "PID; ns; msx40", 500, -265, -220, 500, 0, 6500); + + hImplantIons = new TH2F("hImplantIons", "Implat low gain; x ; y", 400, 0, 1, 400, 0, 1); + hImplantBeta = new TH2F("hImplantBeta", "Implat high gain; x ; y", 400, 0, 1, 400, 0, 1); + + hImplantX = new TH2F("hImplantX", "Implat X; low ; high", 400, 0, 1, 400, 0, 1); + hImplantY = new TH2F("hImplantY", "Implat Y; low ; high", 400, 0, 1, 400, 0, 1); + + + hImplantDyIons = new TH2F("hImplantDyIonsow", "Implat low; sum_a; dy", 400, 0, 80000, 400, 0, 80000); + hImplantDyBeta = new TH2F("hImplantDyBetaigh", "Implat high; sum_a; dy", 400, 0, 8000, 400, 0, 8000); + + hDecay = new TH1F("hDecay", "Decay (beta.time - ion.time) ; [ticks]; counts", 100, 0, 2000); + hDecay_veto = new TH1F("hDecay_veto", "Decay (beta.time - ion.time) [veto]; [ticks]; counts", 100, 0, 2000); + + hFlag = new TH1F("hFlag", "Flag ( 1 = beam, 2 = Ions, 4 = Beta)", 8, 0, 8); + hvetoFlag = new TH1F("hvetoFlag", "vetoFlag ( 1 = front, 4 = rear)", 8, 0, 8); + + eCorr = LoadCorrectionParameters("correction_e.dat"); + + /** + cutFile = new TFile(cutFileName, "READ"); + bool listExist = cutFile->GetListOfKeys()->Contains("cutList"); + if( listExist ) { + cutList = (TObjArray*) cutFile->FindObjectAny("cutList"); + numCut = cutList->GetLast()+1; + printf("----- found %d cuts \n", numCut); + for( int k = 0; k < numCut; k++){ + if( cutList->At(k) != NULL ){ + printf("found a cut at %2d \n", k); + } + } + } + * */ + + printf("============ start \n"); + +} + + +Bool_t peachCake::Process(Long64_t entry){ + + nProcessed++; + + b_multi->GetEntry(entry); + b_detID->GetEntry(entry); + b_e->GetEntry(entry); + b_e_t->GetEntry(entry); + b_cfd->GetEntry(entry); + + energy = TMath::QuietNaN(); + Long64_t startTimeL = 0; + Long64_t startTimeR = 0; + Long64_t stopTime = 0; + double cfdL = TMath::QuietNaN(); + double cfdR = TMath::QuietNaN(); + double cfd2 = TMath::QuietNaN(); + + Long64_t startTime40 = 0; + double cfd40 = TMath::QuietNaN(); + + double a_L[5] = {TMath::QuietNaN()}; ///0 = dy, 1-4 = anode + double a_H[5] = {TMath::QuietNaN()}; + + veto_f = TMath::QuietNaN(); + veto_r = TMath::QuietNaN(); + + veto_f_time = 0; + veto_r_time = 0; + + crossEnergy = 0; /// for analog signal, cross T2 + crossTime = 0; + + for( int i = 0; i < 4 ; i++) { + dyIonsTime[i] = 0; + dyBetaTime[i] = 0; + dyIonsEnergy[i] = 0; + dyBetaEnergy[i] = 0; + } + + int multiDyBeta = 0; + int multiDyIons = 0; + + flag = 0; + vetoFlag = 0; + + //======= Scanning the event + for( int i = 0; i < multi; i++){ + + if( 0 <= detID[i] && detID[i] < 100 ){ // Clover + if( e[i] > 0 ) { + int id = detID[i]; + double eCal = 0; + for( int k = 0; k < int(eCorr[id].size()); k++){ + eCal += eCorr[id][k] * TMath::Power(e[i], k); + } + hCloverID->Fill(detID[i], eCal); + } + } + + if( detID[i] == 200 ) { //dy Beta + if ( multiDyBeta < 4 ){ + dyBetaTime[multiDyBeta] = e_t[i]; + dyBetaEnergy[multiDyBeta] = e[i]; + multiDyBeta ++; + } + } + if( 201 <= detID[i] && detID[i] < 250){ // Implant detector, high gain, Beta + a_H[detID[i] - 200] = e[i]; + } + + if( detID[i] == 250 ) { //dy Ions + if ( multiDyIons < 4 ){ + dyIonsTime[multiDyIons] = e_t[i]; + dyIonsEnergy[multiDyIons] = e[i]; + multiDyIons ++; + } + } + if( 251 <= detID[i] && detID[i] < 260){ // Implant detector, low gain, Ions + a_L[detID[i] - 250] = e[i]; + } + + if( detID[i] == 290 ) { + veto_f = e[i]; + veto_f_time = e_t[i]; + if( veto_f > 0 && (vetoFlag & 1) == 0) vetoFlag += 1; + } + if( detID[i] == 291 ) { + veto_r = e[i]; + veto_r_time = e_t[i]; + if( veto_r > 0 && (vetoFlag & 2) == 0) vetoFlag += 2; + } + + if( detID[i] == 300 ) { // image scint L + if( e[i] > 1000 ) { + startTimeL = e_t[i]; + cfdL = cfd[i]/16384.; + } + } + if( detID[i] == 301 ) { // image scint R + startTimeR = e_t[i]; + cfdR = cfd[i]/16384.; + } + if( detID[i] == 306 ) { // cross T2 analog + if( (vetoFlag & 4) == 0) vetoFlag += 4; + crossEnergy = e[i]; + crossTime = e_t[i]; + } + if( detID[i] == 307 ) { // cross B2 logic + stopTime = e_t[i]; + cfd2 = cfd[i]/16384.; + if( (vetoFlag & 4) == 0) vetoFlag += 4; + } + if( detID[i] == 308 ) energy = e[i]; //MSX40 + //if( detID[i] == 309 ) energy = e[i]; //MSX100 + + if( detID[i] == 310 ) { //MSX40 logic + startTime40 = e_t[i]; + cfd40 = cfd[i]/16384.; + } + } + + //================================== PID + Long64_t startTime = startTimeL; + double startCFD = cfdL; + + TOF = TMath::QuietNaN(); + if( startTime > 0 && stopTime > 0 ){ + if( stopTime > startTime ){ + TOF = double(stopTime - startTime) - startCFD + cfd2; + }else{ + TOF = -double(startTime - stopTime) - startCFD + cfd2 ; + } + TOF = TOF * 8; // to ns + + hPID->Fill(TOF, energy); + if( energy > 0 ) flag += 1; + } + + //================================== Implant + double sumA_L = 0; + double sumA_H = 0; + for( int i = 1; i <= 4 ; i++){ + if( a_L[i] > 0 ) sumA_L += a_L[i]; + if( a_H[i] > 0 ) sumA_H += a_H[i]; + } + + xIons = (a_L[3]+a_L[2])/sumA_L; + yIons = (a_L[1]+a_L[2])/sumA_L; + xBeta = (a_H[3]+a_H[2])/sumA_H; + yBeta = (a_H[1]+a_H[2])/sumA_H; + + + if( (flag & 1) == 0 ) { + hImplantIons->Fill(xIons, yIons); + hImplantBeta->Fill(xBeta, yBeta); + hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]); + hImplantDyBeta->Fill(sumA_H, dyBetaEnergy[0]); + + hImplantX->Fill(xIons, xBeta); + hImplantY->Fill(yIons, yBeta); + } + + if( !TMath::IsNaN(veto_f) && !TMath::IsNaN(veto_r) && crossTime == 0 ){ + //hImplantIons->Fill(xIons, yIons); + //hImplantBeta->Fill(xBeta, yBeta); + //hImplantDyIons->Fill(sumA_L, dyIonsEnergy[0]); + //hImplantDyBeta->Fill(sumA_H, dyBetaEnergy[0]); + // + //hImplantX->Fill(xIons, xBeta); + //hImplantY->Fill(yIons, yBeta); + } + + if( 0 < xIons && xIons < 1 && 0 < yIons && yIons < 1 ) flag += 2; + if( 0 < xBeta && xBeta < 1 && 0 < yBeta && yBeta < 1 ) flag += 4; + + hFlag->Fill(flag); + hvetoFlag->Fill(vetoFlag); + + + //============================== debug + if ( false ) { + + printf("################# %lld, multi:%d\n", entry, multi); + + printf("----- beam Line \n"); + printf(" MSX100 eneergy : %f \n", energy); + printf(" startTime : %llu, CFD : %f \n", startTime, startCFD); + printf(" stopTime : %llu, CFD : %f \n", stopTime, cfd2); + printf(" TOF : %f ns \n", TOF); + + printf("----- Ions, low gain \n"); + for( int i = 1; i <= 4; i++ )printf("a%d : %f\n", i, a_L[i]); + printf("sum A : %f \n", sumA_L); + printf("x : %f , y : %f \n", xIons, yIons); + for ( int i = 0; i < 4; i++ )printf("dy : %u, %llu \n", dyIonsEnergy[i], dyIonsTime[i]); + + printf("----- Beta, high gain \n"); + for( int i = 1; i <= 4; i++ )printf("a%d : %f\n", i, a_H[i]); + printf("sum A : %f \n", sumA_H); + printf("x : %f , y : %f \n", xBeta, yBeta); + for ( int i = 0; i < 4; i++ )printf("dy : %u, %llu \n", dyBetaEnergy[i], dyBetaTime[i]); + + printf("----- Veto \n"); + printf(" cross Time : %llu \n", crossTime); + printf(" cross energy : %d \n", crossEnergy); + printf("veto_f energy : %f \n", veto_f); + printf( "veto_f Time : %llu \n", veto_f_time); + printf("veto_r energy : %f \n", veto_r); + printf( "veto_r Time : %llu \n", veto_r_time); + + printf("============ flag : %d, vetoFlag : %d\n", flag, vetoFlag); + + + } + + //============================= Rejection + if ( flag == 0 ) return kTRUE; + + //============================= Progress + saveFile->cd(); + newTree->Fill(); + + + clock.Stop("timer"); + Double_t time = clock.GetRealTime("timer"); + clock.Start("timer"); + + if ( !shown ) { + if (fmod(time, 10) < 1 ){ + printf( "%12llu[%2d%%]|%3.0f min %5.2f sec | expect:%5.2f min\r", + nProcessed, + TMath::Nint((nProcessed+1)*100./totnumEntry), + TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60., + totnumEntry*time/(nProcessed+1.)/60.); + shown = 1;} + }else{ + if (fmod(time, 10) > 9 ){ + shown = 0; + } + } + + return kTRUE; +} + + +void peachCake::Terminate(){ + + printf("\n===============\n"); + saveFile->cd(); + newTree->Write(); + saveFile->Close(); + + gStyle->SetOptStat("neiou"); + TCanvas * c1 = new TCanvas("c1", "c1", 1500, 1000); + c1->Divide(3,3); + + int padID=1; + + c1->cd(padID); + c1->cd(padID)->SetLogz(); + hPID->Draw("colz"); + + padID ++; + c1->cd(padID); + c1->cd(padID)->SetGrid(); + c1->cd(padID)->SetLogz(); + hCloverID->Draw("colz"); hCloverID->SetNdivisions(-13); + + padID ++; + c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantIons->Draw("colz"); + + padID ++; + c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantBeta->Draw("colz"); + + //padID ++; + //c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantBeta_veto->Draw("colz"); + + padID ++; + c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantDyIons->Draw("colz"); + + padID ++; + c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantDyBeta->Draw("colz"); + + padID ++; + c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantX->Draw("colz"); + + padID ++; + c1->cd(padID); c1->cd(padID)->SetLogz(); hImplantY->Draw("colz"); + + padID ++; + c1->cd(padID); hFlag->Draw(""); + hvetoFlag->SetLineColor(2); hvetoFlag->Draw("same"); + + +} diff --git a/peachCake.h b/peachCake.h new file mode 100644 index 0000000..6ca68f4 --- /dev/null +++ b/peachCake.h @@ -0,0 +1,201 @@ +////////////////////////////////////////////////////////// +// This class has been automatically generated on +// Fri May 13 17:59:58 2022 by ROOT version 6.24/06 +// from TTree tree/root_data/run-0143-00.root +// found on file: root_data/run-0143-00.root +////////////////////////////////////////////////////////// + +#ifndef peachCake_h +#define peachCake_h + +#include +#include +#include +#include +#include +#include "cmath" + +// Header file for the classes stored in the TTree if any. + +class peachCake : public TSelector { +public : + TChain *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 detID[416]; //[multi] + Double_t e[416]; //[multi] + ULong64_t e_t[416]; //[multi] + Int_t cfd[416]; //[multi] + Int_t qdc[416][8]; //[multi] + Int_t multiClover; + Int_t multiBeam; + Int_t runID; + + // List of branches + TBranch *b_event_ID; //! + TBranch *b_multi; //! + TBranch *b_detID; //! + TBranch *b_e; //! + TBranch *b_e_t; //! + TBranch *b_cfd; //! + TBranch *b_qdc; //! + TBranch *b_multiplicity_Clover; //! + TBranch *b_multiplicity_Beam; //! + TBranch *b_runID; //! + + peachCake(TTree * /*tree*/ =0) : fChain(0) { } + virtual ~peachCake() { } + 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(peachCake,0); + + //=========== varibales in new tree + TFile * saveFile; + TTree * newTree; + TString saveFileName; + int totnumEntry; //of the original file + + //tree + Int_t eventID; + UInt_t crossEnergy; + ULong64_t crossTime; + Double_t TOF; + Double_t energy; + + Short_t flag; /// 1 = has beam, 2 = has Ions, 4 = hasBeta, default 0; + + Double_t xIons; + Double_t yIons; + Double_t xBeta; + Double_t yBeta; + UInt_t dyIonsEnergy[4]; + ULong64_t dyIonsTime[4]; + UInt_t dyBetaEnergy[4]; + ULong64_t dyBetaTime[4]; + Double_t veto_f; + ULong64_t veto_f_time; + Double_t veto_r; + ULong64_t veto_r_time; + + Short_t vetoFlag; /// default = 0, 1 = front, 2 = rear, 4 = cross + + //clock + TBenchmark clock; + Bool_t shown; + + Long64_t nProcessed; + +}; + +#endif + +#ifdef peachCake_cxx +void peachCake::Init(TTree *tree){ + + // Set branch addresses and branch pointers + if (!tree) return; + + totnumEntry = tree->GetEntries(); + printf("============= num entry : %d\n", totnumEntry); + + fChain = (TChain *)tree; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("evID", &evID, &b_event_ID); + fChain->SetBranchAddress("multi", &multi, &b_multi); + fChain->SetBranchAddress("detID", detID, &b_detID); + fChain->SetBranchAddress("e", e, &b_e); + fChain->SetBranchAddress("e_t", e_t, &b_e_t); + fChain->SetBranchAddress("cfd", cfd, &b_cfd); + fChain->SetBranchAddress("qdc", qdc, &b_qdc); + fChain->SetBranchAddress("multiClover", &multiClover, &b_multiplicity_Clover); + fChain->SetBranchAddress("multiBeam", &multiBeam, &b_multiplicity_Beam); + fChain->SetBranchAddress("runID", &runID, &b_runID); + + + //============ new tree + saveFileName = "test"; + int oldSeg = -100; + int numFile = fChain->GetListOfFiles()->GetLast()+1; + for( int i = 0; i < numFile; i++){ + TString name = fChain->GetListOfFiles()->At(i)->GetTitle(); + int pos = name.Last('/'); + name = name.Remove(0,pos+1); //this should be run-XXXX-XX.root + if( i == 0 ) { + saveFileName = name; + pos = saveFileName.Last('-'); + saveFileName = saveFileName.Remove(pos); //run-XXXX + } + pos = name.Last('-'); + int segNum = atoi(name.Remove(0, pos+1).Remove(2)); + if( oldSeg == -100 ) oldSeg = segNum; + saveFileName = saveFileName + Form("_%02d",segNum); + if( segNum == oldSeg + 1 && i != numFile -1 ){ + int len = saveFileName.Length(); + saveFileName = saveFileName.Remove(len - 3); + oldSeg = segNum; + } + } + saveFileName = saveFileName + ".root"; + + printf("=========== saveFile : %s \n", saveFileName.Data()); + + saveFile = new TFile(saveFileName, "recreate"); + newTree = new TTree("tree", "tree"); + + newTree->Branch("eventID", &eventID, "eventID/l"); + //newTree->Branch("runID", &runID, "runID/I"); + newTree->Branch("TOF", &TOF, "TOF/D"); + newTree->Branch("energy", &energy, "energy/D"); + newTree->Branch("crossTime", &crossTime, "crossTime/l"); + newTree->Branch("crossEnergy", &crossEnergy, "crossEnergy/i"); + newTree->Branch("flag", &flag, "flag/S"); + newTree->Branch("vetoFlag", &vetoFlag, "vetoFlag/S"); + newTree->Branch("xIons", &xIons, "xIons/D"); + newTree->Branch("yIons", &yIons, "yIons/D"); + newTree->Branch("xBeta", &xBeta, "xBeta/D"); + newTree->Branch("yBeta", &yBeta, "yBeta/D"); + newTree->Branch("dyIonsTime", dyIonsTime, "dyIonsTime[4]/l"); + newTree->Branch("dyBetaTime", dyBetaTime, "dyBetaTime[4]/l"); + newTree->Branch("veto_f", &veto_f, "veto_f/D"); + newTree->Branch("veto_r", &veto_r, "veto_r/D"); + + //==== clock + clock.Reset(); + clock.Start("timer"); + shown = 0; + + nProcessed = 0; + +} + +Bool_t peachCake::Notify(){ + + return kTRUE; +} + +void peachCake::SlaveBegin(TTree * /*tree*/){ + TString option = GetOption(); +} + +void peachCake::SlaveTerminate(){ +} + + +#endif // #ifdef peachCake_cxx diff --git a/pixie_evt b/pixie_evt new file mode 120000 index 0000000..2409bc7 --- /dev/null +++ b/pixie_evt @@ -0,0 +1 @@ +/mnt/data_1TB/e21062/pixie_evt \ No newline at end of file diff --git a/rawdata b/rawdata new file mode 120000 index 0000000..aa88799 --- /dev/null +++ b/rawdata @@ -0,0 +1 @@ +/mnt/data_1TB/e21062/rawdata \ No newline at end of file diff --git a/rootlogon.C b/rootlogon.C new file mode 100644 index 0000000..1c81f74 --- /dev/null +++ b/rootlogon.C @@ -0,0 +1,42 @@ +{ + printf("=================== Root v%s\n", gROOT->GetVersion()); + gStyle->SetOptStat("neiou"); + + + gStyle->SetPalette(1); + gStyle->SetLineWidth(2); + gStyle->SetStatW(0.3); + gStyle->SetStatH(0.3); + gStyle->SetTitleW(0.4); + gStyle->SetTitleH(0.1); + gStyle->SetCanvasBorderMode(0); + gStyle->SetPadBorderMode(0); + gStyle->SetPadColor(0); + gStyle->SetFrameBorderMode(0); + gStyle->SetFrameLineWidth(2); + gStyle->SetCanvasColor(0); + gStyle->SetOptDate(4); + gStyle->GetAttDate()->SetTextFont(62); + gStyle->GetAttDate()->SetTextSize(0.02); + gStyle->GetAttDate()->SetTextAngle(0); + gStyle->GetAttDate()->SetTextAlign(11); + gStyle->GetAttDate()->SetTextColor(1); + gStyle->SetDateX(0); + gStyle->SetDateY(0); + + gROOT->SetStyle("Plain"); // plain histogram style + gStyle->SetOptStat("nemruoi"); // expanded stats box + gStyle->SetPadTickX(1); // tic marks on all axes + gStyle->SetPadTickY(1); // + gStyle->SetOptFit(1111); // the results of the fits + gStyle->SetPadGridX(kTRUE); // draw horizontal and vertical grids + gStyle->SetPadGridY(kTRUE); + gStyle->SetPalette(1,0); // pretty and useful palette + gStyle->SetHistLineWidth(2); // a thicker histogram line + gStyle->SetFrameFillColor(10); // a different frame colour + gStyle->SetTitleFillColor(33); // title colour to highlight it + gStyle->SetTitleW(.76); // Title Width + gStyle->SetTitleH(.07); // Title height + gStyle->SetHistMinimumZero(true); // Suppresses Histogram Zero Supression + +} diff --git a/script.C b/script.C new file mode 100644 index 0000000..9216510 --- /dev/null +++ b/script.C @@ -0,0 +1,10 @@ +{ + TChain * chain = new TChain("tree"); + + chain->Add("root_data/run-0238-[0-4][0-9].root"); + + chain->GetListOfFiles()->Print(); + + chain->Process("peachCake.C+"); + +}