CMS 3D CMS Logo
DTMuonMillepede.cc
Go to the documentation of this file.
1 #include "Alignment/MuonAlignmentAlgorithms/interface/DTMuonMillepede.h"
2 
3 #include <array>
4 #include <iostream>
5 
6 DTMuonMillepede::DTMuonMillepede(std::string path,
7  int n_files,
8  float MaxPt,
9  float MinPt,
10  int nPhihits,
11  int nThetahits,
12  int workingmode,
13  int nMtxSection) {
14  ntuplePath = path;
15  numberOfRootFiles = n_files;
16 
17  ptMax = MaxPt;
18  ptMin = MinPt;
19 
20  nPhiHits = nPhihits;
21  nThetaHits = nThetahits;
22 
23  TDirectory *dirSave = gDirectory;
24 
25  //Interface to Survey information
26  myPG = new ReadPGInfo("./InternalData2009.root");
27 
28  f = new TFile("./LocalMillepedeResults.root", "RECREATE");
29  f->cd();
30 
31  setBranchTree();
32 
33  initNTuples(nMtxSection);
34 
35  calculationMillepede(workingmode);
36 
37  f->Write();
38  f->Close();
39 
40  dirSave->cd();
41 }
42 
43 DTMuonMillepede::~DTMuonMillepede() {}
44 
45 void DTMuonMillepede::calculationMillepede(int workingmode) {
46  //Init Matrizes
47  TMatrixD ****C = new TMatrixD ***[5];
48  TMatrixD ****b = new TMatrixD ***[5];
49 
50  for (int whI = -2; whI < 3; ++whI) {
51  C[whI + 2] = new TMatrixD **[4];
52  b[whI + 2] = new TMatrixD **[4];
53  for (int stI = 1; stI < 5; ++stI) {
54  C[whI + 2][stI - 1] = new TMatrixD *[14];
55  b[whI + 2][stI - 1] = new TMatrixD *[14];
56  for (int seI = 1; seI < 15; ++seI) {
57  if (seI > 12 && stI != 4)
58  continue;
59  if (stI == 4) {
60  C[whI + 2][stI - 1][seI - 1] = new TMatrixD(24, 24);
61  b[whI + 2][stI - 1][seI - 1] = new TMatrixD(24, 1);
62  } else {
63  C[whI + 2][stI - 1][seI - 1] = new TMatrixD(60, 60);
64  b[whI + 2][stI - 1][seI - 1] = new TMatrixD(60, 1);
65  }
66  }
67  }
68  }
69 
70  //Run over the TTree
71  if (workingmode <= 3) {
72  Int_t nentries = (Int_t)tali->GetEntries();
73  for (Int_t i = 0; i < nentries; i++) {
74  tali->GetEntry(i);
75 
76  if (i % 100000 == 0)
77  std::cout << "Analyzing track number " << i << std::endl;
78 
79  //Basic cuts
80  if (pt > ptMax || pt < ptMin)
81  continue;
82 
83  for (int counter = 0; counter < nseg; ++counter) {
84  if (nphihits[counter] < nPhiHits || (nthetahits[counter] < nThetaHits && st[counter] < 4))
85  continue;
86 
87  TMatrixD A(12, 60);
88  TMatrixD R(12, 1);
89  TMatrixD E(12, 12);
90  TMatrixD B(12, 4);
91  TMatrixD I(12, 12);
92  if (st[counter] == 4) {
93  A.ResizeTo(8, 24);
94  R.ResizeTo(8, 1);
95  E.ResizeTo(8, 8);
96  B.ResizeTo(8, 2);
97  I.ResizeTo(8, 8);
98  }
99 
100  for (int counterHi = 0; counterHi < nhits[counter]; counterHi++) {
101  int row = 0;
102  if (sl[counter][counterHi] == 3) {
103  row = 4 + (la[counter][counterHi] - 1);
104  } else if (sl[counter][counterHi] == 2) {
105  row = 8 + (la[counter][counterHi] - 1);
106  } else {
107  row = (la[counter][counterHi] - 1);
108  }
109 
110  float x = xc[counter][counterHi];
111  float y = yc[counter][counterHi];
112  float xp = xcp[counter][counterHi];
113  float yp = ycp[counter][counterHi];
114  float zp = zc[counter][counterHi];
115  float error = ex[counter][counterHi];
116  float dxdz = dxdzSl[counter];
117  float dydz = dydzSl[counter];
118 
119  if (st[counter] != 4) {
120  if (sl[counter][counterHi] == 2) {
121  A(row, row * 5) = -1.0;
122  A(row, row * 5 + 1) = dydz;
123  A(row, row * 5 + 2) = y * dydz;
124  A(row, row * 5 + 3) = -x * dydz;
125  A(row, row * 5 + 4) = -x;
126  R(row, 0) = y - yp;
127  } else {
128  A(row, row * 5 + 0) = -1.0;
129  A(row, row * 5 + 1) = dxdz;
130  A(row, row * 5 + 2) = y * dxdz;
131  A(row, row * 5 + 3) = -x * dxdz;
132  A(row, row * 5 + 4) = y;
133  R(row, 0) = x - xp;
134  }
135  } else {
136  if (sl[counter][counterHi] != 2) {
137  A(row, row * 3 + 0) = -1.0;
138  A(row, row * 3 + 1) = dxdz;
139  A(row, row * 3 + 2) = -x * dxdz;
140  R(row, 0) = x - xp;
141  }
142  }
143 
144  E(row, row) = 1.0 / error;
145  I(row, row) = 1.0;
146 
147  if (sl[counter][counterHi] != 2) {
148  B(row, 0) = 1.0;
149  B(row, 1) = zp;
150  } else {
151  B(row, 2) = 1.0;
152  B(row, 3) = zp;
153  }
154  }
155 
156  TMatrixD Bt = B;
157  Bt.T();
158  TMatrixD At = A;
159  At.T();
160  TMatrixD gamma = (Bt * E * B);
161  gamma.Invert();
162  *(C[wh[counter] + 2][st[counter] - 1][sr[counter] - 1]) += At * E * (B * gamma * Bt * E * A - A);
163  *(b[wh[counter] + 2][st[counter] - 1][sr[counter] - 1]) += At * E * (B * gamma * Bt * E * I - I) * R;
164  }
165  }
166  }
167 
168  if (workingmode == 3)
169  for (int wheel = -2; wheel < 3; ++wheel)
170  for (int station = 1; station < 5; ++station)
171  for (int sector = 1; sector < 15; ++sector) {
172  if (sector > 12 && station != 4)
173  continue;
174 
175  TMatrixD Ctr = *C[wheel + 2][station - 1][sector - 1];
176  TMatrixD btr = *b[wheel + 2][station - 1][sector - 1];
177 
178  TString CtrName = "Ctr_";
179  CtrName += wheel;
180  CtrName += "_";
181  CtrName += station;
182  CtrName += "_";
183  CtrName += sector;
184  Ctr.Write(CtrName);
185 
186  TString btrName = "btr_";
187  btrName += wheel;
188  btrName += "_";
189  btrName += station;
190  btrName += "_";
191  btrName += sector;
192  btr.Write(btrName);
193  }
194 
195  // My Final Calculation and Constraints
196  if (workingmode % 2 == 0) {
197  for (int wheel = -2; wheel < 3; ++wheel)
198  for (int station = 1; station < 5; ++station)
199  for (int sector = 1; sector < 15; ++sector) {
200  if (sector > 12 && station != 4)
201  continue;
202 
203  if (workingmode == 4) {
204  for (int mf = -1; mf <= -1; mf++) {
205  TMatrixD Ch = getMatrixFromFile("Ctr_", wheel, station, sector, mf);
206  TMatrixD bh = getMatrixFromFile("btr_", wheel, station, sector, mf);
207 
208  *(C[wheel + 2][station - 1][sector - 1]) += Ch;
209  *(b[wheel + 2][station - 1][sector - 1]) += bh;
210  }
211  }
212 
213  TMatrixD Ctr = *C[wheel + 2][station - 1][sector - 1];
214  TMatrixD btr = *b[wheel + 2][station - 1][sector - 1];
215 
216  TMatrixD Ccs = getCcsMatrix(wheel, station, sector);
217  TMatrixD bcs = getbcsMatrix(wheel, station, sector);
218 
219  TMatrixD SC = Ctr + Ccs;
220  TMatrixD Sb = btr + bcs;
221 
222  SC.Invert();
223 
224  TMatrixD solution = SC * Sb;
225  for (int icounter = 0; icounter < SC.GetNrows(); icounter++) {
226  for (int jcounter = 0; jcounter < SC.GetNrows(); jcounter++) {
227  cov[icounter][jcounter] = SC(icounter, jcounter);
228  }
229  }
230 
231  int nLayer = 12, nDeg = 5;
232  if (station == 4) {
233  nLayer = 8;
234  nDeg = 3;
235  }
236  for (int counterLayer = 0; counterLayer < nLayer; ++counterLayer)
237  for (int counterDeg = 0; counterDeg < nDeg; ++counterDeg) {
238  if (counterLayer < 4) {
239  slC[counterLayer] = 1;
240  laC[counterLayer] = counterLayer + 1;
241  } else if (counterLayer > 3 && counterLayer < 8) {
242  slC[counterLayer] = 3;
243  laC[counterLayer] = counterLayer - 3;
244  } else {
245  slC[counterLayer] = 2;
246  laC[counterLayer] = counterLayer - 7;
247  }
248 
249  if (counterLayer < 8) {
250  dx[counterLayer] = solution(counterLayer * nDeg, 0);
251  dy[counterLayer] = 0.0;
252  } else {
253  dx[counterLayer] = 0.0;
254  dy[counterLayer] = solution(counterLayer * nDeg, 0);
255  }
256  dz[counterLayer] = solution(counterLayer * nDeg + 1, 0);
257 
258  if (station != 4) {
259  phix[counterLayer] = solution(counterLayer * nDeg + 2, 0);
260  phiy[counterLayer] = solution(counterLayer * nDeg + 3, 0);
261  phiz[counterLayer] = solution(counterLayer * nDeg + 4, 0);
262  } else {
263  phix[counterLayer] = 0.;
264  phiy[counterLayer] = solution(counterLayer * nDeg + 2, 0);
265  phiz[counterLayer] = 0.;
266  }
267  }
268 
269  whC = wheel;
270  stC = station;
271  srC = sector;
272 
273  ttreeOutput->Fill();
274  }
275  }
276 
277  //Final Calculation and Constraints
278  // for(int wheel = -2; wheel < 3; ++wheel) {
279  // for(int station = 1; station < 5; ++station) {
280  // for(int sector = 1; sector < 15; ++sector) {
281  // if(sector > 12 && station != 4) continue;
282  // TMatrixD Ctr = prepareForLagrange(*C[wheel+2][station-1][sector-1]);
283  // TMatrixD btr = prepareForLagrange(*b[wheel+2][station-1][sector-1]);
284  // TMatrixD Cqc = prepareForLagrange(getCqcMatrix(wheel, station, sector));
285  // TMatrixD bqc = prepareForLagrange(getbqcMatrix(wheel, station, sector));
286  // TMatrixD Csurvey = prepareForLagrange(getCsurveyMatrix(wheel, station, sector));
287  // TMatrixD bsurvey = prepareForLagrange(getbsurveyMatrix(wheel, station, sector));
288  // TMatrixD Clag = getLagMatrix(wheel, station, sector);
289  // TMatrixD SC = Ctr + Cqc + Csurvey + Clag;
290  // TMatrixD Sb = btr + bqc + bsurvey;
291  //
292  // SC.Invert();
293 
294  // TMatrixD solution = SC*Sb;
295  // for(int icounter = 0; icounter < SC.GetNrows(); icounter++) {
296  // for(int jcounter = 0; jcounter < SC.GetNrows(); jcounter++) {
297  // cov[icounter][jcounter] = SC(icounter, jcounter);
298  // }
299  // }
300  // whC = wheel;
301  // stC = station;
302  // srC = sector;
303 
304  // for(int c = 0; c < 60; ++c) {
305  // for(int s = 0; s < 60; ++s) {
306  // cov[c][s] = SC(c, s);
307  // }
308  // }
309 
310  // for(int counterLayer = 0; counterLayer < 12; ++counterLayer) {
311  // for(int counterDeg = 0; counterDeg < 5; ++counterDeg) {
312  // if(counterLayer > 7 && station == 4) continue;
313  // if(counterLayer < 4) {
314  // slC[counterLayer] = 1;
315  // laC[counterLayer] = counterLayer+1;
316  // } else if(counterLayer > 3 && counterLayer < 8) {
317  // slC[counterLayer] = 3;
318  // laC[counterLayer] = counterLayer-3;
319  // } else {
320  // slC[counterLayer] = 2;
321  // laC[counterLayer] = counterLayer-7;
322  // }
323  // if(counterLayer < 8) {
324  // dx[counterLayer] = solution(counterLayer*5, 0);
325  // dy[counterLayer] = 0.0;
326  // } else {
327  // dx[counterLayer] = 0.0;
328  // dy[counterLayer] = solution(counterLayer*5, 0);
329  // }
330  // dz[counterLayer] = solution(counterLayer*5+1, 0);
331  // phix[counterLayer] = solution(counterLayer*5+2, 0);
332  // phiy[counterLayer] = solution(counterLayer*5+3, 0);
333  // phiz[counterLayer] = solution(counterLayer*5+4, 0);
334  // }
335  // }
336  // ttreeOutput->Fill();
337  // }
338  // }
339  // }
340  // f->Write();
341 }
342 
343 TMatrixD DTMuonMillepede::getCcsMatrix(int wh, int st, int se) {
344  int size = 60;
345  if (st == 4)
346  size = 24;
347 
348  TMatrixD matrix(size, size);
349 
350  float Error[3][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
351  {0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005},
352  {0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005}};
353 
354  float TollerancyPosition = 0.01;
355  float TollerancyRotation = 0.0001;
356 
357  TMatrixD ConsQC = myPG->giveQCCal(wh, st, se);
358 
359  bool UseQC = true;
360  if (ConsQC.GetNrows() < 70)
361  UseQC = false;
362 
363  int nLayer = 12, nDeg = 5;
364  if (st == 4) {
365  nLayer = 8;
366  nDeg = 3;
367  }
368 
369  for (int la = 0; la < nLayer; la++)
370  for (int dg = 0; dg < nDeg; dg++) {
371  int index = la * nDeg + dg;
372 
373  int rdg = dg + 1;
374  if (la < 8 && rdg == 1)
375  rdg = 0;
376  if (st == 4 && rdg == 3)
377  rdg = 4;
378 
379  double ThisError2 = Error[la / 4][rdg] * Error[la / 4][rdg];
380  if (rdg < 2) {
381  if (UseQC) {
382  ThisError2 += ConsQC(la, 1) * ConsQC(la, 1);
383  } else {
384  ThisError2 += TollerancyPosition * TollerancyPosition;
385  }
386  } else if (rdg == 2) {
387  ThisError2 += TollerancyPosition * TollerancyPosition;
388  } else {
389  ThisError2 += TollerancyRotation * TollerancyRotation;
390  }
391 
392  matrix(index, index) = 1. / ThisError2;
393  }
394 
395  return matrix;
396 }
397 
398 TMatrixD DTMuonMillepede::getbcsMatrix(int wh, int st, int se) {
399  int size = 60;
400  if (st == 4)
401  size = 24;
402 
403  TMatrixD matrix(size, 1);
404 
405  float Error[3][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
406  {0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005},
407  {0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005}};
408 
409  float TollerancyPosition = 0.01;
410  float TollerancyRotation = 0.0001;
411 
412  TMatrixD ConsQC = myPG->giveQCCal(wh, st, se);
413 
414  bool UseQC = true;
415  if (ConsQC.GetNrows() < 70)
416  UseQC = false;
417 
418  TMatrixD Survey = myPG->giveSurvey(wh, st, se);
419 
420  int nLayer = 12, nDeg = 5;
421  if (st == 4) {
422  nLayer = 8;
423  nDeg = 3;
424  }
425 
426  for (int la = 0; la < nLayer; la++)
427  for (int dg = 0; dg < nDeg; dg++) {
428  int index = la * nDeg + dg;
429 
430  int rdg = dg + 1;
431  if (la < 8 && rdg == 1)
432  rdg = 0;
433  if (st == 4 && rdg == 3)
434  rdg = 4;
435 
436  double ThisError2 = Error[la / 4][rdg] * Error[la / 4][rdg];
437  if (rdg < 2) {
438  if (UseQC) {
439  ThisError2 += ConsQC(la, 1) * ConsQC(la, 1);
440  } else {
441  ThisError2 += TollerancyPosition * TollerancyPosition;
442  }
443  } else if (rdg == 2) {
444  ThisError2 += TollerancyPosition * TollerancyPosition;
445  } else {
446  ThisError2 += TollerancyRotation * TollerancyRotation;
447  }
448 
449  float Constraint = 0.;
450  if (la > 3)
451  Constraint += Survey(rdg, la / 4 - 1);
452  if (UseQC && rdg == 0)
453  Constraint += ConsQC(la, 0);
454  ;
455  if (UseQC && rdg == 1)
456  Constraint -= ConsQC(la, 0);
457 
458  matrix(index, 0) = Constraint / ThisError2;
459  }
460 
461  return matrix;
462 }
463 
464 TMatrixD DTMuonMillepede::getMatrixFromFile(const TString &Code, int wh, int st, int se, int mf) {
465  TString MtxFileName = "./LocalMillepedeMatrix_";
466  MtxFileName += mf;
467  MtxFileName += ".root";
468  if (mf == -1)
469  MtxFileName = "./LocalMillepedeMatrix.root";
470 
471  TDirectory *dirSave = gDirectory;
472 
473  TFile *MatrixFile = new TFile(MtxFileName);
474 
475  TString MtxName = Code;
476  MtxName += wh;
477  MtxName += "_";
478  MtxName += st;
479  MtxName += "_";
480  MtxName += se;
481  TMatrixD *ThisMtx = (TMatrixD *)MatrixFile->Get(MtxName);
482 
483  MatrixFile->Close();
484 
485  dirSave->cd();
486 
487  return *ThisMtx;
488 }
489 
490 TMatrixD DTMuonMillepede::getLagMatrix(int wh, int st, int se) {
491  TMatrixD matrix(60 + 6, 60 + 6);
492  if (st == 4)
493  matrix.ResizeTo(40 + 6, 40 + 6);
494 
495  for (int counterDeg = 0; counterDeg < 5; ++counterDeg) {
496  for (int counterLayer = 0; counterLayer < 12; ++counterLayer) {
497  if (st == 4) {
498  matrix(40 + counterDeg, counterLayer * 5 + counterDeg) = 10000.0;
499  matrix(counterLayer * 5 + counterDeg, 40 + counterDeg) = 10000.0;
500  } else {
501  int realCounter = counterDeg + 1;
502  if (counterDeg == 1 && counterLayer < 8) {
503  realCounter = counterDeg - 1;
504  }
505  matrix(counterLayer * 5 + counterDeg, 40 + realCounter) = 10000.0;
506  if ((realCounter == 0 && counterLayer > 7) || (realCounter == 1 && counterLayer < 7))
507  continue;
508  matrix(60 + realCounter, counterLayer * 5 + counterDeg) = 10000.0;
509  }
510  }
511  }
512  return matrix;
513 }
514 
515 TMatrixD DTMuonMillepede::getCqcMatrix(int wh, int st, int se) {
516  TMatrixD surv = myPG->giveQCCal(wh, st, se);
517  TMatrixD sigmaQC(5, 12);
518 
519  TMatrixD matrix(60, 60);
520  if (st == 4)
521  matrix.ResizeTo(40, 40);
522 
523  if (surv.GetNrows() < 7) {
524  for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
525  for (int counterLayer = 0; counterLayer < 12; ++counterLayer) {
526  if (st != 4 && counterLayer > 7)
527  continue;
528  if (counterDeg == 0) {
529  sigmaQC(counterDeg, counterLayer) = 0.05;
530  } else if (counterDeg < 3) {
531  sigmaQC(counterDeg, counterLayer) = 0.05;
532  } else {
533  sigmaQC(counterDeg, counterLayer) = 0.05;
534  }
535  }
536  }
537  } else {
538  float meanvarSL1 =
539  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)) /
540  10000.0;
541  float meanvarSL2 = 0;
542  if (surv.GetNrows() > 9)
543  meanvarSL2 = sqrt(surv(8, 1) * surv(8, 1) + surv(9, 1) * surv(9, 1) + surv(10, 1) * surv(10, 1) +
544  surv(11, 1) * surv(11, 1)) /
545  10000.0;
546  float meanvarSL3 =
547  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)) /
548  10000.0;
549  for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
550  for (int counterLayer = 0; counterLayer < 12; ++counterLayer) {
551  if (st != 4 && counterLayer > 7)
552  continue;
553  float meanerror = 0;
554  if (counterLayer < 4) {
555  meanerror = meanvarSL1;
556  } else if (counterLayer > 3 && counterLayer < 8) {
557  meanerror = meanvarSL2;
558  } else {
559  meanerror = meanvarSL3;
560  }
561  if (counterDeg == 0) {
562  sigmaQC(counterDeg, counterLayer) =
563  sqrt((surv(counterLayer, 1) / 10000.0) * (surv(counterLayer, 1) / 10000.0) + meanerror * meanerror);
564  } else if (counterDeg < 3) {
565  sigmaQC(counterDeg, counterLayer) = 0.05;
566  } else {
567  sigmaQC(counterDeg, counterLayer) = 0.05;
568  }
569  }
570  }
571 
572  std::array<std::array<double, 12>, 12> Eta{};
573  for (size_t counterLayer = 0; counterLayer < Eta.size(); counterLayer++) {
574  if (counterLayer > 7 && st == 4)
575  continue;
576  for (size_t counterLayer2 = 0; counterLayer2 < Eta[counterLayer].size(); counterLayer2++) {
577  if (counterLayer2 > 7 && st == 4)
578  continue;
579  if ((counterLayer2 < 4 && counterLayer < 4) || (counterLayer2 > 3 && counterLayer > 3)) {
580  if (counterLayer == counterLayer2) {
581  Eta[counterLayer][counterLayer2] = 3.0 / (4.0);
582  } else {
583  Eta[counterLayer][counterLayer2] = -1.0 / (4.0);
584  }
585  }
586  }
587  }
588 
589  for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
590  for (size_t counterLayer = 0; counterLayer < 12; counterLayer++) {
591  if (counterLayer > 7 && st == 4)
592  continue;
593  for (size_t counterLayer2 = 0; counterLayer2 < 12; counterLayer2++) {
594  if (counterLayer2 > 7 && st == 4)
595  continue;
596  for (size_t counterLayer3 = 0; counterLayer3 < 12; counterLayer3++) {
597  if (counterLayer3 > 7 && st == 4)
598  continue;
599  matrix(5 * counterLayer2 + counterDeg, 5 * counterLayer3 + counterDeg) +=
600  Eta[counterLayer][counterLayer2] * Eta[counterLayer][counterLayer3] /
601  (sigmaQC(counterDeg, counterLayer) * sigmaQC(counterDeg, counterLayer));
602  }
603  }
604  }
605  }
606  }
607  return matrix;
608 }
609 
610 TMatrixD DTMuonMillepede::getbqcMatrix(int wh, int st, int se) {
611  TMatrixD surv = myPG->giveQCCal(wh, st, se);
612  TMatrixD ResQC(5, 12);
613  TMatrixD sigmaQC(5, 12);
614  TMatrixD matrix(60, 1);
615  if (st == 4)
616  matrix.ResizeTo(40, 1);
617 
618  if (surv.GetNrows() < 7) {
619  for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
620  for (int counterLayer = 0; counterLayer < 12; ++counterLayer) {
621  if (st != 4 && counterLayer > 7)
622  continue;
623  if (counterDeg == 0) {
624  sigmaQC(counterDeg, counterLayer) = 0.05;
625  } else if (counterDeg < 3) {
626  sigmaQC(counterDeg, counterLayer) = 0.05;
627  } else {
628  sigmaQC(counterDeg, counterLayer) = 0.05;
629  }
630  }
631  }
632  } else {
633  for (int counterChamber = 0; counterChamber < 12; ++counterChamber) {
634  ResQC(0, counterChamber) = -surv(counterChamber, 0) / 10000.0;
635  }
636  float meanvarSL1 =
637  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)) /
638  10000.0;
639  float meanvarSL2 = 0;
640  if (surv.GetNrows() > 9) {
641  meanvarSL2 = sqrt(surv(8, 1) * surv(8, 1) + surv(9, 1) * surv(9, 1) + surv(10, 1) * surv(10, 1) +
642  surv(11, 1) * surv(11, 1)) /
643  10000.0;
644  }
645  float meanvarSL3 =
646  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)) /
647  10000.0;
648  for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
649  for (int counterLayer = 0; counterLayer < 12; ++counterLayer) {
650  if (st != 4 && counterLayer > 7)
651  continue;
652  float meanerror = 0;
653  if (counterLayer < 4) {
654  meanerror = meanvarSL1;
655  } else if (counterLayer > 3 && counterLayer < 8) {
656  meanerror = meanvarSL2;
657  } else {
658  meanerror = meanvarSL3;
659  }
660  if (counterDeg == 0) {
661  sigmaQC(counterDeg, counterLayer) =
662  sqrt((surv(counterLayer, 1) / 10000.0) * (surv(counterLayer, 1) / 10000.0) + meanerror * meanerror);
663  } else if (counterDeg < 3) {
664  sigmaQC(counterDeg, counterLayer) = 0.05;
665  } else {
666  sigmaQC(counterDeg, counterLayer) = 0.05;
667  }
668  }
669  }
670  std::array<std::array<double, 12>, 12> Eta{};
671  for (size_t counterLayer = 0; counterLayer < Eta.size(); counterLayer++) {
672  if (counterLayer > 7 && st == 4)
673  continue;
674  for (size_t counterLayer2 = 0; counterLayer2 < Eta[counterLayer].size(); counterLayer2++) {
675  if (counterLayer2 > 7 && st == 4)
676  continue;
677  if ((counterLayer2 < 4 && counterLayer < 4) || (counterLayer2 > 3 && counterLayer > 3)) {
678  if (counterLayer == counterLayer2) {
679  Eta[counterLayer][counterLayer2] = 3.0 / (4.0);
680  } else {
681  Eta[counterLayer][counterLayer2] = -1.0 / (4.0);
682  }
683  } else {
684  Eta[counterLayer][counterLayer2] = 0.0;
685  }
686  }
687  }
688 
689  for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
690  for (size_t counterLayer = 0; counterLayer < 12; counterLayer++) {
691  if (counterLayer > 7 && st == 4)
692  continue;
693  for (size_t counterLayer2 = 0; counterLayer2 < 12; counterLayer2++) {
694  if (counterLayer2 > 7 && st == 4)
695  continue;
696  float mean = 0;
697  if (counterDeg != 0) {
698  if (counterLayer < 4) {
699  mean = (ResQC(counterDeg, 0) + ResQC(counterDeg, 1) + ResQC(counterDeg, 2) + ResQC(counterDeg, 3)) / 4.0;
700  } else if (counterLayer > 3 && counterLayer < 8) {
701  mean = (ResQC(counterDeg, 4) + ResQC(counterDeg, 5) + ResQC(counterDeg, 6) + ResQC(counterDeg, 7)) / 4.0;
702  } else {
703  mean =
704  (ResQC(counterDeg, 8) + ResQC(counterDeg, 9) + ResQC(counterDeg, 10) + ResQC(counterDeg, 11)) / 4.0;
705  }
706  }
707  matrix(5 * counterLayer2 + counterDeg, 0) +=
708  Eta[counterLayer][counterLayer2] * (ResQC(counterDeg, counterLayer) - mean) /
709  (sigmaQC(counterDeg, counterLayer) * sigmaQC(counterDeg, counterLayer));
710  }
711  }
712  }
713  }
714  return matrix;
715 }
716 
717 // TMatrixD DTMuonMillepede::getCsurveyMatrix(int wh, int st, int se) {
718 
719 // int size = 60;
720 // if(st == 4) size = 40;
721 
722 // TMatrixD matrix(size+6, size+6);
723 // //Careful with the sign
724 // float error[] = {0.05, 0.05, 0.05, 0.005, 0.005, 0.005};
725 // for(int counterLayer = 0; counterLayer < size/5; counterLayer++) {
726 // float k = 1;
727 // if(counterLayer < 4) k = -1;
728 // for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
729 // for(int counterLayer2 = 0; counterLayer2 < size/5; counterLayer2++) {
730 // float k2 = 1;
731 // if(counterLayer2 < 4) k2 = -1;
732 // matrix(5*counterLayer+counterDeg, 5*counterLayer2+counterDeg) = k2*k/(16.0*error[counterDeg]*error[counterDeg]);
733 // }
734 // }
735 // }
736 
737 // return matrix;
738 
739 // }
740 
741 TMatrixD DTMuonMillepede::getCsurveyMatrix(int wh, int st, int se) {
742  int size = 60;
743  if (st == 4)
744  size = 40;
745 
746  TMatrixD matrix(size + 6, size + 6);
747  //Careful with the sign
748  float error[2][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
749  {0.05, 0.05, 0.05, 0.005, 0.005, 0.005}};
750  for (int counterLayer = 0; counterLayer < size / 5; counterLayer++) {
751  for (int counterLayer2 = 0; counterLayer2 < size / 5; counterLayer2++) {
752  int sl1 = counterLayer / 4;
753  int sl2 = counterLayer2 / 4;
754  if (sl1 == sl2) {
755  for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
756  int counterDegAux = counterDeg + 1;
757  if (counterLayer < 8 && counterDeg == 1)
758  counterDegAux = 0;
759  int sl = (sl1 + 1) / 2;
760  matrix(5 * counterLayer + counterDeg, 5 * counterLayer2 + counterDeg) =
761  1.0 / (16.0 * error[sl][counterDegAux] * error[sl][counterDegAux]);
762  }
763  }
764  }
765  }
766 
767  return matrix;
768 }
769 
770 // TMatrixD DTMuonMillepede::getbsurveyMatrix(int wh, int st, int se) {
771 
772 // TMatrixD survey = myPG->giveSurvey(wh, st, se);
773 // int size = 60;
774 // if(st == 4) size = 40;
775 
776 // TMatrixD matrix(size+6, 1);
777 // //Careful with the sign
778 // float error[] = {0.05, 0.05, 0.05, 0.005, 0.005, 0.005};
779 // for(int counterLayer = 0; counterLayer < size/5; counterLayer++) {
780 // float k = 1;
781 // if(counterLayer < 4) k = -1;
782 // for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
783 // int counterDegAux = counterDeg+1;
784 // if(counterLayer < 8 && counterDeg == 1) counterDegAux = 0;
785 // matrix(5*counterLayer+counterDeg, 0) = k*survey(counterDegAux, 0)/(16.0*error[counterDeg]*error[counterDeg]);
786 // }
787 // }
788 
789 // return matrix;
790 
791 // }
792 
793 TMatrixD DTMuonMillepede::getbsurveyMatrix(int wh, int st, int se) {
794  TMatrixD survey = myPG->giveSurvey(wh, st, se);
795  int size = 60;
796  if (st == 4)
797  size = 40;
798 
799  TMatrixD matrix(size + 6, 1);
800  //Careful with the sign
801  float error[2][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
802  {0.05, 0.05, 0.05, 0.005, 0.005, 0.005}};
803  for (int counterLayer = 0; counterLayer < size / 5; counterLayer++) {
804  for (int counterDeg = 0; counterDeg < 5; counterDeg++) {
805  int counterDegAux = counterDeg + 1;
806  if (counterLayer < 8 && counterDeg == 1)
807  counterDegAux = 0;
808  int superlayerAux = counterLayer / 4;
809  int sl = (superlayerAux + 1) / 2;
810  matrix(5 * counterLayer + counterDeg, 0) =
811  survey(counterDegAux, superlayerAux) / (4.0 * error[sl][counterDegAux] * error[sl][counterDegAux]);
812  }
813  }
814 
815  return matrix;
816 }
817 
818 TMatrixD DTMuonMillepede::prepareForLagrange(const TMatrixD &m) {
819  TMatrixD updatedMatrix = m;
820  updatedMatrix.ResizeTo(m.GetNrows() + 6, m.GetNcols() + 6);
821  return updatedMatrix;
822 }
823 
824 void DTMuonMillepede::setBranchTree() {
825  ttreeOutput = new TTree("DTLocalMillepede", "DTLocalMillepede");
826 
827  ttreeOutput->Branch("wh", &whC, "wh/I");
828  ttreeOutput->Branch("st", &stC, "st/I");
829  ttreeOutput->Branch("sr", &srC, "sr/I");
830  ttreeOutput->Branch("cov", cov, "cov[60][60]/F");
831  ttreeOutput->Branch("sl", slC, "sl[12]/I");
832  ttreeOutput->Branch("la", laC, "la[12]/I");
833  ttreeOutput->Branch("dx", dx, "dx[12]/F");
834  ttreeOutput->Branch("dy", dy, "dy[12]/F");
835  ttreeOutput->Branch("dz", dz, "dz[12]/F");
836  ttreeOutput->Branch("phix", phix, "phix[12]/F");
837  ttreeOutput->Branch("phiy", phiy, "phiy[12]/F");
838  ttreeOutput->Branch("phiz", phiz, "phiz[12]/F");
839 }
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
leef::Error
edm::ErrorSummaryEntry Error
Definition: LogErrorEventFilter.cc:29
DTMuonMillepede::calculationMillepede
void calculationMillepede(int)
Definition: DTMuonMillepede.cc:45
dydz
float dydz
Definition: beamSpotDipStandalone.cc:59
DTMuonMillepede::getMatrixFromFile
TMatrixD getMatrixFromFile(const TString &Code, int, int, int, int)
Definition: DTMuonMillepede.cc:464
mps_fire.i
i
Definition: mps_fire.py:429
B
Definition: APVGainStruct.h:7
dxdz
float dxdz
Definition: beamSpotDipStandalone.cc:58
castor_dqm_sourceclient_file_cfg.path
path
Definition: castor_dqm_sourceclient_file_cfg.py:37
DTMuonMillepede::getbsurveyMatrix
TMatrixD getbsurveyMatrix(int, int, int)
Definition: DTMuonMillepede.cc:793
DTMuonMillepede::getCsurveyMatrix
TMatrixD getCsurveyMatrix(int, int, int)
Definition: DTMuonMillepede.cc:741
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:48
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
DTMuonMillepede::getbcsMatrix
TMatrixD getbcsMatrix(int, int, int)
Definition: DTMuonMillepede.cc:398
CustomPhysics_cfi.gamma
gamma
Definition: CustomPhysics_cfi.py:17
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
ReadPGInfo::giveQCCal
TMatrixD giveQCCal(int, int, int)
Definition: ReadPGInfo.cc:348
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:84
Exhume::I
const std::complex< double > I
Definition: I.h:8
DTMuonMillepede::getCcsMatrix
TMatrixD getCcsMatrix(int, int, int)
Definition: DTMuonMillepede.cc:343
correctionTermsCaloMet_cff.C
C
Definition: correctionTermsCaloMet_cff.py:34
DTMuonMillepede::getCqcMatrix
TMatrixD getCqcMatrix(int, int, int)
Definition: DTMuonMillepede.cc:515
genfragment_ptgun_cfg.MaxPt
MaxPt
Definition: genfragment_ptgun_cfg.py:12
DTMuonMillepede::getbqcMatrix
TMatrixD getbqcMatrix(int, int, int)
Definition: DTMuonMillepede.cc:610
b
double b
Definition: hdecay.h:120
DTMuonMillepede::prepareForLagrange
TMatrixD prepareForLagrange(const TMatrixD &)
Definition: DTMuonMillepede.cc:818
DTMuonMillepede::DTMuonMillepede
DTMuonMillepede(std::string, int, float, float, int, int, int, int)
Definition: DTMuonMillepede.cc:6
ReadPGInfo::giveSurvey
TMatrixD giveSurvey(int, int, int)
Definition: ReadPGInfo.cc:387
counter
static std::atomic< unsigned int > counter
Definition: SharedResourceNames.cc:17
nano_mu_digi_cff.sector
sector
Definition: nano_mu_digi_cff.py:20
gather_cfg.cout
cout
Definition: gather_cfg.py:144
A
Definition: APVGainStruct.h:7
DTMuonMillepede::cov
float cov[60][60]
Definition: DTMuonMillepede.h:62
electronAnalyzer_cfi.MinPt
MinPt
Definition: electronAnalyzer_cfi.py:36
DTMuonMillepede::getLagMatrix
TMatrixD getLagMatrix(int, int, int)
Definition: DTMuonMillepede.cc:490
Generated for CMSSW Reference Manual by doxygen 1.8.14