From 997a5eac6ea6e57c8c3ed59acdecfd642ce51239 Mon Sep 17 00:00:00 2001 From: Philipp Klaus Date: Tue, 19 Aug 2014 16:03:37 +0200 Subject: [PATCH] s-curves: Faster version of process_hld_rootmulti.c (w/o fitting of the s-curves) --- s-curves/Makefile | 2 + s-curves/process_hld_rootmulti_short.c | 824 +++++++++++++++++++++++++ 2 files changed, 826 insertions(+) create mode 100644 s-curves/process_hld_rootmulti_short.c diff --git a/s-curves/Makefile b/s-curves/Makefile index 846a34a..ec9fb7d 100644 --- a/s-curves/Makefile +++ b/s-curves/Makefile @@ -5,7 +5,9 @@ CLDFLAGS = -std=c++11 `root-config --cflags --glibs` -lboost_system -lboost_file all: process_hld_rootmulti.c process_hld_rootmulti_short.c $(CC) $(CLDFLAGS) process_hld_rootmulti.c -o process_hld_rootmulti + $(CC) $(CLDFLAGS) process_hld_rootmulti_short.c -o process_hld_rootmulti_short clean: rm -f process_hld_rootmulti + rm -f process_hld_rootmulti_short diff --git a/s-curves/process_hld_rootmulti_short.c b/s-curves/process_hld_rootmulti_short.c new file mode 100644 index 0000000..aff41f5 --- /dev/null +++ b/s-curves/process_hld_rootmulti_short.c @@ -0,0 +1,824 @@ +#include + + +#include +#include +#include +#include +#include + +#include +#include +#include "TROOT.h" +#include "TObject.h" +#include "TBrowser.h" +#include "TH1F.h" +#include "TMath.h" +#include "TF1.h" +#include "TCanvas.h" +#include "TRandom3.h" +#include "TApplication.h" + + +void usage(void) { + printf("\nUsage:"); + printf(" process_hld_root \n"); + printf("\n"); + exit (0); +} + +char value[33] = ""; +char *htoi(const char *ptr) { + char ch = *ptr; + int i; + const char* quads[] = {"0000", "0001", "0010", "0011", "0100", "0101", + "0110", "0111", "1000", "1001", "1010", "1011", + "1100", "1101", "1110", "1111"}; + + memset(value, '\0', sizeof(value)); + + while (ch == ' ' || ch == '\t') + ch = *(++ptr); + + for (i = 0; i < 8; i++) { + if (ch >= '0' && ch <= '9') + strncat(value, quads[ch - '0'], 4); + if (ch >= 'A' && ch <= 'F') + strncat(value, quads[10 + ch - 'A'], 4); + if (ch >= 'a' && ch <= 'f') + strncat(value, quads[10 + ch - 'a'], 4); + + ch = *(++ptr); + } + return value; +} + + + +// char output[256] = "/local.1/daq/workdir/scurve/output.txt\0"; + +char output[256] = "/local.1/daq/workdir/scurve/output.txt\0"; + +FILE *fo; +int mode = 0; +int an_mode = 0; +char timestamp[9]; +char external[9]; +char hfnum[9]; +int fnum; +int row; +char hthresh[5]; +char hrun[5]; +int thresh; +int run; +int prev_thresh = 999; +int prev_run = 999; +unsigned char byte; +int part = -1; +int seq = 0; +int ctr; +int pctr; +unsigned char word[5]; +char hword[9]; +int i; +unsigned char rocsize[3]; +char hrocsize[5]; +int size; +int datacounter; +int counts; +char bitword[33]; +char d0bitword[17]; +char d1bitword[17]; +int j; +char str[1153]; +char d0string[577]; +char d1string[577]; + +uint ctoint; +Float_t* pixelprob; +Float_t* pixelprob_arr; +UShort_t* rowcount; +UShort_t threshold; +TTree* scurveTree; +TTree* paraTree[10]; +TString outFileName; +TFile* outFile; +Int_t rowcurr; +bool firstout; +bool filled; +UInt_t pixel; +Float_t mean, sigma, chi2; +Int_t entries; +TCanvas* canv; +TH1F* histos = 0; + +void analyze(void){ + //printf("got: %s %d %d\n",hword,mode, an_mode); + + if (mode == 1){ + if (an_mode == 0){ + counts = 0; + an_mode++; + } + else{ + printf("Error: an_mode not 0 in mode 1!\n"); + exit(0); + } + } + else if (mode == 2){ + switch(an_mode){ + case 0: + break; + case 1: + if (strcmp(hword,"ffffffff") == 0){ + an_mode++; + } + else{ + printf("Error: ffffffff not found (%s)!\n",hword); + exit(0); + } + break; + case 2: + if (strcmp(hword,"f0020001") == 0){ + an_mode++; + } + else{ + printf("Error: f0020001 not found (%s)! Wrong Version or Testmode not activated.\n",hword); + exit(0); + } + break; + case 3: +// if (strcmp(hword,"aaa0aaa0") == 0){ +// an_mode++; +// } +// else{ +// printf("Error: aaa0aaa0 not found (%s)! Wrong ID.\n",hword); +// exit(0); +// } + an_mode++; + break; + case 4: + an_mode++; // ignore status + break; + case 5: + an_mode++; // reserved + break; + case 6: + an_mode++; // reserved + break; + case 7: + memset(external, '\0', sizeof(external)); + strcpy(external,hword); + //printf("Final copied string : %s\n", external); + //printf("From : %s\n", hword); + hthresh[0] = hword[0]; + hthresh[1] = hword[1]; + hthresh[2] = hword[2]; + hthresh[3] = hword[3]; + hthresh[4] = '\0'; + hrun[0] = hword[4]; + hrun[1] = hword[5]; + hrun[2] = hword[6]; + hrun[3] = hword[7]; + hrun[4] = '\0'; + thresh = (int)strtol(hthresh,NULL,16); + run = (int)strtol(hrun,NULL,16); + if (thresh != prev_thresh){ + fprintf(fo,">>>Threshold %d\n", thresh); // --> 1 <-- + +//---------------------------------- + if(firstout) + { + for(int i=0;i<1152*576;i++) + { + rowcurr = (int)(i/1152); + if( rowcount[rowcurr] != 0 ) + { +// printf("%5i %5i %10.6f %10.1i ",i, threshold, pixelprob[i], rowcount[rowcurr]); + pixelprob[i] = 1.*pixelprob[i]/rowcount[rowcurr]; +// printf("%10.6f \n",pixelprob[i]); + } + else + { + pixelprob[i] = 0; + } + } + + scurveTree->Fill(); + + for(int i=0;i<1152*576;i++) {pixelprob[i]=0;} + for(int i=0;i<576;i++) {rowcount[i]=0;} + } + else + { + firstout = true; + } + +// printf("Entries: %i \n", scurveTree->GetEntries()); +//---------------------------------- + + printf("%d ",thresh); + fflush(stdout); + prev_thresh = thresh; + } + if (run != prev_run){ + fprintf(fo," >>>Run %d\n", run); // --> 2 <-- + prev_run = run; + } + an_mode++; + break; + case 8: + memset(hfnum, '\0', sizeof(hfnum)); + strcpy(hfnum,hword); + //printf("Final copied string : %s\n", hfnum); + //printf("From : %s\n", hword); + fnum = (int)strtol(hfnum,NULL,16); + //printf("fnum : %d\n", fnum); + if (run==0){ + an_mode = 0; + } + else{ + memset(str, '\0', sizeof(str)); + memset(d0string, '\0', sizeof(d0string)); + memset(d1string, '\0', sizeof(d1string)); + an_mode++; + } + break; + case 9: + counts++; + if (fnum == 1){ + an_mode = 0; + } + else{ + memset(bitword, '\0', sizeof(bitword)); + std::reverse( hword, &hword[ strlen( hword ) ] ); + strcpy(bitword,htoi(hword)); + strcat(str,bitword); + + if (counts == 36){ + row = fnum-1; + if (row > 572){ + row++; + } + fprintf(fo," >>>row %d external %s framenum %s %s\n",row,external,hfnum,str); // --> 3 <-- +//---------------------------------- + threshold = thresh; + if(row>=576) {printf("%s %i","Wrong Row Number in 'row' =", row); exit(0);} + rowcount[row]++; +// printf("%i %i\n", threshold, row); + for(int strpos=0;strpos %i %i\n",sizeof(string), row); +//---------------------------------- + an_mode = 0; + } + } + break; + default : + printf("Error: an_mode not 1-9 in mode 2!"); + exit(0); + + } + + } + else if (mode == 3){ + an_mode = 0; + } + else{ + printf("Error: mode not 1-3!\n"); + exit(0); + } + + +} + +int main(int argc, char *argv[]) { + printf("-------------\n"); + char* filename; + + if (argc == 2){ + filename = &argv[1][0]; + printf("Using hld file %s\n",filename); + printf("The output %s is NOT used!\n",output); + } + else{ + usage(); + } + + printf("Starting HLD (multi) file analysis!\n"); + + strcat(filename,"*\0"); + glob_t globbuf; + glob(filename, 0, NULL, &globbuf); + size_t num; + FILE *fp; + + fo = fopen(output,"w"); +//---------------------------------- + TApplication theApp("Analysis", &argc, argv); + + + firstout = false; + pixelprob = new Float_t[1152*576]; for(int i=0;i<1152*576;i++) {pixelprob[i]=0;} + rowcount = new UShort_t[576]; for(int i=0;i<576;i++) {rowcount[i]=0;} + + outFileName = "output.root\0"; + outFile = new TFile (outFileName,"RECREATE"); + + scurveTree = new TTree("scurvesall", "scurvesall"); + scurveTree -> Branch("threshold" , &threshold , "threshold/s" , 32000); + scurveTree -> Branch("pixelprob" , pixelprob , "pixelprob[663552]/F" , 32000); + scurveTree -> Branch("frames" , rowcount , "frames[576]/s" , 32000); + +//---------------------------------- + + + for (num = 0; num < globbuf.gl_pathc; num++) + { + char *file = globbuf.gl_pathv[num]; + printf("\n\n --> Using file %s ",file); + if (!(fp = fopen(file, "r"))){ + printf("File not found!\n"); + exit(0); + } + + printf("\nThreshold: "); + + part = -1; + seq = 0; + while ( fread(&byte,1,1,fp) > 0){ + //printf("part: %d\n",part); + switch (part){ + case -1: + seq++; + if (seq == 32) { + seq = 0; + part = 0; + } + break; + case 0: + seq++; + if (seq == 4) { + seq = 0; + part = 1; + } + break; + case 1: + seq++; + if (seq == 4) { + seq = 0; + part = 2; + } + break; + case 2: + seq++; + if (seq == 4) { + seq = 0; + part = 3; + } + break; + case 3: + seq++; + if (seq == 4) { + seq = 0; + part = 4; + } + break; + case 4: + seq++; + if (seq == 4) { + seq = 0; + part = 5; + } + break; + case 5: + seq++; + if (seq == 4) { + seq = 0; + part = 6; + } + break; + case 6: + seq++; + if (seq == 4) { + seq = 0; + part = 7; + } + break; + case 7: + seq++; + if (seq == 4) { + seq = 0; + part = 8; + } + break; + case 8: + seq++; + if (seq == 4) { + ctr = ctr + 1; + ctr = ctr % 2; + seq = 0; + part = 9; + } + break; + case 9: + seq++; + if (seq == 4) { + ctr = ctr + 1; + ctr = ctr % 2; + seq = 0; + part = 10; + } + break; + case 10: + seq++; + if (seq == 4) { + ctr = ctr + 1; + ctr = ctr % 2; + seq = 0; + part = 11; + } + break; + case 11: + seq++; + if (seq == 4) { + ctr = ctr + 1; + ctr = ctr % 2; + seq = 0; + part = 12; + } + break; + case 12: // ROC Header + word[seq] = byte; + seq++; + //printf("byte: %02x",byte); + if (seq == 4) { + word[4] = '\0'; + ctr = ctr + 1; + ctr = ctr % 2; + for (i=0;i<4;i++){ + sprintf(&hword[i*2],"%02x", word[i]); + } + hword[8] = '\0'; + //printf("hword: %s\n",hword); + if (strcmp(hword,"002cd882") != 0){ + if (strcmp(hword,"002cd000") != 0){ + printf("\nWarning: ROC data header neither 002cd882 nor 002cd000 (%s)!\n",hword); + } + //exit(0); + } + + mode = 1; // start mode + analyze(); + rocsize[0] = word[0]; + rocsize[1] = word[1]; + rocsize[2] = '\0'; + for (i=0;i<2;i++){ + sprintf(&hrocsize[i*2],"%02x", rocsize[i]); + } + hrocsize[4] = '\0'; + size = (int)strtol(hrocsize,NULL,16); + //printf("size: %d\n",size); + datacounter = size; + if (size % 2 == 1){ + pctr = 1; + } + else{ + pctr = 0; + } + seq = 0; + part = 13; + + } + break; + case 13: + word[seq] = byte; + seq++; + if (seq == 4) { + word[4] = '\0'; + ctr = ctr + 1; + ctr = ctr % 2; + for (i=0;i<4;i++){ + sprintf(&hword[i*2],"%02x", word[i]); + } + hword[8] = '\0'; + //printf("dword: %s\n",hword); + mode = 2; // data mode + analyze(); + datacounter--; + seq = 0; + part = 13; + if (datacounter == 0){ + if (pctr == 0){ + part = 14; + } + else{ + part = 15; + } + } + } + break; + case 14: + seq++; + if (seq == 8) { + //printf("end\n"); + mode = 3; // end of data block + analyze(); + seq = 0; + if (ctr == 1){ + ctr = 0; + part = 100; + } + else{ + part = 0; + } + } + break; + case 15: + seq++; + if (seq == 16) { + //printf("end\n"); + mode = 3; // end of data block + analyze(); + seq = 0; + if (ctr == 1){ + ctr = 0; + part = 100; + } + else{ + part = 0; + } + } + break; + + case 100: + seq++; + if (seq == 4) { + seq = 0; + part = 0; + } + break; + + } + + } + } +// ------------------------------------------ + for(int i=0;i<1152*576;i++) + { + rowcurr = (int)(i/1152); + if( rowcount[rowcurr] != 0 ) + { + pixelprob[i] = 1.*pixelprob[i]/rowcount[rowcurr]; + } + else + { + pixelprob[i] = 0; + } + } + scurveTree->Fill(); + for(int i=0;i<1152*576;i++) {pixelprob[i]=0;} + for(int i=0;i<576;i++) {rowcount[i]=0;} + + +// Save +// printf("\n%s","-------------\n"); +// printf("Write data: %i Thresholds\n", entries); + scurveTree ->Write("",TObject::kOverwrite); + outFile->Write(); +// outFile->Save(); +// printf("\n%s","-------------\n"); +// printf("%s","-------------\n"); +// printf("%s created\n",Form("%s",outFileName.Data())); +// printf("%s","-------------\n"); + + entries = (int) scurveTree->GetEntries(); + + printf("\n-------------\n"); + printf("Entries: %i\n",entries); + printf("Seperate runs for different banks!\n"); + + scurveTree->SetBranchAddress( "threshold" , &threshold ); + scurveTree->SetBranchAddress( "pixelprob" , pixelprob ); + scurveTree->SetBranchAddress( "frames" , rowcount ); + + TTree *scurveTreeInd[10]; + + int thresh_tmp=255; + int bank=-1; + + for(int i=0;iGetEntry(i); + + if ( thresh_tmp>threshold ) + { + printf("X "); + if( bank!=-1 ) { printf("%i !\n", thresh_tmp); } + bank++; + scurveTreeInd[bank] = new TTree( Form("scurves%i",bank), Form("scurves%i",bank) ); + scurveTreeInd[bank] -> Branch("threshold" , &threshold , "threshold/s" , 32000); + scurveTreeInd[bank] -> Branch("pixelprob" , pixelprob , "pixelprob[663552]/F" , 32000); + scurveTreeInd[bank] -> Branch("frames" , rowcount , "frames[576]/s" , 32000); + + printf("Entries for bank %i found! Process threshold %i ... ", bank, threshold ); + } + else if ( thresh_tmp+1Fill(); + thresh_tmp = threshold; + } + + printf("%i !\n", thresh_tmp); + + printf("-------------\n"); + + for(int i=0;iWrite("",TObject::kOverwrite); + outFile->Write(); +// outFile->Save(); + printf("Data written for bank %i\n",i); + } + printf("-------------\n"); + + pixelprob_arr = new Float_t[entries*1152*576]; + //canv = new TCanvas("canv","canv",600,400); + //canv->cd(); + + TH1F *histos, *histo[3]; + histos = new TH1F( "test", "test" , 255,0.,255.); + histo[0] = new TH1F( "mean" , "mean" , 2550,0.,255.); + histo[1] = new TH1F( "sigma" , "sigma" , 1000,0.,100.); + histo[2] = new TH1F( "chi2" , "chi2" , 1000,0.,1.); + + TF1* erf = new TF1("erf","0.5*(1+TMath::Erf((-x+[0])/[1]))",0,100); + + int chicount; + int meanmin, meanmax; + float meanmean; + float vsigma; + float vtmp=1; + +// printf("Perform fitting of S-Curces:\n"); +// +// for(int ba=0;ba Branch("pixel" , &pixel , "pixel/i" , 32000); +// paraTree[ba] -> Branch("mean" , &mean , "mean/F" , 32000); +// paraTree[ba] -> Branch("sigma" , &sigma , "sigma/F" , 32000); +// paraTree[ba] -> Branch("chi2" , &chi2 , "chi2/F" , 32000); +// +// entries = (int) scurveTreeInd[ba]->GetEntries(); +// +// for(int thr=0;thrGetEntry(thr); +// +// for(int pix=0;pix<1152*576;pix++) +// { +// pixelprob_arr[pix*entries+thr] = pixelprob[pix]; +// } +// } +// +// for(int pix=0;pix<1152*576;pix++) +// { +// if(pix%100==0) { +// //printf("\r Bank %i - %6.1f %%", ba, 100.*pix/(1152*576)); +// fflush(stdout); +// } +// if( TMath::RMS(entries,&pixelprob_arr[pix*entries])!=0 ) +// { +// delete histos; +// histos = new TH1F( "test", "test" , 255,0.,255.); +// +// meanmin=0; +// meanmax=255; +// +// for(int thr=0;thrFill(thr,pixelprob_arr[pix*entries+thr]); +// +// if( pixelprob_arr[pix*entries+thr]!=1 && meanmin==0 ) { meanmin=thr-1; } +// if( pixelprob_arr[pix*entries+thr]==0 && vtmp!=0 ) { meanmax=thr; } +// vtmp = pixelprob_arr[pix*entries+thr]; +// } +// +// vsigma = (meanmax-meanmin)/2.; +// meanmean = (meanmax+meanmin)/2.; +// erf->SetParameters(meanmean,vsigma); +// histos->Fit(erf,"Q"); +// +// pixel = pix; +// mean = erf->GetParameter(0); +// sigma = erf->GetParameter(1)/TMath::Sqrt(2); +// chi2 = erf->GetChisquare()/erf->GetNDF(); +// +// if( meanmeanmax || sigma>75 || sigma<1 ) +// { +// chicount=0; +// //printf("\r----- \n"); +// //printf("%7i (%10i) %10.2f (%10i) %10.2f (%10f) %10.3f\n", pixel, meanmin, mean,meanmax, sigma,vsigma, chi2); +// +// while( (meanmeanmax || sigma>75 || sigma<1) && chicount<10) +// { +// erf->SetParameters(meanmean,1.*(chicount+1)); +// histos->Fit(erf,"QM"); +// pixel = pix; +// mean = erf->GetParameter(0); +// sigma = erf->GetParameter(1)/TMath::Sqrt(2); +// chi2 = erf->GetChisquare()/erf->GetNDF(); +// chicount++; +// +// //printf("%7i (%10i) %10.2f (%10i) %10.2f (%10f) %10.3f\n", pixel, meanmin, mean,meanmax, sigma,vsigma, chi2); +// } +// //histos->Draw(); +// //canv->Update(); +// } +// +// histo[0]->Fill(mean); +// histo[1]->Fill(sigma); +// histo[2]->Fill(chi2); +// +// if(sigma>=75 || chi2>=100) +// { +// //printf("\r----- \n"); +// //printf("\r--> "); +// //printf("%7i (%10i) %10.2f (%10i) %10.2f (%10f) %10.3f\n", pixel, meanmin, mean,meanmax, sigma,vsigma, chi2); +// //histos->Draw(); +// //canv->Update(); +// } +// +// if( sigma<1 ) +// { +// //printf("\r----- \n"); +// //printf("\r--> "); +// //printf("%7i (%10i) %10.2f (%10i) %10.2f (%10f) %10.3f\n", pixel, meanmin, mean,meanmax, sigma,vsigma, chi2); +// //histos->Draw(); +// //canv->Update(); +// } +// +// paraTree[ba]->Fill(); +// +// // if(1) +// // { +// // canv->Update(); +// // } +// +// // printf("\r%9i %5i %5i (%i) : %6.2f \u00b1 %5.2f %6.2f \u00b1 %5.2f %6.2f\n", pix, (int)(pix/1152),pix%1152, int((pix%1152)/288), erf->GetParameter(0), erf->GetParError(0), erf->GetParameter(1)/TMath::Sqrt(2), erf->GetParError(1)/TMath::Sqrt(2), erf->GetChisquare()); +// } +// } +// +// paraTree[ba]->Write("",TObject::kOverwrite); +// outFile->Write(); +// // outFile->Save(); +// } +// +// canv->Divide(3,1); +// canv->cd(1); +// histo[0]->Draw(); +// canv->cd(2); +// histo[1]->Draw(); +// canv->cd(3); +// histo[2]->Draw(); +// +// canv->Update(); + + outFile->Write(); + outFile->Save(); + outFile->Close(); + + delete[] pixelprob; + delete[] rowcount; + delete[] pixelprob_arr; + +// scurveTree->Delete(); +// paraTree->Delete(); + +// delete scurveTree; +// delete paraTree; +// gDirectory->Delete("scurvesall"); + + printf("-------------\n"); + printf("\n%s created\n",Form("%s",outFileName.Data())); + printf("Finished!\n"); + printf("-------------\n"); + gApplication->Terminate(); + //theApp.Run(kTRUE); +// ------------------------------------------ + + + + +} + -- 2.43.0