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