Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:56:48 

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