From 6a3411f7ecbdba3c763152d2b1e204a1f1eb522b Mon Sep 17 00:00:00 2001
From: Philipp Klaus <klaus@physik.uni-frankfurt.de>
Date: Tue, 19 Aug 2014 18:39:05 +0200
Subject: [PATCH] s-curves: latest process_rootmulti_4xbanks_long from Samir

---
 s-curves/Makefile                         |  10 +-
 s-curves/README.txt                       |  33 ++
 s-curves/ana.c                            | 417 ++++++++++++++++++++++
 s-curves/process_rootmulti_4xbanks_long.c |  20 +-
 4 files changed, 472 insertions(+), 8 deletions(-)
 create mode 100644 s-curves/ana.c

diff --git a/s-curves/Makefile b/s-curves/Makefile
index d4ef599..f4180b4 100644
--- a/s-curves/Makefile
+++ b/s-curves/Makefile
@@ -3,11 +3,17 @@ CFLAGS=-c -g -Wall -std=c++11 `root-config --cflags`
 LDFLAGS=`root-config --glibs`
 CLDFLAGS = -std=c++11 `root-config --cflags --glibs` -lboost_system -lboost_filesystem 
 
-all: process_rootmulti_4xbanks.c process_rootmulti_4xbanks_long.c
+SOURCES= process_rootmulti_4xbanks_long.c
+OBJECTS=$(subst .cc,.o,$(SOURCES))
+
+all: process_rootmulti_4xbanks process_rootmulti_4xbanks_long
+
+process_rootmulti_4xbanks: process_rootmulti_4xbanks.c
 	$(CC) $(CLDFLAGS) process_rootmulti_4xbanks.c -o process_rootmulti_4xbanks
+
+process_rootmulti_4xbanks_long: process_rootmulti_4xbanks_long.c
 	$(CC) $(CLDFLAGS) process_rootmulti_4xbanks_long.c -o process_rootmulti_4xbanks_long
 
 clean:
 	rm -f process_rootmulti_4xbanks
 	rm -f process_rootmulti_4xbanks_long
-
diff --git a/s-curves/README.txt b/s-curves/README.txt
index 0e36b24..94a7776 100644
--- a/s-curves/README.txt
+++ b/s-curves/README.txt
@@ -1,6 +1,7 @@
 Tools to analyze s-curve measurements
 =====================================
 
+
 These tools were developed locally on jspc55 and their history
 was then reconstructed to be put into this Git repository.
 
@@ -11,6 +12,38 @@ Those .hld files then have to be processed by the C++ program
 The ROOT data file can then be further analyzed with a simple ROOT macro,
 such as `view_all_averages.c` and `view_rootfile.c`.
 
+
+Howto for the latest process_rootmulti_4xbanks_long.c and ana.c
+---------------------------------------------------------------------
+
+1.) process_rootmulti_4xbanks_long.c
+
+-> Make sure to put the right output folder for output.txt!
+-> perform 'pico process_rootmulti_4xbanks_long.c' and change the corresponding line
+	-> "str + o" (save) and  "str + x" (exit)
+-> perform 'make'
+
+2.) change the naming of the data to be processed, 
+
+-> e.g. 
+	-> file_1.hld
+	-> file_2.hld
+	-> file_3.hld
+-> run process_rootmulti_4xbanks_long
+	-> e.g. './process_rootmulti_4xbanks_long /local.1/probetest/NA61/A02/Scurve/file_'
+
+3.)
+-> run ana.c
+	-> make root available: '. thisroot.sh'
+	-> 'root -l'
+	-> compile ana.c: '.L ana.c+'
+	-> 'rewrite()'
+	-> 'calcpara()'
+	-> 'plotsc()'
+
+4.) move the file you want to not be overwritten to another folder
+
+
 Authors
 -------
 
diff --git a/s-curves/ana.c b/s-curves/ana.c
new file mode 100644
index 0000000..7dd183f
--- /dev/null
+++ b/s-curves/ana.c
@@ -0,0 +1,417 @@
+#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
diff --git a/s-curves/process_rootmulti_4xbanks_long.c b/s-curves/process_rootmulti_4xbanks_long.c
index f09d347..efb8c68 100644
--- a/s-curves/process_rootmulti_4xbanks_long.c
+++ b/s-curves/process_rootmulti_4xbanks_long.c
@@ -53,9 +53,10 @@ char *htoi(const char *ptr) {
 
 
 
-// char output[256] = "/local.1/daq/workdir/scurve/output.txt\0";
+char output[256] = "/local.1/probetest/NA61/software/output.txt\0";
 
-char output[256] = "/local.1/daq/workdir/scurve/output.txt\0";
+// char output[256] = "/local.1/daq/workdir/scurve/output.txt\0";
+// char output[256] = "/u/samir/files/software/proto/scurves/output.txt\0";
 
 FILE *fo;
 int mode = 0;
@@ -296,6 +297,8 @@ void analyze(void){
 }
 
 int main(int argc, char *argv[]) {
+	gDirectory->SaveSelf();
+	
 		printf("-------------\n");
   char* filename;
 
@@ -584,15 +587,15 @@ int main(int argc, char *argv[]) {
 	for(int i=0;i<576;i++) {rowcount[i]=0;}
 	
 	
-// Save
+/*// Save
 // 	printf("\n%s","-------------\n");
 // 	printf("Write data: %i Thresholds\n", entries);
 	scurveTree	->Write("",TObject::kOverwrite);
 	outFile->Write();
-// 	outFile->Save();
+	outFile->Save();
 // 	printf("\n%s","-------------\n");
 // 	printf("%s","-------------\n");
-// 	printf("%s created\n",Form("%s",outFileName.Data()));
+	printf("%s created\n",Form("%s",outFileName.Data()));
 // 	printf("%s","-------------\n");
 	
 	entries = (int) scurveTree->GetEntries();
@@ -784,7 +787,7 @@ int main(int argc, char *argv[]) {
 		}
 		
 		paraTree[ba]->Write("",TObject::kOverwrite);
-		outFile->Write();
+// 		outFile->Write();
 // 		outFile->Save();
 	}
 // 	
@@ -798,6 +801,11 @@ int main(int argc, char *argv[]) {
 // 
 // 	canv->Update();
 	
+// 	outFile->Write();
+// 	outFile->Save();
+	
+	scurvesall->Delete("all");*/
+	
 	outFile->Write();
 	outFile->Save();
 	outFile->Close();
-- 
2.43.0