193 lines
6.5 KiB
C++
193 lines
6.5 KiB
C++
#include "LevenbergMarquardt.h"
|
|
#include "utils.h"
|
|
#include <TCanvas.h>
|
|
#include <TF1.h>
|
|
#include <TFile.h>
|
|
#include <TH1F.h>
|
|
#include <TTree.h>
|
|
|
|
#include <Eigen/Dense>
|
|
#include <iostream>
|
|
|
|
using namespace std;
|
|
using Eigen::MatrixXd;
|
|
using Eigen::VectorXd;
|
|
|
|
void getADC(const char *fin) {
|
|
// read file
|
|
TFile *fRun = new TFile(fin);
|
|
TTree *t = (TTree *)fRun->Get("Tree1");
|
|
|
|
// data numbers
|
|
Long64_t ntot = t->GetEntriesFast();
|
|
|
|
// five X-4 block
|
|
UInt_t bArray[16];
|
|
UInt_t cArray[16];
|
|
UInt_t dArray[16];
|
|
UInt_t eArray[16];
|
|
UInt_t fArray[16];
|
|
|
|
// calibration coefficient
|
|
Float_t CalC[5][8];
|
|
Float_t Calk1[5][8];
|
|
Float_t Calk2[5][8];
|
|
|
|
// read adc data
|
|
t->SetBranchAddress("adc0ch0", &bArray[0]);
|
|
t->SetBranchAddress("adc0ch1", &bArray[1]);
|
|
t->SetBranchAddress("adc0ch2", &bArray[2]);
|
|
t->SetBranchAddress("adc0ch3", &bArray[3]);
|
|
t->SetBranchAddress("adc0ch4", &bArray[4]);
|
|
t->SetBranchAddress("adc0ch5", &bArray[5]);
|
|
t->SetBranchAddress("adc0ch6", &bArray[6]);
|
|
t->SetBranchAddress("adc0ch7", &bArray[7]);
|
|
t->SetBranchAddress("adc0ch8", &bArray[8]);
|
|
t->SetBranchAddress("adc0ch9", &bArray[9]);
|
|
t->SetBranchAddress("adc0ch10", &bArray[10]);
|
|
t->SetBranchAddress("adc0ch11", &bArray[11]);
|
|
t->SetBranchAddress("adc0ch12", &bArray[12]);
|
|
t->SetBranchAddress("adc0ch13", &bArray[13]);
|
|
t->SetBranchAddress("adc0ch14", &bArray[14]);
|
|
t->SetBranchAddress("adc0ch15", &bArray[15]);
|
|
|
|
t->SetBranchAddress("adc0ch16", &cArray[0]);
|
|
t->SetBranchAddress("adc0ch17", &cArray[1]);
|
|
t->SetBranchAddress("adc0ch18", &cArray[2]);
|
|
t->SetBranchAddress("adc0ch19", &cArray[3]);
|
|
t->SetBranchAddress("adc0ch20", &cArray[4]);
|
|
t->SetBranchAddress("adc0ch21", &cArray[5]);
|
|
t->SetBranchAddress("adc0ch22", &cArray[6]);
|
|
t->SetBranchAddress("adc0ch23", &cArray[7]);
|
|
t->SetBranchAddress("adc0ch24", &cArray[8]);
|
|
t->SetBranchAddress("adc0ch25", &cArray[9]);
|
|
t->SetBranchAddress("adc0ch26", &cArray[10]);
|
|
t->SetBranchAddress("adc0ch27", &cArray[11]);
|
|
t->SetBranchAddress("adc0ch28", &cArray[12]);
|
|
t->SetBranchAddress("adc0ch29", &cArray[13]);
|
|
t->SetBranchAddress("adc0ch30", &cArray[14]);
|
|
t->SetBranchAddress("adc0ch31", &cArray[15]);
|
|
|
|
t->SetBranchAddress("adc1ch0", &dArray[0]);
|
|
t->SetBranchAddress("adc1ch1", &dArray[1]);
|
|
t->SetBranchAddress("adc1ch2", &dArray[2]);
|
|
t->SetBranchAddress("adc1ch3", &dArray[3]);
|
|
t->SetBranchAddress("adc1ch4", &dArray[4]);
|
|
t->SetBranchAddress("adc1ch5", &dArray[5]);
|
|
t->SetBranchAddress("adc1ch6", &dArray[6]);
|
|
t->SetBranchAddress("adc1ch7", &dArray[7]);
|
|
t->SetBranchAddress("adc1ch8", &dArray[8]);
|
|
t->SetBranchAddress("adc1ch9", &dArray[9]);
|
|
t->SetBranchAddress("adc1ch10", &dArray[10]);
|
|
t->SetBranchAddress("adc1ch11", &dArray[11]);
|
|
t->SetBranchAddress("adc1ch12", &dArray[12]);
|
|
t->SetBranchAddress("adc1ch13", &dArray[13]);
|
|
t->SetBranchAddress("adc1ch14", &dArray[14]);
|
|
t->SetBranchAddress("adc1ch15", &dArray[15]);
|
|
|
|
t->SetBranchAddress("adc1ch16", &eArray[0]);
|
|
t->SetBranchAddress("adc1ch17", &eArray[1]);
|
|
t->SetBranchAddress("adc1ch18", &eArray[2]);
|
|
t->SetBranchAddress("adc1ch19", &eArray[3]);
|
|
t->SetBranchAddress("adc1ch20", &eArray[4]);
|
|
t->SetBranchAddress("adc1ch21", &eArray[5]);
|
|
t->SetBranchAddress("adc1ch22", &eArray[6]);
|
|
t->SetBranchAddress("adc1ch23", &eArray[7]);
|
|
t->SetBranchAddress("adc1ch24", &eArray[8]);
|
|
t->SetBranchAddress("adc1ch25", &eArray[9]);
|
|
t->SetBranchAddress("adc1ch26", &eArray[10]);
|
|
t->SetBranchAddress("adc1ch27", &eArray[11]);
|
|
t->SetBranchAddress("adc1ch28", &eArray[12]);
|
|
t->SetBranchAddress("adc1ch29", &eArray[13]);
|
|
t->SetBranchAddress("adc1ch30", &eArray[14]);
|
|
t->SetBranchAddress("adc1ch31", &eArray[15]);
|
|
|
|
t->SetBranchAddress("adc2ch0", &fArray[0]);
|
|
t->SetBranchAddress("adc2ch1", &fArray[1]);
|
|
t->SetBranchAddress("adc2ch2", &fArray[2]);
|
|
t->SetBranchAddress("adc2ch3", &fArray[3]);
|
|
t->SetBranchAddress("adc2ch4", &fArray[4]);
|
|
t->SetBranchAddress("adc2ch5", &fArray[5]);
|
|
t->SetBranchAddress("adc2ch6", &fArray[6]);
|
|
t->SetBranchAddress("adc2ch7", &fArray[7]);
|
|
t->SetBranchAddress("adc2ch8", &fArray[8]);
|
|
t->SetBranchAddress("adc2ch9", &fArray[9]);
|
|
t->SetBranchAddress("adc2ch10", &fArray[10]);
|
|
t->SetBranchAddress("adc2ch11", &fArray[11]);
|
|
t->SetBranchAddress("adc2ch12", &fArray[12]);
|
|
t->SetBranchAddress("adc2ch13", &fArray[13]);
|
|
t->SetBranchAddress("adc2ch14", &fArray[14]);
|
|
t->SetBranchAddress("adc2ch15", &fArray[15]);
|
|
|
|
TH1F *Left = new TH1F("Left", "Left", 300, 0, 300);
|
|
TH1F *Right = new TH1F("Right", "Right", 300, 0, 300);
|
|
|
|
for (int i = 0; i < ntot; i++) {
|
|
t->GetEntry(i);
|
|
Left->Fill(bArray[0]);
|
|
Right->Fill(bArray[8]);
|
|
}
|
|
|
|
// use matrix and least square method to get coefficient
|
|
// SigmaClip *clip = new SigmaClip(5);
|
|
int cnt = 0;
|
|
double vX[300];
|
|
for (int i = 2; i < 300; i++)
|
|
if (Left->GetBinContent(i) > 0) {
|
|
vX[cnt] = i;
|
|
cnt += 1;
|
|
}
|
|
|
|
VectorXd mY(cnt);
|
|
MatrixXd mX(cnt, 3);
|
|
MatrixXd mB(3, 1);
|
|
double n;
|
|
double *parma = new double[3];
|
|
double sigma, b0, b1, b2;
|
|
|
|
cnt = 0;
|
|
for (int i = 2; i < 300; i++) {
|
|
n = Left->GetBinContent(i);
|
|
if (n == 0) continue;
|
|
mY(cnt, 0) = log(n);
|
|
mX(cnt, 0) = 1, mX(cnt, 1) = i, mX(cnt, 2) = i * i;
|
|
cnt += 1;
|
|
}
|
|
|
|
mB = (mX.transpose() * mX).inverse() * mX.transpose() * mY;
|
|
b0 = mB(0, 0), b1 = mB(1, 0), b2 = mB(2, 0);
|
|
sigma = -1 / b2;
|
|
parma[1] = b1 * sigma / 2;
|
|
parma[0] = exp(b0 + parma[1] * parma[1] / sigma);
|
|
parma[2] = sqrt(sigma);
|
|
|
|
cout << parma[0] << ", " << parma[1] << ", " << parma[2] << endl;
|
|
|
|
LevenbergMarquardt LM(3, parma, Gaussian, GaussianJacobian);
|
|
for (int i = 2; i < 300; i++) {
|
|
n = Left->GetBinContent(i);
|
|
if (n == 0) continue;
|
|
LM.addData(i, n);
|
|
}
|
|
parma = LM.solve();
|
|
cout << parma[0] << ", " << parma[1] << ", " << parma[2] << endl;
|
|
|
|
// TCanvas *c1 = new TCanvas("c1", "Fitting with Gaussian function");
|
|
// c1->Divide(1, 2);
|
|
|
|
// c1->cd(1);
|
|
// c1->SetGrid();
|
|
// Left->Draw();
|
|
|
|
// c1->cd(2);
|
|
// c1->SetGrid();
|
|
// Right->Draw();
|
|
|
|
// TF1 *GausFunc = new TF1("GausFunc", "[0]*exp(-(x-[1])*(x-[1])/2.0/[2]/[2])", 0, 300);
|
|
// GausFunc->SetParLimits(0, 1000, 2400);
|
|
// GausFunc->SetParLimits(1, 100, 200);
|
|
// GausFunc->SetParLimits(2, 0, 15);
|
|
// Left->Fit(GausFunc, "R");
|
|
// Right->Fit(GausFunc, "R");
|
|
}
|