rtang/XIAEventBuilder
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

snapshot

Browse Source
This commit is contained in:
Ryan Tang 2021-12-24 16:29:50 -05:00
parent baa211d30c
commit ac7246bc13
7 changed files with 97 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
Analyzer.C
View File

@ -29,7 +29,7 @@ vector<vector<double>> eCorr;
//############################################ histogram declaration //############################################ histogram declaration
TH2F * heVID; TH2F * hevID;
TH1F * he[NCRYSTAL]; TH1F * he[NCRYSTAL];
TH2F * hgg[NCRYSTAL][NCRYSTAL]; TH2F * hgg[NCRYSTAL][NCRYSTAL];
@ -37,7 +37,7 @@ TH2F * hgg[NCRYSTAL][NCRYSTAL];
TH2F * hcoin; TH2F * hcoin;
///----- after calibration and BGO veto ///----- after calibration and BGO veto
TH2F * heCalVID; TH2F * heCalvID;
TH1F * heCal[NCRYSTAL]; TH1F * heCal[NCRYSTAL];
TH2F * hcoinBGO; TH2F * hcoinBGO;
@ -50,8 +50,12 @@ void Analyzer::Begin(TTree * tree){
printf("======================== Histograms declaration\n"); printf("======================== Histograms declaration\n");
heVID = new TH2F("heVID", "e vs ID; det ID; e [ch]", NCRYSTAL, 0, NCRYSTAL, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]); hevID = new TH2F("hevID", "e vs ID; det ID; e [ch]", NCRYSTAL, 0, NCRYSTAL, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
heCalVID = new TH2F("heCalVID", Form("eCal vs ID (BGO veto > %.1f); det ID; Energy [keV]", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); heCalvID = new TH2F("heCalvID", Form("eCal vs ID (BGO veto > %.1f); det ID; Energy [keV]", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]);
hevID->SetNdivisions(-409, "X");
heCalvID->SetNdivisions(-409, "X");
for( int i = 0; i < NCRYSTAL; i ++){ for( int i = 0; i < NCRYSTAL; i ++){
he[i] = new TH1F( Form("he%02d", i), Form("e -%02d", i), rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]); he[i] = new TH1F( Form("he%02d", i), Form("e -%02d", i), rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
heCal[i] = new TH1F(Form("heCal%02d", i), Form("eCal -%02d (BGO veto > %.1f); Energy [keV]; count / %d keV", i, BGO_threshold, energyRange[0]), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); heCal[i] = new TH1F(Form("heCal%02d", i), Form("eCal -%02d (BGO veto > %.1f); Energy [keV]; count / %d keV", i, BGO_threshold, energyRange[0]), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]);
@ -110,7 +114,7 @@ Bool_t Analyzer::Process(Long64_t entry){
//if( pileup[detID] == 1 ) continue; //if( pileup[detID] == 1 ) continue;
//======== Fill raw data //======== Fill raw data
heVID->Fill( detID, e[detID]); hevID->Fill( detID, e[detID]);
he[detID]->Fill(e[detID]); he[detID]->Fill(e[detID]);
@ -137,7 +141,7 @@ Bool_t Analyzer::Process(Long64_t entry){
} }
heCalVID->Fill( detID, eCal); heCalvID->Fill( detID, eCal);
heCal[detID]->Fill(eCal); heCal[detID]->Fill(eCal);
for( int detJ = detID +1; detJ < NCRYSTAL; detJ++) { for( int detJ = detID +1; detJ < NCRYSTAL; detJ++) {
@ -165,11 +169,11 @@ void Analyzer::Terminate(){
cCanvas->cd(1); cCanvas->cd(1);
cCanvas->cd(1)->SetLogz(1); cCanvas->cd(1)->SetLogz(1);
heVID->Draw("colz"); hevID->Draw("colz");
cCanvas->cd(2); cCanvas->cd(2);
cCanvas->cd(2)->SetLogz(1); cCanvas->cd(2)->SetLogz(1);
heCalVID->Draw("colz"); heCalvID->Draw("colz");
cCanvas->cd(3); cCanvas->cd(3);
cCanvas->cd(3)->SetLogz(1); cCanvas->cd(3)->SetLogz(1);

88
Monitor_raw.C
View File

@ -11,7 +11,8 @@
//############################################ User setting //############################################ User setting
int rawEnergyRange[2] = {500, 6000}; // in ch int rawEnergyRange[2] = {500, 6000}; // in ch, {min, max}
int energyRange[3] = {2, 40, 2000}; // in keV, {resol, min, max}
TString e_corr = "correction_e.dat"; TString e_corr = "correction_e.dat";
@ -28,8 +29,10 @@ vector<vector<double>> eCorr;
//############################################ histogram declaration //############################################ histogram declaration
TH1F * hTDiff; TH1F * he[NCRYSTAL];
TH2F * hTDiffvEventID; TH1F * heCal[NCRYSTAL];
TH2F * heCalvID;
//############################################ BEGIN //############################################ BEGIN
void Monitor_raw::Begin(TTree * tree){ void Monitor_raw::Begin(TTree * tree){
@ -38,9 +41,27 @@ void Monitor_raw::Begin(TTree * tree){
NumEntries = tree->GetEntries(); NumEntries = tree->GetEntries();
printf("======================== Creating histograms\n"); printf("======================== Creating histograms\n");
for( int i = 0 ; i < NCRYSTAL ; i++){
he[i] = new TH1F (Form("he%02d", i), Form("e%02d", i), rawEnergyRange[1]-rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
heCal[i] = new TH1F (Form("heCal%02d", i), Form("e%02d (Cali.)", i), (energyRange[2]-energyRange[1])/energyRange[0], energyRange[1], energyRange[2]);
switch (i%4){
case 0: he[i]->SetLineColor(2);break;
case 1: he[i]->SetLineColor(kYellow+3);break;
case 2: he[i]->SetLineColor(kGreen+2);break;
case 3: he[i]->SetLineColor(4);break;
}
switch (i%4){
case 0: heCal[i]->SetLineColor(2);break;
case 1: heCal[i]->SetLineColor(kYellow+3);break;
case 2: heCal[i]->SetLineColor(kGreen+2);break;
case 3: heCal[i]->SetLineColor(4);break;
}
}
hTDiff = new TH1F("hTDiff", "time different between this and next event", 2000, -1000, 1000); heCalvID = new TH2F("heCalvID", "ID vs Energy (Cali.); ID; Energy", NCRYSTAL, 0, NCRYSTAL, (energyRange[2]-energyRange[1])/energyRange[0], energyRange[1], energyRange[2]);
hTDiffvEventID = new TH2F("hTDiffvEventID", "time different between this and next event; eventID; TDiff", 50, 0, maxEvent, 2000, -1000, 1000);
printf("======================== end of histograms creation.\n"); printf("======================== end of histograms creation.\n");
@ -61,7 +82,6 @@ Bool_t Monitor_raw::Process(Long64_t entry){
if (ProcessedEntries>NumEntries*Frac-1) { if (ProcessedEntries>NumEntries*Frac-1) {
TString msg; msg.Form("%llu", NumEntries/1000); TString msg; msg.Form("%llu", NumEntries/1000);
int len = msg.Sizeof(); int len = msg.Sizeof();
//printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\033[A\n",
printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n", printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n",
Frac*100, len, ProcessedEntries/1000,NumEntries/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac); Frac*100, len, ProcessedEntries/1000,NumEntries/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac);
StpWatch.Start(kFALSE); StpWatch.Start(kFALSE);
@ -73,19 +93,22 @@ Bool_t Monitor_raw::Process(Long64_t entry){
Terminate(); Terminate();
} }
b_ID->GetEntry(entry);
b_energy->GetEntry(entry); b_energy->GetEntry(entry);
b_time_stamp->GetEntry(entry); b_timestamp->GetEntry(entry);
ULong64_t t0 = t; int detID = mapping[ID];
if( entry < (Long64_t) NumEntries ) b_time_stamp->GetEntry(entry+1);; if( 0 <= detID && detID < 100 ){
he[detID]->Fill(e);
//double eCal = ApplyCorrection(eCorr, detID, e);
double eCal = eCorr[detID][0]+eCorr[detID][1]*e;
heCal[detID]->Fill(eCal);
heCalvID->Fill(detID, eCal);
}
ULong64_t t1 = t;
int tDiff = (int) t1 - t0;
hTDiff->Fill( tDiff );
hTDiffvEventID->Fill( entry, tDiff);
return kTRUE; return kTRUE;
} }
@ -96,18 +119,35 @@ void Monitor_raw::Terminate(){
printf("============================== finishing.\n"); printf("============================== finishing.\n");
gROOT->cd(); gROOT->cd();
TCanvas * cc = new TCanvas("cc", "cc", 2000, 1000); int nCrystalPerClover = 4;
int nClover = NCRYSTAL / nCrystalPerClover;
TCanvas * cc = new TCanvas("cc", "cc", 2000, 2000);
if( cc->GetShowEventStatus() == 0 ) cc->ToggleEventStatus(); if( cc->GetShowEventStatus() == 0 ) cc->ToggleEventStatus();
cc->Divide(1, 9, 0);
cc->Divide(2,1);
cc->cd(1); for (Int_t i = 0; i < nClover; i++) {
hTDiff->Draw(); int canvasID = i + 1;
cc->cd(canvasID);
cc->cd(canvasID)->SetGrid();
cc->cd(canvasID)->SetTickx(2);
cc->cd(canvasID)->SetTicky(2);
cc->cd(canvasID)->SetBottomMargin(0.06);
cc->cd(canvasID)->SetLogy();
for( Int_t j = 0; j < nCrystalPerClover; j++){
int hID = nCrystalPerClover*i+ j ;
heCal[hID]->Draw("same");
//he[hID]->Draw("same");
}
}
cc->SetCrosshair(1);
cc->cd(2); TCanvas * c1 = new TCanvas("c1", "c1", 1000, 1000);
cc->cd(2)->SetGrid(); if( c1->GetShowEventStatus() == 0 ) c1->ToggleEventStatus();
hTDiffvEventID->SetMarkerStyle(3); c1->SetLogz();
hTDiffvEventID->Draw(""); c1->SetGridx();
heCalvID->SetNdivisions(-409, "X");
heCalvID->Draw("colz");
} }

12
Monitor_raw.h
View File

@ -19,15 +19,15 @@ public :
// Declaration of leaf types // Declaration of leaf types
Long64_t evID; Long64_t evID;
UShort_t detID; UShort_t ID;
UShort_t e; UShort_t e;
ULong64_t t; ULong64_t t;
// List of branches // List of branches
TBranch *b_data_ID; //! TBranch *b_data_ID; //!
TBranch *b_det_ID; //! TBranch *b_ID; //!
TBranch *b_energy; //! TBranch *b_energy; //!
TBranch *b_time_stamp; //! TBranch *b_timestamp; //!
Monitor_raw(TTree * /*tree*/ =0) : fChain(0) { } Monitor_raw(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~Monitor_raw() { } virtual ~Monitor_raw() { }
@ -59,9 +59,9 @@ void Monitor_raw::Init(TTree *tree)
fChain->SetMakeClass(1); fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_data_ID); fChain->SetBranchAddress("evID", &evID, &b_data_ID);
fChain->SetBranchAddress("detID", &detID, &b_det_ID); fChain->SetBranchAddress("id", &ID, &b_ID);
fChain->SetBranchAddress("e", &e, &b_energy); fChain->SetBranchAddress("e", &e, &b_energy);
fChain->SetBranchAddress("t", &t, &b_time_stamp); fChain->SetBranchAddress("t", &t, &b_timestamp);
} }
Bool_t Monitor_raw::Notify() Bool_t Monitor_raw::Notify()

7
armory/evt2hist.cpp
View File

@ -168,10 +168,11 @@ int main(int argn, char **argv) {
///} ///}
///data.Print(); ///data.Print();
int detID = map[data.id]; int detID = mapping[data.id];
if( 0 <= detID && detID < 100 ){ if( 0 <= detID && detID < 100 ){
if( corrFile != ""){ if( corrFile != ""){
double eCal = ApplyCorrection(corr, detID, data.energy); //double eCal = ApplyCorrection(corr, detID, data.energy);
double eCal = corr[detID][0] + corr[detID][1]*data.energy;
he[detID]->Fill(eCal); he[detID]->Fill(eCal);
}else{ }else{
he[detID]->Fill(data.energy); he[detID]->Fill(data.energy);
@ -230,7 +231,7 @@ int main(int argn, char **argv) {
fclose(inFile); fclose(inFile);
printf("\n============= reasched end of file\n"); printf("\n============= reasched end of file\n");
printf("Crtl+C to end program.\n"); printf("\nCrtl+C to end program.\n");
app->Run(); app->Run();

21
armory/evt2root.cpp
View File

@ -63,10 +63,9 @@ public:
//############################################# //#############################################
int main(int argn, char **argv) { int main(int argn, char **argv) {
if (argn != 2 && argn != 3 ) { if (argn != 2 ) {
printf("Usage :\n"); printf("Usage :\n");
printf("%s [evt File] [timeWindow] \n", argv[0]); printf("%s [evt File] \n", argv[0]);
printf(" timeWindow : number of tick, 1 tick = 10 ns. default = 100 \n");
return 1; return 1;
} }
@ -75,7 +74,6 @@ int main(int argn, char **argv) {
outFileName.Remove(inFileName.First('.')); outFileName.Remove(inFileName.First('.'));
outFileName.Append("_raw.root"); outFileName.Append("_raw.root");
long int inFilePos = 0;
TBenchmark gClock; TBenchmark gClock;
gClock.Reset(); gClock.Reset();
gClock.Start("timer"); gClock.Start("timer");
@ -96,6 +94,7 @@ int main(int argn, char **argv) {
fseek(inFile, 0L, SEEK_END); fseek(inFile, 0L, SEEK_END);
long int inFileSize = ftell(inFile); long int inFileSize = ftell(inFile);
rewind(inFile); ///back to the File begining rewind(inFile); ///back to the File begining
long int inFilePos = 0;
printf(" in file: %s\n", inFileName.Data()); printf(" in file: %s\n", inFileName.Data());
printf("out file: %s\n", outFileName.Data()); printf("out file: %s\n", outFileName.Data());
@ -106,21 +105,20 @@ int main(int argn, char **argv) {
TTree * tree = new TTree("tree", "tree"); TTree * tree = new TTree("tree", "tree");
tree->Branch("evID", &measureID, "data_ID/L"); tree->Branch("evID", &measureID, "data_ID/L");
tree->Branch("detID", &data.id, "det_ID/s"); tree->Branch("id", &data.id, "ID/s");
tree->Branch("e", &data.energy, "energy/s"); tree->Branch("e", &data.energy, "energy/s");
tree->Branch("t", &data.time, "time_stamp/l"); tree->Branch("t", &data.time, "timestamp/l");
//tree->Branch("tdiff", &data.timeDiff, "time_Diff/L"); //tree->Branch("tdiff", &data.timeDiff, "time_Diff/L");
//=======TODO online event building
unsigned int header[4]; //read 4 header, unsigned int = 4 byte = 32 bits. unsigned int header[4]; //read 4 header, unsigned int = 4 byte = 32 bits.
unsigned long long nWords = 0; unsigned long long nWords = 0;
ULong64_t timeLast = 0; //ULong64_t timeLast = 0;
//=============== Read File //=============== Read File
/// while ( ! feof(inFile) ){ /// while ( ! feof(inFile) ){
while ( inFilePos <= inFileSize ){ while( inFilePos < inFileSize || feof(inFile) ){
fread(header, sizeof(header), 1, inFile); fread(header, sizeof(header), 1, inFile);
inFilePos = ftell(inFile); inFilePos = ftell(inFile);
@ -166,13 +164,12 @@ int main(int argn, char **argv) {
//event stats, print status every 10000 events //event stats, print status every 10000 events
if ( measureID % 10000 == 0 ) { if ( measureID % 10000 == 0 ) {
inFilePos = ftell(inFile);
float tempf = (float)inFileSize/(1024.*1024.*1024.); float tempf = (float)inFileSize/(1024.*1024.*1024.);
gClock.Stop("timer"); gClock.Stop("timer");
double time = gClock.GetRealTime("timer"); double time = gClock.GetRealTime("timer");
gClock.Start("timer"); gClock.Start("timer");
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\033[A\r", 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\033[A\r",
measureID, (100*inFilePos/inFileSize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); measureID +1 , (100*inFilePos/inFileSize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
} }
//cern fill tree //cern fill tree
@ -187,7 +184,7 @@ int main(int argn, char **argv) {
gClock.Start("timer"); gClock.Start("timer");
float tempf = (float)inFileSize/(1024.*1024.*1024.); float tempf = (float)inFileSize/(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", 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",
measureID, (100*inFilePos/inFileSize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); measureID+1, (100*inFilePos/inFileSize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
fclose(inFile); fclose(inFile);

2
armory/pixie2root.cpp
View File

@ -307,7 +307,7 @@ int main(int argc, char **argv) {
///========== need a mapping, can reduce the array size, speed up. ///========== need a mapping, can reduce the array size, speed up.
int ch = map[subevt[sevtmult].id]; int ch = mapping[subevt[sevtmult].id];
if ( 0 <= ch && ch < NCRYSTAL ){ if ( 0 <= ch && ch < NCRYSTAL ){
energy[ch] = subevt[sevtmult].energy; energy[ch] = subevt[sevtmult].energy;
etimestamp[ch] = subevt[sevtmult].time; etimestamp[ch] = subevt[sevtmult].time;

2
mapping.h
View File

@ -12,7 +12,7 @@ Other : 200 - 299
#define NOTHER 52 #define NOTHER 52
// 0 1 2 3 4 5 6 7 8 9 // 0 1 2 3 4 5 6 7 8 9
int map[130] = { 0, 1, 2, 3, 100, 4, 5, 6, 7, 101, // 0 int mapping[130] = { 0, 1, 2, 3, 100, 4, 5, 6, 7, 101, // 0
8, 9, 10, 11, 102, -1, 12, 13, 14, 15, // 10 8, 9, 10, 11, 102, -1, 12, 13, 14, 15, // 10
103, 16, 17, 18, 19, 104, -1, -1, -1, -1, // 20 103, 16, 17, 18, 19, 104, -1, -1, -1, -1, // 20
-1, -1, 20, 21, 22, 23, 105, 24, 25, 26, // 30 -1, -1, 20, 21, 22, 23, 105, 24, 25, 26, // 30