CMS 3D CMS Logo
List of all members | Public Member Functions | Private Attributes
DTMuonMillepede Class Reference

#include <DTMuonMillepede.h>

Inheritance diagram for DTMuonMillepede:
DTMuonLocalAlignment

Public Member Functions

void calculationMillepede (int)
 
 DTMuonMillepede (std::string, int, float, float, int, int, int, int)
 
TMatrixD getbcsMatrix (int, int, int)
 
TMatrixD getbqcMatrix (int, int, int)
 
TMatrixD getbsurveyMatrix (int, int, int)
 
TMatrixD getCcsMatrix (int, int, int)
 
TMatrixD getCqcMatrix (int, int, int)
 
TMatrixD getCsurveyMatrix (int, int, int)
 
TMatrixD getLagMatrix (int, int, int)
 
TMatrixD getMatrixFromFile (const TString &Code, int, int, int, int)
 
TMatrixD prepareForLagrange (const TMatrixD &)
 
void setBranchTree ()
 
 ~DTMuonMillepede ()
 
- Public Member Functions inherited from DTMuonLocalAlignment
 DTMuonLocalAlignment ()
 
void initNTuples (int)
 
void setBranchAddressTree ()
 
 ~DTMuonLocalAlignment ()
 

Private Attributes

float cov [60][60]
 
float dx [12]
 
float dy [12]
 
float dz [12]
 
TFile * f
 
int laC [12]
 
ReadPGInfo * myPG
 
int nPhiHits
 
int nThetaHits
 
float phix [12]
 
float phiy [12]
 
float phiz [12]
 
float ptMax
 
float ptMin
 
int slC [12]
 
int srC
 
int stC
 
TTree * ttreeOutput
 
int whC
 

Additional Inherited Members

- Public Attributes inherited from DTMuonLocalAlignment
float charge
 
float dxdzSl [5]
 
float dxdzSlSL1 [5]
 
float dxdzSlSL3 [5]
 
float dydzSl [5]
 
float edxdzSl [5]
 
float edxdzSlSL1 [5]
 
float edxdzSlSL3 [5]
 
float edydzSl [5]
 
float eta
 
float ex [5][14]
 
float excp [5][14]
 
float exdxdzSl [5]
 
float exdxdzSlSL1 [5]
 
float exdxdzSlSL3 [5]
 
float exSl [5]
 
float exSlSL1 [5]
 
float exSlSL3 [5]
 
float eycp [5][14]
 
float eydydzSl [5]
 
float eySl [5]
 
TFile * f
 
int la [5][14]
 
int nhits [5]
 
int nphihits [5]
 
int nseg
 
int nthetahits [5]
 
std::string ntuplePath
 
int numberOfRootFiles
 
float p
 
float phi
 
float pt
 
int sl [5][14]
 
int sr [5]
 
int st [5]
 
TChain * tali
 
int wh [5]
 
float xc [5][14]
 
float xcp [5][14]
 
float xSl [5]
 
float xSL1SL3 [5]
 
float xSL3SL1 [5]
 
float xSlSL1 [5]
 
float xSlSL3 [5]
 
float yc [5][14]
 
float ycp [5][14]
 
float ySl [5]
 
float zc [5][14]
 

Detailed Description

Date
2010/02/25 11:33:32
Revision
1.2
Author
Luca Scodellaro Luca..nosp@m.Scod.nosp@m.ellar.nosp@m.o@ce.nosp@m.rn.ch

Definition at line 19 of file DTMuonMillepede.h.

Constructor & Destructor Documentation

DTMuonMillepede::DTMuonMillepede ( std::string  path,
int  n_files,
float  MaxPt,
float  MinPt,
int  nPhihits,
int  nThetahits,
int  workingmode,
int  nMtxSection 
)

Definition at line 5 of file DTMuonMillepede.cc.

References calculationMillepede(), f, DTMuonLocalAlignment::initNTuples(), genfragment_ptgun_cfg::MaxPt, electronAnalyzer_cfi::MinPt, myPG, nPhiHits, nThetaHits, DTMuonLocalAlignment::ntuplePath, DTMuonLocalAlignment::numberOfRootFiles, castor_dqm_sourceclient_file_cfg::path, ptMax, ptMin, and setBranchTree().

