Back to home page

Project CMSSW displayed by LXR
 
 

    

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   //Interface to Survey information
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   //Init Matrizes
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   //Run over the TTree
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       //Basic cuts
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   // My Final Calculation and Constraints
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   //Final Calculation and Constraints
0277   //   for(int wheel = -2; wheel < 3; ++wheel) {
0278   //     for(int station = 1; station < 5; ++station) {
0279   //       for(int sector = 1; sector < 15; ++sector) {
0280   //         if(sector > 12 && station != 4) continue;
0281   //         TMatrixD Ctr = prepareForLagrange(*C[wheel+2][station-1][sector-1]);
0282   //         TMatrixD btr = prepareForLagrange(*b[wheel+2][station-1][sector-1]);
0283   //         TMatrixD Cqc = prepareForLagrange(getCqcMatrix(wheel, station, sector));
0284   //         TMatrixD bqc = prepareForLagrange(getbqcMatrix(wheel, station, sector));
0285   //         TMatrixD Csurvey = prepareForLagrange(getCsurveyMatrix(wheel, station, sector));
0286   //         TMatrixD bsurvey = prepareForLagrange(getbsurveyMatrix(wheel, station, sector));
0287   //         TMatrixD Clag = getLagMatrix(wheel, station, sector);
0288   //         TMatrixD SC = Ctr + Cqc + Csurvey + Clag;
0289   //         TMatrixD Sb = btr + bqc + bsurvey;
0290   //
0291   //         SC.Invert();
0292 
0293   //         TMatrixD solution = SC*Sb;
0294   //         for(int icounter = 0; icounter < SC.GetNrows(); icounter++) {
0295   //           for(int jcounter = 0; jcounter < SC.GetNrows(); jcounter++) {
0296   //             cov[icounter][jcounter] = SC(icounter, jcounter);
0297   //           }
0298   //         }
0299   //         whC = wheel;
0300   //         stC = station;
0301   //         srC = sector;
0302 
0303   //         for(int c = 0; c < 60; ++c) {
0304   //           for(int s = 0; s < 60; ++s) {
0305   //             cov[c][s] = SC(c, s);
0306   //           }
0307   //         }
0308 
0309   //         for(int counterLayer = 0; counterLayer < 12; ++counterLayer) {
0310   //           for(int counterDeg = 0; counterDeg < 5; ++counterDeg) {
0311   //             if(counterLayer > 7 && station == 4) continue;
0312   //             if(counterLayer < 4) {
0313   //               slC[counterLayer] = 1;
0314   //               laC[counterLayer] = counterLayer+1;
0315   //             } else if(counterLayer > 3 && counterLayer < 8) {
0316   //               slC[counterLayer] = 3;
0317   //               laC[counterLayer] = counterLayer-3;
0318   //             } else {
0319   //               slC[counterLayer] = 2;
0320   //               laC[counterLayer] = counterLayer-7;
0321   //             }
0322   //             if(counterLayer < 8) {
0323   //               dx[counterLayer] = solution(counterLayer*5, 0);
0324   //               dy[counterLayer] = 0.0;
0325   //             } else {
0326   //               dx[counterLayer] = 0.0;
0327   //               dy[counterLayer] = solution(counterLayer*5, 0);
0328   //             }
0329   //             dz[counterLayer] = solution(counterLayer*5+1, 0);
0330   //             phix[counterLayer] = solution(counterLayer*5+2, 0);
0331   //             phiy[counterLayer] = solution(counterLayer*5+3, 0);
0332   //             phiz[counterLayer] = solution(counterLayer*5+4, 0);
0333   //           }
0334   //         }
0335   //         ttreeOutput->Fill();
0336   //       }
0337   //     }
0338   //   }
0339   //  f->Write();
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 // TMatrixD DTMuonMillepede::getCsurveyMatrix(int wh, int st, int se) {
0720 
0721 //   int size = 60;
0722 //   if(st == 4) size = 40;
0723 
0724 //   TMatrixD matrix(size+6, size+6);
0725 //   //Careful with the sign
0726 //   float error[] = {0.05, 0.05, 0.05, 0.005, 0.005, 0.005};
0727 //   for(int counterLayer = 0; counterLayer < size/5; counterLayer++) {
0728 //     float k = 1;
0729 //     if(counterLayer < 4) k = -1;
0730 //     for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
0731 //       for(int counterLayer2 = 0; counterLayer2 < size/5; counterLayer2++) {
0732 //         float k2 = 1;
0733 //         if(counterLayer2 < 4) k2 = -1;
0734 //         matrix(5*counterLayer+counterDeg, 5*counterLayer2+counterDeg) =  k2*k/(16.0*error[counterDeg]*error[counterDeg]);
0735 //       }
0736 //     }
0737 //   }
0738 
0739 //   return matrix;
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   //Careful with the sign
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 // TMatrixD DTMuonMillepede::getbsurveyMatrix(int wh, int st, int se) {
0773 
0774 //   TMatrixD survey = myPG->giveSurvey(wh, st, se);
0775 //   int size = 60;
0776 //   if(st == 4) size = 40;
0777 
0778 //   TMatrixD matrix(size+6, 1);
0779 //   //Careful with the sign
0780 //   float error[] = {0.05, 0.05, 0.05, 0.005, 0.005, 0.005};
0781 //   for(int counterLayer = 0; counterLayer < size/5; counterLayer++) {
0782 //     float k = 1;
0783 //     if(counterLayer < 4) k = -1;
0784 //     for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
0785 //       int counterDegAux = counterDeg+1;
0786 //       if(counterLayer < 8 && counterDeg == 1) counterDegAux = 0;
0787 //       matrix(5*counterLayer+counterDeg, 0) =  k*survey(counterDegAux, 0)/(16.0*error[counterDeg]*error[counterDeg]);
0788 //     }
0789 //   }
0790 
0791 //   return matrix;
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   //Careful with the sign
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!