File indexing completed on 2023-03-17 11:23:16
0001 #include <iostream>
0002 #include <stdlib.h>
0003 #include <cmath>
0004 #include "TROOT.h"
0005 #include "TFile.h"
0006 #include "TRint.h"
0007 #include "TH1.h"
0008 #include "TF1.h"
0009 #include "TF2.h"
0010 #include "TH2.h"
0011 #include "TCanvas.h"
0012 #include "TMinuit.h"
0013 #include "TSystem.h"
0014 #include "TGraphErrors.h"
0015 #include "TMatrixD.h"
0016 #include "TMatrixDSym.h"
0017 #include "NtupleHelper.h"
0018 #include "global.h"
0019
0020
0021
0022 const char par_name[dim][20] = {
0023 "z0 ", "sigma ", "x0 ", "y0 ", "dxdz ", "dydz ", "sigmaBeam ", "c0 ", "c1 ", "emittance ", "betastar "};
0024
0025 Double_t params[dim], errparams[dim];
0026 Double_t sfpar[dim], errsfpar[dim];
0027
0028 constexpr Double_t step[dim] = {1.e-5, 1.e-5, 1.e-3, 1.e-3, 1.e-3, 1.e-3, 1.e-5, 1.e-5, 1.e-5, 1.e-5, 1.e-5};
0029 zData zdata;
0030 int tmpNtrks_ = 0;
0031 int fnthite = 0;
0032 TMatrixD ftmp;
0033
0034 Double_t fd0cut = 4.0;
0035
0036 int sequence[11] = {500, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000};
0037 Double_t sequenceD[11] = {500, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000};
0038 Double_t zeros[11] = {0.};
0039
0040 int iloop_ = 0;
0041
0042 int Nsequence_ = 11;
0043 Double_t res_d0phi_X[11] = {0.};
0044 Double_t res_d0phi_Xerr[11] = {0.};
0045 Double_t res_d0phi_Y[11] = {0.};
0046 Double_t res_d0phi_Yerr[11] = {0.};
0047 Double_t res_d0phi_dXdz[11] = {0.};
0048 Double_t res_d0phi_dXdzerr[11] = {0.};
0049 Double_t res_d0phi_dYdz[11] = {0.};
0050 Double_t res_d0phi_dYdzerr[11] = {0.};
0051
0052 Double_t res_Z[11] = {0.};
0053 Double_t res_sigmaZ[11] = {0.};
0054 Double_t res_Zerr[11] = {0.};
0055 Double_t res_sigmaZerr[11] = {0.};
0056 Double_t res_Z_lh[11] = {0.};
0057 Double_t res_sigmaZ_lh[11] = {0.};
0058 Double_t res_Zerr_lh[11] = {0.};
0059 Double_t res_sigmaZerr_lh[11] = {0.};
0060
0061 TGraphErrors *g_X;
0062 TGraphErrors *g_Y;
0063 TGraphErrors *g_dXdz;
0064 TGraphErrors *g_dYdz;
0065 TGraphErrors *g_Z;
0066 TGraphErrors *g_sigmaZ;
0067 TGraphErrors *g_Z_lh;
0068 TGraphErrors *g_sigmaZ_lh;
0069
0070 Double_t zdis(Double_t z, Double_t sigma, Double_t *parms) {
0071
0072
0073
0074
0075
0076 Double_t sig = sqrt(sigma * sigma + parms[Par_Sigma] * parms[Par_Sigma]);
0077 Double_t result = (exp(-((z - parms[Par_Z0]) * (z - parms[Par_Z0])) / (2.0 * sig * sig))) / (sig * sqrt2pi);
0078 return result;
0079 }
0080
0081 void zfcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *params, Int_t iflag) {
0082
0083
0084
0085 f = 0.0;
0086 for (zDataConstIter iter = zdata.begin(); iter != zdata.end(); ++iter) {
0087 f = log(zdis(iter->Z, iter->SigZ, params)) + f;
0088 }
0089 f = -2.0 * f;
0090 return;
0091 }
0092
0093 float betaf(float betastar, float emmitance, float z, float z0) {
0094 float x = sqrt(emmitance * (betastar + (((z - z0) * (z - z0)) / betastar)));
0095 return x;
0096 }
0097
0098 Double_t ddis(Double_t z, Double_t sigma, Double_t d, Double_t sigmad, Double_t *parms) {
0099
0100
0101
0102
0103 Double_t sig = betaf(parms[Par_beta], parms[Par_eps], z, parms[Par_Z0]);
0104 sig = sqrt(sig * sig + sigmad * sigmad);
0105 Double_t result = (exp(-(d * d) / (2.0 * sig * sig))) / (sig * sqrt2pi);
0106 return result;
0107 }
0108
0109 Double_t allddis(Double_t z, Double_t d, Double_t sigmad, Double_t phi0, Double_t *parms) {
0110
0111
0112
0113
0114 Double_t sig =
0115 sqrt(parms[Par_Sigbeam] * parms[Par_Sigbeam] + 10 * sigmad * sigmad);
0116 Double_t dprime =
0117 d - (parms[Par_x0] + z * parms[Par_dxdz]) * sin(phi0) + (parms[Par_y0] + z * parms[Par_dydz]) * cos(phi0);
0118 Double_t result = (exp(-(dprime * dprime) / (2.0 * sig * sig))) / (sig * sqrt2pi);
0119 return result;
0120 }
0121
0122 Double_t allddisBeta(Double_t z, Double_t sigma, Double_t d, Double_t sigmad, Double_t phi0, Double_t *parms) {
0123
0124
0125
0126
0127 Double_t sigmabeam =
0128 sqrt(parms[Par_eps] * (parms[Par_beta] + (((z - parms[Par_Z0]) * (z - parms[Par_Z0])) / parms[Par_beta])));
0129
0130 Double_t sig = sqrt(sigmabeam * sigmabeam + sigmad * sigmad);
0131
0132 Double_t dprime =
0133 d - (parms[Par_x0] + z * parms[Par_dxdz]) * sin(phi0) + (parms[Par_y0] + z * parms[Par_dydz]) * cos(phi0);
0134 Double_t result = (exp(-(dprime * dprime) / (2.0 * sig * sig))) / (sig * sqrt2pi);
0135 return result;
0136 }
0137
0138 Double_t allddis2(Double_t z, Double_t sigma, Double_t d, Double_t pt, Double_t phi0, Double_t *parms) {
0139
0140
0141
0142 Double_t sigmad = parms[Par_c0] + parms[Par_c1] / pt;
0143 Double_t sig = sqrt(parms[Par_Sigbeam] * parms[Par_Sigbeam] + sigmad * sigmad);
0144 Double_t dprime =
0145 d - (parms[Par_x0] + z * parms[Par_dxdz]) * sin(phi0) + (parms[Par_y0] + z * parms[Par_dydz]) * cos(phi0);
0146 Double_t result = (exp(-(dprime * dprime) / (2.0 * sig * sig))) / (sig * sqrt2pi);
0147 return result;
0148 }
0149
0150 void dfcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *params, Int_t iflag) {
0151
0152
0153
0154 f = 0.0;
0155 for (zDataConstIter iter = zdata.begin(); iter != zdata.end(); ++iter) {
0156
0157
0158 f = log(allddis(iter->Z, iter->D, iter->SigD, iter->Phi, params)) + f;
0159 }
0160 f = -2.0 * f;
0161 return;
0162 }
0163
0164 void dfcnbeta(Int_t &npar, Double_t *gin, Double_t &f, Double_t *params, Int_t iflag) {
0165
0166
0167
0168 f = 0.0;
0169 for (zDataConstIter iter = zdata.begin(); iter != zdata.end(); ++iter) {
0170 f = log(allddisBeta(iter->Z, iter->SigZ, iter->D, iter->SigD, iter->Phi, params)) + f;
0171 }
0172 f = -2.0 * f;
0173 return;
0174 }
0175
0176 void dfcn2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *params, Int_t iflag) {
0177
0178
0179
0180 f = 0.0;
0181 for (zDataConstIter iter = zdata.begin(); iter != zdata.end(); ++iter) {
0182 f = log(allddis2(iter->Z, iter->SigZ, iter->D, iter->Pt, iter->Phi, params) * zdis(iter->Z, iter->SigZ, params)) +
0183 f;
0184 }
0185 f = -2.0 * f;
0186 return;
0187 }
0188
0189 void cfcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *params, Int_t iflag) {
0190
0191
0192
0193 f = 0.0;
0194 for (zDataConstIter iter = zdata.begin(); iter != zdata.end(); ++iter) {
0195 f = log(allddis(iter->Z, iter->D, iter->SigD, iter->Phi, params) * zdis(iter->Z, iter->SigZ, params)) + f;
0196 }
0197 f = -2.0 * f;
0198 return;
0199 }
0200
0201 void fit(TMatrixD &x, TMatrixDSym &V) {
0202 TMatrixDSym Vint(4);
0203 TMatrixD b(4, 1);
0204 Double_t weightsum = 0;
0205 tmpNtrks_ = 0;
0206
0207 Vint.Zero();
0208 b.Zero();
0209 TMatrixD g(4, 1);
0210 TMatrixDSym temp(4);
0211 for (zDataConstIter i = zdata.begin(); i != zdata.end(); ++i) {
0212
0213 if (i->weight2 == 0)
0214 continue;
0215
0216 g(0, 0) = sin(i->Phi);
0217 g(1, 0) = -cos(i->Phi);
0218 g(2, 0) = i->Z * g(0, 0);
0219 g(3, 0) = i->Z * g(1, 0);
0220
0221 temp.Zero();
0222 for (int j = 0; j < 4; ++j) {
0223 for (int k = j; k < 4; ++k) {
0224 temp(j, k) += g(j, 0) * g(k, 0);
0225 }
0226 }
0227 double sigmabeam2 = 0.002 * 0.002;
0228
0229 double sigma2 = sigmabeam2 + (i->SigD) * (i->SigD);
0230
0231 TMatrixD ftmptrans(1, 4);
0232 ftmptrans = ftmptrans.Transpose(ftmp);
0233 TMatrixD dcor = ftmptrans * g;
0234
0235 bool pass = true;
0236 if ((fnthite > 0) && (std::abs(i->D - dcor(0, 0)) > fd0cut))
0237 pass = false;
0238
0239
0240
0241
0242 if (pass) {
0243 Vint += (temp * (1 / sigma2));
0244 b += ((i->D / sigma2) * g);
0245 weightsum += sqrt(i->weight2);
0246 tmpNtrks_++;
0247 }
0248 }
0249 Double_t determinant;
0250 V = Vint.InvertFast(&determinant);
0251 x = V * b;
0252 ftmp = x;
0253
0254 std::cout << "number of tracks used in this iteration = " << tmpNtrks_ << std::endl;
0255 std::cout << "Sum of all weights:" << weightsum << "\n";
0256 std::cout << "x0 :" << x(0, 0) << " +- " << sqrt(V(0, 0)) << " cm \n";
0257 std::cout << "y0 :" << x(1, 0) << " +- " << sqrt(V(1, 1)) << " cm \n";
0258 std::cout << "x slope :" << x(2, 0) << " +- " << sqrt(V(2, 2)) << " cm \n";
0259 std::cout << "y slope :" << x(3, 0) << " +- " << sqrt(V(3, 3)) << " cm \n";
0260
0261 res_d0phi_X[iloop_] = x(0, 0);
0262 res_d0phi_Xerr[iloop_] = sqrt(V(0, 0));
0263 res_d0phi_Y[iloop_] = x(1, 0);
0264 res_d0phi_Yerr[iloop_] = sqrt(V(1, 1));
0265 res_d0phi_dXdz[iloop_] = x(2, 0);
0266 res_d0phi_dXdzerr[iloop_] = sqrt(V(2, 2));
0267 res_d0phi_dYdz[iloop_] = x(3, 0);
0268 res_d0phi_dYdzerr[iloop_] = sqrt(V(3, 3));
0269 }
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292 Double_t cutOnChi2(const TMatrixD &x, Double_t chi2cut) {
0293 double sigmabeam2 = 0.002 * 0.002;
0294 TMatrixD g(1, 4);
0295 Double_t weightsum = 0;
0296 for (zDataIter i = zdata.begin(); i != zdata.end(); ++i) {
0297 g(0, 0) = sin(i->Phi);
0298 g(0, 1) = -cos(i->Phi);
0299 g(0, 2) = i->Z * g(0, 0);
0300 g(0, 3) = i->Z * g(0, 1);
0301 TMatrixD dcor = g * x;
0302
0303 Double_t chi2 = (i->D - dcor(0, 0)) * (i->D - dcor(0, 0)) / (sigmabeam2 + (i->SigD) * (i->SigD));
0304 if (chi2 > chi2cut) {
0305 i->weight2 = 0;
0306 } else {
0307 i->weight2 = 1;
0308 weightsum += 1;
0309 }
0310 }
0311
0312 return weightsum;
0313 }
0314
0315 void fitAll() {
0316 std::cout << "starting fits" << std::endl;
0317
0318
0319 TH1F *h1z = new TH1F("h1z", "z distribution", 100, -50., 50.);
0320 for (zDataConstIter i = zdata.begin(); i != zdata.end(); ++i) {
0321 h1z->Fill(i->Z, i->SigZ);
0322 }
0323 h1z->Sumw2();
0324 h1z->Fit("gaus");
0325
0326
0327 TF1 *fgaus = h1z->GetFunction("gaus");
0328
0329 res_Z[iloop_] = fgaus->GetParameter(1);
0330 res_sigmaZ[iloop_] = fgaus->GetParameter(2);
0331 res_Zerr[iloop_] = fgaus->GetParError(1);
0332 res_sigmaZerr[iloop_] = fgaus->GetParError(2);
0333
0334 std::cout << "======== End of Chi2 Fit for Z ========\n" << std::endl;
0335
0336
0337 TMinuit *gmMinuit = new TMinuit(2);
0338 gmMinuit->SetFCN(zfcn);
0339 std::cout << "SetFCN done" << std::endl;
0340 int ierflg = 0;
0341 sfpar[Par_Sigma] = 7.55;
0342 sfpar[Par_Z0] = 0.;
0343 errsfpar[Par_Sigma] = 0.;
0344 errsfpar[Par_Z0] = 0.;
0345
0346
0347
0348
0349
0350 gmMinuit->SetFCN(zfcn);
0351 for (int i = 0; i < 2; i++) {
0352 gmMinuit->mnparm(i, par_name[i], sfpar[i], step[i], 0, 0, ierflg);
0353
0354 }
0355
0356
0357
0358 gmMinuit->Migrad();
0359
0360
0361
0362 for (int i = 0; i < 2; i++) {
0363 gmMinuit->GetParameter(i, sfpar[i], errsfpar[i]);
0364 }
0365
0366 res_Z_lh[iloop_] = sfpar[Par_Z0];
0367 res_sigmaZ_lh[iloop_] = sfpar[Par_Sigma];
0368 res_Zerr_lh[iloop_] = errsfpar[Par_Z0];
0369 res_sigmaZerr_lh[iloop_] = errsfpar[Par_Sigma];
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382 std::cout << "======== End of Maximum LH Fit for Z ========\n" << std::endl;
0383
0384
0385
0386 tmpNtrks_ = 0;
0387
0388 TMatrixD xmat(4, 1);
0389 TMatrixDSym Verr(4);
0390 fit(xmat, Verr);
0391 fnthite++;
0392
0393
0394
0395 int fminNtrks = 100;
0396 double fconvergence = 0.9;
0397
0398 while (tmpNtrks_ > fconvergence * zdata.size()) {
0399
0400
0401
0402
0403
0404
0405
0406 fit(xmat, Verr);
0407 fd0cut /= 1.5;
0408
0409 if (tmpNtrks_ > fconvergence * zdata.size() && tmpNtrks_ > fminNtrks)
0410 fnthite++;
0411 }
0412 std::cout << " number of iterations: " << fnthite << std::endl;
0413
0414 std::cout << "======== End of d0-phi Fit ========\n" << std::endl;
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446
0447
0448
0449
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459
0460
0461
0462
0463
0464
0465
0466
0467
0468
0469
0470
0471
0472
0473
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495
0496
0497
0498
0499
0500
0501
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522
0523
0524
0525
0526
0527
0528
0529
0530 }
0531
0532 int main(int argc, char **argv) {
0533
0534
0535
0536 std::cout << "name: " << argv[1] << std::endl;
0537
0538 ftmp.ResizeTo(4, 1);
0539 ftmp.Zero();
0540
0541 int cut_on_events = 0;
0542 bool run_sequence = false;
0543
0544 if (argc >= 3) {
0545 TString tmpst(argv[2]);
0546 cut_on_events = tmpst.Atoi();
0547 if (cut_on_events == -1) {
0548 run_sequence = true;
0549 std::cout << " run a sequence of fits for 0.5k, 1k, 2k, 5k, 10k" << std::endl;
0550 } else {
0551 std::cout << " maximum number of tracks = " << cut_on_events << std::endl;
0552 }
0553 }
0554
0555 TString in_filename(argv[1]);
0556 TString tmp_string = in_filename;
0557 TString out_filename = tmp_string.Replace(in_filename.Length() - 5, in_filename.Length(), "_results.root");
0558 if (argc == 4) {
0559 out_filename = TString(argv[3]);
0560 }
0561
0562
0563
0564
0565 TFile *outroot = new TFile(out_filename, "RECREATE");
0566
0567 NtupleHelper *t = new NtupleHelper(in_filename);
0568 std::cout << " ntuple initialized" << std::endl;
0569
0570
0571
0572
0573 t->Book();
0574 std::cout << " ntuple booked" << std::endl;
0575
0576 if (run_sequence) {
0577 for (int iloop = 0; iloop < Nsequence_; iloop++) {
0578 std::cout << " loop over # " << sequence[iloop] << " tracks." << std::endl;
0579 zdata = t->Loop(sequence[iloop]);
0580 fitAll();
0581 std::cout << " loop has finished." << std::endl;
0582 iloop_++;
0583 }
0584
0585 outroot->cd();
0586 g_X = new TGraphErrors(Nsequence_, sequenceD, res_d0phi_X, zeros, res_d0phi_Xerr);
0587 g_Y = new TGraphErrors(Nsequence_, sequenceD, res_d0phi_Y, zeros, res_d0phi_Yerr);
0588 g_dXdz = new TGraphErrors(Nsequence_, sequenceD, res_d0phi_dXdz, zeros, res_d0phi_dXdzerr);
0589 g_dYdz = new TGraphErrors(Nsequence_, sequenceD, res_d0phi_dYdz, zeros, res_d0phi_dYdzerr);
0590
0591 g_Z = new TGraphErrors(Nsequence_, sequenceD, res_Z, zeros, res_Zerr);
0592 g_sigmaZ = new TGraphErrors(Nsequence_, sequenceD, res_sigmaZ, zeros, res_sigmaZerr);
0593
0594 g_Z_lh = new TGraphErrors(Nsequence_, sequenceD, res_Z_lh, zeros, res_Zerr_lh);
0595 g_sigmaZ_lh = new TGraphErrors(Nsequence_, sequenceD, res_sigmaZ_lh, zeros, res_sigmaZerr_lh);
0596
0597 g_X->SetName("beamX");
0598 g_Y->SetName("beamY");
0599 g_dXdz->SetName("beamdXdz");
0600 g_dYdz->SetName("beamdYdz");
0601 g_Z->SetName("beamZ");
0602 g_sigmaZ->SetName("beamSigmaZ");
0603 g_Z_lh->SetName("beamZ_likelihood");
0604 g_sigmaZ_lh->SetName("beamSigmaZ_likelihood");
0605
0606
0607
0608
0609 g_X->Write();
0610 g_Y->Write();
0611 g_dXdz->Write();
0612 g_dYdz->Write();
0613 g_Z->Write();
0614 g_sigmaZ->Write();
0615 g_Z_lh->Write();
0616 g_sigmaZ_lh->Write();
0617
0618 outroot->Close();
0619
0620 } else {
0621 zdata = t->Loop(cut_on_events);
0622 fitAll();
0623 std::cout << " finished loop over events" << std::endl;
0624 }
0625
0626
0627
0628 return 0;
0629 }