#include #include #include #include #include #include #include "TSystem.h" #include "TObject.h" #include "TFile.h" #include "TTree.h" #include "TString.h" #include "TMath.h" #include "TGraph.h" #include "TLatex.h" #include "TBenchmark.h" #include "TH1F.h" #include "TApplication.h" #include "TCanvas.h" #include "TClonesArray.h" #include "../mapping.h" #include "../armory/AnalysisLibrary.h" #define MAX_CRATES 2 #define MAX_BOARDS_PER_CRATE 13 #define MAX_CHANNELS_PER_BOARD 16 #define BOARD_START 2 //#############################TODO // 1) multiple file // 2) Change to GUI // 4) eventBuilding class measurment{ 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; 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]; UShort_t id; Int_t detID; ULong64_t eventID; UShort_t trace[1024]; measurment(){ Clear(); }; ~measurment(){}; void Clear(){ ch = 0; slot = 0; crate = 0; headerLength = 0; eventLength = 0; pileup = false; time = 0; cfd = 0; energy = 0; trace_length = 0; trace_out_of_range = 0; id = 0; detID = -1; eventID = 0; trailing = 0; leading = 0; gap = 0; baseline = 0; for( int i = 0; i < 8; i++) QDCsum[i] = -1; for( int i = 0; i < 1024; i++) trace[i] = 0; } void Print(){ printf("============== eventID : %llu\n", eventID); printf("Crate: %d, Slot: %d, Ch: %d | id: %d = detID : %d \n", crate, slot, ch, id, detID); 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); 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); } printf(" QDCsum : \n"); for( int i = 0; i < 8; i++) printf(" %-10d\n", QDCsum[i]); } if( eventLength > headerLength ){ printf(" trace:\n"); for( int i = 0 ; i < trace_length ; i++)printf("%3d| %-10d\n",i, trace[i]); } } }; //############################################# // main //############################################### int main(int argn, char **argv) { if (argn < 2 || argn > 6 ) { printf("Usage :\n"); printf("%s [evt File] [E corr] [raw E threshold] [Save Hist] [Save Root]\n", argv[0]); printf(" [E corr] : correction file for gamma energy \n"); printf(" [raw E threshold] : min raw E \n"); printf(" [save Hist] : 1/0 \n"); printf(" [save Root] : 1/0 \n"); return 1; } TString inFileName = argv[1]; TString corrFile = ""; std::vector> eCorr; if( argn >= 3 ){ corrFile = argv[2]; eCorr.clear(); eCorr = LoadCorrectionParameters(corrFile); } int rawEnergyThreshold = 100; if( argn >= 4 ) rawEnergyThreshold = atoi(argv[3]); bool isSaveHist = false; ///save gamma hist for calibration if( argn >= 5 ) isSaveHist = atoi(argv[5]); bool isSaveRoot = false; ///save data into root if( argn >= 6 ) isSaveRoot = atoi(argv[6]); long int inFilePos; TBenchmark gClock; gClock.Reset(); gClock.Start("timer"); ULong64_t measureID = -1; measurment data; printf("====================================\n"); FILE * inFile = fopen(inFileName, "r"); if( inFile == NULL ){ printf("Cannot read file : %s \n", inFileName.Data()); return -404; } 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("--------------------------------\n"); TFile * fFile = NULL; TTree * tree = NULL; if( isSaveRoot ){ fFile = new TFile("temp.root", "RECREATE"); tree = new TTree("tree", "tree"); tree->Branch("detID", &data.detID, "detID/s"); tree->Branch("e", &data.energy, "energy/s"); tree->Branch("e_t", &data.time, "timestamp/l"); 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"); } //================ get file size fseek(inFile, 0L, SEEK_END); long int inFileSize = ftell(inFile); rewind(inFile); ///back to the File begining unsigned long long fpos = 0; //================ Historgrams TH1F * he[NCRYSTAL]; for( int i = 0 ; i < NCRYSTAL; i++){ he[i] = new TH1F(Form("he%02d", i), Form("e-%2d", i), 2000, 0, 2000); 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; } } 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(1, 9, 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); //=============== Read File unsigned int header[4]; //read 4 header, unsigned int = 4 byte = 32 bits. while ( ! feof(inFile) ){ fread(header, sizeof(header), 1, inFile); measureID ++; /// see the Pixie-16 user manual, Table4-2 data.eventID = measureID; 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 = header[2] >> 16 ; data.energy = (header[3] & 0xFFFF ) /2; // I don;t know why it has to "rebin" data.trace_length = (header[3] >> 16) & 0x7FFF; data.trace_out_of_range = header[3] >> 31; data.id = data.crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data.slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + data.ch; data.detID = mapping[data.id]; ///======== read QDCsum if( data.headerLength >= 4 ){ unsigned int extraHeader[data.headerLength-4]; fread(extraHeader, sizeof(extraHeader), 1, inFile); if( data.headerLength > 12){ data.trailing = extraHeader[0]; data.leading = extraHeader[1]; data.gap = extraHeader[2]; data.baseline = extraHeader[3]; } 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]; } } ///====== read trace if( data.eventLength > data.headerLength ){ unsigned int traceBlock[data.trace_length / 2]; fread(traceBlock, sizeof(traceBlock), 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 ; } } ///if( measureID < 10 ) { /// printf("----------------------event Length: %u, fpos: %llu byte (%lld words)\n", data.eventLength, fpos, fpos/4); /// for(int i = 0; i < 4; i++) printf(" %x\n", header[i]); /// data.Print(); ///} //=== jump to next measurement. obsolete, we read the whole block ///if( data.eventLength > 4 ) { /// if( fseek(inFile, sizeof(int) * (data.eventLength-data.headerLength), SEEK_CUR) != 0 ) break;; ///} fpos = ftell(inFile); //==== Fill Histogram if( 0 <= data.detID && data.detID < 100 && data.energy > rawEnergyThreshold ){ if( corrFile != ""){ ///========= apply correction int order = (int) eCorr[data.detID].size(); double eCal = 0; for( int k = 0; k < order ; k++){ eCal += eCorr[data.detID][k] * TMath::Power(data.energy, k); } he[data.detID]->Fill(eCal); }else{ he[data.detID]->Fill(data.energy); } } //===== 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( isSaveRoot ){ fFile->cd(); tree->Fill(); } //==== event stats, print status every 10000 events if ( measureID % 10000 == 0 ) { /// update file size, slow down? fseek(inFile, 0L, SEEK_END); inFileSize = ftell(inFile); fseek(inFile, fpos, SEEK_SET); inFilePos = ftell(inFile); 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", measureID, inFilePos/(1024.*1024.*1024.), inFileSize/1024./1024./1024, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); } //==== 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(); ///} 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(); inFilePos = ftell(inFile); 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", measureID, inFilePos/(1024.*1024.*1024.), inFileSize/1024./1024./1024, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); fclose(inFile); printf("\n\n\n============= reached end of file\n"); if( isSaveHist ) { printf(" save gamma histograms : \033[1;3mhist.root\033[m\n"); TFile * fHist = new TFile("hist.root", "RECREATE"); for( int i = 0; i < NCRYSTAL; i++) he[i]->Write("", TObject::kOverwrite); fHist->Close(); } if( isSaveRoot){ printf(" save into Root : \033[1;3mtemp.root\033[m\n"); fFile->cd(); tree->Write(); fFile->Close(); } printf("Crtl+C to end program.\n"); app->Run(); }