12  {
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 }
DTMuonMillepede::calculationMillepede
void calculationMillepede(int)
Definition: DTMuonMillepede.cc:44
castor_dqm_sourceclient_file_cfg.path
path
Definition: castor_dqm_sourceclient_file_cfg.py:37
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:48
genfragment_ptgun_cfg.MaxPt
MaxPt
Definition: genfragment_ptgun_cfg.py:12
electronAnalyzer_cfi.MinPt
MinPt
Definition: electronAnalyzer_cfi.py:36
DTMuonMillepede::~DTMuonMillepede ( )

Definition at line 42 of file DTMuonMillepede.cc.

42 {}

Member Function Documentation

void DTMuonMillepede::calculationMillepede ( int  workingmode)

Definition at line 44 of file DTMuonMillepede.cc.

References MaterialEffects_cfi::A, TtFullHadDaughter::B, b, gen::C, counter, gather_cfg::cout, cov, dx, BeamSpotFakeParameters_cfi::dxdz, DTMuonLocalAlignment::dxdzSl, dy, BeamSpotFakeParameters_cfi::dydz, DTMuonLocalAlignment::dydzSl, dz, relativeConstraints::error, DTMuonLocalAlignment::ex, CustomPhysics_cfi::gamma, getbcsMatrix(), getCcsMatrix(), getMatrixFromFile(), Exhume::I, mps_fire::i, DTMuonLocalAlignment::la, laC, DTMuonLocalAlignment::nhits, DTMuonLocalAlignment::nphihits, nPhiHits, DTMuonLocalAlignment::nseg, DTMuonLocalAlignment::nthetahits, nThetaHits, phix, phiy, phiz, DTMuonLocalAlignment::pt, ptMax, ptMin, dttmaxenums::R, DTMuonLocalAlignment::sl, slC, DTMuonLocalAlignment::sr, srC, DTMuonLocalAlignment::st, relativeConstraints::station, stC, DTMuonLocalAlignment::tali, ttreeOutput, DTMuonLocalAlignment::wh, whC, makeMuonMisalignmentScenario::wheel, x, DTMuonLocalAlignment::xc, DTMuonLocalAlignment::xcp, y, DTMuonLocalAlignment::yc, DTMuonLocalAlignment::ycp, and DTMuonLocalAlignment::zc.

Referenced by DTMuonMillepede().

44  {
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 }
DTMuonMillepede::getMatrixFromFile
TMatrixD getMatrixFromFile(const TString &Code, int, int, int, int)
Definition: DTMuonMillepede.cc:463
mps_fire.i
i
Definition: mps_fire.py:341
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
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
MaterialEffects_cfi.A
A
Definition: MaterialEffects_cfi.py:11
b
double b
Definition: hdecay.h:118
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
TMatrixD DTMuonMillepede::getbcsMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 397 of file DTMuonMillepede.cc.

References ReadPGInfo::giveQCCal(), ReadPGInfo::giveSurvey(), DTMuonLocalAlignment::la, makeMuonMisalignmentScenario::matrix, myPG, and findQualityFiles::size.

Referenced by calculationMillepede().

397  {
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 }
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
leef::Error
edm::ErrorSummaryEntry Error
Definition: LogErrorEventFilter.cc:29
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:48
ReadPGInfo::giveQCCal
TMatrixD giveQCCal(int, int, int)
Definition: ReadPGInfo.cc:348
ReadPGInfo::giveSurvey
TMatrixD giveSurvey(int, int, int)
Definition: ReadPGInfo.cc:387
TMatrixD DTMuonMillepede::getbqcMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 611 of file DTMuonMillepede.cc.

References ReadPGInfo::giveQCCal(), makeMuonMisalignmentScenario::matrix, SiStripPI::mean, myPG, and mathSSE::sqrt().

611  {
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 }
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:48
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
ReadPGInfo::giveQCCal
TMatrixD giveQCCal(int, int, int)
Definition: ReadPGInfo.cc:348
TMatrixD DTMuonMillepede::getbsurveyMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 795 of file DTMuonMillepede.cc.

References relativeConstraints::error, ReadPGInfo::giveSurvey(), makeMuonMisalignmentScenario::matrix, myPG, findQualityFiles::size, and DTMuonLocalAlignment::sl.

795  {
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 }
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:48
ReadPGInfo::giveSurvey
TMatrixD giveSurvey(int, int, int)
Definition: ReadPGInfo.cc:387
TMatrixD DTMuonMillepede::getCcsMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 342 of file DTMuonMillepede.cc.

References ReadPGInfo::giveQCCal(), DTMuonLocalAlignment::la, makeMuonMisalignmentScenario::matrix, myPG, and findQualityFiles::size.

