initial commit
This commit is contained in:
commit
e3aa858ff0
8
.gitignore
vendored
Normal file
8
.gitignore
vendored
Normal file
|
|
@ -0,0 +1,8 @@
|
|||
*.pcm
|
||||
*.d
|
||||
*.so
|
||||
*.root
|
||||
|
||||
legacy_code
|
||||
nscl2pixie
|
||||
nscl2pixie_haha
|
||||
267
Cali_gamma_histogram.C
Normal file
267
Cali_gamma_histogram.C
Normal file
|
|
@ -0,0 +1,267 @@
|
|||
#include <TFile.h>
|
||||
#include <TTree.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TROOT.h>
|
||||
#include <TSystem.h>
|
||||
#include <TStyle.h>
|
||||
#include <TProfile.h>
|
||||
#include <TH2F.h>
|
||||
#include <TH1F.h>
|
||||
#include <TF1.h>
|
||||
#include <TMath.h>
|
||||
#include <TGraph.h>
|
||||
#include <TLine.h>
|
||||
#include <TSpectrum.h>
|
||||
#include "armory/AnalysisLibrary.h"
|
||||
|
||||
vector<double> 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<double> 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<double> energy;
|
||||
vector<double> refEnergy;
|
||||
vector<double> 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<double> 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<double> k1 = fitEnergy;
|
||||
vector<double> k2 = refEnergy;
|
||||
//===== when nPeaks != refEnergy.size(), need to matching the two vector size by checking the r-squared.
|
||||
if( option != -1 ){
|
||||
vector<vector<double>> 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<double> 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;
|
||||
|
||||
}
|
||||
246
DecayFinder.C
Normal file
246
DecayFinder.C
Normal file
|
|
@ -0,0 +1,246 @@
|
|||
#define DecayFinder_cxx
|
||||
|
||||
#include "DecayFinder.h"
|
||||
#include <TH2.h>
|
||||
#include <TH1.h>
|
||||
#include <TMath.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TStyle.h>
|
||||
#include <TCutG.h>
|
||||
#include <TObjArray.h>
|
||||
#include <TRandom.h>
|
||||
|
||||
//===================== 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");
|
||||
|
||||
}
|
||||
136
DecayFinder.h
Normal file
136
DecayFinder.h
Normal file
|
|
@ -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 <TROOT.h>
|
||||
#include <TChain.h>
|
||||
#include <TFile.h>
|
||||
#include <TSelector.h>
|
||||
#include <TBenchmark.h>
|
||||
|
||||
// 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
|
||||
54
PIDCutChecker.C
Normal file
54
PIDCutChecker.C
Normal file
|
|
@ -0,0 +1,54 @@
|
|||
#include <TH2F.h>
|
||||
#include <TFile.h>
|
||||
#include <TChain.h>
|
||||
#include <TROOT.h>
|
||||
#include <TStyle.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TCutG.h>
|
||||
#include <TString.h>
|
||||
#include <TObjArray.h>
|
||||
|
||||
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();
|
||||
|
||||
}
|
||||
90
PIDCutCreator.C
Normal file
90
PIDCutCreator.C
Normal file
|
|
@ -0,0 +1,90 @@
|
|||
#include <TH2F.h>
|
||||
#include <TFile.h>
|
||||
#include <TChain.h>
|
||||
#include <TROOT.h>
|
||||
#include <TStyle.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TCutG.h>
|
||||
#include <TString.h>
|
||||
#include <TObjArray.h>
|
||||
|
||||
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");
|
||||
|
||||
}
|
||||
296
armory/AnalysisLibrary.h
Normal file
296
armory/AnalysisLibrary.h
Normal file
|
|
@ -0,0 +1,296 @@
|
|||
#ifndef Analysis_Library_h
|
||||
#define Analysis_Library_h
|
||||
|
||||
#include <TF1.h>
|
||||
#include <TGraph.h>
|
||||
#include <TSpectrum.h>
|
||||
#include <TMath.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <vector>
|
||||
#include <string>
|
||||
|
||||
std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, int shift = 0){
|
||||
|
||||
std::vector<std::string> output;
|
||||
|
||||
size_t pos;
|
||||
do{
|
||||
pos = tempLine.find(splitter); // fine splitter
|
||||
if( pos == 0 ){ //check if it is splitter again
|
||||
tempLine = tempLine.substr(pos+1);
|
||||
continue;
|
||||
}
|
||||
|
||||
std::string secStr;
|
||||
if( pos == std::string::npos ){
|
||||
secStr = tempLine;
|
||||
}else{
|
||||
secStr = tempLine.substr(0, pos+shift);
|
||||
tempLine = tempLine.substr(pos+shift);
|
||||
}
|
||||
|
||||
//check if secStr is begin with space
|
||||
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<std::vector<double>> combination(std::vector<double> arr, int r){
|
||||
|
||||
std::vector<std::vector<double>> output;
|
||||
|
||||
int n = arr.size();
|
||||
std::vector<int> 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<double> 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<double> 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<double> dataX, std::vector<double> 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<std::vector<double>> FindMatchingPair(std::vector<double> enX, std::vector<double> enY){
|
||||
|
||||
//output[0] = fitEnergy;
|
||||
//output[1] = refEnergy;
|
||||
|
||||
int nX = enX.size();
|
||||
int nY = enY.size();
|
||||
|
||||
std::vector<double> fitEnergy;
|
||||
std::vector<double> refEnergy;
|
||||
|
||||
if( nX > nY ){
|
||||
|
||||
std::vector<std::vector<double>> 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<std::vector<double>> 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<std::vector<double>> haha;
|
||||
haha.push_back(fitEnergy);
|
||||
haha.push_back(refEnergy);
|
||||
|
||||
return haha;
|
||||
|
||||
}
|
||||
|
||||
|
||||
std::vector<std::vector<double>> LoadCorrectionParameters(TString corrFile, bool show=false){
|
||||
|
||||
printf(" load correction parameters : %s", corrFile.Data());
|
||||
std::ifstream file;
|
||||
file.open(corrFile.Data());
|
||||
|
||||
std::vector<std::vector<double>> corr;
|
||||
corr.clear();
|
||||
|
||||
std::vector<double> 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<std::string> 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<double> temp = {0, 1};
|
||||
for( int i = 0; i < 36; i++){
|
||||
corr.push_back(temp);
|
||||
}
|
||||
}
|
||||
|
||||
return corr;
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
195
armory/Analyzer_Utili.c
Normal file
195
armory/Analyzer_Utili.c
Normal file
|
|
@ -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<double> 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<double> peaks = fitAuto(he[i], BG, threshold, sigmaMax, peakDensity);
|
||||
vector<vector<double>> output = FindMatchingPair(peaks, refEnergy);
|
||||
|
||||
vector<double> haha1 = output[0];
|
||||
vector<double> 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);
|
||||
}
|
||||
|
||||
}
|
||||
2787
armory/AutoFit.C
Normal file
2787
armory/AutoFit.C
Normal file
File diff suppressed because it is too large
Load Diff
110
armory/DataBlock.h
Normal file
110
armory/DataBlock.h
Normal file
|
|
@ -0,0 +1,110 @@
|
|||
|
||||
#ifndef DATABLOCK_H
|
||||
#define DATABLOCK_H
|
||||
|
||||
#include <stdio.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#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
|
||||
BIN
armory/EventBuilder_evt
Executable file
BIN
armory/EventBuilder_evt
Executable file
Binary file not shown.
270
armory/EventBuilder_evt.cpp
Normal file
270
armory/EventBuilder_evt.cpp
Normal file
|
|
@ -0,0 +1,270 @@
|
|||
/**********************************************************/
|
||||
/* */
|
||||
/* Modified by Ryan From */
|
||||
/* */
|
||||
/* PXI SCAN CODE -- J.M. Allmond (ORNL) -- July 2016 */
|
||||
/* */
|
||||
/**********************************************************/
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
#include <stdbool.h>
|
||||
|
||||
#include "TFile.h"
|
||||
#include "TTree.h"
|
||||
#include "TMath.h"
|
||||
#include "TBenchmark.h"
|
||||
|
||||
#define 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;
|
||||
}
|
||||
|
||||
|
||||
206
armory/EventBuilder_seperated.cpp
Normal file
206
armory/EventBuilder_seperated.cpp
Normal file
|
|
@ -0,0 +1,206 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <math.h>
|
||||
#include <stdbool.h>
|
||||
|
||||
#include "TFile.h"
|
||||
#include "TTree.h"
|
||||
#include "TMath.h"
|
||||
#include "TBenchmark.h"
|
||||
#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] <timeWindows> <SaveFileName>\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);
|
||||
|
||||
}
|
||||
111
armory/LoadCorrection.h
Normal file
111
armory/LoadCorrection.h
Normal file
|
|
@ -0,0 +1,111 @@
|
|||
#ifndef Analysis_Library_h
|
||||
#define Analysis_Library_h
|
||||
|
||||
#include <TF1.h>
|
||||
#include <TGraph.h>
|
||||
#include <TSpectrum.h>
|
||||
#include <TMath.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <vector>
|
||||
#include <string>
|
||||
|
||||
std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, int shift = 0){
|
||||
|
||||
std::vector<std::string> output;
|
||||
|
||||
size_t pos;
|
||||
do{
|
||||
pos = tempLine.find(splitter); // fine splitter
|
||||
if( pos == 0 ){ //check if it is splitter again
|
||||
tempLine = tempLine.substr(pos+1);
|
||||
continue;
|
||||
}
|
||||
|
||||
std::string secStr;
|
||||
if( pos == std::string::npos ){
|
||||
secStr = tempLine;
|
||||
}else{
|
||||
secStr = tempLine.substr(0, pos+shift);
|
||||
tempLine = tempLine.substr(pos+shift);
|
||||
}
|
||||
|
||||
//check if secStr is begin with space
|
||||
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<std::vector<double>> LoadCorrectionParameters(TString corrFile, bool show=false){
|
||||
|
||||
printf(" load correction parameters : %s", corrFile.Data());
|
||||
std::ifstream file;
|
||||
file.open(corrFile.Data());
|
||||
|
||||
std::vector<std::vector<double>> corr;
|
||||
corr.clear();
|
||||
|
||||
std::vector<double> 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<std::string> 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<double> temp = {0, 1};
|
||||
for( int i = 0; i < 36; i++){
|
||||
corr.push_back(temp);
|
||||
}
|
||||
}
|
||||
|
||||
return corr;
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
BIN
armory/MergeEVT
Executable file
BIN
armory/MergeEVT
Executable file
Binary file not shown.
135
armory/MergeEVT.cpp
Normal file
135
armory/MergeEVT.cpp
Normal file
|
|
@ -0,0 +1,135 @@
|
|||
#include <stdio.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include "TFile.h"
|
||||
#include "TTree.h"
|
||||
#include "TString.h"
|
||||
#include "TMath.h"
|
||||
#include "TBenchmark.h"
|
||||
#include <vector>
|
||||
|
||||
#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);
|
||||
|
||||
}
|
||||
BIN
armory/ev22txt
Executable file
BIN
armory/ev22txt
Executable file
Binary file not shown.
81
armory/ev22txt.cpp
Normal file
81
armory/ev22txt.cpp
Normal file
|
|
@ -0,0 +1,81 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <iostream>
|
||||
#include <string.h>
|
||||
|
||||
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;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
BIN
armory/evt2hist
Executable file
BIN
armory/evt2hist
Executable file
Binary file not shown.
344
armory/evt2hist.cpp
Normal file
344
armory/evt2hist.cpp
Normal file
|
|
@ -0,0 +1,344 @@
|
|||
#include <stdio.h>
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include <vector>
|
||||
#include <thread>
|
||||
|
||||
#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<std::vector<double>> 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();
|
||||
|
||||
}
|
||||
373
armory/evtReader.h
Normal file
373
armory/evtReader.h
Normal file
|
|
@ -0,0 +1,373 @@
|
|||
#ifndef EVTREADER_H
|
||||
#define EVTREADER_H
|
||||
|
||||
#include <stdio.h> /// for FILE
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#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
|
||||
25
armory/makefile
Normal file
25
armory/makefile
Normal file
|
|
@ -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
|
||||
87
convertPixie.sh
Executable file
87
convertPixie.sh
Executable file
|
|
@ -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
|
||||
|
||||
|
||||
|
||||
|
||||
52
correction_Clover.dat
Normal file
52
correction_Clover.dat
Normal file
|
|
@ -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
|
||||
1
correction_e.dat
Symbolic link
1
correction_e.dat
Symbolic link
|
|
@ -0,0 +1 @@
|
|||
degai.cal
|
||||
64
degai.cal
Normal file
64
degai.cal
Normal file
|
|
@ -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
|
||||
56
mapping.h
Normal file
56
mapping.h
Normal file
|
|
@ -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
|
||||
288
nscl2pixie.c
Normal file
288
nscl2pixie.c
Normal file
|
|
@ -0,0 +1,288 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
/*
|
||||
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);
|
||||
|
||||
}
|
||||
380
peachCake.C
Normal file
380
peachCake.C
Normal file
|
|
@ -0,0 +1,380 @@
|
|||
#define peachCake_cxx
|
||||
|
||||
#include "peachCake.h"
|
||||
#include <TH2.h>
|
||||
#include <TH1.h>
|
||||
#include <TMath.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TStyle.h>
|
||||
#include <TCutG.h>
|
||||
|
||||
#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<vector<double>> 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");
|
||||
|
||||
|
||||
}
|
||||
201
peachCake.h
Normal file
201
peachCake.h
Normal file
|
|
@ -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 <TROOT.h>
|
||||
#include <TChain.h>
|
||||
#include <TFile.h>
|
||||
#include <TSelector.h>
|
||||
#include <TBenchmark.h>
|
||||
#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
|
||||
42
rootlogon.C
Normal file
42
rootlogon.C
Normal file
|
|
@ -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
|
||||
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user