From 85a2dca8b5ab4fdb2aebb5d1a3a860575a774390 Mon Sep 17 00:00:00 2001
From: "Muzalevsky I.A" <muzalevsky@exolina2.jinr.ru>
Date: Tue, 10 Jan 2017 20:09:34 +0300
Subject: [PATCH] script for analizing data with CFD method modified

---
 convertTektronix/IntegralCFD.c | 53 ++++++++++++++++++++--------------
 1 file changed, 31 insertions(+), 22 deletions(-)

diff --git a/convertTektronix/IntegralCFD.c b/convertTektronix/IntegralCFD.c
index b46d420..b7e60da 100644
--- a/convertTektronix/IntegralCFD.c
+++ b/convertTektronix/IntegralCFD.c
@@ -1,9 +1,10 @@
 void IntegralCFD(){
-	
+	gSystem->Load("../libData.so");
+
 	const Int_t Tdelta=12;
 	const Double_t coeff=0.5;
 	Int_t j,i,nentry,NumM,imax,imin;
-	Double_t A[1000]; /// амплитуды с сырых данных (отрицательные)
+	Double_t A1[1000],A2[1000],A3[1000],A4[1000]; /// амплитуды с сырых данных (отрицательные)
 	Double_t T[1000];
 	Double_t Am[1000]; // амплитуды с перевёрнутых пиков
 	Double_t Tim[1000];
@@ -14,23 +15,32 @@ void IntegralCFD(){
 	//TF1 *func = new TF1("fit",p,0,20e-09,1); // функция фитирования
 	TF1 *fit1 = new TF1("fit1","[0]");
 	fit1->SetParName(0,"Rlevel");
-	fit1->SetRange(0,2e-08);
+	fit1->SetRange(0,20);
 
 	TF1 *fit2 = new TF1("fit2","[0]"); ////////// создание функций фитирования
-	fit2->SetRange(8e-08,1e-07);
+	fit2->SetRange(80,100);
 	fit2->SetParName(0,"Llevel");
 
 	TF1 *fit3 = new TF1("fit3","[0]*x+[1]");
 	fit3->SetParName(0,"LinPar");
 
 	TFile *f = new TFile("exp2.root","READ");
-	TTree *tree1 = (TTree*)f->Get("tree");
+	TTree *tree1 = (TTree*)f->Get("rtree");
 //////	
-	TBranch *branch1 = tree1->GetBranch("A");
-	branch1->SetAddress(A);
+	TBranch *branch1 = tree1->GetBranch("A1");
+	branch1->SetAddress(A1);
+
+	TBranch *branch2 = tree1->GetBranch("A2");
+	branch2->SetAddress(A2);
+
+	TBranch *branch3 = tree1->GetBranch("A3");
+	branch3->SetAddress(A3);
+
+	TBranch *branch4 = tree1->GetBranch("A4");
+	branch4->SetAddress(A4);
 	
-	TBranch *branch2 = tree1->GetBranch("T");
-	branch2->SetAddress(T);
+	TBranch *branch5 = tree1->GetBranch("T");
+	branch5->SetAddress(T);
 //////
 	TFile *f1 = new TFile("exp2Integral.root","RECREATE"); // файл для записи
 	TTree *theTree = new TTree("theTree","peak");
@@ -40,7 +50,7 @@ void IntegralCFD(){
 	theTree->Branch("zeroLEVEL1",&zeroLEVEL1,"zeroLEVEL1/D");   
 	theTree->Branch("zCFD",&zCFD,"zCFD/D"); 
 	theTree->Branch("tm",&tm,"tm/D");   	//// tm - положение пика(максимума амплитуды) записывается время(T[i]) при котором этот максимум был
-	theTree->Branch("t1_15",&t1_15,"t1	_15/D"); 		/// время когда амплитуда достигла 20 процентов от максимума( смотрим слева направо)
+	theTree->Branch("t1_15",&t1_15,"t1_15/D"); 		/// время когда амплитуда достигла 20 процентов от максимума( смотрим слева направо)
 	theTree->Branch("t2_15",&t2_15,"t2_15/D"); 
 	theTree->Branch("m",&m,"m/D"); // максимальная амплитуда
 	theTree->Branch("m1_15",&m1_15,"m1_15/D"); // амплитуда в первой точке превышения 20%
@@ -64,10 +74,10 @@ void IntegralCFD(){
 		for(i=0;i<1000;i++) {Acfd[i]=0;}
 		for(i=0;i<1000;i++)// нахожу интеграл, полный, точку максимума для фрейма номер j
 		{  			
-			Am[i] = -A[i];
+			Am[i] = -A1[i];
 			Tim[i] = T[i];
 			if( i>Tdelta ){
-				Acfd[i] = A[i]*coeff;
+				Acfd[i] = A1[i]*coeff;
 				Acfd[i] = Acfd[i] + Am[i-Tdelta];
 			}
 			if((Am[i]>m) || (Am[i]==m)){m = Am[i];tm = T[i];}
@@ -75,8 +85,8 @@ void IntegralCFD(){
 		}
 		maxCFD = 0; minCFD = 0;
 		for(i=0;i<1000;i++){ // поиск максимума и минимума у преобразованного сигнала, для поиска пересечения с нулём
-			if(Acfd[i]>maxCFD) {maxCFD = Acfd[i];imax=i;t1=Tim[i];}	
-			if(Acfd[i]<minCFD) {minCFD = Acfd[i];imin=i;t2=Tim[i];}	
+			if(Acfd[i]>maxCFD) {maxCFD = Acfd[i];imax=i;}	
+			if(Acfd[i]<minCFD) {minCFD = Acfd[i];imin=i;}	
 		}
 
 		/*TGraph *gr3 = new TGraph(1000,Tim,Acfd);
@@ -84,11 +94,11 @@ void IntegralCFD(){
 		c1->Update();		
 		cout<<endl<<imin<<"  "<<imax<<endl;	*/	
 
-		NumM = (tm - 5e-013)*(1e+10); // номер точки (0,1000) с максимальной амплитудой сигнала
+		NumM = (tm - 5e-04)*(10); // номер точки (0,1000) с максимальной амплитудой сигнала
 		Integrall=0;
 		for(i=0;i<1000;i++)	// нахожу интеграл, c границами
 		{  			
-			if( (Tim[i]>(tm-3e-9)) && (Tim[i]<(tm+17e-9)) ) {Integrall = Integrall + Am[i];}
+			if( (Tim[i]>(tm-3)) && (Tim[i]<(tm+17)) ) {Integrall = Integrall + Am[i];}
 		}
 		//cout<<endl<<"integral= "<<Integrall<<endl;
 ////////////////// 		
@@ -102,7 +112,7 @@ void IntegralCFD(){
 		t2_15=0;m2_15 = 0; mL_15 = 0;			
 		for(i=NumM;i>0;i--)
 		{
-			if((Am[i]-zeroLEVEL1)<0.1*(m-zeroLEVEL1)) {t1_15 = Tim[i+1];m1_15 = A[i+1];break;}						
+			if((Am[i]-zeroLEVEL1)<0.1*(m-zeroLEVEL1)) {t1_15 = Tim[i+1];m1_15 = Am[i+1];break;}						
 		}
 		for(i=NumM;i>0;i--)
 		{
@@ -114,25 +124,24 @@ void IntegralCFD(){
 		}				
 		for(i=NumM;i<1000;i++)
 		{
-			if((Am[i]-zeroLEVEL1)<0.15*(m-zeroLEVEL1)) {t2_15 = Tim[i-1];m2_15 = A[i-1];break;}						
+			if((Am[i]-zeroLEVEL1)<0.15*(m-zeroLEVEL1)) {t2_15 = Tim[i-1];m2_15 = Am[i-1];break;}						
 		}
 		for(i=NumM;i<1000;i++)
 		{
-				if((Am[i]-zeroLEVEL1)>0.15*(m-zeroLEVEL1)) {tL_15 = Tim[i-1];mL_15 = A[i-1];}						
+				if((Am[i]-zeroLEVEL1)>0.15*(m-zeroLEVEL1)) {tL_15 = Tim[i-1];mL_15 = Am[i-1];}						
 		}
 
 		deltaZero = 100;
 		for(i=imin;i<imax;i++){   // поиск нуля в сигнале CFD
-			if(abs(A[i])<deltaZero) {deltaZero = abs(A[i]); Tcfd = T[i];}	
+			if(abs(Acfd[i])<deltaZero) {deltaZero = abs(Acfd[i]); Tcfd = T[i];}	
 		}
 		
-		
 		fit3->SetRange(t1_15,t1_90);
 		gr1->Fit("fit3","R","goff");
 		c1->Update();
 		ChiSquare = fit3->GetChisquare();
 		LinPar = fit3->GetParameter(0);
-
+		
 		theTree->Fill();
 	}	
 //////////	
-- 
2.18.1