|
|
|
#include <DTMuonMillepede.h>
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 (TString Code, int, int, int, int) |
| TMatrixD | prepareForLagrange (const TMatrixD &) |
| void | setBranchTree () |
| ~DTMuonMillepede () | |
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 |
| 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(), step1_ZMM_7Tev::MinPt, myPG, nPhiHits, nThetaHits, DTMuonLocalAlignment::ntuplePath, DTMuonLocalAlignment::numberOfRootFiles, getHLTPrescaleColumns::path, ptMax, ptMin, and setBranchTree().
{
ntuplePath = path;
numberOfRootFiles = n_files;
ptMax = MaxPt;
ptMin = MinPt;
nPhiHits = nPhihits;
nThetaHits = nThetahits;
TDirectory * dirSave = gDirectory;
//Interface to Survey information
myPG = new ReadPGInfo("./InternalData2009.root");
f = new TFile("./LocalMillepedeResults.root", "RECREATE");
f->cd();
setBranchTree();
initNTuples(nMtxSection);
calculationMillepede(workingmode);
f->Write();
f->Close();
dirSave->cd();
}
| DTMuonMillepede::~DTMuonMillepede | ( | ) |
Definition at line 41 of file DTMuonMillepede.cc.
{}
| void DTMuonMillepede::calculationMillepede | ( | int | workingmode | ) |
Definition at line 45 of file DTMuonMillepede.cc.
References funct::A, b, funct::C, gather_cfg::cout, cov, dx, DTMuonLocalAlignment::dxdzSl, dy, DTMuonLocalAlignment::dydzSl, dz, error, DTMuonLocalAlignment::ex, getbcsMatrix(), getCcsMatrix(), getMatrixFromFile(), Exhume::I, i, DTMuonLocalAlignment::la, laC, DTMuonLocalAlignment::nhits, DTMuonLocalAlignment::nphihits, nPhiHits, DTMuonLocalAlignment::nseg, nThetaHits, DTMuonLocalAlignment::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().
{
//Init Matrizes
TMatrixD ****C = new TMatrixD ***[5];
TMatrixD ****b = new TMatrixD ***[5];
for(int whI = -2; whI < 3; ++whI) {
C[whI+2] = new TMatrixD **[4];
b[whI+2] = new TMatrixD **[4];
for(int stI = 1; stI < 5; ++stI) {
C[whI+2][stI-1] = new TMatrixD * [14];
b[whI+2][stI-1] = new TMatrixD * [14];
for(int seI = 1; seI < 15; ++seI) {
if(seI > 12 && stI != 4) continue;
if(stI == 4) {
C[whI+2][stI-1][seI-1] = new TMatrixD(24,24);
b[whI+2][stI-1][seI-1] = new TMatrixD(24,1);
} else {
C[whI+2][stI-1][seI-1] = new TMatrixD(60,60);
b[whI+2][stI-1][seI-1] = new TMatrixD(60,1);
}
}
}
}
//Run over the TTree
if (workingmode<=3) {
Int_t nentries = (Int_t)tali->GetEntries();
for (Int_t i=0;i<nentries;i++) {
tali->GetEntry(i);
if (i%100000==0) std::cout << "Analyzing track number " << i << std::endl;
//Basic cuts
if(pt > ptMax || pt < ptMin) continue;
for(int counter = 0; counter < nseg; ++counter) {
if (nphihits[counter]<nPhiHits || (nthetahits[counter]<nThetaHits && st[counter]<4)) continue;
TMatrixD A(12,60);
TMatrixD R(12,1);
TMatrixD E(12,12);
TMatrixD B(12, 4);
TMatrixD I(12, 12);
if(st[counter] == 4) {
A.ResizeTo(8,24);
R.ResizeTo(8,1);
E.ResizeTo(8,8);
B.ResizeTo(8,2);
I.ResizeTo(8,8);
}
for(int counterHi = 0; counterHi < nhits[counter]; counterHi++) {
int row = 0;
if(sl[counter][counterHi] == 3) {
row = 4 + (la[counter][counterHi]-1);
} else if(sl[counter][counterHi] == 2) {
row = 8 + (la[counter][counterHi]-1);
} else {
row = (la[counter][counterHi]-1);
}
float x = xc[counter][counterHi];
float y = yc[counter][counterHi];
float xp = xcp[counter][counterHi];
float yp = ycp[counter][counterHi];
float zp = zc[counter][counterHi];
float error = ex[counter][counterHi];
float dxdz = dxdzSl[counter];
float dydz = dydzSl[counter];
if(st[counter]!=4) {
if(sl[counter][counterHi] == 2) {
A(row, row*5) = -1.0;
A(row, row*5+1) = dydz;
A(row, row*5+2) = y*dydz;
A(row, row*5+3) = -x*dydz;
A(row, row*5+4) = -x;
R(row, 0) = y-yp;
} else {
A(row, row*5+0) = -1.0;
A(row, row*5+1) = dxdz;
A(row, row*5+2) = y*dxdz;
A(row, row*5+3) = -x*dxdz;
A(row, row*5+4) = y;
R(row, 0) = x-xp;
}
} else {
if(sl[counter][counterHi]!=2) {
A(row, row*3+0) = -1.0;
A(row, row*3+1) = dxdz;
A(row, row*3+2) = -x*dxdz;
R(row, 0) = x-xp;
}
}
E(row, row) = 1.0/error;
I(row, row) = 1.0;
if(sl[counter][counterHi] != 2) {
B(row, 0) = 1.0;
B(row, 1) = zp;
} else {
B(row, 2) = 1.0;
B(row, 3) = zp;
}
}
TMatrixD Bt = B; Bt.T();
TMatrixD At = A; At.T();
TMatrixD gamma = (Bt*E*B); gamma.Invert();
*(C[wh[counter]+2][st[counter]-1][sr[counter]-1]) += At*E*(B*gamma*Bt*E*A-A);
*(b[wh[counter]+2][st[counter]-1][sr[counter]-1]) += At*E*(B*gamma*Bt*E*I-I)*R;
}
}
}
if (workingmode==3)
for(int wheel = -2; wheel < 3; ++wheel)
for(int station = 1; station < 5; ++station)
for(int sector = 1; sector < 15; ++sector) {
if(sector > 12 && station != 4) continue;
TMatrixD Ctr = *C[wheel+2][station-1][sector-1];
TMatrixD btr = *b[wheel+2][station-1][sector-1];
TString CtrName = "Ctr_"; CtrName += wheel; CtrName += "_"; CtrName += station;
CtrName += "_"; CtrName += sector;
Ctr.Write(CtrName);
TString btrName = "btr_"; btrName += wheel; btrName += "_"; btrName += station;
btrName += "_"; btrName += sector;
btr.Write(btrName);
}
// My Final Calculation and Constraints
if (workingmode%2==0) {
for(int wheel = -2; wheel < 3; ++wheel)
for(int station = 1; station < 5; ++station)
for(int sector = 1; sector < 15; ++sector) {
if(sector > 12 && station != 4) continue;
if (workingmode==4) {
for (int mf = -1; mf<=-1; mf++) {
TMatrixD Ch = getMatrixFromFile("Ctr_",wheel,station,sector,mf);
TMatrixD bh = getMatrixFromFile("btr_",wheel,station,sector,mf);
*(C[wheel+2][station-1][sector-1]) += Ch;
*(b[wheel+2][station-1][sector-1]) += bh;
}
}
TMatrixD Ctr = *C[wheel+2][station-1][sector-1];
TMatrixD btr = *b[wheel+2][station-1][sector-1];
TMatrixD Ccs = getCcsMatrix(wheel, station, sector);
TMatrixD bcs = getbcsMatrix(wheel, station, sector);
TMatrixD SC = Ctr + Ccs;
TMatrixD Sb = btr + bcs;
SC.Invert();
TMatrixD solution = SC*Sb;
for(int icounter = 0; icounter < SC.GetNrows(); icounter++) {
for(int jcounter = 0; jcounter < SC.GetNrows(); jcounter++) {
cov[icounter][jcounter] = SC(icounter, jcounter);
}
}
int nLayer = 12, nDeg = 5;
if (station==4) { nLayer = 8; nDeg = 3; }
for(int counterLayer = 0; counterLayer < nLayer; ++counterLayer)
for(int counterDeg = 0; counterDeg < nDeg; ++counterDeg) {
if(counterLayer < 4) {
slC[counterLayer] = 1;
laC[counterLayer] = counterLayer+1;
} else if(counterLayer > 3 && counterLayer < 8) {
slC[counterLayer] = 3;
laC[counterLayer] = counterLayer-3;
} else {
slC[counterLayer] = 2;
laC[counterLayer] = counterLayer-7;
}
if(counterLayer < 8) {
dx[counterLayer] = solution(counterLayer*nDeg, 0);
dy[counterLayer] = 0.0;
} else {
dx[counterLayer] = 0.0;
dy[counterLayer] = solution(counterLayer*nDeg, 0);
}
dz[counterLayer] = solution(counterLayer*nDeg+1, 0);
if (station!=4) {
phix[counterLayer] = solution(counterLayer*nDeg+2, 0);
phiy[counterLayer] = solution(counterLayer*nDeg+3, 0);
phiz[counterLayer] = solution(counterLayer*nDeg+4, 0);
} else {
phix[counterLayer] = 0.;
phiy[counterLayer] = solution(counterLayer*nDeg+2, 0);
phiz[counterLayer] = 0.;
}
}
whC = wheel; stC = station; srC = sector;
ttreeOutput->Fill();
}
}
//Final Calculation and Constraints
// for(int wheel = -2; wheel < 3; ++wheel) {
// for(int station = 1; station < 5; ++station) {
// for(int sector = 1; sector < 15; ++sector) {
// if(sector > 12 && station != 4) continue;
// TMatrixD Ctr = prepareForLagrange(*C[wheel+2][station-1][sector-1]);
// TMatrixD btr = prepareForLagrange(*b[wheel+2][station-1][sector-1]);
// TMatrixD Cqc = prepareForLagrange(getCqcMatrix(wheel, station, sector));
// TMatrixD bqc = prepareForLagrange(getbqcMatrix(wheel, station, sector));
// TMatrixD Csurvey = prepareForLagrange(getCsurveyMatrix(wheel, station, sector));
// TMatrixD bsurvey = prepareForLagrange(getbsurveyMatrix(wheel, station, sector));
// TMatrixD Clag = getLagMatrix(wheel, station, sector);
// TMatrixD SC = Ctr + Cqc + Csurvey + Clag;
// TMatrixD Sb = btr + bqc + bsurvey;
//
// SC.Invert();
// TMatrixD solution = SC*Sb;
// for(int icounter = 0; icounter < SC.GetNrows(); icounter++) {
// for(int jcounter = 0; jcounter < SC.GetNrows(); jcounter++) {
// cov[icounter][jcounter] = SC(icounter, jcounter);
// }
// }
// whC = wheel;
// stC = station;
// srC = sector;
// for(int c = 0; c < 60; ++c) {
// for(int s = 0; s < 60; ++s) {
// cov[c][s] = SC(c, s);
// }
// }
// for(int counterLayer = 0; counterLayer < 12; ++counterLayer) {
// for(int counterDeg = 0; counterDeg < 5; ++counterDeg) {
// if(counterLayer > 7 && station == 4) continue;
// if(counterLayer < 4) {
// slC[counterLayer] = 1;
// laC[counterLayer] = counterLayer+1;
// } else if(counterLayer > 3 && counterLayer < 8) {
// slC[counterLayer] = 3;
// laC[counterLayer] = counterLayer-3;
// } else {
// slC[counterLayer] = 2;
// laC[counterLayer] = counterLayer-7;
// }
// if(counterLayer < 8) {
// dx[counterLayer] = solution(counterLayer*5, 0);
// dy[counterLayer] = 0.0;
// } else {
// dx[counterLayer] = 0.0;
// dy[counterLayer] = solution(counterLayer*5, 0);
// }
// dz[counterLayer] = solution(counterLayer*5+1, 0);
// phix[counterLayer] = solution(counterLayer*5+2, 0);
// phiy[counterLayer] = solution(counterLayer*5+3, 0);
// phiz[counterLayer] = solution(counterLayer*5+4, 0);
// }
// }
// ttreeOutput->Fill();
// }
// }
// }
// f->Write();
}
| TMatrixD DTMuonMillepede::getbcsMatrix | ( | int | wh, |
| int | st, | ||
| int | se | ||
| ) |
Definition at line 392 of file DTMuonMillepede.cc.
References ReadPGInfo::giveQCCal(), ReadPGInfo::giveSurvey(), getHLTprescales::index, DTMuonLocalAlignment::la, makeMuonMisalignmentScenario::matrix, myPG, and findQualityFiles::size.
Referenced by calculationMillepede().
{
int size = 60;
if(st == 4) size = 24;
TMatrixD matrix(size, 1);
float Error[3][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
{0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005},
{0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005}};
float TollerancyPosition = 0.01;
float TollerancyRotation = 0.0001;
TMatrixD ConsQC = myPG->giveQCCal(wh, st, se);
bool UseQC = true;
if (ConsQC.GetNrows()<70) UseQC = false;
TMatrixD Survey = myPG->giveSurvey(wh, st, se);
int nLayer = 12, nDeg = 5;
if (st==4) { nLayer = 8; nDeg = 3; }
for (int la = 0; la<nLayer; la++)
for (int dg = 0; dg<nDeg; dg++) {
int index = la*nDeg + dg;
int rdg = dg + 1;
if (la<8 && rdg==1) rdg = 0;
if (st==4 && rdg==3) rdg = 4;
double ThisError2 = Error[la/4][rdg]*Error[la/4][rdg];
if (rdg<2) {
if (UseQC) { ThisError2 += ConsQC(la, 1)*ConsQC(la, 1);
} else { ThisError2 += TollerancyPosition*TollerancyPosition; }
} else if (rdg==2) { ThisError2 += TollerancyPosition*TollerancyPosition;
} else { ThisError2 += TollerancyRotation*TollerancyRotation; }
float Constraint = 0.;
if (la>3)
Constraint += Survey(rdg, la/4-1);
if (UseQC && rdg==0) Constraint += ConsQC(la, 0);;
if (UseQC && rdg==1) Constraint -= ConsQC(la, 0);
matrix(index, 0) = Constraint/ThisError2;
}
return matrix;
}
| 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().
{
TMatrixD surv = myPG->giveQCCal(wh, st, se);
TMatrixD ResQC(5, 12);
TMatrixD sigmaQC(5, 12);
TMatrixD matrix(60, 1);
if(st == 4) matrix.ResizeTo(40, 1);
if(surv.GetNrows() < 7) {
for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
for(int counterLayer = 0; counterLayer < 12; ++counterLayer) {
if(st != 4 && counterLayer > 7) continue;
if(counterDeg == 0) {
sigmaQC(counterDeg, counterLayer) = 0.05;
} else if(counterDeg < 3) {
sigmaQC(counterDeg, counterLayer) = 0.05;
} else {
sigmaQC(counterDeg, counterLayer) = 0.05;
}
}
}
} else {
for(int counterChamber = 0; counterChamber < 12; ++counterChamber) {
ResQC(0, counterChamber) = -surv(counterChamber, 0)/10000.0;
}
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;
float meanvarSL2 = 0;
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;
}
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;
for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
for(int counterLayer = 0; counterLayer < 12; ++counterLayer) {
if(st != 4 && counterLayer > 7) continue;
float meanerror = 0;
if(counterLayer < 4) {
meanerror = meanvarSL1;
} else if(counterLayer > 3 && counterLayer < 8){
meanerror = meanvarSL2;
} else {
meanerror = meanvarSL3;
}
if(counterDeg == 0) {
sigmaQC(counterDeg, counterLayer) = sqrt((surv(counterLayer, 1)/10000.0)*(surv(counterLayer, 1)/10000.0)+meanerror*meanerror);
} else if(counterDeg < 3) {
sigmaQC(counterDeg, counterLayer) = 0.05;
} else {
sigmaQC(counterDeg, counterLayer) = 0.05;
}
}
}
double **Eta = new double *[12];
for(int counterLayer = 0; counterLayer < 12; counterLayer++) {
if(counterLayer > 7 && st == 4) continue;
Eta[counterLayer] = new double [5];
for(int counterLayer2 = 0; counterLayer2 < 12; counterLayer2++) {
if(counterLayer > 7 && st == 4) continue;
if((counterLayer2 < 4 && counterLayer < 4) || (counterLayer2 > 3 && counterLayer > 3)) {
if(counterLayer == counterLayer2) {
Eta[counterLayer][counterLayer2] = 3.0/(4.0);
} else {
Eta[counterLayer][counterLayer2] = -1.0/(4.0);
}
} else {
Eta[counterLayer][counterLayer2] = 0.0;
}
}
}
for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
for(int counterLayer = 0; counterLayer < 12; counterLayer++) {
if(counterLayer > 7 && st == 4) continue;
for(int counterLayer2 = 0; counterLayer2 < 12; counterLayer2++) {
if(counterLayer2 > 7 && st == 4) continue;
float mean = 0;
if(counterDeg != 0) {
if(counterLayer < 4) {
mean = (ResQC(counterDeg, 0)+ResQC(counterDeg, 1)+ResQC(counterDeg, 2)+ResQC(counterDeg, 3))/4.0;
} else if(counterLayer > 3 && counterLayer < 8) {
mean = (ResQC(counterDeg, 4)+ResQC(counterDeg, 5)+ResQC(counterDeg, 6)+ResQC(counterDeg, 7))/4.0;
} else {
mean = (ResQC(counterDeg, 8)+ResQC(counterDeg, 9)+ResQC(counterDeg, 10)+ResQC(counterDeg, 11))/4.0;
}
}
matrix(5*counterLayer2+counterDeg, 0) += Eta[counterLayer][counterLayer2]*(ResQC(counterDeg,counterLayer)-mean)/(sigmaQC(counterDeg,counterLayer)*sigmaQC(counterDeg,counterLayer));
}
}
}
}
return matrix;
}
| TMatrixD DTMuonMillepede::getbsurveyMatrix | ( | int | wh, |
| int | st, | ||
| int | se | ||
| ) |
Definition at line 749 of file DTMuonMillepede.cc.
References error, ReadPGInfo::giveSurvey(), makeMuonMisalignmentScenario::matrix, myPG, findQualityFiles::size, and DTMuonLocalAlignment::sl.
{
TMatrixD survey = myPG->giveSurvey(wh, st, se);
int size = 60;
if(st == 4) size = 40;
TMatrixD matrix(size+6, 1);
//Careful with the sign
float error[2][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
{0.05, 0.05, 0.05, 0.005, 0.005, 0.005}};
for(int counterLayer = 0; counterLayer < size/5; counterLayer++) {
for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
int counterDegAux = counterDeg+1;
if(counterLayer < 8 && counterDeg == 1) counterDegAux = 0;
int superlayerAux = counterLayer/4;
int sl = (superlayerAux + 1)/2;
matrix(5*counterLayer+counterDeg, 0) = survey(counterDegAux, superlayerAux)/(4.0*error[sl][counterDegAux]*error[sl][counterDegAux]);
}
}
return matrix;
}
| TMatrixD DTMuonMillepede::getCcsMatrix | ( | int | wh, |
| int | st, | ||
| int | se | ||
| ) |
Definition at line 344 of file DTMuonMillepede.cc.
References ReadPGInfo::giveQCCal(), getHLTprescales::index, DTMuonLocalAlignment::la, makeMuonMisalignmentScenario::matrix, myPG, and findQualityFiles::size.
Referenced by calculationMillepede().
{
int size = 60;
if(st == 4) size = 24;
TMatrixD matrix(size, size);
float Error[3][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
{0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005},
{0.005, 0.005, 0.005, 0.00005, 0.00005, 0.00005}};
float TollerancyPosition = 0.01;
float TollerancyRotation = 0.0001;
TMatrixD ConsQC = myPG->giveQCCal(wh, st, se);
bool UseQC = true;
if (ConsQC.GetNrows()<70) UseQC = false;
int nLayer = 12, nDeg = 5;
if (st==4) { nLayer = 8; nDeg = 3; }
for (int la = 0; la<nLayer; la++)
for (int dg = 0; dg<nDeg; dg++) {
int index = la*nDeg + dg;
int rdg = dg + 1;
if (la<8 && rdg==1) rdg = 0;
if (st==4 && rdg==3) rdg = 4;
double ThisError2 = Error[la/4][rdg]*Error[la/4][rdg];
if (rdg<2) {
if (UseQC) { ThisError2 += ConsQC(la, 1)*ConsQC(la, 1);
} else { ThisError2 += TollerancyPosition*TollerancyPosition; }
} else if (rdg==2) { ThisError2 += TollerancyPosition*TollerancyPosition;
} else { ThisError2 += TollerancyRotation*TollerancyRotation; }
matrix(index, index) = 1./ThisError2;
}
return matrix;
}
| 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().
{
TMatrixD surv = myPG->giveQCCal(wh, st, se);
TMatrixD sigmaQC(5, 12);
TMatrixD matrix(60, 60);
if(st == 4) matrix.ResizeTo(40, 40);
if(surv.GetNrows() < 7) {
for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
for(int counterLayer = 0; counterLayer < 12; ++counterLayer) {
if(st != 4 && counterLayer > 7) continue;
if(counterDeg == 0) {
sigmaQC(counterDeg, counterLayer) = 0.05;
} else if(counterDeg < 3) {
sigmaQC(counterDeg, counterLayer) = 0.05;
} else {
sigmaQC(counterDeg, counterLayer) = 0.05;
}
}
}
} else {
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;
float meanvarSL2 = 0;
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;
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;
for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
for(int counterLayer = 0; counterLayer < 12; ++counterLayer) {
if(st != 4 && counterLayer > 7) continue;
float meanerror = 0;
if(counterLayer < 4) {
meanerror = meanvarSL1;
} else if(counterLayer > 3 && counterLayer < 8){
meanerror = meanvarSL2;
} else {
meanerror = meanvarSL3;
}
if(counterDeg == 0) {
sigmaQC(counterDeg, counterLayer) = sqrt((surv(counterLayer, 1)/10000.0)*(surv(counterLayer, 1)/10000.0)+meanerror*meanerror);
} else if(counterDeg < 3) {
sigmaQC(counterDeg, counterLayer) = 0.05;
} else {
sigmaQC(counterDeg, counterLayer) = 0.05;
}
}
}
double **Eta = new double *[12];
for(int counterLayer = 0; counterLayer < 12; counterLayer++) {
if(counterLayer > 7 && st == 4) continue;
Eta[counterLayer] = new double [5];
for(int counterLayer2 = 0; counterLayer2 < 12; counterLayer2++) {
if(counterLayer > 7 && st == 4) continue;
if((counterLayer2 < 4 && counterLayer < 4) || (counterLayer2 > 3 && counterLayer > 3)) {
if(counterLayer == counterLayer2) {
Eta[counterLayer][counterLayer2] = 3.0/(4.0);
} else {
Eta[counterLayer][counterLayer2] = -1.0/(4.0);
}
} else {
Eta[counterLayer][counterLayer2] = 0.0;
}
}
}
for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
for(int counterLayer = 0; counterLayer < 12; counterLayer++) {
if(counterLayer > 7 && st == 4) continue;
for(int counterLayer2 = 0; counterLayer2 < 12; counterLayer2++) {
if(counterLayer2 > 7 && st == 4) continue;
for(int counterLayer3 = 0; counterLayer3 < 12; counterLayer3++) {
if(counterLayer3 > 7 && st == 4) continue;
matrix(5*counterLayer2+counterDeg, 5*counterLayer3+counterDeg) +=
Eta[counterLayer][counterLayer2]*Eta[counterLayer][counterLayer3]/(sigmaQC(counterDeg,counterLayer)*sigmaQC(counterDeg,counterLayer));
}
}
}
}
}
return matrix;
}
| TMatrixD DTMuonMillepede::getCsurveyMatrix | ( | int | wh, |
| int | st, | ||
| int | se | ||
| ) |
Definition at line 694 of file DTMuonMillepede.cc.
References error, makeMuonMisalignmentScenario::matrix, findQualityFiles::size, and DTMuonLocalAlignment::sl.
{
int size = 60;
if(st == 4) size = 40;
TMatrixD matrix(size+6, size+6);
//Careful with the sign
float error[2][6] = {{0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001},
{0.05, 0.05, 0.05, 0.005, 0.005, 0.005}};
for(int counterLayer = 0; counterLayer < size/5; counterLayer++) {
for(int counterLayer2 = 0; counterLayer2 < size/5; counterLayer2++) {
int sl1 = counterLayer/4;
int sl2 = counterLayer2/4;
if (sl1==sl2) {
for(int counterDeg = 0; counterDeg < 5; counterDeg++) {
int counterDegAux = counterDeg+1;
if(counterLayer < 8 && counterDeg == 1) counterDegAux = 0;
int sl = (sl1 + 1)/2;
matrix(5*counterLayer+counterDeg, 5*counterLayer2+counterDeg) = 1.0/(16.0*error[sl][counterDegAux]*error[sl][counterDegAux]);
}
}
}
}
return matrix;
}
| TMatrixD DTMuonMillepede::getLagMatrix | ( | int | wh, |
| int | st, | ||
| int | se | ||
| ) |
Definition at line 466 of file DTMuonMillepede.cc.
References makeMuonMisalignmentScenario::matrix.
{
TMatrixD matrix(60+6,60+6);
if(st == 4) matrix.ResizeTo(40+6, 40+6);
for(int counterDeg = 0; counterDeg < 5; ++counterDeg) {
for(int counterLayer = 0; counterLayer < 12; ++counterLayer) {
if(st == 4) {
matrix(40+counterDeg, counterLayer*5+counterDeg) = 10000.0;
matrix(counterLayer*5+counterDeg, 40+counterDeg) = 10000.0;
} else {
int realCounter = counterDeg+1;
if(counterDeg == 1 && counterLayer < 8) {
realCounter = counterDeg-1;
}
matrix(counterLayer*5+counterDeg,40+realCounter) = 10000.0;
if( (realCounter == 0 && counterLayer > 7) ||
(realCounter == 1 && counterLayer < 7)) continue;
matrix(60+realCounter,counterLayer*5+counterDeg) = 10000.0;
}
}
}
return matrix;
}
| TMatrixD DTMuonMillepede::getMatrixFromFile | ( | 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().
{
TString MtxFileName = "./LocalMillepedeMatrix_"; MtxFileName += mf; MtxFileName += ".root";
if (mf==-1) MtxFileName = "./LocalMillepedeMatrix.root";
TDirectory *dirSave = gDirectory;
TFile *MatrixFile = new TFile(MtxFileName);
TString MtxName = Code; MtxName += wh; MtxName += "_"; MtxName += st; MtxName += "_"; MtxName += se;
TMatrixD *ThisMtx = (TMatrixD *)MatrixFile->Get(MtxName);
MatrixFile->Close();
dirSave->cd();
return *ThisMtx;
}
| TMatrixD DTMuonMillepede::prepareForLagrange | ( | const TMatrixD & | m | ) |
| 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().
{
ttreeOutput = new TTree("DTLocalMillepede", "DTLocalMillepede");
ttreeOutput->Branch("wh", &whC, "wh/I");
ttreeOutput->Branch("st", &stC, "st/I");
ttreeOutput->Branch("sr", &srC, "sr/I");
ttreeOutput->Branch("cov", cov, "cov[60][60]/F");
ttreeOutput->Branch("sl", slC, "sl[12]/I");
ttreeOutput->Branch("la", laC, "la[12]/I");
ttreeOutput->Branch("dx", dx, "dx[12]/F");
ttreeOutput->Branch("dy", dy, "dy[12]/F");
ttreeOutput->Branch("dz", dz, "dz[12]/F");
ttreeOutput->Branch("phix", phix, "phix[12]/F");
ttreeOutput->Branch("phiy", phiy, "phiy[12]/F");
ttreeOutput->Branch("phiz", phiz, "phiz[12]/F");
}
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] |
Reimplemented from DTMuonLocalAlignment.
Definition at line 54 of file DTMuonMillepede.h.
Referenced by 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().
1.7.3