Referenced by calculationMillepede().

342  {
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 }
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
leef::Error
edm::ErrorSummaryEntry Error
Definition: LogErrorEventFilter.cc:29
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:48
ReadPGInfo::giveQCCal
TMatrixD giveQCCal(int, int, int)
Definition: ReadPGInfo.cc:348
TMatrixD DTMuonMillepede::getCqcMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 514 of file DTMuonMillepede.cc.

References ReadPGInfo::giveQCCal(), makeMuonMisalignmentScenario::matrix, myPG, and mathSSE::sqrt().

514  {
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 }
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:48
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
ReadPGInfo::giveQCCal
TMatrixD giveQCCal(int, int, int)
Definition: ReadPGInfo.cc:348
TMatrixD DTMuonMillepede::getCsurveyMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 743 of file DTMuonMillepede.cc.

References relativeConstraints::error, makeMuonMisalignmentScenario::matrix, findQualityFiles::size, and DTMuonLocalAlignment::sl.

743  {
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 }
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
TMatrixD DTMuonMillepede::getLagMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 489 of file DTMuonMillepede.cc.

References makeMuonMisalignmentScenario::matrix.

489  {
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 }
TMatrixD DTMuonMillepede::getMatrixFromFile ( const TString &  Code,
int  wh,
int  st,
int  se,
int  mf 
)

Definition at line 463 of file DTMuonMillepede.cc.

References DTMuonLocalAlignment::st, and DTMuonLocalAlignment::wh.

Referenced by calculationMillepede().

463  {
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 }
TMatrixD DTMuonMillepede::prepareForLagrange ( const TMatrixD &  m)

Definition at line 820 of file DTMuonMillepede.cc.

References visualization-live-secondInstance_cfg::m.

820  {
821  TMatrixD updatedMatrix = m;
822  updatedMatrix.ResizeTo(m.GetNrows() + 6, m.GetNcols() + 6);
823  return updatedMatrix;
824 }
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:60
void DTMuonMillepede::setBranchTree ( )

Definition at line 826 of file DTMuonMillepede.cc.

References cov, dx, dy, dz, laC, phix, phiy, phiz, slC, srC, stC, ttreeOutput, and whC.

Referenced by DTMuonMillepede().

826  {
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 }
DTMuonMillepede::cov
float cov[60][60]
Definition: DTMuonMillepede.h:62

Member Data Documentation

float DTMuonMillepede::cov[60][60]
private

Definition at line 62 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::dx[12]
private

Definition at line 61 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::dy[12]
private

Definition at line 61 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::dz[12]
private

Definition at line 61 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), Electron.Electron::cutBasedId(), ntupleDataFormat.Track::dzPull(), and setBranchTree().

TFile* DTMuonMillepede::f
private

Definition at line 50 of file DTMuonMillepede.h.

Referenced by Vispa.Views.RootCanvasView.RootCanvasView::createGraph(), Vispa.Views.RootCanvasView.RootCanvasView::createLegoPlot(), DTMuonMillepede(), and ztail.Decoder::initial_synchronize().

int DTMuonMillepede::laC[12]
private

Definition at line 60 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

ReadPGInfo* DTMuonMillepede::myPG
private

Definition at line 48 of file DTMuonMillepede.h.

Referenced by DTMuonMillepede(), getbcsMatrix(), getbqcMatrix(), getbsurveyMatrix(), getCcsMatrix(), and getCqcMatrix().

int DTMuonMillepede::nPhiHits
private

Definition at line 55 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and DTMuonMillepede().

int DTMuonMillepede::nThetaHits
private

Definition at line 55 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and DTMuonMillepede().

float DTMuonMillepede::phix[12]
private

Definition at line 61 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::phiy[12]
private

Definition at line 61 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::phiz[12]
private

Definition at line 61 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::ptMax
private

Definition at line 53 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and DTMuonMillepede().

float DTMuonMillepede::ptMin
private

Definition at line 53 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and DTMuonMillepede().

int DTMuonMillepede::slC[12]
private

Definition at line 60 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

int DTMuonMillepede::srC
private

Definition at line 59 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

int DTMuonMillepede::stC
private

Definition at line 59 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

TTree* DTMuonMillepede::ttreeOutput
private

Definition at line 51 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

int DTMuonMillepede::whC
private

Definition at line 59 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

Generated for CMSSW Reference Manual by doxygen 1.8.11