From 3c16621bc9cc45c24ad8d182db104bb997395e28 Mon Sep 17 00:00:00 2001
From: "Muzalevsky I.A" <muzalevsky@exolina2.jinr.ru>
Date: Mon, 26 Dec 2016 17:10:54 +0300
Subject: [PATCH] added script for fitting data

---
 convertTektronix/IntegralFULL.c | 110 ++++++++++++++++++++++++++++++++
 1 file changed, 110 insertions(+)
 create mode 100644 convertTektronix/IntegralFULL.c

diff --git a/convertTektronix/IntegralFULL.c b/convertTektronix/IntegralFULL.c
new file mode 100644
index 0000000..7139501
--- /dev/null
+++ b/convertTektronix/IntegralFULL.c
@@ -0,0 +1,110 @@
+void IntegralFULL(){
+
+	Int_t j,i,nentry,NumM;
+	Double_t A[1000]; /// амплитуды с сырых данных (отрицательные)
+	Double_t T[1000];
+	Double_t Am[1000]; // амплитуды с перевёрнутых пиков
+	Double_t Tim[1000];
+	Double_t sum,zeroLEVEL1,zeroLEVEL2,tm,t1_15,m,m1_15,t2_15,m2_15,tL_15,mL_15,Integrall,LinPar,ChiSquare,t_90;
+	TString var;
+
+	//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);
+
+	TF1 *fit2 = new TF1("fit2","[0]"); ////////// создание функций фитирования
+	fit2->SetRange(8e-08,1e-07);
+	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");
+//////	
+	TBranch *branch1 = tree1->GetBranch("A");
+	branch1->SetAddress(A);
+	
+	TBranch *branch2 = tree1->GetBranch("T");
+	branch2->SetAddress(T);
+//////
+	TFile *f1 = new TFile("exp2Integral.root","RECREATE"); // файл для записи
+	TTree *theTree = new TTree("theTree","peak");
+	theTree->Branch("Am",Am,"Am[1000]/D");
+	theTree->Branch("Tim",Tim,"Tim[1000]/D");
+	theTree->Branch("sum",&sum,"sum/D");
+	theTree->Branch("zeroLEVEL1",&zeroLEVEL1,"zeroLEVEL1/D");   
+	theTree->Branch("zeroLEVEL2",&zeroLEVEL2,"zeroLEVEL2/D"); 
+	theTree->Branch("tm",&tm,"tm/D");   	//// tm - положение пика(максимума амплитуды) записывается время(T[i]) при котором этот максимум был
+	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%
+	theTree->Branch("m2_15",&m2_15,"m2_15/D"); // амплитуда во второй точке превышения 20%
+	theTree->Branch("tL_15",&tL_15,"tL_15/D"); // время последнего превышения 20%
+	theTree->Branch("mL_15",&mL_15,"mL_15/D"); //амплитуда в последнем превышении 20%
+	theTree->Branch("Integrall",&Integrall,"Integrall/D"); //интеграл в диапозоне 20 нс
+	theTree->Branch("LinPar",&LinPar,"LinPar/D"); //параметр фита прямой границы сигнала
+	theTree->Branch("ChiSquare",&ChiSquare,"ChiSquare/D"); //хи квадрат
+///////// заполнение дерева
+	TCanvas *c1 = new TCanvas("c","zero line"); 
+	j=0;
+	for(j=0;j<tree1->GetEntries();j++)
+	{
+		tree1->GetEntry(j);
+		sum=0; // обнуление интеграла
+		m=0;  	// и максимума
+		for(i=0;i<1000;i++){  			// нахожу интеграл, полный, точку максимума для фрейма номер j
+			Am[i] = -A[i];
+			Tim[i] = T[i];
+			if((Am[i]>m) || (Am[i]==m)){m = Am[i];tm = T[i];}
+			sum = sum + Am[i];
+		}
+		NumM = (tm - 5e-013)*(1e+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];}
+		}
+////////////////// делаю только для первых 12 фреймов
+			{		
+			TGraph *gr1 = new TGraph(1000,Tim,Am); // поиск нуля
+			gr1->Draw("Al*");
+			gr1->Fit("fit1","R");
+			zeroLEVEL1  = fit1->GetParameter(0);			
+			c1->Update();
+			gr1->Fit("fit2","R");
+			zeroLEVEL2  = fit2->GetParameter(0);			
+			c1->Update();			
+			m1_15=0;
+			t1_15=0; //поиск начальной границы поднятия на 20%
+			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;}						
+			}
+			for(i=NumM;i>0;i--)
+			{
+				if((Am[i]-zeroLEVEL1)<0.9*(m-zeroLEVEL1)) {t_90 = Tim[i+1];break;}						
+			}			
+			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;}						
+			}
+			for(i=NumM;i<1000;i++)
+			{
+				if((Am[i]-zeroLEVEL1)>0.15*(m-zeroLEVEL1)) {tL_15 = Tim[i-1];mL_15 = A[i-1];}						
+			}
+			fit3->SetRange(t1_15,t_90);
+			gr1->Fit("fit3","R");
+			c1->Update();
+			ChiSquare = fit3->GetChisquare();
+			LinPar = fit3->GetParameter(0);
+		}
+		theTree->Fill();
+	}
+//////////	
+	theTree->Write();
+	f1->Close();
+	f->Close();
+}
-- 
2.18.1