IntegralCFD.c

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 A1[1000],A2[1000],A3[1000],A4[1000]; /// амплитуды с сырых данных (отрицательные)
	Double_t T[1000];
	Double_t Am[1000]; // амплитуды с перевёрнутых пиков
	Double_t Tim[1000];
	Double_t Acfd[1000];
	Double_t sum,zeroLEVEL1,zeroLEVEL2,tm,t1_15,m,m1_15,t2_15,m2_15,tL_15,mL_15,Integrall,LinPar,ChiSquare,t1_90,t2_90,minCFD,maxCFD,zCFD,Tcfd,deltaZero;
	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,20);

	TF1 *fit2 = new TF1("fit2","[0]"); ////////// создание функций фитирования
	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("rtree");
//////	
	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 *branch5 = tree1->GetBranch("T");
	branch5->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("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("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"); //хи квадрат
	theTree->Branch("t2_90",&t2_90,"t2_90/D"); // 90% от максимума после максимума
	theTree->Branch("Acfd",Acfd,"Acfd[1000]/D"); // am cdf
	theTree->Branch("Tcfd",&Tcfd,"Tcfd/D"); // момент пересечения нуля в CFD
///////// заполнение дерева
	TCanvas *c1 = new TCanvas("c","zero line"); 
	j=0;
	for(j=0;j<tree1->GetEntries();j++) // j - номер события , i - внутри!
	{
		tree1->GetEntry(j);
		sum=0; // обнуление интеграла
		m=0;  	// и максимума
		for(i=0;i<1000;i++) {Acfd[i]=0;}
		for(i=0;i<1000;i++)// нахожу интеграл, полный, точку максимума для фрейма номер j
		{  			
			Am[i] = -A1[i];
			Tim[i] = T[i];
			if( i>Tdelta ){
				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];}
			sum = sum + Am[i];
		}
		maxCFD = 0; minCFD = 0;
		for(i=0;i<1000;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);
		gr3->Draw("Al*");
		c1->Update();		
		cout<<endl<<imin<<"  "<<imax<<endl;	*/	

		NumM = (tm - 5e-04)*(10); // номер точки (0,1000) с максимальной амплитудой сигнала
		Integrall=0;
		for(i=0;i<1000;i++)	// нахожу интеграл, c границами
		{  			
			if( (Tim[i]>(tm-3)) && (Tim[i]<(tm+17)) ) {Integrall = Integrall + Am[i];}
		}
		//cout<<endl<<"integral= "<<Integrall<<endl;
////////////////// 		
		TGraph *gr1 = new TGraph(1000,Tim,Am); // поиск нуля
		gr1->Draw("Al*");
		gr1->Fit("fit1","R","goff");
		zeroLEVEL1  = fit1->GetParameter(0);			
		c1->Update();					
		m1_15=0;
		t1_15=0; //поиск начальной границы поднятия на 10%
		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 = Am[i+1];break;}						
		}
		for(i=NumM;i>0;i--)
		{
			if((Am[i]-zeroLEVEL1)<0.9*(m-zeroLEVEL1)) {t1_90 = Tim[i+1];break;}						
		}
		for(i=NumM;i<1000;i++)
		{
			if((Am[i]-zeroLEVEL1)<0.9*(m-zeroLEVEL1)) {t2_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 = 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 = Am[i-1];}						
		}

		deltaZero = 100;
		for(i=imin;i<imax;i++){   // поиск нуля в сигнале CFD
			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();
	}	
//////////	
	theTree->Write();
	f1->Close();
	f->Close();
}