--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.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 "TGraph.h"
+
+//####################################################################
+
+void rewrite() {
+ cout<<"-------------"<<endl;
+ cout<<"REWRITE FILE:"<<endl;
+ UShort_t threshold;
+ Int_t* pixel = new Int_t [1152*576];
+ Float_t* pixelprob = new Float_t [1152*576];
+ Float_t* pixelprobout = new Float_t [1152*576];
+ UShort_t* frames = new UShort_t [576];
+
+ for(int i=0;i<1152*576;i++) { pixel[i]=i; }
+
+// File IN
+ TString RootFile = "output.root";
+ TFile* f = new TFile(RootFile);
+ TTree* tree = (TTree*)f->Get( "scurvesall" );
+
+ tree -> SetBranchAddress ("threshold" , &threshold );
+ tree -> SetBranchAddress ("pixelprob" , pixelprob );
+ tree -> SetBranchAddress ("frames" , frames );
+
+// File OUT
+ TString RootFileOut = "outputsc.root";
+ TFile* fout = new TFile(RootFileOut,"RECREATE");
+ TTree* treeout = new TTree( "scurvesall" , "scurvesall" );
+
+ treeout -> Branch("threshold" , &threshold , "threshold/s" , 32000);
+// treeout -> Branch("pixel" , pixel , "pixel[663552]/I" , 32000);
+ treeout -> Branch("pixelprob" , pixelprobout , "pixelprob[663552]/F" , 32000);
+// --------------
+ Int_t nthr = tree->GetEntries();
+ if (nthr % 4 != 0){
+ printf("Error: %d not divisible by 4!", nthr);
+ exit(1);
+ }
+ nthr = nthr/4;
+
+ for(int i=0;i<nthr;i++)
+ {
+ for(int j=0;j<4;j++)
+ {
+ tree->GetEntry(i+j*nthr);
+
+ for(int pix=0;pix<1152*576;pix++)
+ {
+ if( pix%1152>=1152/4*(j) && pix%1152<1152/4*(j+1) )
+ {
+ pixelprobout[pix] = pixelprob[pix];
+ }
+ }
+ }
+
+ cout<<"\rThreshold: "<<threshold<<" ( " <<i+1<<" of "<<nthr<<" ) ";
+ cout.flush();
+ treeout->Fill();
+ }
+ cout<<"\r ...done! "<<endl;
+ cout<<"-------------"<<endl;
+ printf("\r Save File!\n");
+ treeout -> Write("",TObject::kOverwrite);
+ fout->Save();
+ fout->Close();
+ f->Close();
+ cout<<"-------------"<<endl;
+}
+
+//####################################################################
+
+void plotsc() {
+
+ UShort_t threshold;
+// Int_t* pixel = new Int_t [1152*576];
+ Float_t* pixelprob = new Float_t [1152*576];
+
+// File IN
+ TString RootFile = "outputsc.root";
+ TFile* f = new TFile(RootFile);
+ TTree* tree = (TTree*)f->Get( "scurvesall" );
+
+ tree -> SetBranchAddress ("threshold" , &threshold );
+// tree -> SetBranchAddress ("pixel" , pixel );
+ tree -> SetBranchAddress ("pixelprob" , pixelprob );
+// --------------
+ Int_t nthr = tree->GetEntries();
+
+ TH1F* plot[1152*576];
+ for(int i=0;i<1152*576;i++)
+ {
+ plot[i] = new TH1F( Form("plot%i",i), Form("plot%i",i), nthr,0,255);
+ }
+
+ for(int i=0;i<nthr;i++)
+ {
+ tree->GetEntry(i);
+
+ for(int pix=0;pix<1152*576;pix++)
+ {
+ plot[pix] ->Fill(threshold,pixelprob[pix]);
+
+ if(pix==1152) { cout<<threshold<<"\t"<<pix<<"\t"<<pixelprob[1152]<<"\t"<<endl; }
+ }
+ }
+
+ TCanvas* cm2 = new TCanvas("a","<", 50,100,900,600);
+ cm2->Divide(4,1);
+
+ int row, column;
+ int matrix;
+ bool FIRST[4]= {0,0,0,0};
+
+ for(int pix=0;pix<1152*576;pix++)
+ {
+
+ column = pix%1152;
+ row = (int)(pix/1152);
+
+ matrix = int(column/288);
+
+ if(1) // Choose your conditions for which pixel you want to plot the scurves
+ {
+
+ if( row!=0 && row!=573 )
+ {
+ cm2->cd(matrix+1);
+
+ // if(row<288)
+ // {
+ plot[pix]->SetLineColor(2);
+ // }
+ // else
+ // {
+ // plot[pix]->SetLineColor(3);
+ // }
+
+ if( FIRST[matrix] )
+ {
+ plot[pix]->Draw("same");
+ }
+ else
+ {
+ plot[pix]->Draw();
+ FIRST[matrix] = 1;
+ }
+ }
+ }
+ }
+
+// f->Close();
+}
+
+//####################################################################
+
+void calcpara() {
+
+ UShort_t threshold;
+ Int_t pixel, row, column;
+// Int_t* pixel = new Int_t [1152*576];
+ Float_t* pixelprob = new Float_t [1152*576];
+ Float_t mean, sigma, chi2;
+
+// File IN
+ TString RootFile = "outputsc.root";
+ TFile* f = new TFile(RootFile);
+ TTree* tree = (TTree*)f->Get( "scurvesall" );
+
+ tree -> SetBranchAddress ("threshold" , &threshold );
+// tree -> SetBranchAddress ("pixel" , pixel );
+ tree -> SetBranchAddress ("pixelprob" , pixelprob );
+
+// File OUT
+ TString RootFileOut = "outputpara.root";
+ TFile* fout = new TFile(RootFileOut,"RECREATE");
+ TTree* paraTree[4];
+
+ for(int ba=0;ba<4;ba++)
+ {
+ paraTree[ba] = new TTree( Form("para%i",ba), Form("para%i",ba));
+ paraTree[ba] -> Branch("pixel" , &pixel , "pixel/i" , 32000);
+ paraTree[ba] -> Branch("row" , &row , "row/i" , 32000);
+ paraTree[ba] -> Branch("column" , &column , "column/F" , 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);
+ }
+// --------------
+ Int_t entries = tree->GetEntries();
+
+ printf("\n-------------\n");
+ printf("Found thresholds: %i\n",entries);
+ printf("-------------\n");
+
+ Float_t* pixelprob_arr = new Float_t[entries*1152*576];
+ Int_t* threshold_arr = new Int_t[entries];
+
+ TH1F *histos, *histo[12];
+ histos = new TH1F( "test", "test" , entries,0.,255.);
+ histo[0] = new TH1F( "mean0" , "mean0" , 2550,0.,255.);
+ histo[1] = new TH1F( "mean1" , "mean1" , 2550,0.,255.);
+ histo[2] = new TH1F( "mean2" , "mean2" , 2550,0.,255.);
+ histo[3] = new TH1F( "mean3" , "mean3" , 2550,0.,255.);
+ histo[4] = new TH1F( "sigma0", "sigma0" , 1000,0.,100.);
+ histo[5] = new TH1F( "sigma1", "sigma1" , 1000,0.,100.);
+ histo[6] = new TH1F( "sigma2", "sigma2" , 1000,0.,100.);
+ histo[7] = new TH1F( "sigma3", "sigma3" , 1000,0.,100.);
+ histo[8] = new TH1F( "chi20" , "chi20" , 1000,0.,1.);
+ histo[9] = new TH1F( "chi21" , "chi21" , 1000,0.,1.);
+ histo[10] = new TH1F( "chi22" , "chi22" , 1000,0.,1.);
+ histo[11] = new TH1F( "chi23" , "chi23" , 1000,0.,1.);
+
+ TF1* erf = new TF1("erf","0.5*(1+TMath::Erf((-x+[0])/[1]))",0,100);
+
+ int chicount=0;
+ int meanmin, meanmax;
+ float meanmean;
+ float vsigma;
+ float vtmp=1;
+// int column;
+// int row;
+ int matrix;
+ float MAXX[3] = {0, 0, 0};
+ double MAXY[3] = {0, 0, 0};
+
+ TCanvas* canv = new TCanvas("sc","sc", 50,100,900,600);
+ canv->Divide(3,1);
+
+ printf("Perform fitting of S-Curces:\n");
+
+ for(int thr=0;thr<entries;thr++)
+ {
+ tree->GetEntry(thr);
+
+ for(int pix=0;pix<1152*576;pix++)
+ {
+ pixelprob_arr[pix*entries+thr] = pixelprob[pix];
+ }
+
+ threshold_arr[thr] = (int)threshold;
+
+ }
+
+ for(int pix=0;pix<1152*576;pix++)
+ {
+// if(pix%100==0) {
+// printf("\r - %6.1f %% - %5i -", 100.*pix/(1152*576), pix);
+// // canv->Update();
+// fflush(stdout);
+// }
+
+ if( TMath::RMS(entries,&pixelprob_arr[pix*entries])!=0 )
+ {
+ delete histos;
+ histos = new TH1F( "test", "test" , entries,0.,255.);
+
+ meanmin=0;
+ meanmax=255;
+
+ for(int thr=0;thr<entries;thr++)
+ {
+ histos->Fill(threshold_arr[thr],pixelprob_arr[pix*entries+thr]);
+
+ if( pixelprob_arr[pix*entries+thr]!=1 && meanmin==0 ) { meanmin=threshold_arr[thr]-1; }
+ if( pixelprob_arr[pix*entries+thr]==0 && vtmp!=0 ) { meanmax=threshold_arr[thr]; }
+ vtmp = pixelprob_arr[pix*entries+thr];
+ }
+
+ vsigma = (meanmax-meanmin)/2.;
+ meanmean = (meanmax+meanmin)/2.;
+ erf->SetParameters(meanmean,vsigma);
+
+ pixel = pix;
+ column = pix%1152;
+ row = (int)(pix/1152);
+ matrix = int(column/288);
+
+ if( row!=0 && row!=573 )
+ {
+ histos->Fit(erf,"QN");
+
+ mean = erf->GetParameter(0);
+ sigma = erf->GetParameter(1)/TMath::Sqrt(2);
+ chi2 = erf->GetChisquare()/erf->GetNDF();
+
+ chicount=0;
+
+ 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<0.0001) && chicount<10)
+ {
+ erf->SetParameters(meanmean,1.*(chicount+1));
+ histos->Fit(erf,"QMN");
+ 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);
+ }
+ }
+
+ histo[0+matrix]->Fill(mean);
+ histo[4+matrix]->Fill(sigma);
+ histo[8+matrix]->Fill(chi2);
+
+ if( mean > MAXX[0] ) { MAXX[0]=mean+1; }
+ if( sigma > MAXX[1] ) { MAXX[1]=sigma+1; }
+ if( chi2 > MAXX[2] ) { MAXX[2]=chi2+0.01; }
+
+ if(pix%100==0)
+ {
+ printf("\r Pixel: %5i - %6.1f %% - fit attempts: %5i - %6.2f %6.2f %6.5f", pix, 100.*pix/(1152*576), chicount, mean, sigma, chi2); fflush(stdout);
+
+ if(pix%5000==0)
+ {
+
+ if( histo[0]->GetBinContent(histo[0]->GetMaximumBin()) > MAXY[0] ) { MAXY[0] = histo[0]->GetBinContent(histo[0]->GetMaximumBin())+1; }
+ if( histo[1]->GetBinContent(histo[1]->GetMaximumBin()) > MAXY[0] ) { MAXY[0] = histo[1]->GetBinContent(histo[1]->GetMaximumBin())+1; }
+ if( histo[2]->GetBinContent(histo[2]->GetMaximumBin()) > MAXY[0] ) { MAXY[0] = histo[2]->GetBinContent(histo[2]->GetMaximumBin())+1; }
+ if( histo[3]->GetBinContent(histo[3]->GetMaximumBin()) > MAXY[0] ) { MAXY[0] = histo[3]->GetBinContent(histo[3]->GetMaximumBin())+1; }
+
+ if( histo[4]->GetBinContent(histo[4]->GetMaximumBin()) > MAXY[1] ) { MAXY[1] = histo[4]->GetBinContent(histo[4]->GetMaximumBin())+1; }
+ if( histo[5]->GetBinContent(histo[5]->GetMaximumBin()) > MAXY[1] ) { MAXY[1] = histo[5]->GetBinContent(histo[5]->GetMaximumBin())+1; }
+ if( histo[6]->GetBinContent(histo[6]->GetMaximumBin()) > MAXY[1] ) { MAXY[1] = histo[6]->GetBinContent(histo[6]->GetMaximumBin())+1; }
+ if( histo[7]->GetBinContent(histo[7]->GetMaximumBin()) > MAXY[1] ) { MAXY[1] = histo[7]->GetBinContent(histo[7]->GetMaximumBin())+1; }
+
+ if( histo[8]->GetBinContent(histo[8]->GetMaximumBin()) > MAXY[2] ) { MAXY[2] = histo[8]->GetBinContent(histo[8]->GetMaximumBin())+1; }
+ if( histo[9]->GetBinContent(histo[9]->GetMaximumBin()) > MAXY[2] ) { MAXY[2] = histo[9]->GetBinContent(histo[9]->GetMaximumBin())+1; }
+ if( histo[10]->GetBinContent(histo[10]->GetMaximumBin()) > MAXY[2] ) { MAXY[2] = histo[10]->GetBinContent(histo[10]->GetMaximumBin())+1; }
+ if( histo[11]->GetBinContent(histo[11]->GetMaximumBin()) > MAXY[2] ) { MAXY[2] = histo[11]->GetBinContent(histo[11]->GetMaximumBin())+1; }
+
+ canv->cd(1);
+ histo[0]->SetLineColor(1);
+ histo[0]->Draw();
+ histo[0]->SetAxisRange(0,MAXX[0]);
+ histo[0]->SetMaximum(MAXY[0]);
+ histo[1]->SetLineColor(2);
+ histo[1]->Draw("same");
+ histo[2]->SetLineColor(3);
+ histo[2]->Draw("same");
+ histo[3]->SetLineColor(4);
+ histo[3]->Draw("same");
+
+ canv->cd(2);
+ histo[4]->SetLineColor(1);
+ histo[4]->Draw();
+ histo[4]->SetAxisRange(0,MAXX[1]);
+ histo[4]->SetMaximum(MAXY[1]);
+ histo[5]->SetLineColor(2);
+ histo[5]->Draw("same");
+ histo[6]->SetLineColor(3);
+ histo[6]->Draw("same");
+ histo[7]->SetLineColor(4);
+ histo[7]->Draw("same");
+
+ canv->cd(3);
+ histo[8]->SetLineColor(1);
+ histo[8]->Draw();
+ histo[8]->SetAxisRange(0,MAXX[2]);
+ histo[8]->SetMaximum(MAXY[2]);
+ histo[9]->SetLineColor(2);
+ histo[9]->Draw("same");
+ histo[10]->SetLineColor(3);
+ histo[10]->Draw("same");
+ histo[11]->SetLineColor(4);
+ histo[11]->Draw("same");
+
+ canv->Update();
+ }
+ }
+
+ paraTree[matrix]->Fill();
+ }
+ }
+ }
+
+
+ for(int ba=0;ba<4;ba++)
+ {
+ paraTree[ba]->Write("",TObject::kOverwrite);
+ }
+
+ fout->Save();
+ fout->Close();
+ f->Close();
+ cout<<"\n-------------"<<endl;
+}
+
+//####################################################################
+
+void ana() {
+
+// rewrite();
+ plotsc();
+// calcpara();
+}
+
+//####################################################################
\ No newline at end of file
jspc29 / mvd_soft.git / commitdiff