LaserAlignment/src/LASBarrelAlgorithm.cc

0001
0002
0003 #include "Alignment/LaserAlignment/interface/LASBarrelAlgorithm.h"
0004
0005 // this is ugly but we need it for Minuit
0006 // until a better solution is at hand
0007 LASGlobalData<LASCoordinateSet>* aMeasuredCoordinates;
0008 LASGlobalData<LASCoordinateSet>* aNominalCoordinates;
0009
0010 ///
0011 ///
0012 ///
0013 LASBarrelAlgorithm::LASBarrelAlgorithm() { minuit = new TMinuit(52); }
0014
0015 ///
0016 /// The minimization of the equation system for the barrel.
0017 /// For documentation, please refer to the TkLasATModel TWiki page:
0018 ///   TWiki > CMS Web > CMSTrackerLaserAlignmenSystem > TkLasBarrelAlgorithm > TkLasATModel
0019 ///
0020 LASBarrelAlignmentParameterSet LASBarrelAlgorithm::CalculateParameters(
0021     LASGlobalData<LASCoordinateSet>& measuredCoordinates, LASGlobalData<LASCoordinateSet>& nominalCoordinates) {
0022   std::cout << " [LASBarrelAlgorithm::CalculateParameters] -- Starting." << std::endl;
0023
0024   ///////////////////////////////////////////////////////////////////////////////////////////////////
0025   // for testing..
0026   //ReadMisalignmentFromFile( "misalign-var.txt", measuredCoordinates, nominalCoordinates );
0027   ///////////////////////////////////////////////////////////////////////////////////////////////////
0028
0029   // statics for minuit
0030   aMeasuredCoordinates = &measuredCoordinates;
0031   aNominalCoordinates = &nominalCoordinates;
0032
0033   // minimizer and variables for it
0034   minuit->SetFCN(fcn);
0035   double arglist[10];
0036   int _ierflg = 0;
0037
0038   // toggle minuit blabla
0039   arglist[0] = -1;
0040   minuit->mnexcm("SET PRI", arglist, 1, _ierflg);
0041
0042   // set par errors
0043   arglist[0] = 1;
0044   minuit->mnexcm("SET ERR", arglist, 1, _ierflg);
0045
0046   //
0047   // define 52 parameters
0048   //
0049
0050   // start values: to be evacuated to cfg
0051   static float _vstart[52] = {
0052       0.00, 0.00, 0.0,  0.0,  0.0,  0.0,               // subdet for TIB+
0053       0.00, 0.00, 0.0,  0.0,  0.0,  0.0,               // subdet for TIB-
0054       0.00, 0.00, 0.0,  0.0,  0.0,  0.0,               // subdet for TOB+
0055       0.00, 0.00, 0.0,  0.0,  0.0,  0.0,               // subdet for TOB-
0056       0.00, 0.00, 0.0,  0.0,  0.0,  0.0,               // subdet for TEC+
0057       0.00, 0.00, 0.0,  0.0,  0.0,  0.0,               // subdet for TEC-
0058       0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,  // beams 0-3
0059       0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00   // beams 4-7
0060   };
0061
0062   ///////////////////////////////////////////////////////////////////////////////////////////////////
0063   // ReadStartParametersFromFile( "startParameters.txt", _vstart ); // debug
0064   ///////////////////////////////////////////////////////////////////////////////////////////////////
0065
0066   // step sizes: to be tuned, to be evacuated to cfg
0067   static float _vstep[52] = {
0068       0.001, 0.001, 0.1,   0.1,   0.1,   0.1,                  // subdet for TIB+
0069       0.001, 0.001, 0.1,   0.1,   0.1,   0.1,                  // subdet for TIB-
0070       0.001, 0.001, 0.1,   0.1,   0.1,   0.1,                  // subdet for TOB+
0071       0.001, 0.001, 0.1,   0.1,   0.1,   0.1,                  // subdet for TOB-
0072       0.001, 0.001, 0.1,   0.1,   0.1,   0.1,                  // subdet for TEC+
0073       0.001, 0.001, 0.1,   0.1,   0.1,   0.1,                  // subdet for TEC-
0074       0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001,  // beams 0-3
0075       0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001   // beams 4-7
0076   };
0077
0078   // subdetector parameters for TIB+:
0079
0080   // rotation around z of first end face
0081   minuit->mnparm(0, "subRot1TIB+", _vstart[0], _vstep[0], 0, 0, _ierflg);
0082   // rotation around z of second end face
0083   minuit->mnparm(1, "subRot2TIB+", _vstart[1], _vstep[1], 0, 0, _ierflg);
0084   // translation along x of first end face
0085   minuit->mnparm(2, "subTransX1TIB+", _vstart[2], _vstep[2], 0, 0, _ierflg);
0086   // translation along x of second end face
0087   minuit->mnparm(3, "subTransX2TIB+", _vstart[3], _vstep[3], 0, 0, _ierflg);
0088   // translation along y of first end face
0089   minuit->mnparm(4, "subTransY1TIB+", _vstart[4], _vstep[4], 0, 0, _ierflg);
0090   // translation along y of second  end face
0091   minuit->mnparm(5, "subTransY2TIB+", _vstart[5], _vstep[5], 0, 0, _ierflg);
0092
0093   // subdetector parameters for TIB-:
0094
0095   // rotation around z of first end face
0096   minuit->mnparm(6, "subRot1TIB-", _vstart[6], _vstep[6], 0, 0, _ierflg);
0097   // rotation around z of second end face
0098   minuit->mnparm(7, "subRot2TIB-", _vstart[7], _vstep[7], 0, 0, _ierflg);
0099   // translation along x of first end face
0100   minuit->mnparm(8, "subTransX1TIB-", _vstart[8], _vstep[8], 0, 0, _ierflg);
0101   // translation along x of second end face
0102   minuit->mnparm(9, "subTransX2TIB-", _vstart[9], _vstep[9], 0, 0, _ierflg);
0103   // translation along y of first end face
0104   minuit->mnparm(10, "subTransY1TIB-", _vstart[10], _vstep[10], 0, 0, _ierflg);
0105   // translation along y of second  end face
0106   minuit->mnparm(11, "subTransY2TIB-", _vstart[11], _vstep[11], 0, 0, _ierflg);
0107
0108   // subdetector parameters for TOB+:
0109
0110   // rotation around z of first end face
0111   minuit->mnparm(12, "subRot1TOB+", _vstart[12], _vstep[12], 0, 0, _ierflg);
0112   // rotation around z of second end face
0113   minuit->mnparm(13, "subRot2TOB+", _vstart[13], _vstep[13], 0, 0, _ierflg);
0114   // translation along x of first end face
0115   minuit->mnparm(14, "subTransX1TOB+", _vstart[14], _vstep[14], 0, 0, _ierflg);
0116   // translation along x of second end face
0117   minuit->mnparm(15, "subTransX2TOB+", _vstart[15], _vstep[15], 0, 0, _ierflg);
0118   // translation along y of first end face
0119   minuit->mnparm(16, "subTransY1TOB+", _vstart[16], _vstep[16], 0, 0, _ierflg);
0120   // translation along y of second  end face
0121   minuit->mnparm(17, "subTransY2TOB+", _vstart[17], _vstep[17], 0, 0, _ierflg);
0122
0123   // subdetector parameters for TOB-:
0124
0125   // rotation around z of first end face
0126   minuit->mnparm(18, "subRot1TOB-", _vstart[18], _vstep[18], 0, 0, _ierflg);
0127   // rotation around z of second end face
0128   minuit->mnparm(19, "subRot2TOB-", _vstart[19], _vstep[19], 0, 0, _ierflg);
0129   // translation along x of first end face
0130   minuit->mnparm(20, "subTransX1TOB-", _vstart[20], _vstep[20], 0, 0, _ierflg);
0131   // translation along x of second end face
0132   minuit->mnparm(21, "subTransX2TOB-", _vstart[21], _vstep[21], 0, 0, _ierflg);
0133   // translation along y of first end face
0134   minuit->mnparm(22, "subTransY1TOB-", _vstart[22], _vstep[22], 0, 0, _ierflg);
0135   // translation along y of second  end face
0136   minuit->mnparm(23, "subTransY2TOB-", _vstart[23], _vstep[23], 0, 0, _ierflg);
0137
0138   // subdetector parameters for TEC+:
0139
0140   // rotation around z of first end face
0141   minuit->mnparm(24, "subRot1TEC+", _vstart[24], _vstep[24], 0, 0, _ierflg);
0142   // rotation around z of second end face
0143   minuit->mnparm(25, "subRot2TEC+", _vstart[25], _vstep[25], 0, 0, _ierflg);
0144   // translation along x of first end face
0145   minuit->mnparm(26, "subTransX1TEC+", _vstart[26], _vstep[26], 0, 0, _ierflg);
0146   // translation along x of second end face
0147   minuit->mnparm(27, "subTransX2TEC+", _vstart[27], _vstep[27], 0, 0, _ierflg);
0148   // translation along y of first end face
0149   minuit->mnparm(28, "subTransY1TEC+", _vstart[28], _vstep[28], 0, 0, _ierflg);
0150   // translation along y of second  end face
0151   minuit->mnparm(29, "subTransY2TEC+", _vstart[29], _vstep[29], 0, 0, _ierflg);
0152
0153   // subdetector parameters for TEC-:
0154
0155   // rotation around z of first end face
0156   minuit->mnparm(30, "subRot1TEC-", _vstart[30], _vstep[30], 0, 0, _ierflg);
0157   // rotation around z of second end face
0158   minuit->mnparm(31, "subRot2TEC-", _vstart[31], _vstep[31], 0, 0, _ierflg);
0159   // translation along x of first end face
0160   minuit->mnparm(32, "subTransX1TEC-", _vstart[32], _vstep[32], 0, 0, _ierflg);
0161   // translation along x of second end face
0162   minuit->mnparm(33, "subTransX2TEC-", _vstart[33], _vstep[33], 0, 0, _ierflg);
0163   // translation along y of first end face
0164   minuit->mnparm(34, "subTransY1TEC-", _vstart[34], _vstep[34], 0, 0, _ierflg);
0165   // translation along y of second  end face
0166   minuit->mnparm(35, "subTransY2TEC-", _vstart[35], _vstep[35], 0, 0, _ierflg);
0167
0168   // beam parameters (+-z pairs, duplicated for beams 0-7):
0169
0170   // rotation around z at zt1
0171   minuit->mnparm(36, "beamRot1Beam0", _vstart[36], _vstep[36], 0, 0, _ierflg);
0172   // rotation around z at zt2
0173   minuit->mnparm(37, "beamRot2Beam0", _vstart[37], _vstep[37], 0, 0, _ierflg);
0174
0175   // rotation around z at zt1
0176   minuit->mnparm(38, "beamRot1Beam1", _vstart[38], _vstep[38], 0, 0, _ierflg);
0177   // rotation around z at zt2
0178   minuit->mnparm(39, "beamRot2Beam1", _vstart[39], _vstep[39], 0, 0, _ierflg);
0179
0180   // rotation around z at zt1
0181   minuit->mnparm(40, "beamRot1Beam2", _vstart[40], _vstep[40], 0, 0, _ierflg);
0182   // rotation around z at zt2
0183   minuit->mnparm(41, "beamRot2Beam2", _vstart[41], _vstep[41], 0, 0, _ierflg);
0184
0185   // rotation around z at zt1
0186   minuit->mnparm(42, "beamRot1Beam3", _vstart[42], _vstep[42], 0, 0, _ierflg);
0187   // rotation around z at zt2
0188   minuit->mnparm(43, "beamRot2Beam3", _vstart[43], _vstep[43], 0, 0, _ierflg);
0189
0190   // rotation around z at zt1
0191   minuit->mnparm(44, "beamRot1Beam4", _vstart[44], _vstep[44], 0, 0, _ierflg);
0192   // rotation around z at zt2
0193   minuit->mnparm(45, "beamRot2Beam4", _vstart[45], _vstep[45], 0, 0, _ierflg);
0194
0195   // rotation around z at zt1
0196   minuit->mnparm(46, "beamRot1Beam5", _vstart[46], _vstep[46], 0, 0, _ierflg);
0197   // rotation around z at zt2
0198   minuit->mnparm(47, "beamRot2Beam5", _vstart[47], _vstep[47], 0, 0, _ierflg);
0199
0200   // rotation around z at zt1
0201   minuit->mnparm(48, "beamRot1Beam6", _vstart[48], _vstep[48], 0, 0, _ierflg);
0202   // rotation around z at zt2
0203   minuit->mnparm(49, "beamRot2Beam6", _vstart[49], _vstep[49], 0, 0, _ierflg);
0204
0205   // rotation around z at zt1
0206   minuit->mnparm(50, "beamRot1Beam7", _vstart[50], _vstep[50], 0, 0, _ierflg);
0207   // rotation around z at zt2
0208   minuit->mnparm(51, "beamRot2Beam7", _vstart[51], _vstep[51], 0, 0, _ierflg);
0209
0210   // we fix the respective outer disks 9 of each endcap
0211   // as a reference system (pars 25,27,29 & 30,32,34)
0212   // note: minuit numbering is fortran style...
0213   arglist[0] = 26;
0214   arglist[1] = 28;
0215   arglist[2] = 30;
0216   //  minuit->mnexcm( "FIX", arglist ,3, _ierflg ); // TEC+
0217   arglist[0] = 31;
0218   arglist[1] = 33;
0219   arglist[2] = 35;
0220   //  minuit->mnexcm( "FIX", arglist ,3, _ierflg ); // TEC-
0221
0222   ///////////////////////////////////////////////////////////////////////////////////////////////////
0223   // DEBUG: FIX BEAM PARAMETERS /////////////////////////////////////////////////////////////////////
0224   double parlist[16];
0225   for (int par = 37; par <= 52; ++par)
0226     parlist[par - 37] = par;
0227   minuit->mnexcm("FIX", parlist, 16, _ierflg);
0228   ///////////////////////////////////////////////////////////////////////////////////////////////////
0229
0230   ///////////////////////////////////////////////////////////////////////////////////////////////////
0231   // DEBUG: FIX ALGN PARAMETERS /////////////////////////////////////////////////////////////////////
0232   //   double parlist[36];
0233   //   for( int par = 1; par <= 36; ++par ) parlist[par-1] = par;
0234   //   minuit->mnexcm( "FIX", parlist ,36, _ierflg );
0235   ///////////////////////////////////////////////////////////////////////////////////////////////////
0236
0237   // now ready for minimization step
0238   arglist[0] = 10000;
0239   arglist[1] = 0.1;
0240   minuit->mnexcm("MIGRAD", arglist, 2, _ierflg);  // minimizer
0241   //  minuit->mnexcm( "MINOS", arglist , 1, _ierflg ); // error recalculation
0242
0243   // now fill the result vector.
0244   // turned out that the parameter numbering is stupid, change this later..
0245   LASBarrelAlignmentParameterSet theResult;
0246   double par = 0., parError = 0.;
0247
0248   // TEC+ rot
0249   minuit->GetParameter(24, par, parError);
0250   theResult.GetParameter(0, 0, 0).first = par;
0251   theResult.GetParameter(0, 0, 0).second = parError;
0252   minuit->GetParameter(25, par, parError);
0253   theResult.GetParameter(0, 1, 0).first = par;
0254   theResult.GetParameter(0, 1, 0).second = parError;
0255   // TEC+ x
0256   minuit->GetParameter(26, par, parError);
0257   theResult.GetParameter(0, 0, 1).first = par;
0258   theResult.GetParameter(0, 0, 1).second = parError;
0259   minuit->GetParameter(27, par, parError);
0260   theResult.GetParameter(0, 1, 1).first = par;
0261   theResult.GetParameter(0, 1, 1).second = parError;
0262   // TEC+ x
0263   minuit->GetParameter(28, par, parError);
0264   theResult.GetParameter(0, 0, 2).first = par;
0265   theResult.GetParameter(0, 0, 2).second = parError;
0266   minuit->GetParameter(29, par, parError);
0267   theResult.GetParameter(0, 1, 2).first = par;
0268   theResult.GetParameter(0, 1, 2).second = parError;
0269
0270   // TEC- rot
0271   minuit->GetParameter(30, par, parError);
0272   theResult.GetParameter(1, 0, 0).first = par;
0273   theResult.GetParameter(1, 0, 0).second = parError;
0274   minuit->GetParameter(31, par, parError);
0275   theResult.GetParameter(1, 1, 0).first = par;
0276   theResult.GetParameter(1, 1, 0).second = parError;
0277   // TEC- x
0278   minuit->GetParameter(32, par, parError);
0279   theResult.GetParameter(1, 0, 1).first = par;
0280   theResult.GetParameter(1, 0, 1).second = parError;
0281   minuit->GetParameter(33, par, parError);
0282   theResult.GetParameter(1, 1, 1).first = par;
0283   theResult.GetParameter(1, 1, 1).second = parError;
0284   // TEC- x
0285   minuit->GetParameter(34, par, parError);
0286   theResult.GetParameter(1, 0, 2).first = par;
0287   theResult.GetParameter(1, 0, 2).second = parError;
0288   minuit->GetParameter(35, par, parError);
0289   theResult.GetParameter(1, 1, 2).first = par;
0290   theResult.GetParameter(1, 1, 2).second = parError;
0291
0292   // TIB+ rot
0293   minuit->GetParameter(0, par, parError);
0294   theResult.GetParameter(2, 0, 0).first = par;
0295   theResult.GetParameter(2, 0, 0).second = parError;
0296   minuit->GetParameter(1, par, parError);
0297   theResult.GetParameter(2, 1, 0).first = par;
0298   theResult.GetParameter(2, 1, 0).second = parError;
0299   // TIB+ x
0300   minuit->GetParameter(2, par, parError);
0301   theResult.GetParameter(2, 0, 1).first = par;
0302   theResult.GetParameter(2, 0, 1).second = parError;
0303   minuit->GetParameter(3, par, parError);
0304   theResult.GetParameter(2, 1, 1).first = par;
0305   theResult.GetParameter(2, 1, 1).second = parError;
0306   // TIB+ x
0307   minuit->GetParameter(4, par, parError);
0308   theResult.GetParameter(2, 0, 2).first = par;
0309   theResult.GetParameter(2, 0, 2).second = parError;
0310   minuit->GetParameter(5, par, parError);
0311   theResult.GetParameter(2, 1, 2).first = par;
0312   theResult.GetParameter(2, 1, 2).second = parError;
0313
0314   // TIB- rot
0315   minuit->GetParameter(6, par, parError);
0316   theResult.GetParameter(3, 0, 0).first = par;
0317   theResult.GetParameter(3, 0, 0).second = parError;
0318   minuit->GetParameter(7, par, parError);
0319   theResult.GetParameter(3, 1, 0).first = par;
0320   theResult.GetParameter(3, 1, 0).second = parError;
0321   // TIB- x
0322   minuit->GetParameter(8, par, parError);
0323   theResult.GetParameter(3, 0, 1).first = par;
0324   theResult.GetParameter(3, 0, 1).second = parError;
0325   minuit->GetParameter(9, par, parError);
0326   theResult.GetParameter(3, 1, 1).first = par;
0327   theResult.GetParameter(3, 1, 1).second = parError;
0328   // TIB- x
0329   minuit->GetParameter(10, par, parError);
0330   theResult.GetParameter(3, 0, 2).first = par;
0331   theResult.GetParameter(3, 0, 2).second = parError;
0332   minuit->GetParameter(11, par, parError);
0333   theResult.GetParameter(3, 1, 2).first = par;
0334   theResult.GetParameter(3, 1, 2).second = parError;
0335
0336   // TOB+ rot
0337   minuit->GetParameter(12, par, parError);
0338   theResult.GetParameter(4, 0, 0).first = par;
0339   theResult.GetParameter(4, 0, 0).second = parError;
0340   minuit->GetParameter(13, par, parError);
0341   theResult.GetParameter(4, 1, 0).first = par;
0342   theResult.GetParameter(4, 1, 0).second = parError;
0343   // TOB+ x
0344   minuit->GetParameter(14, par, parError);
0345   theResult.GetParameter(4, 0, 1).first = par;
0346   theResult.GetParameter(4, 0, 1).second = parError;
0347   minuit->GetParameter(15, par, parError);
0348   theResult.GetParameter(4, 1, 1).first = par;
0349   theResult.GetParameter(4, 1, 1).second = parError;
0350   // TOB+ x
0351   minuit->GetParameter(16, par, parError);
0352   theResult.GetParameter(4, 0, 2).first = par;
0353   theResult.GetParameter(4, 0, 2).second = parError;
0354   minuit->GetParameter(17, par, parError);
0355   theResult.GetParameter(4, 1, 2).first = par;
0356   theResult.GetParameter(4, 1, 2).second = parError;
0357
0358   // TOB- rot
0359   minuit->GetParameter(18, par, parError);
0360   theResult.GetParameter(5, 0, 0).first = par;
0361   theResult.GetParameter(5, 0, 0).second = parError;
0362   minuit->GetParameter(19, par, parError);
0363   theResult.GetParameter(5, 1, 0).first = par;
0364   theResult.GetParameter(5, 1, 0).second = parError;
0365   // TOB- x
0366   minuit->GetParameter(20, par, parError);
0367   theResult.GetParameter(5, 0, 1).first = par;
0368   theResult.GetParameter(5, 0, 1).second = parError;
0369   minuit->GetParameter(21, par, parError);
0370   theResult.GetParameter(5, 1, 1).first = par;
0371   theResult.GetParameter(5, 1, 1).second = parError;
0372   // TOB- x
0373   minuit->GetParameter(22, par, parError);
0374   theResult.GetParameter(5, 0, 2).first = par;
0375   theResult.GetParameter(5, 0, 2).second = parError;
0376   minuit->GetParameter(23, par, parError);
0377   theResult.GetParameter(5, 1, 2).first = par;
0378   theResult.GetParameter(5, 1, 2).second = parError;
0379
0380   std::cout << " [LASBarrelAlgorithm::CalculateParameters] -- Done." << std::endl;
0381
0382   return theResult;
0383 }
0384
0385 ///
0386 /// minuit chisquare func
0387 ///
0388 void fcn(int& npar, double* gin, double& f, double* par, int iflag) {
0389   double chisquare = 0.;
0390
0391   // the loop object and its variables
0392   LASGlobalLoop moduleLoop;
0393   int det, beam, pos, disk;
0394
0395   /////////////////////////////////////////////////////////////////////////////
0396   // ADJUST THIS ALSO IN LASGeometryUpdater
0397   /////////////////////////////////////////////////////////////////////////////
0398
0399   // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
0400   // parameters are: det, side(0=+/1=-), z(lower/upper). TECs have no side (use side = 0)
0401   std::vector<std::vector<std::vector<double> > > endFaceZPositions(
0402       4, std::vector<std::vector<double> >(2, std::vector<double>(2, 0.)));
0403   endFaceZPositions.at(0).at(0).at(0) = 1322.5;   // TEC+, *, disk1 ///
0404   endFaceZPositions.at(0).at(0).at(1) = 2667.5;   // TEC+, *, disk9 /// SIDE INFORMATION
0405   endFaceZPositions.at(1).at(0).at(0) = -2667;    // TEC-, *, disk1 /// MEANINGLESS FOR TEC -> USE .at(0)!
0406   endFaceZPositions.at(1).at(0).at(1) = -1322.5;  // TEC-, *, disk9 ///
0407   endFaceZPositions.at(2).at(1).at(0) = -700.;    // TIB,  -, small z
0408   endFaceZPositions.at(2).at(1).at(1) = -300.;    // TIB,  -, large z
0409   endFaceZPositions.at(2).at(0).at(0) = 300.;     // TIB,  +, small z
0410   endFaceZPositions.at(2).at(0).at(1) = 700.;     // TIB,  +, large z
0411   endFaceZPositions.at(3).at(1).at(0) = -1090.;   // TOB,  -, small z
0412   endFaceZPositions.at(3).at(1).at(1) = -300.;    // TOB,  -, large z
0413   endFaceZPositions.at(3).at(0).at(0) = 300.;     // TOB,  +, small z
0414   endFaceZPositions.at(3).at(0).at(1) = 1090.;    // TOB,  +, large z
0415
0416   // the z positions of the virtual planes at which the beam parameters are measured
0417   std::vector<double> disk9EndFaceZPositions(2, 0.);
0418   disk9EndFaceZPositions.at(0) = -2667.5;  // TEC- disk9
0419   disk9EndFaceZPositions.at(1) = 2667.5;   // TEC+ disk9
0420
0421   // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
0422   double detReducedZ[2] = {0., 0.};
0423   // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
0424   double beamReducedZ[2] = {0., 0.};
0425
0426   // calculate residual for TIBTOB
0427   det = 2;
0428   beam = 0;
0429   pos = 0;
0430   do {
0431     // define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB-
0432     const int theSide = pos < 3 ? 0 : 1;
0433
0434     // this is the path the beam has to travel radially after being reflected
0435     // by the AT mirrors (TIB:50mm, TOB:36mm) -> used for beam parameters
0436     const double radialOffset = det == 2 ? 50. : 36.;
0437
0438     // reduced module's z position with respect to the subdetector endfaces
0439     detReducedZ[0] =
0440         aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetZ() - endFaceZPositions.at(det).at(theSide).at(0);
0441     detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0442     detReducedZ[1] =
0443         endFaceZPositions.at(det).at(theSide).at(1) - aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetZ();
0444     detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0445
0446     // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
0447     beamReducedZ[0] =
0448         (aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset) - disk9EndFaceZPositions.at(0);
0449     beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0450     beamReducedZ[1] =
0451         disk9EndFaceZPositions.at(1) - (aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetZ() - radialOffset);
0452     beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0453
0454     // phi residual for this module as measured
0455     const double measuredResidual = aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetPhi() -  //&
0456                                     aNominalCoordinates->GetTIBTOBEntry(det, beam, pos).GetPhi();
0457
0458     // shortcuts for speed
0459     const double currentPhi = aNominalCoordinates->GetTIBTOBEntry(det, beam, pos).GetPhi();
0460     const double currentR = aNominalCoordinates->GetTIBTOBEntry(det, beam, pos).GetR();
0461
0462     // phi residual for this module calculated from the parameters and nominal coordinates:
0463     // this is the sum over the contributions from all parameters
0464     double calculatedResidual = 0.;
0465
0466     // note that the contributions ym_{i,j,k} given in the tables in TkLasATModel TWiki
0467     // are defined as R*phi, so here they are divided by the R_j factors (we minimize delta phi)
0468
0469     // unfortunately, minuit keeps parameters in a 1-dim array,
0470     // so we need to address the correct par[] for the 4 cases TIB+/TIB-/TOB+/TOB-
0471     int indexBase = 0;
0472     if (det == 2) {  // TIB
0473       if (theSide == 0)
0474         indexBase = 0;  // TIB+
0475       if (theSide == 1)
0476         indexBase = 6;  // TIB-
0477     }
0478     if (det == 3) {  // TOB
0479       if (theSide == 0)
0480         indexBase = 12;  // TOB+
0481       if (theSide == 1)
0482         indexBase = 18;  // TOB-
0483     }
0484
0485     // par[0] ("subRot1"): rotation around z of first end face
0486     calculatedResidual += detReducedZ[1] * par[indexBase + 0];
0487
0488     // par[1] ("subRot2"): rotation around z of second end face
0489     calculatedResidual += detReducedZ[0] * par[indexBase + 1];
0490
0491     // par[2] ("subTransX1"): translation along x of first end face
0492     calculatedResidual += detReducedZ[1] * sin(currentPhi) / currentR * par[indexBase + 2];
0493
0494     // par[3] ("subTransX2"): translation along x of second end face
0495     calculatedResidual += detReducedZ[0] * sin(currentPhi) / currentR * par[indexBase + 3];
0496
0497     // par[4] ("subTransY1"): translation along y of first end face
0498     calculatedResidual += -1. * detReducedZ[1] * cos(currentPhi) / currentR * par[indexBase + 4];
0499
0500     // par[5] ("subTransY2"): translation along y of second end face
0501     calculatedResidual += -1. * detReducedZ[0] * cos(currentPhi) / currentR * par[indexBase + 5];
0502
0503     // now come the 8*2 beam parameters, calculate the respective parameter index base first (-> which beam)
0504     indexBase = 36 + beam * 2;
0505
0506     // (there's no TIB/TOB/+/- distinction here for the beams)
0507
0508     // ("beamRot1"): rotation around z at zt1
0509     calculatedResidual += beamReducedZ[1] * par[indexBase];
0510
0511     // ("beamRot2"): rotation around z at zt2
0512     calculatedResidual += beamReducedZ[0] * par[indexBase + 1];
0513
0514     // now calculate the chisquare
0515     chisquare += pow(measuredResidual - calculatedResidual, 2) /
0516                  pow(aMeasuredCoordinates->GetTIBTOBEntry(det, beam, pos).GetPhiError(), 2);
0517
0518   } while (moduleLoop.TIBTOBLoop(det, beam, pos));
0519
0520   // calculate residual for TEC AT
0521   det = 0;
0522   beam = 0;
0523   disk = 0;
0524   do {
0525     // define the side: TECs sides already disentangled by the "det" index, so fix this to zero
0526     const int theSide = 0;
0527
0528     // reduced module's z position with respect to the subdetector endfaces
0529     detReducedZ[0] =
0530         aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetZ() - endFaceZPositions.at(det).at(theSide).at(0);
0531     detReducedZ[0] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0532     detReducedZ[1] =
0533         endFaceZPositions.at(det).at(theSide).at(1) - aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetZ();
0534     detReducedZ[1] /= (endFaceZPositions.at(det).at(theSide).at(1) - endFaceZPositions.at(det).at(theSide).at(0));
0535
0536     // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
0537     beamReducedZ[0] = aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetZ() - disk9EndFaceZPositions.at(0);
0538     beamReducedZ[0] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0539     beamReducedZ[1] = disk9EndFaceZPositions.at(1) - aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetZ();
0540     beamReducedZ[1] /= (disk9EndFaceZPositions.at(1) - disk9EndFaceZPositions.at(0));
0541
0542     // phi residual for this module as measured
0543     const double measuredResidual = aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetPhi() -  //&
0544                                     aNominalCoordinates->GetTEC2TECEntry(det, beam, disk).GetPhi();
0545
0546     // shortcuts for speed
0547     const double currentPhi = aNominalCoordinates->GetTEC2TECEntry(det, beam, disk).GetPhi();
0548     const double currentR = aNominalCoordinates->GetTEC2TECEntry(det, beam, disk).GetR();
0549
0550     // phi residual for this module calculated from the parameters and nominal coordinates:
0551     // this is the sum over the contributions from all parameters
0552     double calculatedResidual = 0.;
0553
0554     // note that the contributions ym_{i,j,k} given in the tables in TkLasATModel TWiki
0555     // are defined as R*phi, so here they are divided by the R_j factors (we minimize delta phi)
0556
0557     // there's also a distinction between TEC+/- parameters in situ (det==0 ? <TEC+> : <TEC->)
0558
0559     // par[0] ("subRot1"): rotation around z of first end face
0560     calculatedResidual += detReducedZ[1] * (det == 0 ? par[24] : par[30]);
0561
0562     // par[1] ("subRot2"): rotation around z of second end face
0563     calculatedResidual += detReducedZ[0] * (det == 0 ? par[25] : par[31]);
0564
0565     // par[2] ("subTransX1"): translation along x of first end face
0566     calculatedResidual += detReducedZ[1] * sin(currentPhi) * (det == 0 ? par[26] : par[32]) / currentR;
0567
0568     // par[3] ("subTransX2"): translation along x of second end face
0569     calculatedResidual += detReducedZ[0] * sin(currentPhi) * (det == 0 ? par[27] : par[33]) / currentR;
0570
0571     // par[4] ("subTransY1"): translation along y of first end face
0572     calculatedResidual += -1. * detReducedZ[1] * cos(currentPhi) * (det == 0 ? par[28] : par[34]) / currentR;
0573
0574     // par[5] ("subTransY2"): translation along y of second end face
0575     calculatedResidual += -1. * detReducedZ[0] * cos(currentPhi) * (det == 0 ? par[29] : par[35]) / currentR;
0576
0577     // now come the 8*2 beam parameters; calculate the respective parameter index base first (-> which beam)
0578     const unsigned int indexBase = 36 + beam * 2;
0579
0580     // there's no TEC+/- distinction here
0581
0582     // par[6] ("beamRot1"): rotation around z at zt1
0583     calculatedResidual += beamReducedZ[1] * par[indexBase];
0584
0585     // par[7] ("beamRot2"): rotation around z at zt2
0586     calculatedResidual += beamReducedZ[0] * par[indexBase + 1];
0587
0588     // now calculate the chisquare
0589     // TODO: check for phi != 0 !!!
0590     chisquare += pow(measuredResidual - calculatedResidual, 2) /
0591                  pow(aMeasuredCoordinates->GetTEC2TECEntry(det, beam, disk).GetPhiError(), 2);
0592
0593   } while (moduleLoop.TEC2TECLoop(det, beam, disk));
0594
0595   // return the chisquare by ref
0596   f = chisquare;
0597 }
0598
0599 ///
0600 /// Print resulting parameters to stdout
0601 /// and to a file - for debugging only
0602 ///
0603 void LASBarrelAlgorithm::Dump(void) {
0604   if (!minuit) {
0605     std::cerr << " [LASBarrelAlgorithm::Dump] ** WARNING: minimizer object uninitialized." << std::endl;
0606     return;
0607   }
0608
0609   std::cout << std::endl << " [LASBarrelAlgorithm::Dump] -- Parameter dump: " << std::endl;
0610
0611   const int subdetParMap[6] = {24, 30, 0, 6, 12, 18};  // map to one-dim array
0612   const std::string subdetNames[6] = {" TEC+  ", " TEC-  ", " TIB+  ", " TIB-  ", " TOB+  ", " TOB-  "};
0613   double value, error;
0614
0615   std::cout << " Detector parameters: " << std::endl;
0616   std::cout << " -------------------" << std::endl;
0617   std::cout << " Values:     PHI1         X1          Y1         PHI2         X2          Y2   " << std::endl;
0618   for (int subdet = 0; subdet < 6; ++subdet) {
0619     std::cout << subdetNames[subdet];
0620     for (int par = subdetParMap[subdet]; par <= subdetParMap[subdet] + 4; par += 2) {
0621       minuit->GetParameter(par, value, error);
0622       std::cout << std::setw(12) << std::setprecision(6) << std::fixed << value;
0623     }
0624     for (int par = subdetParMap[subdet] + 1; par <= subdetParMap[subdet] + 5; par += 2) {
0625       minuit->GetParameter(par, value, error);
0626       std::cout << std::setw(12) << std::setprecision(6) << std::fixed << value;
0627     }
0628     std::cout << std::endl;
0629   }
0630
0631   std::cout << " Errors:     PHI1         X1          Y1         PHI2         X2          Y2   " << std::endl;
0632   for (int subdet = 0; subdet < 6; ++subdet) {
0633     std::cout << subdetNames[subdet];
0634     for (int par = subdetParMap[subdet]; par <= subdetParMap[subdet] + 4; par += 2) {
0635       minuit->GetParameter(par, value, error);
0636       std::cout << std::setw(12) << std::setprecision(6) << std::fixed << error;
0637     }
0638     for (int par = subdetParMap[subdet] + 1; par <= subdetParMap[subdet] + 5; par += 2) {
0639       minuit->GetParameter(par, value, error);
0640       std::cout << std::setw(12) << std::setprecision(6) << std::fixed << error;
0641     }
0642     std::cout << std::endl;
0643   }
0644
0645   std::cout << std::endl;
0646   std::cout << " Beam parameters: " << std::endl;
0647   std::cout << " ---------------" << std::endl;
0648   std::cout << " Values:  PHI1        PHI2" << std::endl;
0649   for (int beam = 0; beam < 8; ++beam) {
0650     std::cout << " " << beam << "  ";
0651     for (int z = 0; z < 2; ++z) {
0652       minuit->GetParameter(36 + 2 * beam + z, value, error);
0653       std::cout << std::setw(12) << std::setprecision(6) << std::fixed << value;
0654     }
0655     std::cout << std::endl;
0656   }
0657
0658   std::cout << " Errors:  PHI1        PHI2" << std::endl;
0659   for (int beam = 0; beam < 8; ++beam) {
0660     std::cout << " " << beam << "  ";
0661     for (int z = 0; z < 2; ++z) {
0662       minuit->GetParameter(36 + 2 * beam + z, value, error);
0663       std::cout << std::setw(12) << std::setprecision(6) << std::fixed << error;
0664     }
0665     std::cout << std::endl;
0666   }
0667
0668   // det parameters once again without leading column (for easy read-in), into a file
0669   std::ofstream file("/afs/cern.ch/user/o/olzem/public/parameters_det.txt");
0670   for (int subdet = 0; subdet < 6; ++subdet) {
0671     for (int par = subdetParMap[subdet]; par <= subdetParMap[subdet] + 4; par += 2) {
0672       minuit->GetParameter(par, value, error);
0673       file << std::setw(12) << std::setprecision(6) << std::fixed << value;
0674     }
0675     for (int par = subdetParMap[subdet] + 1; par <= subdetParMap[subdet] + 5; par += 2) {
0676       minuit->GetParameter(par, value, error);
0677       file << std::setw(12) << std::setprecision(6) << std::fixed << value;
0678     }
0679     file << std::endl;
0680   }
0681   file.close();
0682
0683   // same for beam parameters
0684   file.open("/afs/cern.ch/user/o/olzem/public/parameters_beam.txt");
0685   for (int beam = 0; beam < 8; ++beam) {
0686     for (int z = 0; z < 2; ++z) {
0687       minuit->GetParameter(36 + 2 * beam + z, value, error);
0688       file << std::setw(12) << std::setprecision(6) << std::fixed << value;
0689     }
0690     file << std::endl;
0691   }
0692   file.close();
0693
0694   std::cout << " [LASBarrelAlgorithm::Dump] -- End parameter dump." << std::endl;
0695   std::cout << std::endl;
0696 }
0697
0698 ///
0699 /// allows to push in a simple simulated misalignment for quick internal testing purposes;
0700 /// overwrites LASGlobalData<LASCoordinateSet>& measuredCoordinates;
0701 /// call at beginning of LASBarrelAlgorithm::CalculateParameters method
0702 ///
0703 /// one line per module,
0704 /// format for TEC:              det ring beam disk phi phiErr
0705 /// format for TEC(at) & TIBTOB: det beam   z  "-1" phi phiErr
0706 ///
0707 void LASBarrelAlgorithm::ReadMisalignmentFromFile(const char* filename,
0708                                                   LASGlobalData<LASCoordinateSet>& measuredCoordinates,
0709                                                   LASGlobalData<LASCoordinateSet>& nominalCoordinates) {
0710   std::ifstream file(filename);
0711   if (file.bad()) {
0712     std::cerr << " [LASBarrelAlgorithm::ReadMisalignmentFromFile] ** ERROR: cannot open file \"" << filename << "\"."
0713               << std::endl;
0714     return;
0715   }
0716
0717   std::cerr
0718       << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0719       << std::endl;
0720   std::cerr
0721       << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0722       << std::endl;
0723   std::cerr
0724       << " [LASBarrelAlgorithm::ReadMisalignmentFromFile] ** WARNING: you are reading a fake measurement from a file!"
0725       << std::endl;
0726   std::cerr
0727       << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0728       << std::endl;
0729   std::cerr
0730       << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0731       << std::endl;
0732
0733   // the measured coordinates will finally be overwritten;
0734   // first, set them to the nominal values
0735   measuredCoordinates = nominalCoordinates;
0736
0737   // and put large errors on all values;
0738   {
0739     LASGlobalLoop moduleLoop;
0740     int det, ring, beam, disk, pos;
0741
0742     det = 0;
0743     ring = 0;
0744     beam = 0;
0745     disk = 0;
0746     do {
0747       measuredCoordinates.GetTECEntry(det, ring, beam, disk).SetPhiError(1000.);
0748     } while (moduleLoop.TECLoop(det, ring, beam, disk));
0749
0750     det = 2;
0751     beam = 0;
0752     pos = 0;
0753     do {
0754       measuredCoordinates.GetTIBTOBEntry(det, beam, pos).SetPhiError(1000.);
0755     } while (moduleLoop.TIBTOBLoop(det, beam, pos));
0756
0757     det = 0;
0758     beam = 0;
0759     disk = 0;
0760     do {
0761       measuredCoordinates.GetTEC2TECEntry(det, beam, disk).SetPhiError(1000.);
0762     } while (moduleLoop.TEC2TECLoop(det, beam, disk));
0763   }
0764
0765   // buffers for read-in
0766   int det, beam, z, ring;
0767   double phi, phiError;
0768
0769   while (!file.eof()) {
0770     file >> det;
0771     if (file.eof())
0772       break;  // do not read the last line twice, do not fill trash if file empty
0773
0774     file >> beam;
0775     file >> z;
0776     file >> ring;
0777     file >> phi;
0778     file >> phiError;
0779
0780     if (det > 1) {  // TIB/TOB
0781       measuredCoordinates.GetTIBTOBEntry(det, beam, z).SetPhi(phi);
0782       measuredCoordinates.GetTIBTOBEntry(det, beam, z).SetPhiError(phiError);
0783     } else {            // TEC or TEC(at)
0784       if (ring > -1) {  // TEC
0785         measuredCoordinates.GetTECEntry(det, ring, beam, z).SetPhi(phi);
0786         measuredCoordinates.GetTECEntry(det, ring, beam, z).SetPhiError(phiError);
0787       } else {  // TEC(at)
0788         measuredCoordinates.GetTEC2TECEntry(det, beam, z).SetPhi(phi);
0789         measuredCoordinates.GetTEC2TECEntry(det, beam, z).SetPhiError(phiError);
0790       }
0791     }
0792   }
0793
0794   file.close();
0795 }
0796
0797 ///
0798 /// this function is here only for debugging, don't use.
0799 /// file format:
0800 /// <phi1> <x1> <y1> <phi2> <x2> <y2> // for TEC*
0801 ///   "     "    "     "     "    "   // TEC-
0802 /// .. then for TIB+, TIB-, TOB+, TOB-
0803 /// index 1 if for lower z, 2 for higher z
0804 ///
0805 void LASBarrelAlgorithm::ReadStartParametersFromFile(const char* filename, float values[52]) {
0806   std::ifstream file(filename);
0807   if (file.bad()) {
0808     std::cerr << " [LASBarrelAlgorithm::ReadStartParametersFromFile] ** ERROR: cannot open file \"" << filename << "\"."
0809               << std::endl;
0810     return;
0811   }
0812
0813   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0814                "@@@@@@@@@@@@"
0815             << std::endl;
0816   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0817                "@@@@@@@@@@@@"
0818             << std::endl;
0819   std::cerr << " [LASBarrelAlgorithm::ReadStartParametersFrom File] ** WARNING: you are reading parameter start values "
0820                "from a file!"
0821             << std::endl;
0822   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0823                "@@@@@@@@@@@@"
0824             << std::endl;
0825   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
0826                "@@@@@@@@@@@@"
0827             << std::endl;
0828
0829   // map to the minuit par array
0830   const int subdetParMap[6] = {24, 30, 0, 6, 12, 18};
0831
0832   for (int det = 0; det < 6; ++det) {
0833     file >> values[subdetParMap[det]];      // phi1
0834     file >> values[subdetParMap[det] + 2];  // x1
0835     file >> values[subdetParMap[det] + 4];  // y1
0836     file >> values[subdetParMap[det] + 1];  // phi2
0837     file >> values[subdetParMap[det] + 3];  // x2
0838     file >> values[subdetParMap[det] + 5];  // y2
0839   }
0840 }