--- /dev/null
+#include <stdio.h>
+
+
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <glob.h>
+#include <algorithm>
+
+#include <TTree.h>
+#include <TFile.h>
+#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 <hld multifile base name>\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<sizeof(str)-1;strpos++)
+ {
+ ctoint = str[strpos]-'0';
+ pixelprob[row*1152+strpos] += ctoint;
+ if( !(ctoint==0 || ctoint==1) ) {printf("%s %i","Wrong Pixel Signal in 'string' =", ctoint); exit(0);}
+// printf("%i ",(int)pixelprob[row*1152+strpos]);
+ }
+// printf("--> %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;i<entries;i++)
+ {
+ scurveTree->GetEntry(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+1<threshold)
+ {
+ printf("Y ");
+ printf("%i ... ", threshold );
+ }
+
+ scurveTreeInd[bank]->Fill();
+ thresh_tmp = threshold;
+ }
+
+ printf("%i !\n", thresh_tmp);
+
+ printf("-------------\n");
+
+ for(int i=0;i<bank+1;i++)
+ {
+ scurveTreeInd[i]->Write("",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<bank+1;ba++)
+// {
+// paraTree[ba] = new TTree( Form("para%i",ba), Form("para%i",ba));
+// paraTree[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;thr<entries;thr++)
+// {
+// scurveTreeInd[ba]->GetEntry(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;thr<entries;thr++)
+// {
+// histos->Fill(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( mean<meanmin || mean>meanmax || 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( (mean<meanmin || mean>meanmax || 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);
+// ------------------------------------------
+
+
+
+
+}
+
jspc29 / mvd_soft.git / commitdiff