CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 20 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 6 of file DTMuonMillepede.cc.

References calculationMillepede(), f, DTMuonLocalAlignment::initNTuples(), HLT_25ns14e33_v1_cff::MaxPt, step1_ZMM_7Tev::MinPt, myPG, nPhiHits, nThetaHits, DTMuonLocalAlignment::ntuplePath, DTMuonLocalAlignment::numberOfRootFiles, cmsHarvester::path, ptMax, ptMin, and setBranchTree().

9  {
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 }
DTMuonMillepede::calculationMillepede
void calculationMillepede(int)
Definition: DTMuonMillepede.cc:45
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:52
cmsHarvester.path
tuple path
else: Piece not in the list, fine.
Definition: cmsHarvester.py:1591
step1_ZMM_7Tev.MinPt
tuple MinPt
Definition: step1_ZMM_7Tev.py:113
HLT_25ns14e33_v1_cff.MaxPt
tuple MaxPt
Definition: HLT_25ns14e33_v1_cff.py:15976
DTMuonMillepede::~DTMuonMillepede ( )

Definition at line 41 of file DTMuonMillepede.cc.

41 {}

Member Function Documentation

void DTMuonMillepede::calculationMillepede ( int  workingmode)

Definition at line 45 of file DTMuonMillepede.cc.

References HLT_25ns14e33_v1_cff::A, b, funct::C, counter, gather_cfg::cout, cov, dx, DTMuonLocalAlignment::dxdzSl, dy, DTMuonLocalAlignment::dydzSl, dz, relativeConstraints::error, DTMuonLocalAlignment::ex, getbcsMatrix(), getCcsMatrix(), getMatrixFromFile(), Exhume::I, 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, x, DTMuonLocalAlignment::xc, DTMuonLocalAlignment::xcp, detailsBasic3DVector::y, DTMuonLocalAlignment::yc, DTMuonLocalAlignment::ycp, and DTMuonLocalAlignment::zc.

Referenced by DTMuonMillepede().

45  {
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 }
HLT_25ns14e33_v1_cff.A
tuple A
Definition: HLT_25ns14e33_v1_cff.py:396
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::getbcsMatrix
TMatrixD getbcsMatrix(int, int, int)
Definition: DTMuonMillepede.cc:392
double_binary
Definition: DDStreamer.cc:227
Exhume::I
const std::complex< double > I
Definition: I.h:8
DTMuonMillepede::getCcsMatrix
TMatrixD getCcsMatrix(int, int, int)
Definition: DTMuonMillepede.cc:344
b
double b
Definition: hdecay.h:120
counter
static std::atomic< unsigned int > counter
Definition: SharedResourceNames.cc:15
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
x
Definition: DDAxes.h:10
DTMuonMillepede::cov
float cov[60][60]
Definition: DTMuonMillepede.h:66
TMatrixD DTMuonMillepede::getbcsMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 392 of file DTMuonMillepede.cc.

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

Referenced by calculationMillepede().

392  {
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 }
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:52
ReadPGInfo::giveQCCal
TMatrixD giveQCCal(int, int, int)
Definition: ReadPGInfo.cc:803
ReadPGInfo::giveSurvey
TMatrixD giveSurvey(int, int, int)
Definition: ReadPGInfo.cc:840
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
TMatrixD DTMuonMillepede::getbqcMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 576 of file DTMuonMillepede.cc.

References reco::tau::disc::Eta(), ReadPGInfo::giveQCCal(), makeMuonMisalignmentScenario::matrix, timingPdfMaker::mean, myPG, and mathSSE::sqrt().

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

Definition at line 749 of file DTMuonMillepede.cc.

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

749  {
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 }
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:52
ReadPGInfo::giveSurvey
TMatrixD giveSurvey(int, int, int)
Definition: ReadPGInfo.cc:840
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
TMatrixD DTMuonMillepede::getCcsMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 344 of file DTMuonMillepede.cc.

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

Referenced by calculationMillepede().

344  {
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 }
DTMuonMillepede::myPG
ReadPGInfo * myPG
Definition: DTMuonMillepede.h:52
ReadPGInfo::giveQCCal
TMatrixD giveQCCal(int, int, int)
Definition: ReadPGInfo.cc:803
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
TMatrixD DTMuonMillepede::getCqcMatrix ( int  wh,
int  st,
int  se 
)

Definition at line 492 of file DTMuonMillepede.cc.

References reco::tau::disc::Eta(), ReadPGInfo::giveQCCal(), makeMuonMisalignmentScenario::matrix, myPG, and mathSSE::sqrt().

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

Definition at line 694 of file DTMuonMillepede.cc.

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

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

Definition at line 466 of file DTMuonMillepede.cc.

References makeMuonMisalignmentScenario::matrix.

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

Definition at line 446 of file DTMuonMillepede.cc.

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

Referenced by calculationMillepede().

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

Definition at line 774 of file DTMuonMillepede.cc.

References contentValuesFiles::m.

774  {
775 
776  TMatrixD updatedMatrix = m;
777  updatedMatrix.ResizeTo(m.GetNrows()+6, m.GetNcols()+6);
778  return updatedMatrix;
779 
780 }
void DTMuonMillepede::setBranchTree ( )

Definition at line 783 of file DTMuonMillepede.cc.

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

Referenced by DTMuonMillepede().

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

Member Data Documentation

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

Definition at line 66 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::dx[12]
private

Definition at line 65 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::dy[12]
private

Definition at line 65 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::dz[12]
private

Definition at line 65 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

TFile* DTMuonMillepede::f
private

Definition at line 54 of file DTMuonMillepede.h.

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

int DTMuonMillepede::laC[12]
private

Definition at line 64 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

ReadPGInfo* DTMuonMillepede::myPG
private

Definition at line 52 of file DTMuonMillepede.h.

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

int DTMuonMillepede::nPhiHits
private

Definition at line 59 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and DTMuonMillepede().

int DTMuonMillepede::nThetaHits
private

Definition at line 59 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and DTMuonMillepede().

float DTMuonMillepede::phix[12]
private

Definition at line 65 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::phiy[12]
private

Definition at line 65 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::phiz[12]
private

Definition at line 65 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

float DTMuonMillepede::ptMax
private

Definition at line 57 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and DTMuonMillepede().

float DTMuonMillepede::ptMin
private

Definition at line 57 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and DTMuonMillepede().

int DTMuonMillepede::slC[12]
private

Definition at line 64 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

int DTMuonMillepede::srC
private

Definition at line 63 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

int DTMuonMillepede::stC
private

Definition at line 63 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

TTree* DTMuonMillepede::ttreeOutput
private

Definition at line 55 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

int DTMuonMillepede::whC
private

Definition at line 63 of file DTMuonMillepede.h.

Referenced by calculationMillepede(), and setBranchTree().

Generated for CMSSW Reference Manual by doxygen 1.8.5