File indexing completed on 2023-03-17 10:39:49
0001 #include "Alignment/MuonAlignmentAlgorithms/interface/DTMuonMillepede.h"
0002
0003 #include <iostream>
0004
0005 DTMuonMillepede::DTMuonMillepede(std::string path,
0006 int n_files,
0007 float MaxPt,
0008 float MinPt,
0009 int nPhihits,
0010 int nThetahits,
0011 int workingmode,
0012 int nMtxSection) {
0013 ntuplePath = path;
0014 numberOfRootFiles = n_files;
0015
0016 ptMax = MaxPt;
0017 ptMin = MinPt;
0018
0019 nPhiHits = nPhihits;
0020 nThetaHits = nThetahits;
0021
0022 TDirectory *dirSave = gDirectory;
0023
0024
0025 myPG = new ReadPGInfo("./InternalData2009.root");
0026
0027 f = new TFile("./LocalMillepedeResults.root", "RECREATE");
0028 f->cd();
0029
0030 setBranchTree();
0031
0032 initNTuples(nMtxSection);
0033
0034 calculationMillepede(workingmode);
0035
0036 f->Write();
0037 f->Close();
0038
0039 dirSave->cd();
0040 }
0041
0042 DTMuonMillepede::~DTMuonMillepede() {}
0043
0044 void DTMuonMillepede::calculationMillepede(int workingmode) {
0045
0046 TMatrixD ****C = new TMatrixD ***[5];
0047 TMatrixD ****b = new TMatrixD ***[5];
0048
0049 for (int whI = -2; whI < 3; ++whI) {
0050 C[whI + 2] = new TMatrixD **[4];
0051 b[whI + 2] = new TMatrixD **[4];
0052 for (int stI = 1; stI < 5; ++stI) {
0053 C[whI + 2][stI - 1] = new TMatrixD *[14];
0054 b[whI + 2][stI - 1] = new TMatrixD *[14];
0055 for (int seI = 1; seI < 15; ++seI) {
0056 if (seI > 12 && stI != 4)
0057 continue;
0058 if (stI == 4) {
0059 C[whI + 2][stI - 1][seI - 1] = new TMatrixD(24, 24);
0060 b[whI + 2][stI - 1][seI - 1] = new TMatrixD(24, 1);
0061 } else {
0062 C[whI + 2][stI - 1][seI - 1] = new TMatrixD(60, 60);
0063 b[whI + 2][stI - 1][seI - 1] = new TMatrixD(60, 1);
0064 }
0065 }
0066 }
0067 }
0068
0069
0070 if (workingmode <= 3) {
0071 Int_t nentries = (Int_t)tali->GetEntries();
0072 for (Int_t i = 0; i < nentries; i++) {
0073 tali->GetEntry(i);
0074
0075 if (i % 100000 == 0)
0076 std::cout << "Analyzing track number " << i << std::endl;
0077
0078
0079 if (pt > ptMax || pt < ptMin)
0080 continue;
0081
0082 for (int counter = 0; counter < nseg; ++counter) {
0083 if (nphihits[counter] < nPhiHits || (nthetahits[counter] < nThetaHits && st[counter] < 4))
0084 continue;
0085
0086 TMatrixD A(12, 60);
0087 TMatrixD R(12, 1);
0088 TMatrixD E(12, 12);
0089 TMatrixD B(12, 4);
0090 TMatrixD I(12, 12);
0091 if (st[counter] == 4) {
0092 A.ResizeTo(8, 24);
0093 R.ResizeTo(8, 1);
0094 E.ResizeTo(8, 8);
0095 B.ResizeTo(8, 2);
0096 I.ResizeTo(8, 8);
0097 }
0098
0099 for (int counterHi = 0; counterHi < nhits[counter]; counterHi++) {
0100 int row = 0;
0101 if (sl[counter][counterHi] == 3) {
0102 row = 4 + (la[counter][counterHi] - 1);
0103 } else if (sl[counter][counterHi] == 2) {
0104 row = 8 + (la[counter][counterHi] - 1);
0105 } else {
0106 row = (la[counter][counterHi] - 1);
0107 }
0108
0109 float x = xc[counter][counterHi];
0110 float y = yc[counter][counterHi];
0111 float xp = xcp[counter][counterHi];
0112 float yp = ycp[counter][counterHi];
0113 float zp = zc[counter][counterHi];
0114 float error = ex[counter][counterHi];
0115 float dxdz = dxdzSl[counter];
0116 float dydz = dydzSl[counter];
0117
0118 if (st[counter] != 4) {
0119 if (sl[counter][counterHi] == 2) {
0120 A(row, row * 5) = -1.0;
0121 A(row, row * 5 + 1) = dydz;
0122 A(row, row * 5 + 2) = y * dydz;
0123 A(row, row * 5 + 3) = -x * dydz;
0124 A(row, row * 5 + 4) = -x;
0125 R(row, 0) = y - yp;
0126 } else {
0127 A(row, row * 5 + 0) = -1.0;
0128 A(row, row * 5 + 1) = dxdz;
0129 A(row, row * 5 + 2) = y * dxdz;
0130 A(row, row * 5 + 3) = -x * dxdz;
0131 A(row, row * 5 + 4) = y;
0132 R(row, 0) = x - xp;
0133 }
0134 } else {
0135 if (sl[counter][counterHi] != 2) {
0136 A(row, row * 3 + 0) = -1.0;
0137 A(row, row * 3 + 1) = dxdz;
0138 A(row, row * 3 + 2) = -x * dxdz;
0139 R(row, 0) = x - xp;
0140 }
0141 }
0142
0143 E(row, row) = 1.0 / error;
0144 I(row, row) = 1.0;
0145
0146 if (sl[counter][counterHi] != 2) {
0147 B(row, 0) = 1.0;
0148 B(row, 1) = zp;
0149 } else {
0150 B(row, 2) = 1.0;
0151 B(row, 3) = zp;
0152 }
0153 }
0154
0155 TMatrixD Bt = B;
0156 Bt.T();
0157 TMatrixD At = A;
0158 At.T();
0159 TMatrixD gamma = (Bt * E * B);
0160 gamma.Invert();
0161 *(C[wh[counter] + 2][st[counter] - 1][sr[counter] - 1]) += At * E * (B * gamma * Bt * E * A - A);
0162 *(b[wh[counter] + 2][st[counter] - 1][sr[counter] - 1]) += At * E * (B * gamma * Bt * E * I - I) * R;
0163 }
0164 }
0165 }
0166
0167 if (workingmode == 3)
0168 for (int wheel = -2; wheel < 3; ++wheel)
0169 for (int station = 1; station < 5; ++station)
0170 for (int sector = 1; sector < 15; ++sector) {
0171 if (sector > 12 && station != 4)
0172 continue;
0173
0174 TMatrixD Ctr = *C[wheel + 2][station - 1][sector - 1];
0175 TMatrixD btr = *b[wheel + 2][station - 1][sector - 1];
0176
0177 TString CtrName = "Ctr_";
0178 CtrName += wheel;
0179 CtrName += "_";
0180 CtrName += station;
0181 CtrName += "_";
0182 CtrName += sector;
0183 Ctr.Write(CtrName);
0184
0185 TString btrName = "btr_";
0186 btrName += wheel;
0187 btrName += "_";
0188 btrName += station;
0189 btrName += "_";
0190 btrName += sector;
0191 btr.Write(btrName);
0192 }
0193
0194
0195 if (workingmode % 2 == 0) {
0196 for (int wheel = -2; wheel < 3; ++wheel)
0197 for (int station = 1; station < 5; ++station)
0198 for (int sector = 1; sector < 15; ++sector) {
0199 if (sector > 12 && station != 4)
0200 continue;
0201
0202 if (workingmode == 4) {
0203 for (int mf = -1; mf <= -1; mf++) {
0204 TMatrixD Ch = getMatrixFromFile("Ctr_", wheel, station, sector, mf);
0205 TMatrixD bh = getMatrixFromFile("btr_", wheel, station, sector, mf);
0206
0207 *(C[wheel + 2][station - 1][sector - 1]) += Ch;
0208 *(b[wheel + 2][station - 1][sector - 1]) += bh;
0209 }
0210 }
0211
0212 TMatrixD Ctr = *C[wheel + 2][station - 1][sector - 1];
0213 TMatrixD btr = *b[wheel + 2][station - 1][sector - 1];
0214
0215 TMatrixD Ccs = getCcsMatrix(wheel, station, sector);
0216 TMatrixD bcs = getbcsMatrix(wheel, station, sector);
0217
0218 TMatrixD SC = Ctr + Ccs;
0219 TMatrixD Sb = btr + bcs;
0220
0221 SC.Invert();
0222
0223 TMatrixD solution = SC * Sb;
0224 for (int icounter = 0; icounter < SC.GetNrows(); icounter++) {
0225 for (int jcounter = 0; jcounter < SC.GetNrows(); jcounter++) {
0226 cov[icounter][jcounter] = SC(icounter, jcounter);
0227 }
0228 }
0229
0230 int nLayer = 12, nDeg = 5;
0231 if (station == 4) {
0232 nLayer = 8;
0233 nDeg = 3;
0234 }
0235 for (int counterLayer = 0; counterLayer < nLayer; ++counterLayer)
0236 for (int counterDeg = 0; counterDeg < nDeg; ++counterDeg) {
0237 if (counterLayer < 4) {
0238 slC[counterLayer] = 1;
0239 laC[counterLayer] = counterLayer + 1;
0240 } else if (counterLayer > 3 && counterLayer < 8) {
0241 slC[counterLayer] = 3;
0242 laC[counterLayer] = counterLayer - 3;
0243 } else {
0244 slC[counterLayer] = 2;
0245 laC[counterLayer] = counterLayer - 7;
0246 }
0247
0248 if (counterLayer < 8) {
0249 dx[counterLayer] = solution(counterLayer * nDeg, 0);
0250 dy[counterLayer] = 0.0;
0251 } else {
0252 dx[counterLayer] = 0.0;
0253 dy[counterLayer] = solution(counterLayer * nDeg, 0);
0254 }
0255 dz[counterLayer] = solution(counterLayer * nDeg + 1, 0);
0256
0257 if (station != 4) {
0258 phix[counterLayer] = solution(counterLayer * nDeg + 2, 0);
0259 phiy[counterLayer] = solution(counterLayer * nDeg + 3, 0);
0260 phiz[counterLayer] = solution(counterLayer * nDeg + 4, 0);
0261 } else {
0262 phix[counterLayer] = 0.;
0263 phiy[counterLayer] = solution(counterLayer * nDeg + 2, 0);
0264 phiz[counterLayer] = 0.;
0265 }
0266 }
0267
0268 whC = wheel;
0269 stC = station;
0270 srC = sector;
0271
0272 ttreeOutput->Fill();
0273 }
0274 }
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340 }
0341
0342 TMatrixD DTMuonMillepede::getCcsMatrix(int wh, int st, int se) {
0343 int size = 60;
0344 if (st == 4)
0345 size = 24;
0346
0347 TMatrixD matrix(size, size);
0348
0349 float Error[3][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
0350 {0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005},
0351 {0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005}};
0352
0353 float TollerancyPosition = 0.01;
0354 float TollerancyRotation = 0.0001;
0355
0356 TMatrixD ConsQC = myPG->giveQCCal(wh, st, se);
0357
0358 bool UseQC = true;
0359 if (ConsQC.GetNrows() < 70)
0360 UseQC = false;
0361
0362 int nLayer = 12, nDeg = 5;
0363 if (st == 4) {
0364 nLayer = 8;
0365 nDeg = 3;
0366 }
0367
0368 for (int la = 0; la < nLayer; la++)
0369 for (int dg = 0; dg < nDeg; dg++) {
0370 int index = la * nDeg + dg;
0371
0372 int rdg = dg + 1;
0373 if (la < 8 && rdg == 1)
0374 rdg = 0;
0375 if (st == 4 && rdg == 3)
0376 rdg = 4;
0377
0378 double ThisError2 = Error[la / 4][rdg] * Error[la / 4][rdg];
0379 if (rdg < 2) {
0380 if (UseQC) {
0381 ThisError2 += ConsQC(la, 1) * ConsQC(la, 1);
0382 } else {
0383 ThisError2 += TollerancyPosition * TollerancyPosition;
0384 }
0385 } else if (rdg == 2) {
0386 ThisError2 += TollerancyPosition * TollerancyPosition;
0387 } else {
0388 ThisError2 += TollerancyRotation * TollerancyRotation;
0389 }
0390
0391 matrix(index, index) = 1. / ThisError2;
0392 }
0393
0394 return matrix;
0395 }
0396
0397 TMatrixD DTMuonMillepede::getbcsMatrix(int wh, int st, int se) {
0398 int size = 60;
0399 if (st == 4)
0400 size = 24;
0401
0402 TMatrixD matrix(size, 1);
0403
0404 float Error[3][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
0405 {0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005},
0406 {0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005}};
0407
0408 float TollerancyPosition = 0.01;
0409 float TollerancyRotation = 0.0001;
0410
0411 TMatrixD ConsQC = myPG->giveQCCal(wh, st, se);
0412
0413 bool UseQC = true;
0414 if (ConsQC.GetNrows() < 70)
0415 UseQC = false;
0416
0417 TMatrixD Survey = myPG->giveSurvey(wh, st, se);
0418
0419 int nLayer = 12, nDeg = 5;
0420 if (st == 4) {
0421 nLayer = 8;
0422 nDeg = 3;
0423 }
0424
0425 for (int la = 0; la < nLayer; la++)
0426 for (int dg = 0; dg < nDeg; dg++) {
0427 int index = la * nDeg + dg;
0428
0429 int rdg = dg + 1;
0430 if (la < 8 && rdg == 1)
0431 rdg = 0;
0432 if (st == 4 && rdg == 3)
0433 rdg = 4;
0434
0435 double ThisError2 = Error[la / 4][rdg] * Error[la / 4][rdg];
0436 if (rdg < 2) {
0437 if (UseQC) {
0438 ThisError2 += ConsQC(la, 1) * ConsQC(la, 1);
0439 } else {
0440 ThisError2 += TollerancyPosition * TollerancyPosition;
0441 }
0442 } else if (rdg == 2) {
0443 ThisError2 += TollerancyPosition * TollerancyPosition;
0444 } else {
0445 ThisError2 += TollerancyRotation * TollerancyRotation;
0446 }
0447
0448 float Constraint = 0.;
0449 if (la > 3)
0450 Constraint += Survey(rdg, la / 4 - 1);
0451 if (UseQC && rdg == 0)
0452 Constraint += ConsQC(la, 0);
0453 ;
0454 if (UseQC && rdg == 1)
0455 Constraint -= ConsQC(la, 0);
0456
0457 matrix(index, 0) = Constraint / ThisError2;
0458 }
0459
0460 return matrix;
0461 }
0462
0463 TMatrixD DTMuonMillepede::getMatrixFromFile(const TString &Code, int wh, int st, int se, int mf) {
0464 TString MtxFileName = "./LocalMillepedeMatrix_";
0465 MtxFileName += mf;
0466 MtxFileName += ".root";
0467 if (mf == -1)
0468 MtxFileName = "./LocalMillepedeMatrix.root";
0469
0470 TDirectory *dirSave = gDirectory;
0471
0472 TFile *MatrixFile = new TFile(MtxFileName);
0473
0474 TString MtxName = Code;
0475 MtxName += wh;
0476 MtxName += "_";
0477 MtxName += st;
0478 MtxName += "_";
0479 MtxName += se;
0480 TMatrixD *ThisMtx = (TMatrixD *)MatrixFile->Get(MtxName);
0481
0482 MatrixFile->Close();
0483
0484 dirSave->cd();
0485
0486 return *ThisMtx;
0487 }
0488
0489 TMatrixD DTMuonMillepede::getLagMatrix(int wh, int st, int se) {
0490 TMatrixD matrix(60 + 6, 60 + 6);
0491 if (st == 4)
0492 matrix.ResizeTo(40 + 6, 40 + 6);
0493
0494 for (int counterDeg = 0; counterDeg < 5; ++counterDeg) {
0495 for (int counterLayer = 0; counterLayer < 12; ++counterLayer) {
0496 if (st == 4) {
0497 matrix(40 + counterDeg, counterLayer * 5 + counterDeg) = 10000.0;
0498 matrix(counterLayer * 5 + counterDeg, 40 + counterDeg) = 10000.0;
0499 } else {
0500 int realCounter = counterDeg + 1;
0501 if (counterDeg == 1 && counterLayer < 8) {
0502 realCounter = counterDeg - 1;
0503 }
0504 matrix(counterLayer * 5 + counterDeg, 40 + realCounter) = 10000.0;
0505 if ((realCounter == 0 && counterLayer > 7) || (realCounter == 1 && counterLayer < 7))
0506 continue;
0507 matrix(60 + realCounter, counterLayer * 5 + counterDeg) = 10000.0;
0508 }
0509 }
0510 }
0511 return matrix;
0512 }
0513
0514 TMatrixD DTMuonMillepede::getCqcMatrix(int wh, int st, int se) {
0515 TMatrixD surv = myPG->giveQCCal(wh, st, se);
0516 TMatrixD sigmaQC(5, 12);
0517
0518 TMatrixD matrix(60, 60);
0519 if (st == 4)
0520 matrix.ResizeTo(40, 40);
0521
0522 if (surv.GetNrows() < 7) {
0523 for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
0524 for (int counterLayer = 0; counterLayer < 12; ++counterLayer) {
0525 if (st != 4 && counterLayer > 7)
0526 continue;
0527 if (counterDeg == 0) {
0528 sigmaQC(counterDeg, counterLayer) = 0.05;
0529 } else if (counterDeg < 3) {
0530 sigmaQC(counterDeg, counterLayer) = 0.05;
0531 } else {
0532 sigmaQC(counterDeg, counterLayer) = 0.05;
0533 }
0534 }
0535 }
0536 } else {
0537 float meanvarSL1 =
0538 sqrt(surv(0, 1) * surv(0, 1) + surv(1, 1) * surv(1, 1) + surv(2, 1) * surv(2, 1) + surv(3, 1) * surv(3, 1)) /
0539 10000.0;
0540 float meanvarSL2 = 0;
0541 if (surv.GetNrows() > 9)
0542 meanvarSL2 = sqrt(surv(8, 1) * surv(8, 1) + surv(9, 1) * surv(9, 1) + surv(10, 1) * surv(10, 1) +
0543 surv(11, 1) * surv(11, 1)) /
0544 10000.0;
0545 float meanvarSL3 =
0546 sqrt(surv(4, 1) * surv(4, 1) + surv(5, 1) * surv(5, 1) + surv(6, 1) * surv(6, 1) + surv(7, 1) * surv(7, 1)) /
0547 10000.0;
0548 for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
0549 for (int counterLayer = 0; counterLayer < 12; ++counterLayer) {
0550 if (st != 4 && counterLayer > 7)
0551 continue;
0552 float meanerror = 0;
0553 if (counterLayer < 4) {
0554 meanerror = meanvarSL1;
0555 } else if (counterLayer > 3 && counterLayer < 8) {
0556 meanerror = meanvarSL2;
0557 } else {
0558 meanerror = meanvarSL3;
0559 }
0560 if (counterDeg == 0) {
0561 sigmaQC(counterDeg, counterLayer) =
0562 sqrt((surv(counterLayer, 1) / 10000.0) * (surv(counterLayer, 1) / 10000.0) + meanerror * meanerror);
0563 } else if (counterDeg < 3) {
0564 sigmaQC(counterDeg, counterLayer) = 0.05;
0565 } else {
0566 sigmaQC(counterDeg, counterLayer) = 0.05;
0567 }
0568 }
0569 }
0570 double **Eta = new double *[12];
0571 for (int counterLayer = 0; counterLayer < 12; counterLayer++) {
0572 if (counterLayer > 7 && st == 4)
0573 continue;
0574 Eta[counterLayer] = new double[5];
0575 for (int counterLayer2 = 0; counterLayer2 < 12; counterLayer2++) {
0576 if (counterLayer > 7 && st == 4)
0577 continue;
0578 if ((counterLayer2 < 4 && counterLayer < 4) || (counterLayer2 > 3 && counterLayer > 3)) {
0579 if (counterLayer == counterLayer2) {
0580 Eta[counterLayer][counterLayer2] = 3.0 / (4.0);
0581 } else {
0582 Eta[counterLayer][counterLayer2] = -1.0 / (4.0);
0583 }
0584 } else {
0585 Eta[counterLayer][counterLayer2] = 0.0;
0586 }
0587 }
0588 }
0589
0590 for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
0591 for (int counterLayer = 0; counterLayer < 12; counterLayer++) {
0592 if (counterLayer > 7 && st == 4)
0593 continue;
0594 for (int counterLayer2 = 0; counterLayer2 < 12; counterLayer2++) {
0595 if (counterLayer2 > 7 && st == 4)
0596 continue;
0597 for (int counterLayer3 = 0; counterLayer3 < 12; counterLayer3++) {
0598 if (counterLayer3 > 7 && st == 4)
0599 continue;
0600 matrix(5 * counterLayer2 + counterDeg, 5 * counterLayer3 + counterDeg) +=
0601 Eta[counterLayer][counterLayer2] * Eta[counterLayer][counterLayer3] /
0602 (sigmaQC(counterDeg, counterLayer) * sigmaQC(counterDeg, counterLayer));
0603 }
0604 }
0605 }
0606 }
0607 }
0608 return matrix;
0609 }
0610
0611 TMatrixD DTMuonMillepede::getbqcMatrix(int wh, int st, int se) {
0612 TMatrixD surv = myPG->giveQCCal(wh, st, se);
0613 TMatrixD ResQC(5, 12);
0614 TMatrixD sigmaQC(5, 12);
0615 TMatrixD matrix(60, 1);
0616 if (st == 4)
0617 matrix.ResizeTo(40, 1);
0618
0619 if (surv.GetNrows() < 7) {
0620 for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
0621 for (int counterLayer = 0; counterLayer < 12; ++counterLayer) {
0622 if (st != 4 && counterLayer > 7)
0623 continue;
0624 if (counterDeg == 0) {
0625 sigmaQC(counterDeg, counterLayer) = 0.05;
0626 } else if (counterDeg < 3) {
0627 sigmaQC(counterDeg, counterLayer) = 0.05;
0628 } else {
0629 sigmaQC(counterDeg, counterLayer) = 0.05;
0630 }
0631 }
0632 }
0633 } else {
0634 for (int counterChamber = 0; counterChamber < 12; ++counterChamber) {
0635 ResQC(0, counterChamber) = -surv(counterChamber, 0) / 10000.0;
0636 }
0637 float meanvarSL1 =
0638 sqrt(surv(0, 1) * surv(0, 1) + surv(1, 1) * surv(1, 1) + surv(2, 1) * surv(2, 1) + surv(3, 1) * surv(3, 1)) /
0639 10000.0;
0640 float meanvarSL2 = 0;
0641 if (surv.GetNrows() > 9) {
0642 meanvarSL2 = sqrt(surv(8, 1) * surv(8, 1) + surv(9, 1) * surv(9, 1) + surv(10, 1) * surv(10, 1) +
0643 surv(11, 1) * surv(11, 1)) /
0644 10000.0;
0645 }
0646 float meanvarSL3 =
0647 sqrt(surv(4, 1) * surv(4, 1) + surv(5, 1) * surv(5, 1) + surv(6, 1) * surv(6, 1) + surv(7, 1) * surv(7, 1)) /
0648 10000.0;
0649 for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
0650 for (int counterLayer = 0; counterLayer < 12; ++counterLayer) {
0651 if (st != 4 && counterLayer > 7)
0652 continue;
0653 float meanerror = 0;
0654 if (counterLayer < 4) {
0655 meanerror = meanvarSL1;
0656 } else if (counterLayer > 3 && counterLayer < 8) {
0657 meanerror = meanvarSL2;
0658 } else {
0659 meanerror = meanvarSL3;
0660 }
0661 if (counterDeg == 0) {
0662 sigmaQC(counterDeg, counterLayer) =
0663 sqrt((surv(counterLayer, 1) / 10000.0) * (surv(counterLayer, 1) / 10000.0) + meanerror * meanerror);
0664 } else if (counterDeg < 3) {
0665 sigmaQC(counterDeg, counterLayer) = 0.05;
0666 } else {
0667 sigmaQC(counterDeg, counterLayer) = 0.05;
0668 }
0669 }
0670 }
0671 double **Eta = new double *[12];
0672 for (int counterLayer = 0; counterLayer < 12; counterLayer++) {
0673 if (counterLayer > 7 && st == 4)
0674 continue;
0675 Eta[counterLayer] = new double[5];
0676 for (int counterLayer2 = 0; counterLayer2 < 12; counterLayer2++) {
0677 if (counterLayer > 7 && st == 4)
0678 continue;
0679 if ((counterLayer2 < 4 && counterLayer < 4) || (counterLayer2 > 3 && counterLayer > 3)) {
0680 if (counterLayer == counterLayer2) {
0681 Eta[counterLayer][counterLayer2] = 3.0 / (4.0);
0682 } else {
0683 Eta[counterLayer][counterLayer2] = -1.0 / (4.0);
0684 }
0685 } else {
0686 Eta[counterLayer][counterLayer2] = 0.0;
0687 }
0688 }
0689 }
0690
0691 for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
0692 for (int counterLayer = 0; counterLayer < 12; counterLayer++) {
0693 if (counterLayer > 7 && st == 4)
0694 continue;
0695 for (int counterLayer2 = 0; counterLayer2 < 12; counterLayer2++) {
0696 if (counterLayer2 > 7 && st == 4)
0697 continue;
0698 float mean = 0;
0699 if (counterDeg != 0) {
0700 if (counterLayer < 4) {
0701 mean = (ResQC(counterDeg, 0) + ResQC(counterDeg, 1) + ResQC(counterDeg, 2) + ResQC(counterDeg, 3)) / 4.0;
0702 } else if (counterLayer > 3 && counterLayer < 8) {
0703 mean = (ResQC(counterDeg, 4) + ResQC(counterDeg, 5) + ResQC(counterDeg, 6) + ResQC(counterDeg, 7)) / 4.0;
0704 } else {
0705 mean =
0706 (ResQC(counterDeg, 8) + ResQC(counterDeg, 9) + ResQC(counterDeg, 10) + ResQC(counterDeg, 11)) / 4.0;
0707 }
0708 }
0709 matrix(5 * counterLayer2 + counterDeg, 0) +=
0710 Eta[counterLayer][counterLayer2] * (ResQC(counterDeg, counterLayer) - mean) /
0711 (sigmaQC(counterDeg, counterLayer) * sigmaQC(counterDeg, counterLayer));
0712 }
0713 }
0714 }
0715 }
0716 return matrix;
0717 }
0718
0719
0720
0721
0722
0723
0724
0725
0726
0727
0728
0729
0730
0731
0732
0733
0734
0735
0736
0737
0738
0739
0740
0741
0742
0743 TMatrixD DTMuonMillepede::getCsurveyMatrix(int wh, int st, int se) {
0744 int size = 60;
0745 if (st == 4)
0746 size = 40;
0747
0748 TMatrixD matrix(size + 6, size + 6);
0749
0750 float error[2][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
0751 {0.05, 0.05, 0.05, 0.005, 0.005, 0.005}};
0752 for (int counterLayer = 0; counterLayer < size / 5; counterLayer++) {
0753 for (int counterLayer2 = 0; counterLayer2 < size / 5; counterLayer2++) {
0754 int sl1 = counterLayer / 4;
0755 int sl2 = counterLayer2 / 4;
0756 if (sl1 == sl2) {
0757 for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
0758 int counterDegAux = counterDeg + 1;
0759 if (counterLayer < 8 && counterDeg == 1)
0760 counterDegAux = 0;
0761 int sl = (sl1 + 1) / 2;
0762 matrix(5 * counterLayer + counterDeg, 5 * counterLayer2 + counterDeg) =
0763 1.0 / (16.0 * error[sl][counterDegAux] * error[sl][counterDegAux]);
0764 }
0765 }
0766 }
0767 }
0768
0769 return matrix;
0770 }
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782
0783
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794
0795 TMatrixD DTMuonMillepede::getbsurveyMatrix(int wh, int st, int se) {
0796 TMatrixD survey = myPG->giveSurvey(wh, st, se);
0797 int size = 60;
0798 if (st == 4)
0799 size = 40;
0800
0801 TMatrixD matrix(size + 6, 1);
0802
0803 float error[2][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
0804 {0.05, 0.05, 0.05, 0.005, 0.005, 0.005}};
0805 for (int counterLayer = 0; counterLayer < size / 5; counterLayer++) {
0806 for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
0807 int counterDegAux = counterDeg + 1;
0808 if (counterLayer < 8 && counterDeg == 1)
0809 counterDegAux = 0;
0810 int superlayerAux = counterLayer / 4;
0811 int sl = (superlayerAux + 1) / 2;
0812 matrix(5 * counterLayer + counterDeg, 0) =
0813 survey(counterDegAux, superlayerAux) / (4.0 * error[sl][counterDegAux] * error[sl][counterDegAux]);
0814 }
0815 }
0816
0817 return matrix;
0818 }
0819
0820 TMatrixD DTMuonMillepede::prepareForLagrange(const TMatrixD &m) {
0821 TMatrixD updatedMatrix = m;
0822 updatedMatrix.ResizeTo(m.GetNrows() + 6, m.GetNcols() + 6);
0823 return updatedMatrix;
0824 }
0825
0826 void DTMuonMillepede::setBranchTree() {
0827 ttreeOutput = new TTree("DTLocalMillepede", "DTLocalMillepede");
0828
0829 ttreeOutput->Branch("wh", &whC, "wh/I");
0830 ttreeOutput->Branch("st", &stC, "st/I");
0831 ttreeOutput->Branch("sr", &srC, "sr/I");
0832 ttreeOutput->Branch("cov", cov, "cov[60][60]/F");
0833 ttreeOutput->Branch("sl", slC, "sl[12]/I");
0834 ttreeOutput->Branch("la", laC, "la[12]/I");
0835 ttreeOutput->Branch("dx", dx, "dx[12]/F");
0836 ttreeOutput->Branch("dy", dy, "dy[12]/F");
0837 ttreeOutput->Branch("dz", dz, "dz[12]/F");
0838 ttreeOutput->Branch("phix", phix, "phix[12]/F");
0839 ttreeOutput->Branch("phiy", phiy, "phiy[12]/F");
0840 ttreeOutput->Branch("phiz", phiz, "phiz[12]/F");
0841 }