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