diff --git a/Projects/S034/Calibration/CalibVDrift/CalibVDrift.txt b/Projects/S034/Calibration/CalibVDrift/CalibVDrift.txt new file mode 100644 index 0000000000000000000000000000000000000000..457c0e9a3ff45163c661555b948c5dc5055efe47 --- /dev/null +++ b/Projects/S034/Calibration/CalibVDrift/CalibVDrift.txt @@ -0,0 +1,36 @@ +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 diff --git a/Projects/S034/Calibration/CalibVDrift/MakeHisto.cxx b/Projects/S034/Calibration/CalibVDrift/MakeHisto.cxx new file mode 100644 index 0000000000000000000000000000000000000000..bd9d997aed03995980fc77f3494b1c4f0a12d9d4 --- /dev/null +++ b/Projects/S034/Calibration/CalibVDrift/MakeHisto.cxx @@ -0,0 +1,13 @@ +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(); +} diff --git a/Projects/S034/Calibration/CalibVDrift/TPad_distrib.root b/Projects/S034/Calibration/CalibVDrift/TPad_distrib.root new file mode 100755 index 0000000000000000000000000000000000000000..c0429336e064e31c189bc550167aa868772b9d8e Binary files /dev/null and b/Projects/S034/Calibration/CalibVDrift/TPad_distrib.root differ diff --git a/Projects/S034/Calibration/CalibVDrift/cal.cxx b/Projects/S034/Calibration/CalibVDrift/cal.cxx new file mode 100644 index 0000000000000000000000000000000000000000..4f3ca38dd024b5c0b4a5f2aef3b6b2634cc5f3b8 --- /dev/null +++ b/Projects/S034/Calibration/CalibVDrift/cal.cxx @@ -0,0 +1,212 @@ +/* #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; +} diff --git a/Projects/S034/Calibration/FDC0_Calib.txt b/Projects/S034/Calibration/FDC0_Calib.txt new file mode 100644 index 0000000000000000000000000000000000000000..b42eab0460664741e7a105ad67c8818ffb2e75ef --- /dev/null +++ b/Projects/S034/Calibration/FDC0_Calib.txt @@ -0,0 +1,8 @@ +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 diff --git a/Projects/S034/Calibration/FDC2_Calib.txt b/Projects/S034/Calibration/FDC2_Calib.txt new file mode 100644 index 0000000000000000000000000000000000000000..7c8b9c77ca43222f37849bb2c4e214cccd9e3b2e --- /dev/null +++ b/Projects/S034/Calibration/FDC2_Calib.txt @@ -0,0 +1,14 @@ +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