Commit 6fca7423 authored by Adrien Matta's avatar Adrien Matta
Browse files

* Adding calibration file and macro to S034

parent 4a2cc574
Pipeline #120194 passed with stages
in 9 minutes and 27 seconds
Minos_R1_VDRIFT 0.03508
Minos_R1_OFFSET -1202.54
Minos_R2_VDRIFT 0.03508
Minos_R2_OFFSET -1201.37
Minos_R3_VDRIFT 0.03508
Minos_R3_OFFSET -1428.81
Minos_R4_VDRIFT 0.03508
Minos_R4_OFFSET -1389.57
Minos_R5_VDRIFT 0.03508
Minos_R5_OFFSET -1202.95
Minos_R6_VDRIFT 0.03508
Minos_R6_OFFSET -1387.21
Minos_R7_VDRIFT 0.03508
Minos_R7_OFFSET -1355.77
Minos_R8_VDRIFT 0.03508
Minos_R8_OFFSET -1387.43
Minos_R9_VDRIFT 0.03508
Minos_R9_OFFSET -1220.77
Minos_R10_VDRIFT 0.03508
Minos_R10_OFFSET -1202.43
Minos_R11_VDRIFT 0.03508
Minos_R11_OFFSET -1207.41
Minos_R12_VDRIFT 0.03508
Minos_R12_OFFSET -1194.95
Minos_R13_VDRIFT 0.03508
Minos_R13_OFFSET -1211.2
Minos_R14_VDRIFT 0.03508
Minos_R14_OFFSET -1200.51
Minos_R15_VDRIFT 0.03508
Minos_R15_OFFSET -1201.83
Minos_R16_VDRIFT 0.03508
Minos_R16_OFFSET -1200.02
Minos_R17_VDRIFT 0.03508
Minos_R17_OFFSET -1199.43
Minos_R18_VDRIFT 0.03508
Minos_R18_OFFSET -1214.62
void MakeHisto(){
auto file = new TFile("../../../Outputs/Analysis/s034_test.root","READ");
auto tree = (TTree*) file->FindObjectAny("PhysicsTree");
auto out = new TFile("TPad_distrib_new.root","RECREATE");
for(unsigned int i = 0 ; i < 18 ; i++){
cout << i << endl;
tree->Draw(Form("T_Pad>>TPad_ring%d(2000)",i+1),Form("Ring_Pad==%d",i+1));
// auto h = (TH1F*) gDirectory->FindObjectAny(Form("TPad_ring%d",i+1));
}
out->Write();
out->Close();
}
/* #include "TGraphErrors.h" */
TGraphErrors* GetT(int ring);
TGraphErrors* GetZ(int ring);
TGraphErrors* Calibrate(const TGraphErrors* g, double offset, double vdrift);
void Scale(const TGraphErrors* ref,TGraphErrors* mod);
double chi2(TGraphErrors* g1,TGraphErrors* g2);
double Chi2(const double* parameter);
void NumericalMinimization(const char * minName ,const char *algoName);
TGraphErrors* z;// simulated Z distribution (ref)
TGraphErrors* t;// measured T distribution
TGraphErrors* c;// calibrated Z distribution from t
vector<double> vdrift;
vector<double> evdrift;
vector<double> offset;
vector<double> eoffset;
////////////////////////////////////////////////////////////////////////////////
void cal(){
TCanvas* cv = new TCanvas();
cv->Divide(4,5);
ofstream outfile;
outfile.open("CalibVDrift.txt");
for(unsigned int i = 0 ; i < 18; i++){
cout << "Processing Ring " << i << endl;
cv->cd(i+1);
t = GetT(i+1);
z = GetZ(i+1);
NumericalMinimization("Minuit2","Combined"); //SCAN, MIGRAD, COMBINED
z->Draw("ap");
z->SetLineColor(kRed);
c->SetLineColor(kBlack);
c->Draw("p");
outfile << "Minos_R" << i+1 << "_VDRIFT " << vdrift[i] << endl;
outfile << "Minos_R" << i+1 << "_OFFSET " << offset[i] << endl;
}
outfile.close();
TCanvas* c2 = new TCanvas();
c2->Divide(2,1);
c2->cd(1);
vector<double> x = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18};
vector<double> ex;// = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
ex.resize(18,0);
TGraphErrors* VDrift = new TGraphErrors(vdrift.size(),&x[0],&vdrift[0],&ex[0],&evdrift[0]);
TGraphErrors* Offset= new TGraphErrors(offset.size(),&x[0],&offset[0],&ex[0],&eoffset[0]);
VDrift->Draw("ap");
c2->cd(2);
Offset->Draw("ap");
}
////////////////////////////////////////////////////////////////////////////////
TGraphErrors* GetT(int ring){
static TFile* Tfile = new TFile("TPad_distrib_new.root");
TString name = Form("TPad_ring%d",ring);
TH1F* h = (TH1F*) Tfile->FindObjectAny(name);
h->Scale(1./h->Integral());
unsigned int size = h->GetNbinsX();
vector<double> x ;
vector<double> y ;
vector<double> ex;
vector<double> ey;
for(unsigned int i = 0 ; i < size ; i++){
x.push_back(h->GetBinCenter(i));
y.push_back(h->GetBinContent(i));
ex.push_back(0);
ey.push_back(h->GetBinError(i));
}
TGraphErrors* g = new TGraphErrors(size,&x[0],&y[0],&ex[0],&ey[0]);
return g;
}
////////////////////////////////////////////////////////////////////////////////
TGraphErrors* GetZ(int ring){
static TFile* Tfile = new TFile("ZPad_distrib.root");
TString name = Form("ZPad_Ring%d",ring);
TH1F* h = (TH1F*) Tfile->FindObjectAny(name);
h->Scale(1./h->Integral());
unsigned int size = h->GetNbinsX();
vector<double> x ;
vector<double> y ;
vector<double> ex;
vector<double> ey;
for(unsigned int i = 0 ; i < size ; i++){
x.push_back(h->GetBinCenter(i));
y.push_back(h->GetBinContent(i));
ex.push_back(0);
ey.push_back(h->GetBinError(i));
}
TGraphErrors* g = new TGraphErrors(size,&x[0],&y[0],&ex[0],&ey[0]);
g->SetMarkerColor(kAzure+7);
g->SetLineColor(kAzure+7);
return g;
}
////////////////////////////////////////////////////////////////////////////////
TGraphErrors* Calibrate(const TGraphErrors* g, double offset, double vdrift){
TGraphErrors* res = new TGraphErrors(*g);
unsigned int size = g->GetN();
double x, y;
for(unsigned int i = 0 ; i < size ; i++){
g->GetPoint(i, x, y);
res->SetPoint(i,(x*30+offset)*vdrift,y);
}
res->SetLineColor(kGreen-3);
res->SetMarkerColor(kGreen-3);
return res;
}
////////////////////////////////////////////////////////////////////////////////
void Scale(const TGraphErrors* ref,TGraphErrors* mod){
unsigned int size = ref->GetN();
double scale=0;
int point=0;
double x, y;
for(unsigned int i = 0 ; i < size ; i++){
ref->GetPoint(i,x,y);
double xref = x;
double yref = y;
double ymod = mod->Eval(xref);
if(yref>0.0002){
scale+= ymod/yref;
point++;
}
}
double val = scale/point;
size = mod->GetN();
for(unsigned int i = 0 ; i < size ; i++){
mod->GetPoint(i,x,y);
mod->SetPoint(i,x,y/val);
}
}
////////////////////////////////////////////////////////////////////////////////
double Chi2(const double* parameter){
c = Calibrate(t,parameter[0],parameter[1]);
Scale(z,c);
return chi2(z,c);
}
////////////////////////////////////////////////////////////////////////////////
double chi2(TGraphErrors* g1,TGraphErrors* g2){
double chi2 = 0;
unsigned int size = g1->GetN();
double x, y;
for(unsigned int i = 0 ; i < size ; i++){
g1->GetPoint(i, x, y);
double xg1 = x;
double yg1 = y;
double yg2 = g2->Eval(xg1);
if(yg2 && yg1)
chi2+=(yg1-yg2)*(yg1-yg2)/(g1->GetErrorY(i));
}
return chi2;
}
////////////////////////////////////////////////////////////////////////////////
void NumericalMinimization(const char * minName ,const char *algoName){
ROOT::Math::Minimizer* min =
ROOT::Math::Factory::CreateMinimizer(minName, algoName);
// set tolerance , etc...
min->SetMaxFunctionCalls(1000000); // for Minuit/Minuit2
/* min->SetPrecision(1e-4); */
//min->SetTolerance(0.001);
min->SetValidError(true);
// min->SetPrintLevel(1);
// create funciton wrapper for minimizer
// a IMultiGenFunction type
ROOT::Math::Functor f(&Chi2,2);
min->SetFunction(f);
min->Clear();
double* parameter = new double[2];
parameter[0] = -1200;
parameter[1] = 0.03508 ;
//parameter[1] = 0.0337261;
// Set the free variables to be minimized!
min->SetLimitedVariable(0,"Offset",parameter[0],1,-2000,-1000);
//min->SetFixedVariable(0,"Offset",parameter[0]);
//min->SetLimitedVariable(1,"VDrift",parameter[1],1e-5,0.02,0.04);
// Set a fixed VDrift
min->SetFixedVariable(1,"VDrift",parameter[1]);
// do the minimization
min->Minimize();
const double *xs = min->X();
std::cout << "Offset = " << xs[0] << " +/- " << min->Errors()[0] << std::endl;
std::cout << "VDrift = " << xs[1] << " +/- " << min->Errors()[1] << std::endl;
std::cout << "Chi2= " << Chi2(xs) << std::endl;
vdrift.push_back(xs[1]);
evdrift.push_back(min->Errors()[1]);
offset.push_back(xs[0]);
eoffset.push_back(min->Errors()[0]);
return 0;
}
FDC0_L0 2.5 0.0966973 1769.21
FDC0_L1 2.5 0.096531 1769.04
FDC0_L2 2.5 0.0961997 1769.75
FDC0_L3 2.5 0.0971254 1768.75
FDC0_L4 2.5 0.0964606 1768.51
FDC0_L5 2.5 0.096712 1768.22
FDC0_L6 2.5 0.096629 1768.89
FDC0_L7 2.5 0.0969487 1768.65
FDC2_L0 10 0.0184488 1533.56
FDC2_L1 10 0.0184124 1531.96
FDC2_L2 10 0.0181578 1532.93
FDC2_L3 10 0.0181458 1531.38
FDC2_L4 10 0.0180886 1532.03
FDC2_L5 10 0.0181936 1530.02
FDC2_L6 10 0.0182746 1530.47
FDC2_L7 10 0.0182248 1530.27
FDC2_L8 10 0.0181578 1532.7
FDC2_L9 10 0.0182491 1530.83
FDC2_L10 10 0.018172 1531.31
FDC2_L11 10 0.0183112 1530.59
FDC2_L12 10 0.0183293 1531.59
FDC2_L13 10 0.0184108 1529.31
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment