From 66958a6d2d111eb26bcf1152a433b83edb47ae13 Mon Sep 17 00:00:00 2001 From: Ryan Tang Date: Tue, 7 Dec 2021 18:33:21 -0500 Subject: [PATCH] basic of xia2root can do correction --- .gitignore | 6 + corr.txt | 64 +++++++ makefile | 10 ++ xia2ev2_nopart.cpp | 418 +++++++++++++++++++++++++++++++++++++++++++++ xia2root.cpp | 375 ++++++++++++++++++++++++++++++++++++++++ 5 files changed, 873 insertions(+) create mode 100644 .gitignore create mode 100644 corr.txt create mode 100644 makefile create mode 100644 xia2ev2_nopart.cpp create mode 100644 xia2root.cpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..04b07e5 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +*.root +*.evt +*.ev2 +*.log +xia2root +xia2ev2* diff --git a/corr.txt b/corr.txt new file mode 100644 index 0000000..1646e28 --- /dev/null +++ b/corr.txt @@ -0,0 +1,64 @@ +0 1 0 0 +0 2 0 0 +1000 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 +0 1 0 0 diff --git a/makefile b/makefile new file mode 100644 index 0000000..3e260d4 --- /dev/null +++ b/makefile @@ -0,0 +1,10 @@ +CC=g++ + +all: xia2root +#all: xia2root xia2ev2 + +xia2root: xia2root.cpp + $(CC) xia2root.cpp -o xia2root `root-config --cflags --glibs` + +xia2ev2: xia2ev2_nopart.cpp + $(CC) xia2ev2_nopart.cpp -o xia2ev2_nopart diff --git a/xia2ev2_nopart.cpp b/xia2ev2_nopart.cpp new file mode 100644 index 0000000..c148dd3 --- /dev/null +++ b/xia2ev2_nopart.cpp @@ -0,0 +1,418 @@ +#include +#include +#include +//for 64bit compiler: +#include +#include + +char outdata[70],outloga[70]; + +//int *outenergy=NULL; +int outsize=1; +//char *outcount=NULL; +//char *outdetnum=NULL; +char *outcount=new char[1]; +char *outdetnum=new char[1]; +//outenergy=new int[1]; +int *outenergy=new int[1]; +//outcount[0]=1; +int correctflag = 0,gateflag = 0,CntAll = 0,CntGate = 0; +long maxtime=100; //max number of 10 ns tics for an event +long maxeventime=0,mineventime=400000000,CntEvts[64]; +int multiplicity[10],startdet=15,PileUp[64]; +int evN=0; +int evStart = 0; +int evEnd = 1; +unsigned int *bufsam=NULL; +long *bufsiz=NULL; +long negativedelevtimecount=0,timezero=-30000; +int size1,size2,bsize,bsam=2,bread=1,chan,slot,chann,outevtcount=0; +int pu,energy,evtimehi,pileupcount=0,zerocount=0; +unsigned short cfd,cfdold; +int detectorcount=0,mul=0; +short detenergy[256],detid[256]; +unsigned long maskpu = 2147483648,evtimelo; +long long deltime=0LL,timeinevent=0LL; +unsigned long long evtime,multiplier=4294967296LL,oldtime=0LL; +long long zero=0LL; +int channel[65][100],timelo[65],timehi[65]; +float correctable[65][100]; +float fracx,correct,top,bottom,tdc,percent,averagetime=0; +int ipnt,iadc,itdc; +int pntcnt[65]; + +void InitializeVariables(char **args); +void WriteOutEvent(char **args,FILE *outdat); +void BuildEvent(); + +int main(int argn,char **argv) { + if ((argn<3) || (argn>8)) { + std::cerr<<"Usage:\n "<6) { + evStart=atoi(argv[6]); + maxtime=atoi(argv[4]); + startdet=atoi(argv[5]); + correctflag = 1; + if (evStart == 0) evEnd = 1; + else { + std::cerr<<"You must specify an end event number if the start event number is not 0. Your start event number is "< int for 64 bit + fread(bufsiz,sizeof(int),bread,infile); + //std::cout<=evStart) { + + if (bufsam[0] == 30) { + chan = (bufsam[2]) & (15); + slot = ((bufsam[2]) & (240))/16; + chann = (slot - 2)*16 + chan + 1; + pu = ((bufsam[2]) & (maskpu))/maskpu; + energy = ((bufsam[5]) & 65535); + evtimehi = ((bufsam[4]) & 65535); + evtimelo = bufsam[3]; + evtime = evtimelo + multiplier*evtimehi; + cfd = bufsam[4]/65536; + deltime = evtime - oldtime; + if(deltime < zero) negativedelevtimecount++; + if (energy==0) zerocount++; + if(pu>0)pileupcount++; + if ((pu > 0)&&(chann>0)&&(chann<65)) PileUp[chann-1]++; + + + // ignore pileups + timeinevent=timeinevent+deltime; + oldtime = evtime; + + if(timeinevent > maxtime) { + //This detector starts another event. First write out the previous event. + + WriteOutEvent(argv,outdat); + + // Now store this signal which starts a new event + timeinevent = 1; + timezero = -30000; + detectorcount = 0; + mul = 0; + averagetime = 0; + + BuildEvent(); + }else { + //This is just another detector in the current event + //Now we have to put the event just read into the buffer + + BuildEvent(); + + } //end if(timeinevent + + } //end if bufsam[0]=30 + + } //end if evN >= evStart + + } //end while evN < evEnd + + delete [] bufsiz; + delete [] bufsam; + + if(feof(infile) && evN 0) std::cout<<"PU("<0) fprintf(outlog,"Corrected time walk using file %s \n\n",argv[3]); + if (strcmp(argv[2],"NULL") != 0) fprintf(outlog,"wrote %i events to %s using xia2ev2_part with particle detector number %i\n\n",outevtcount,outdata,startdet); + + fprintf(outlog,"negativedelevtimecount = %i, pileupcount = %i zerocount = %i\n",negativedelevtimecount,pileupcount,zerocount); + + fprintf(outlog,"max time (10 ns) in event = %i min time betweeen events = %i\n\n",mineventime,maxtime); + + + for(int is=0;is<10;is=is+2) fprintf(outlog,"multiplicity %i: counts = %i; multiplicity %i: counts = %i\n",is+1,multiplicity[is],is+2,multiplicity[is+1]); + fprintf(outlog,"\n"); + + for(int is=0; is<64;is++) { + if (CntEvts[is] > 0) fprintf(outlog,"ADC(%i): Counts = %i Pileups = %i\n",is+1,CntEvts[is],PileUp[is]); + } + + percent = 0.; + if (CntAll > 0) percent = 100. * float(CntGate)/ float(CntAll); + fprintf(outlog,"Total coincidence events = %i Number in time gates = %i\n Percent in time gate = %5.2f\n",CntAll,CntGate,percent); + + + if (strcmp(argv[2],"NULL") != 0) fclose(outdat); + return 0; +} + + + +void InitializeVariables(char **args) { + + //initialize correctable + + for (ipnt = 0; ipnt<65; ipnt++) pntcnt[ipnt]=0; + for (int adc=0; adc<65; adc++){ + timelo[adc] = 0; timehi[adc]=32767; + for (ipnt=0; ipnt<100; ipnt++) { + correctable[adc][ipnt] = -1.; + channel[adc][ipnt]=-1; + } + } + + for (int is=0; is<10; is++) multiplicity[is]=0; + for (int is=0; is<64;is++) PileUp[is]=0; + for (int is=0; is<64; is++) CntEvts[is]=0; + + + if (correctflag > 0) { +// Read in correction table + std::cerr<<"correction table file = "<> iadc; + *in >> offset; + *in >> timelo[iadc]; + *in >> timehi[iadc]; + while (iadc >= 0) { + if (startdet == -2) std::cout<<"iadc ="<> channel[iadc][ipnt]; +// std::cout<> iadc; + *in >> offset; + *in >> timelo[iadc]; + *in >> timehi[iadc]; + + } //end if (iadc + ipnt=0; + } //end while (iadc + + + if (startdet == -2) return; + + } //end if (correctflag > 0 + +} //end void InitializaVariables + +void WriteOutEvent(char **args,FILE *outdat) { + + if(timeinevent < mineventime) mineventime = timeinevent; + if (mul>8) mul=9; + multiplicity[mul]++; + +if (strcmp(args[2],"NULL") != 0) { + detectorcount = detectorcount & 255; +// timezero will be > 0 if startdet fired in this event +// if(detectorcount > 0 && timezero > -30000) { + if(detectorcount > 0 && timezero == -30000) { + CntAll ++; +// printf("\n[%i] ",detectorcount); + for (int icnt=1; icnt 0) { //Is there a correct table for this ADC? + ipnt=0; + +// Linear interpolation between points in correct.tab + while (ipnt < pntcnt[iadc] && (detenergy[icnt-1] >= channel[iadc][ipnt])) ipnt++; + top = correctable[iadc][ipnt] - correctable[iadc][ipnt-1]; + bottom = channel[iadc][ipnt] - channel[iadc][ipnt-1]; + fracx = top/bottom; + correct = correctable[iadc][ipnt-1] + fracx*(detenergy[icnt] - channel[iadc][ipnt-1]); + tdc = detenergy[icnt]; + if (detid[icnt-1] != startdet) tdc = 200 - correct + tdc; + if (detid[icnt-1] == startdet) tdc = correct + tdc; + if (tdc <1) tdc=1; if (tdc > 32767) tdc=32767; + detenergy[icnt] = tdc + 0.5; + averagetime += tdc; + if (detid[icnt-1] == startdet) timezero = (long) (tdc+0.5); + } //end if (pntcnt + + else if (timezero == -30000) { + detenergy[icnt] = detenergy[icnt] + 99; + averagetime += detenergy[icnt]; + } + + + } //end for (int icnt= + +// Adjust for timezero + averagetime = 2*averagetime/detectorcount; + for (int icnt=1; icnt timehi[iadc])) gateflag = 1; + } + +// write out event if in gate + + if (gateflag == 0) { + CntGate ++; + outevtcount++; + outcount[0] = detectorcount; + fwrite(outcount,sizeof(char),outsize,outdat); + for (int icnt=0; icnt 0 && timezero > -30000) + } //endif strcmp +// timeinevent = 0; +// timezero=-30000; + +} //end void WriteOutEvent( + +void BuildEvent() { + + if(timeinevent > maxeventime) maxeventime = timeinevent; + mul++; +//first put energy in list +//suppress zero energies and those over range to not confuse gnuscope + if(energy>10 && energy < 32000) { + detid[detectorcount]=chann; + detenergy[detectorcount]=energy; + detectorcount++; + if (chann == startdet) timezero = timeinevent + 99; +//now put timeinevent in list as chann + 64 + detid[detectorcount]=chann+64; + detenergy[detectorcount]=timeinevent; + detectorcount++; + } // end if (energy>10 .. + +} // end BuildEvent diff --git a/xia2root.cpp b/xia2root.cpp new file mode 100644 index 0000000..20037fd --- /dev/null +++ b/xia2root.cpp @@ -0,0 +1,375 @@ +#include +#include +#include +#include +#include +#include "TFile.h" +#include "TTree.h" +#include "TString.h" +#include "TMath.h" +#include + +#define NUMDET 64 /// number of detector +#define STARTDETID 15 + +std::vector SplitStr(std::string tempLine, std::string splitter, int shift = 0){ + + std::vector output; + + size_t pos; + do{ + pos = tempLine.find(splitter); /// fine splitter + if( pos == 0 ){ ///check if it is splitter again + tempLine = tempLine.substr(pos+1); + continue; + } + + std::string secStr; + if( pos == std::string::npos ){ + secStr = tempLine; + }else{ + secStr = tempLine.substr(0, pos+shift); + tempLine = tempLine.substr(pos+shift); + } + + ///check if secStr is begin with space + while( secStr.substr(0, 1) == " "){ + secStr = secStr.substr(1); + }; + + ///check if secStr is end with space + while( secStr.back() == ' '){ + secStr = secStr.substr(0, secStr.size()-1); + } + + output.push_back(secStr); + //printf(" |%s---\n", secStr.c_str()); + + }while(pos != std::string::npos ); + + return output; +} + +std::vector> LoadCorrectionParameters(TString corrFile){ + + printf("==================== load correction parameters : %s", corrFile.Data()); + std::ifstream file; + file.open(corrFile.Data()); + + std::vector> corr; + corr.clear(); + + std::vector detCorr; + detCorr.clear(); + + if( file.is_open() ){ + + while( file.good() ){ + + std::string line; + getline(file, line); + + if( line.substr(0,1) == "#" ) continue; + if( line.substr(0,2) == "//" ) continue; + if( line.size() == 0 ) continue; + + std::vector temp = SplitStr(line, " "); + + detCorr.clear(); + for( int i = 0; i < (int) temp.size() ; i++){ + detCorr.push_back(std::stod(temp[i])); + } + corr.push_back(detCorr); + } + + file.close(); + + printf(".... done\n"); + 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("%6.2f, ", corr[i][j]); + } + printf("%6.2f\n", corr[i][len-1]); + } + + }else{ + printf(".... fail\n"); + } + + return corr; +} + + +//################################################################################### +//################ ########################### +//################ main ########################### +//################ ########################### +//################################################################################### +int main(int argn,char **argv) { + if ( argn == 1 ) { + printf("Usage: \n"); + printf("%s file_in.evt raw_Opt timeWidow correctionFile\n", argv[0]); + printf(" | | |\n"); + printf(" | | + correction file, row for det, col for order of correction\n"); + printf(" | |\n"); + printf(" | + when build event, event build window, 1 = 10 ns, default 100\n"); + printf(" + default 0 = raw, 1 = event build \n"); + /// std::cerr<<"Usage:\n "<= 3 ) rawOpt = atoi(argv[2]); + + int timeWindow = 100; + if ( argn >= 4 ) timeWindow = atoi(argv[3]); + + TString corrFileName = ""; + bool hasCorr = false; + if ( argn == 5 ) { + corrFileName = argv[4]; + hasCorr = true; + } + + FILE *infile=fopen(argv[1],"r"); + if (infile==NULL) { + printf("cannot open file : %s \n", argv[1]); + ///std::cerr<<"Problem opening "<> corr; + if( hasCorr ) corr = LoadCorrectionParameters(corrFileName); + + int chan,slot,chann; + int pu; /// pile up + int energy; + double cEnergy; + unsigned long long evtime; + unsigned short cfd; + + int pileupcount = 0; + int zerocount = 0; + int PileUp[64]; + + const unsigned long maskpu = 2147483648; + const unsigned long multiplier = 4294967296LL; + + double energyA[NUMDET]; + double cEnergyA[NUMDET]; + unsigned long long timeA[NUMDET]; + int puA[NUMDET]; + long long diffTimeA[NUMDET]; + unsigned short cfdA[NUMDET]; + int multi = 0; /// multipicilty in an event + int detMulti[NUMDET]; /// multiplicity in a detector in an event + + TFile * outFile = new TFile(outFileName, "RECREATE"); + outFile->cd(); + TTree * tree = new TTree("tree", "tree"); + + tree->Branch("eventID", &eventID, "event_number/l"); + + if ( rawOpt == 0 ){ /// when save raw data + tree->Branch("chan", &chan, "chan/I"); + tree->Branch("slot", &slot, "slot/I"); + tree->Branch("chann", &chann, "channel number/I"); + tree->Branch("pu", &pu, "pile-up/I"); + tree->Branch("energy", &energy, "energy/I"); + if( hasCorr) tree->Branch("cEnergy", &cEnergy, "corrected_energy/D"); + tree->Branch("time", &evtime, "timestamp/l"); + tree->Branch("cfd", &cfd, "cfd/s"); + }else{ /// when build event by time-window + tree->Branch("energy", energyA, Form("energy[%d]/D", NUMDET)); + if( hasCorr) tree->Branch("cEnergy", cEnergyA, Form("corrected_energy[%d]/D", NUMDET)); + tree->Branch("time", timeA, Form("timestamp[%d]/l", NUMDET)); + tree->Branch("dtime", diffTimeA, Form("diff_time[%d]/L", NUMDET)); + tree->Branch("pu", puA, Form("pile_up[%d]/I", NUMDET)); + tree->Branch("cfd", cfdA, Form("cfd[%d]/I", NUMDET)); + tree->Branch("multi", &multi, "multiplicity/I"); + tree->Branch("detMulti", detMulti, Form("det_multiplicity[%d]/I", NUMDET)); + } + + ///change this for 64bit compiler long *bufsam=NULL; + + //clear energy and time array + for( int i = 0; i < NUMDET; i++){ + energyA[i] = TMath::QuietNaN(); + cEnergyA[i] = TMath::QuietNaN(); + timeA[i] = 0; + diffTimeA[i] = -999; + cfdA[i] = 0; + puA[i] = -1; + detMulti[i] = 0; + } + multi = 0; + + unsigned long long startTime = 0; + long long diffTime = 0; + + int bread = 1; + int bsam = 2; + long * bufsiz=new long[bsam]; + unsigned int *bufsam = NULL; + + printf("============ Start looping events | build event ? %s", rawOpt == 0 ? "No" : "Yes"); + if( rawOpt == 1 ) { + printf(" time window : +- %d click\n", timeWindow); + }else{ + printf("\n"); + } + while ( !feof(infile) ) { + + // get buffer size + ///change long -> int for 64 bit + fread(bufsiz,sizeof(int),bread,infile); /// read 1 int (4 byte) from infile and save to bufsize + int bsize = bufsiz[0] -4 ; + if (feof(infile)) break; + blockNum ++; + + ///change for 64bit bufsam=new long[bsize/4+1]; + bufsam = new unsigned int[bsize/4+1]; + fread((char*)bufsam, 1, bsize, infile); /// read bsize of 1 byte from infile and save to char + ///printf("============ bsize : %d \n", bsize); + ///for( int i = 0; i < bsize; i++) printf("%d, ", bufsam[i]); + ///printf("\n"); + + if (bufsam[0] == 30) { + block30Num ++; + chan = (bufsam[2]) & (15); + slot = ((bufsam[2]) & (240))/16; + chann = (slot - 2)*16 + chan + 1; + pu = ((bufsam[2]) & (maskpu))/maskpu; + energy = ((bufsam[5]) & 65535); + unsigned long long evtimehi = ((bufsam[4]) & 65535); + unsigned long long evtimelo = bufsam[3]; + evtime = evtimelo + multiplier*evtimehi; + cfd = bufsam[4]/65536; + + if ( energy == 0 ) zerocount++; + if ( pu > 0 ) pileupcount++; + if ((pu > 0 ) && ( chann > 0 ) && ( chann < 65 )) PileUp[chann-1]++; + + if( blockNum % 100000 == 0 ) printf("."); + ///if( blockNum % 100000 == 0 ) printf("%llu \n", blockNum); + ///if( block30Num < 50) printf("b30: %10llu, chan: %d, slot: %d, chann: %2d, pu: %2d, energy: %5d, evtime: %13llu, cfd: %d\n", block30Num, chan, slot, chann, pu, energy, evtime, cfd); + + /// energy correction + if ( hasCorr ){ + cEnergy = 0; + int order = (int) corr[chann-1].size(); + for( int i = 0; i < order ; i++){ + cEnergy += corr[chann-1][i] * TMath::Power((double)energy, i); + } + } + + if ( rawOpt == 0 ) { + eventID++; + outFile->cd(); + tree->Fill(); + }else{ /// build event + + if ( startTime == 0 ) startTime = evtime; + + diffTime = evtime - startTime; + + if( -timeWindow < diffTime && diffTime < timeWindow ){ + + if( !TMath::IsNaN(energyA[chann-1]) ) detMulti[chann-1] ++; + energyA[chann-1] = energy; + cEnergyA[chann-1] = cEnergy; + timeA[chann-1] = evtime; + diffTimeA[chann-1] = diffTime; + puA[chann-1] = pu; + detMulti[chann-1]++; + multi++; + + + }else{ + /// fill tree + eventID++; + outFile->cd(); + tree->Fill(); + /// clear energy and time array + multi = 0; + for( int i = 0; i < NUMDET; i++){ + energyA[i] = TMath::QuietNaN(); + cEnergyA[i] = TMath::QuietNaN(); + timeA[i] = 0; + diffTimeA[i] = -999; + puA[i] = -1; + detMulti[i] = 0; + cfdA[i] = 0; + } + + /// fill the 1st data of a new event + startTime = evtime; + energyA[chann-1] = energy; + cEnergyA[chann-1] = cEnergy; + timeA[chann-1] = evtime; + diffTimeA[chann-1] = 0; + puA[chann-1] = pu; + detMulti[chann-1]++; + multi++; + } + } + + } ///end if bufsam[0]=30 + + } + + delete [] bufsiz; + delete [] bufsam; + fclose(infile); + + printf("\n============ end of event loop, totoal block read: %llu \n", blockNum); + + eventID++; + outFile->cd(); + tree->Write(); + outFile->Close(); + + + //========================= Print summary + printf("============================================\n"); + ///printf(" number of block: %llu\n", blockNum); + printf(" number of type 30 block: %llu\n", block30Num); + printf(" event built: %llu\n", eventID); + printf("============================================\n"); + + return 0; +} +