MuonAlignmentAlgorithms/src/DTMuonSLToSL.cc

0001 #include "Alignment/MuonAlignmentAlgorithms/interface/DTMuonSLToSL.h"
0002
0003 DTMuonSLToSL::DTMuonSLToSL(std::string path, int n_files, float MaxPt, float MinPt, TFile *f_) {
0004   ntuplePath = path;
0005   numberOfRootFiles = n_files;
0006
0007   f = f_;
0008
0009   ptMax = MaxPt;
0010   ptMin = MinPt;
0011
0012   setBranchTree();
0013
0014   initNTuples(0);
0015
0016   calculationSLToSL();
0017 }
0018
0019 DTMuonSLToSL::~DTMuonSLToSL() {}
0020
0021 void DTMuonSLToSL::calculationSLToSL() {
0022   TMatrixD ****C13 = new TMatrixD ***[5];
0023   TMatrixD ****b13 = new TMatrixD ***[5];
0024   TMatrixD ****C31 = new TMatrixD ***[5];
0025   TMatrixD ****b31 = new TMatrixD ***[5];
0026
0027   for (int whI = -2; whI < 3; ++whI) {
0028     C13[whI + 2] = new TMatrixD **[4];
0029     b13[whI + 2] = new TMatrixD **[4];
0030     C31[whI + 2] = new TMatrixD **[4];
0031     b31[whI + 2] = new TMatrixD **[4];
0032     for (int stI = 1; stI < 5; ++stI) {
0033       C13[whI + 2][stI - 1] = new TMatrixD *[14];
0034       b13[whI + 2][stI - 1] = new TMatrixD *[14];
0035       C31[whI + 2][stI - 1] = new TMatrixD *[14];
0036       b31[whI + 2][stI - 1] = new TMatrixD *[14];
0037       for (int seI = 1; seI < 15; ++seI) {
0038         if (seI > 12 && stI != 4)
0039           continue;
0040         C13[whI + 2][stI - 1][seI - 1] = new TMatrixD(3, 3);
0041         b13[whI + 2][stI - 1][seI - 1] = new TMatrixD(3, 1);
0042         C31[whI + 2][stI - 1][seI - 1] = new TMatrixD(3, 3);
0043         b31[whI + 2][stI - 1][seI - 1] = new TMatrixD(3, 1);
0044       }
0045     }
0046   }
0047
0048   //Run over the TTree
0049   Int_t nentries = (Int_t)tali->GetEntries();
0050   for (Int_t i = 0; i < nentries; i++) {
0051     tali->GetEntry(i);
0052     //Basic cuts
0053     if (pt > ptMax || pt < ptMin)
0054       continue;
0055
0056     bool repeatedHits = false;
0057     for (int counter = 0; counter < nseg; ++counter) {
0058       //Make sure there are no repeated hits
0059       for (int counterHi = 0; counterHi < nhits[counter]; counterHi++) {
0060         for (int counterHj = 0; counterHj < nhits[counter]; counterHj++) {
0061           if (counterHi == counterHj)
0062             continue;
0063           if (zc[counter][counterHi] == zc[counter][counterHj]) {
0064             repeatedHits = true;
0065           }
0066         }
0067       }
0068       if (repeatedHits == true)
0069         continue;
0070
0071       float x_13 = xSlSL3[counter];
0072       float xp_13 = xSL1SL3[counter];
0073       float x_31 = xSlSL1[counter];
0074       float xp_31 = xSL3SL1[counter];
0075       //float tanphi = dxdzSl[counter];
0076       float tanphi_13 = dxdzSlSL1[counter];
0077       float tanphi_31 = dxdzSlSL3[counter];
0078       int wheel = wh[counter];
0079       int station = st[counter];
0080       int sector = sr[counter];
0081
0082       if (fabs(x_13 - xp_13) < 3 && fabs(x_31 - xp_31) && fabs(tanphi_13 - tanphi_31) < 0.06) {
0083         *(C13[wheel + 2][station - 1][sector - 1]) += returnCSLMatrix(x_13, xp_13, tanphi_13);
0084         *(b13[wheel + 2][station - 1][sector - 1]) += returnbSLMatrix(x_13, xp_13, tanphi_13);
0085
0086         *(C31[wheel + 2][station - 1][sector - 1]) += returnCSLMatrix(x_31, xp_31, tanphi_31);
0087         *(b31[wheel + 2][station - 1][sector - 1]) += returnbSLMatrix(x_31, xp_31, tanphi_31);
0088       }
0089     }
0090   }
0091
0092   for (int wheel = -2; wheel < 3; ++wheel) {
0093     for (int station = 1; station < 5; ++station) {
0094       for (int sector = 1; sector < 15; ++sector) {
0095         if (sector > 12 && station != 4)
0096           continue;
0097         TMatrixD solution13(3, 1);
0098         TMatrixD solution31(3, 1);
0099         TMatrixD C31_copy = *(C31[wheel + 2][station - 1][sector - 1]);
0100         TMatrixD C13_copy = *(C13[wheel + 2][station - 1][sector - 1]);
0101         TMatrixD b31_copy = *(b31[wheel + 2][station - 1][sector - 1]);
0102         TMatrixD b13_copy = *(b13[wheel + 2][station - 1][sector - 1]);
0103
0104         C31_copy.Invert();
0105         C13_copy.Invert();
0106         solution13 = C13_copy * b13_copy;
0107         solution31 = C31_copy * b31_copy;
0108         whC = wheel;
0109         stC = station;
0110         srC = sector;
0111         dx = solution13(0, 0);
0112         dz = solution13(1, 0);
0113         phiy = solution13(2, 0);
0114         for (int c = 0; c < 3; ++c) {
0115           for (int s = 0; s < 3; ++s) {
0116             cov[c][s] = C13_copy(c, s);
0117           }
0118         }
0119         ttreeOutput->Fill();
0120       }
0121     }
0122   }
0123   f->Write();
0124 }
0125
0126 TMatrixD DTMuonSLToSL::returnCSLMatrix(float x, float xp, float tanphi) {
0127   TMatrixD matrix(3, 3);
0128
0129   matrix(0, 0) = 1.0;
0130   matrix(1, 0) = tanphi;
0131   matrix(0, 1) = tanphi;
0132   matrix(1, 1) = tanphi * tanphi;
0133   matrix(0, 2) = tanphi * xp;
0134   matrix(2, 0) = tanphi * xp;
0135   matrix(2, 2) = tanphi * tanphi * xp * xp;
0136   matrix(2, 1) = tanphi * tanphi * xp;
0137   matrix(1, 2) = tanphi * tanphi * xp;
0138
0139   return matrix;
0140 }
0141
0142 TMatrixD DTMuonSLToSL::returnbSLMatrix(float x, float xp, float tanphi) {
0143   TMatrixD matrix(3, 1);
0144
0145   matrix(0, 0) = -(x - xp);
0146   matrix(1, 0) = -(x - xp) * tanphi;
0147   matrix(2, 0) = -(x - xp) * tanphi * xp;
0148
0149   return matrix;
0150 }
0151
0152 void DTMuonSLToSL::setBranchTree() {
0153   ttreeOutput = new TTree("DTSLToSLResult", "DTSLToSLResult");
0154
0155   ttreeOutput->Branch("wh", &whC, "wh/F");
0156   ttreeOutput->Branch("st", &stC, "st/F");
0157   ttreeOutput->Branch("sr", &srC, "sr/F");
0158   ttreeOutput->Branch("dx", &dx, "dx/F");
0159   ttreeOutput->Branch("dz", &dz, "dz/F");
0160   ttreeOutput->Branch("phiy", &phiy, "phiy/F");
0161   ttreeOutput->Branch("cov", cov, "cov[3][3]/F");
0162 }