From 3fb4f66816ee388c3c76e106713f730d5b2e97cd Mon Sep 17 00:00:00 2001 From: "Muzalevsky I.A" Date: Mon, 9 Jan 2017 15:24:43 +0300 Subject: [PATCH] script for CFD method added --- convertTektronix/IntegralCFD.c | 142 +++++++++++++++++++++++++++++++++ 1 file changed, 142 insertions(+) create mode 100644 convertTektronix/IntegralCFD.c diff --git a/convertTektronix/IntegralCFD.c b/convertTektronix/IntegralCFD.c new file mode 100644 index 0000000..b46d420 --- /dev/null +++ b/convertTektronix/IntegralCFD.c @@ -0,0 +1,142 @@ +void IntegralCFD(){ + + 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 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,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("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;jGetEntries();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] = -A[i]; + Tim[i] = T[i]; + if( i>Tdelta ){ + Acfd[i] = A[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;t1=Tim[i];} + if(Acfd[i]Draw("Al*"); + c1->Update(); + cout<(tm-3e-9)) && (Tim[i]<(tm+17e-9)) ) {Integrall = Integrall + Am[i];} + } + //cout<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 = A[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 = 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];} + } + + deltaZero = 100; + for(i=imin;iSetRange(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(); +} -- 2.18.1