CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Alignment/LaserAlignment/src/LASBarrelAlgorithm.cc

Go to the documentation of this file.
00001 
00002 
00003 #include "Alignment/LaserAlignment/interface/LASBarrelAlgorithm.h"
00004 
00005 
00006 // this is ugly but we need it for Minuit
00007 // until a better solution is at hand
00008 LASGlobalData<LASCoordinateSet>* aMeasuredCoordinates;
00009 LASGlobalData<LASCoordinateSet>* aNominalCoordinates;
00010 
00011 
00012 
00013 
00014 
00018 LASBarrelAlgorithm::LASBarrelAlgorithm() {
00019 
00020   minuit = new TMinuit( 52 );
00021   
00022 }
00023 
00024 
00025 
00026 
00027 
00033 LASBarrelAlignmentParameterSet LASBarrelAlgorithm::CalculateParameters( LASGlobalData<LASCoordinateSet>& measuredCoordinates,
00034                                                                         LASGlobalData<LASCoordinateSet>& nominalCoordinates   ) {
00035   
00036   std::cout << " [LASBarrelAlgorithm::CalculateParameters] -- Starting." << std::endl;
00037 
00038 
00040   // for testing..
00041   //ReadMisalignmentFromFile( "misalign-var.txt", measuredCoordinates, nominalCoordinates );
00043     
00044 
00045   // statics for minuit
00046   aMeasuredCoordinates = &measuredCoordinates;
00047   aNominalCoordinates  = &nominalCoordinates;
00048 
00049   // minimizer and variables for it
00050   minuit->SetFCN( fcn );
00051   double arglist[10];
00052   int _ierflg = 0;
00053 
00054   // toggle minuit blabla
00055   arglist[0] = -1;
00056   minuit->mnexcm( "SET PRI", arglist, 1, _ierflg );
00057 
00058   // set par errors
00059   arglist[0] = 1;
00060   minuit->mnexcm( "SET ERR", arglist ,1, _ierflg );
00061 
00062 
00063   //
00064   // define 52 parameters
00065   //
00066 
00067   // start values: to be evacuated to cfg
00068   static float _vstart[52] = {
00069     0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TIB+
00070     0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TIB-
00071     0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TOB+
00072     0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TOB-
00073     0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TEC+
00074     0.00, 0.00, 0.0, 0.0, 0.0, 0.0, // subdet for TEC-
00075     0.00, 0.00,  0.00, 0.00,  0.00, 0.00,  0.00, 0.00, // beams 0-3
00076     0.00, 0.00,  0.00, 0.00,  0.00, 0.00,  0.00, 0.00  // beams 4-7
00077   };
00078 
00079 
00081   // ReadStartParametersFromFile( "startParameters.txt", _vstart ); // debug
00083 
00084 
00085   // step sizes: to be tuned, to be evacuated to cfg
00086   static float _vstep[52] = { 
00087     0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TIB+
00088     0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TIB-
00089     0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TOB+
00090     0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TOB-
00091     0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TEC+
00092     0.001, 0.001, 0.1, 0.1, 0.1, 0.1, // subdet for TEC-
00093     0.001, 0.001, 0.001, 0.001,  0.001, 0.001,  0.001, 0.001, // beams 0-3
00094     0.001, 0.001, 0.001, 0.001,  0.001, 0.001,  0.001, 0.001  // beams 4-7
00095   };
00096 
00097 
00098   // subdetector parameters for TIB+:
00099 
00100   // rotation around z of first end face 
00101   minuit->mnparm(  0, "subRot1TIB+",      _vstart[0],  _vstep[0],  0, 0, _ierflg );
00102   // rotation around z of second end face
00103   minuit->mnparm(  1, "subRot2TIB+",      _vstart[1],  _vstep[1],  0, 0, _ierflg );
00104   // translation along x of first end face
00105   minuit->mnparm(  2, "subTransX1TIB+",   _vstart[2],  _vstep[2],  0, 0, _ierflg );
00106   // translation along x of second end face
00107   minuit->mnparm(  3, "subTransX2TIB+",   _vstart[3],  _vstep[3],  0, 0, _ierflg );
00108   // translation along y of first end face
00109   minuit->mnparm(  4, "subTransY1TIB+",   _vstart[4],  _vstep[4],  0, 0, _ierflg );
00110   // translation along y of second  end face
00111   minuit->mnparm(  5, "subTransY2TIB+",   _vstart[5],  _vstep[5],  0, 0, _ierflg );
00112 
00113   // subdetector parameters for TIB-:
00114 
00115   // rotation around z of first end face 
00116   minuit->mnparm(  6, "subRot1TIB-",      _vstart[6],  _vstep[6],  0, 0, _ierflg );
00117   // rotation around z of second end face
00118   minuit->mnparm(  7, "subRot2TIB-",      _vstart[7],  _vstep[7],  0, 0, _ierflg );
00119   // translation along x of first end face
00120   minuit->mnparm(  8, "subTransX1TIB-",   _vstart[8],  _vstep[8],  0, 0, _ierflg );
00121   // translation along x of second end face
00122   minuit->mnparm(  9, "subTransX2TIB-",   _vstart[9],  _vstep[9],  0, 0, _ierflg );
00123   // translation along y of first end face
00124   minuit->mnparm( 10, "subTransY1TIB-",   _vstart[10], _vstep[10], 0, 0, _ierflg );
00125   // translation along y of second  end face
00126   minuit->mnparm( 11, "subTransY2TIB-",   _vstart[11], _vstep[11], 0, 0, _ierflg );
00127 
00128   // subdetector parameters for TOB+:
00129 
00130   // rotation around z of first end face 
00131   minuit->mnparm( 12, "subRot1TOB+",      _vstart[12], _vstep[12], 0, 0, _ierflg );
00132   // rotation around z of second end face
00133   minuit->mnparm( 13, "subRot2TOB+",      _vstart[13], _vstep[13], 0, 0, _ierflg );
00134   // translation along x of first end face
00135   minuit->mnparm( 14, "subTransX1TOB+",   _vstart[14], _vstep[14], 0, 0, _ierflg );
00136   // translation along x of second end face
00137   minuit->mnparm( 15, "subTransX2TOB+",   _vstart[15], _vstep[15], 0, 0, _ierflg );
00138   // translation along y of first end face
00139   minuit->mnparm( 16, "subTransY1TOB+",   _vstart[16], _vstep[16], 0, 0, _ierflg );
00140   // translation along y of second  end face
00141   minuit->mnparm( 17, "subTransY2TOB+",   _vstart[17], _vstep[17], 0, 0, _ierflg );
00142 
00143   // subdetector parameters for TOB-:
00144 
00145   // rotation around z of first end face 
00146   minuit->mnparm( 18, "subRot1TOB-",      _vstart[18], _vstep[18], 0, 0, _ierflg );
00147   // rotation around z of second end face
00148   minuit->mnparm( 19, "subRot2TOB-",      _vstart[19], _vstep[19], 0, 0, _ierflg );
00149   // translation along x of first end face
00150   minuit->mnparm( 20, "subTransX1TOB-",   _vstart[20], _vstep[20], 0, 0, _ierflg );
00151   // translation along x of second end face
00152   minuit->mnparm( 21, "subTransX2TOB-",   _vstart[21], _vstep[21], 0, 0, _ierflg );
00153   // translation along y of first end face
00154   minuit->mnparm( 22, "subTransY1TOB-",   _vstart[22], _vstep[22], 0, 0, _ierflg );
00155   // translation along y of second  end face
00156   minuit->mnparm( 23, "subTransY2TOB-",   _vstart[23], _vstep[23], 0, 0, _ierflg );
00157 
00158   // subdetector parameters for TEC+:
00159 
00160   // rotation around z of first end face 
00161   minuit->mnparm( 24, "subRot1TEC+",      _vstart[24], _vstep[24], 0, 0, _ierflg );
00162   // rotation around z of second end face
00163   minuit->mnparm( 25, "subRot2TEC+",      _vstart[25], _vstep[25], 0, 0, _ierflg );
00164   // translation along x of first end face
00165   minuit->mnparm( 26, "subTransX1TEC+",   _vstart[26], _vstep[26], 0, 0, _ierflg );
00166   // translation along x of second end face
00167   minuit->mnparm( 27, "subTransX2TEC+",   _vstart[27], _vstep[27], 0, 0, _ierflg );
00168   // translation along y of first end face
00169   minuit->mnparm( 28, "subTransY1TEC+",   _vstart[28], _vstep[28], 0, 0, _ierflg );
00170   // translation along y of second  end face
00171   minuit->mnparm( 29, "subTransY2TEC+",   _vstart[29], _vstep[29], 0, 0, _ierflg );
00172 
00173   // subdetector parameters for TEC-:
00174 
00175   // rotation around z of first end face 
00176   minuit->mnparm( 30, "subRot1TEC-",      _vstart[30], _vstep[30], 0, 0, _ierflg );
00177   // rotation around z of second end face
00178   minuit->mnparm( 31, "subRot2TEC-",      _vstart[31], _vstep[31], 0, 0, _ierflg );
00179   // translation along x of first end face
00180   minuit->mnparm( 32, "subTransX1TEC-",   _vstart[32], _vstep[32], 0, 0, _ierflg );
00181   // translation along x of second end face
00182   minuit->mnparm( 33, "subTransX2TEC-",   _vstart[33], _vstep[33], 0, 0, _ierflg );
00183   // translation along y of first end face
00184   minuit->mnparm( 34, "subTransY1TEC-",   _vstart[34], _vstep[34], 0, 0, _ierflg );
00185   // translation along y of second  end face
00186   minuit->mnparm( 35, "subTransY2TEC-",   _vstart[35], _vstep[35], 0, 0, _ierflg );
00187 
00188 
00189   // beam parameters (+-z pairs, duplicated for beams 0-7):
00190 
00191   // rotation around z at zt1
00192   minuit->mnparm( 36, "beamRot1Beam0",    _vstart[36],  _vstep[36], 0, 0, _ierflg );
00193   // rotation around z at zt2
00194   minuit->mnparm( 37, "beamRot2Beam0",    _vstart[37],  _vstep[37], 0, 0, _ierflg );
00195   
00196   // rotation around z at zt1
00197   minuit->mnparm( 38, "beamRot1Beam1",    _vstart[38],  _vstep[38], 0, 0, _ierflg );
00198   // rotation around z at zt2
00199   minuit->mnparm( 39, "beamRot2Beam1",    _vstart[39],  _vstep[39], 0, 0, _ierflg );
00200 
00201   // rotation around z at zt1
00202   minuit->mnparm( 40, "beamRot1Beam2",    _vstart[40],  _vstep[40], 0, 0, _ierflg );
00203   // rotation around z at zt2
00204   minuit->mnparm( 41, "beamRot2Beam2",    _vstart[41],  _vstep[41], 0, 0, _ierflg );
00205 
00206   // rotation around z at zt1
00207   minuit->mnparm( 42, "beamRot1Beam3",    _vstart[42],  _vstep[42], 0, 0, _ierflg );
00208   // rotation around z at zt2
00209   minuit->mnparm( 43, "beamRot2Beam3",    _vstart[43],  _vstep[43], 0, 0, _ierflg );
00210 
00211   // rotation around z at zt1
00212   minuit->mnparm( 44, "beamRot1Beam4",    _vstart[44],  _vstep[44], 0, 0, _ierflg );
00213   // rotation around z at zt2
00214   minuit->mnparm( 45, "beamRot2Beam4",    _vstart[45],  _vstep[45], 0, 0, _ierflg );
00215   
00216   // rotation around z at zt1
00217   minuit->mnparm( 46, "beamRot1Beam5",    _vstart[46],  _vstep[46], 0, 0, _ierflg );
00218   // rotation around z at zt2
00219   minuit->mnparm( 47, "beamRot2Beam5",    _vstart[47],  _vstep[47], 0, 0, _ierflg );
00220 
00221   // rotation around z at zt1
00222   minuit->mnparm( 48, "beamRot1Beam6",    _vstart[48],  _vstep[48], 0, 0, _ierflg );
00223   // rotation around z at zt2
00224   minuit->mnparm( 49, "beamRot2Beam6",    _vstart[49],  _vstep[49], 0, 0, _ierflg );
00225 
00226   // rotation around z at zt1
00227   minuit->mnparm( 50, "beamRot1Beam7",    _vstart[50],  _vstep[50], 0, 0, _ierflg );
00228   // rotation around z at zt2
00229   minuit->mnparm( 51, "beamRot2Beam7",    _vstart[51],  _vstep[51], 0, 0, _ierflg );
00230 
00231 
00232   // we fix the respective outer disks 9 of each endcap
00233   // as a reference system (pars 25,27,29 & 30,32,34)
00234   // note: minuit numbering is fortran style...
00235   arglist[0] = 26; arglist[1] = 28; arglist[2] = 30;
00236   //  minuit->mnexcm( "FIX", arglist ,3, _ierflg ); // TEC+
00237   arglist[0] = 31; arglist[1] = 33; arglist[2] = 35;
00238   //  minuit->mnexcm( "FIX", arglist ,3, _ierflg ); // TEC-
00239 
00240 
00241 
00242 
00244   // DEBUG: FIX BEAM PARAMETERS /////////////////////////////////////////////////////////////////////
00245     double parlist[16];
00246     for( int par = 37; par <= 52; ++par ) parlist[par-37] = par;
00247     minuit->mnexcm( "FIX", parlist ,16, _ierflg );
00249 
00251   // DEBUG: FIX ALGN PARAMETERS /////////////////////////////////////////////////////////////////////
00252 //   double parlist[36];
00253 //   for( int par = 1; par <= 36; ++par ) parlist[par-1] = par;
00254 //   minuit->mnexcm( "FIX", parlist ,36, _ierflg );
00256 
00257 
00258 
00259   // now ready for minimization step
00260   arglist[0] = 10000;
00261   arglist[1] = 0.1;
00262   minuit->mnexcm( "MIGRAD", arglist , 2, _ierflg ); // minimizer
00263   //  minuit->mnexcm( "MINOS", arglist , 1, _ierflg ); // error recalculation
00264 
00265   // now fill the result vector.
00266   // turned out that the parameter numbering is stupid, change this later..
00267   LASBarrelAlignmentParameterSet theResult;
00268   double par = 0., parError = 0.;
00269 
00270   // TEC+ rot
00271   minuit->GetParameter( 24, par, parError ); theResult.GetParameter( 0, 0, 0 ).first = par; theResult.GetParameter( 0, 0, 0 ).second = parError;
00272   minuit->GetParameter( 25, par, parError ); theResult.GetParameter( 0, 1, 0 ).first = par; theResult.GetParameter( 0, 1, 0 ).second = parError;
00273   // TEC+ x
00274   minuit->GetParameter( 26, par, parError ); theResult.GetParameter( 0, 0, 1 ).first = par; theResult.GetParameter( 0, 0, 1 ).second = parError;
00275   minuit->GetParameter( 27, par, parError ); theResult.GetParameter( 0, 1, 1 ).first = par; theResult.GetParameter( 0, 1, 1 ).second = parError;
00276   // TEC+ x
00277   minuit->GetParameter( 28, par, parError ); theResult.GetParameter( 0, 0, 2 ).first = par; theResult.GetParameter( 0, 0, 2 ).second = parError;
00278   minuit->GetParameter( 29, par, parError ); theResult.GetParameter( 0, 1, 2 ).first = par; theResult.GetParameter( 0, 1, 2 ).second = parError;
00279 
00280   // TEC- rot
00281   minuit->GetParameter( 30, par, parError ); theResult.GetParameter( 1, 0, 0 ).first = par; theResult.GetParameter( 1, 0, 0 ).second = parError;
00282   minuit->GetParameter( 31, par, parError ); theResult.GetParameter( 1, 1, 0 ).first = par; theResult.GetParameter( 1, 1, 0 ).second = parError;
00283   // TEC- x
00284   minuit->GetParameter( 32, par, parError ); theResult.GetParameter( 1, 0, 1 ).first = par; theResult.GetParameter( 1, 0, 1 ).second = parError;
00285   minuit->GetParameter( 33, par, parError ); theResult.GetParameter( 1, 1, 1 ).first = par; theResult.GetParameter( 1, 1, 1 ).second = parError;
00286   // TEC- x
00287   minuit->GetParameter( 34, par, parError ); theResult.GetParameter( 1, 0, 2 ).first = par; theResult.GetParameter( 1, 0, 2 ).second = parError;
00288   minuit->GetParameter( 35, par, parError ); theResult.GetParameter( 1, 1, 2 ).first = par; theResult.GetParameter( 1, 1, 2 ).second = parError;
00289 
00290   // TIB+ rot
00291   minuit->GetParameter(  0, par, parError ); theResult.GetParameter( 2, 0, 0 ).first = par; theResult.GetParameter( 2, 0, 0 ).second = parError;
00292   minuit->GetParameter(  1, par, parError ); theResult.GetParameter( 2, 1, 0 ).first = par; theResult.GetParameter( 2, 1, 0 ).second = parError;
00293   // TIB+ x
00294   minuit->GetParameter(  2, par, parError ); theResult.GetParameter( 2, 0, 1 ).first = par; theResult.GetParameter( 2, 0, 1 ).second = parError;
00295   minuit->GetParameter(  3, par, parError ); theResult.GetParameter( 2, 1, 1 ).first = par; theResult.GetParameter( 2, 1, 1 ).second = parError;
00296   // TIB+ x
00297   minuit->GetParameter(  4, par, parError ); theResult.GetParameter( 2, 0, 2 ).first = par; theResult.GetParameter( 2, 0, 2 ).second = parError;
00298   minuit->GetParameter(  5, par, parError ); theResult.GetParameter( 2, 1, 2 ).first = par; theResult.GetParameter( 2, 1, 2 ).second = parError;
00299 
00300   // TIB- rot
00301   minuit->GetParameter(  6, par, parError ); theResult.GetParameter( 3, 0, 0 ).first = par; theResult.GetParameter( 3, 0, 0 ).second = parError;
00302   minuit->GetParameter(  7, par, parError ); theResult.GetParameter( 3, 1, 0 ).first = par; theResult.GetParameter( 3, 1, 0 ).second = parError;
00303   // TIB- x
00304   minuit->GetParameter(  8, par, parError ); theResult.GetParameter( 3, 0, 1 ).first = par; theResult.GetParameter( 3, 0, 1 ).second = parError;
00305   minuit->GetParameter(  9, par, parError ); theResult.GetParameter( 3, 1, 1 ).first = par; theResult.GetParameter( 3, 1, 1 ).second = parError;
00306   // TIB- x
00307   minuit->GetParameter( 10, par, parError ); theResult.GetParameter( 3, 0, 2 ).first = par; theResult.GetParameter( 3, 0, 2 ).second = parError;
00308   minuit->GetParameter( 11, par, parError ); theResult.GetParameter( 3, 1, 2 ).first = par; theResult.GetParameter( 3, 1, 2 ).second = parError;
00309 
00310   // TOB+ rot
00311   minuit->GetParameter( 12, par, parError ); theResult.GetParameter( 4, 0, 0 ).first = par; theResult.GetParameter( 4, 0, 0 ).second = parError;
00312   minuit->GetParameter( 13, par, parError ); theResult.GetParameter( 4, 1, 0 ).first = par; theResult.GetParameter( 4, 1, 0 ).second = parError;
00313   // TOB+ x
00314   minuit->GetParameter( 14, par, parError ); theResult.GetParameter( 4, 0, 1 ).first = par; theResult.GetParameter( 4, 0, 1 ).second = parError;
00315   minuit->GetParameter( 15, par, parError ); theResult.GetParameter( 4, 1, 1 ).first = par; theResult.GetParameter( 4, 1, 1 ).second = parError;
00316   // TOB+ x
00317   minuit->GetParameter( 16, par, parError ); theResult.GetParameter( 4, 0, 2 ).first = par; theResult.GetParameter( 4, 0, 2 ).second = parError;
00318   minuit->GetParameter( 17, par, parError ); theResult.GetParameter( 4, 1, 2 ).first = par; theResult.GetParameter( 4, 1, 2 ).second = parError;
00319 
00320   // TOB- rot
00321   minuit->GetParameter( 18, par, parError ); theResult.GetParameter( 5, 0, 0 ).first = par; theResult.GetParameter( 5, 0, 0 ).second = parError;
00322   minuit->GetParameter( 19, par, parError ); theResult.GetParameter( 5, 1, 0 ).first = par; theResult.GetParameter( 5, 1, 0 ).second = parError;
00323   // TOB- x
00324   minuit->GetParameter( 20, par, parError ); theResult.GetParameter( 5, 0, 1 ).first = par; theResult.GetParameter( 5, 0, 1 ).second = parError;
00325   minuit->GetParameter( 21, par, parError ); theResult.GetParameter( 5, 1, 1 ).first = par; theResult.GetParameter( 5, 1, 1 ).second = parError;
00326   // TOB- x
00327   minuit->GetParameter( 22, par, parError ); theResult.GetParameter( 5, 0, 2 ).first = par; theResult.GetParameter( 5, 0, 2 ).second = parError;
00328   minuit->GetParameter( 23, par, parError ); theResult.GetParameter( 5, 1, 2 ).first = par; theResult.GetParameter( 5, 1, 2 ).second = parError;
00329 
00330   std::cout << " [LASBarrelAlgorithm::CalculateParameters] -- Done." << std::endl;
00331 
00332   return theResult;
00333 
00334 }
00335 
00336 
00337 
00338 
00339 
00343 void fcn( int &npar, double *gin, double &f, double *par, int iflag )  {
00344 
00345   double chisquare = 0.;
00346 
00347   // the loop object and its variables
00348   LASGlobalLoop moduleLoop;
00349   int det, beam, pos, disk;
00350 
00352   // ADJUST THIS ALSO IN LASGeometryUpdater
00354 
00355   // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
00356   // parameters are: det, side(0=+/1=-), z(lower/upper). TECs have no side (use side = 0)
00357   std::vector<std::vector<std::vector<double> > > endFaceZPositions( 4, std::vector<std::vector<double> >( 2, std::vector<double>( 2, 0. ) ) );
00358   endFaceZPositions.at( 0 ).at( 0 ).at( 0 ) = 1322.5;  // TEC+, *, disk1 ///
00359   endFaceZPositions.at( 0 ).at( 0 ).at( 1 ) = 2667.5;  // TEC+, *, disk9 /// SIDE INFORMATION
00360   endFaceZPositions.at( 1 ).at( 0 ).at( 0 ) = -2667; // TEC-, *, disk1 /// MEANINGLESS FOR TEC -> USE .at(0)!
00361   endFaceZPositions.at( 1 ).at( 0 ).at( 1 ) = -1322.5; // TEC-, *, disk9 ///
00362   endFaceZPositions.at( 2 ).at( 1 ).at( 0 ) = -700.;  // TIB,  -, small z
00363   endFaceZPositions.at( 2 ).at( 1 ).at( 1 ) = -300.;  // TIB,  -, large z
00364   endFaceZPositions.at( 2 ).at( 0 ).at( 0 ) = 300.;   // TIB,  +, small z
00365   endFaceZPositions.at( 2 ).at( 0 ).at( 1 ) = 700.;   // TIB,  +, large z
00366   endFaceZPositions.at( 3 ).at( 1 ).at( 0 ) = -1090.; // TOB,  -, small z
00367   endFaceZPositions.at( 3 ).at( 1 ).at( 1 ) = -300.;  // TOB,  -, large z
00368   endFaceZPositions.at( 3 ).at( 0 ).at( 0 ) = 300.;   // TOB,  +, small z
00369   endFaceZPositions.at( 3 ).at( 0 ).at( 1 ) = 1090.;  // TOB,  +, large z
00370 
00371   // the z positions of the virtual planes at which the beam parameters are measured
00372   std::vector<double> disk9EndFaceZPositions( 2, 0. );
00373   disk9EndFaceZPositions.at( 0 ) = -2667.5; // TEC- disk9
00374   disk9EndFaceZPositions.at( 1 ) =  2667.5; // TEC+ disk9
00375   
00376   // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
00377   double detReducedZ[2] = { 0., 0. };
00378   // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
00379   double beamReducedZ[2] = { 0., 0. };
00380 
00381   // calculate residual for TIBTOB
00382   det = 2; beam = 0; pos = 0;
00383   do {
00384 
00385     // define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB-
00386     const int theSide = pos<3 ? 0 : 1;
00387 
00388     // this is the path the beam has to travel radially after being reflected 
00389     // by the AT mirrors (TIB:50mm, TOB:36mm) -> used for beam parameters
00390     const double radialOffset = det==2 ? 50. : 36.;
00391 
00392     // reduced module's z position with respect to the subdetector endfaces
00393     detReducedZ[0] = aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ() - endFaceZPositions.at( det ).at( theSide ).at( 0 );
00394     detReducedZ[0] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00395     detReducedZ[1] = endFaceZPositions.at( det ).at( theSide ).at( 1 ) - aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ();
00396     detReducedZ[1] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00397 
00398     // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
00399     beamReducedZ[0] = ( aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ() - radialOffset ) - disk9EndFaceZPositions.at( 0 );
00400     beamReducedZ[0] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00401     beamReducedZ[1] = disk9EndFaceZPositions.at( 1 ) - ( aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ() - radialOffset );
00402     beamReducedZ[1] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00403 
00404 
00405     // phi residual for this module as measured
00406     const double measuredResidual = aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhi() - //&
00407       aNominalCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhi();
00408 
00409     // shortcuts for speed
00410     const double currentPhi = aNominalCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhi();
00411     const double currentR   = aNominalCoordinates->GetTIBTOBEntry( det, beam, pos ).GetR();
00412 
00413     // phi residual for this module calculated from the parameters and nominal coordinates:
00414     // this is the sum over the contributions from all parameters
00415     double calculatedResidual = 0.;
00416 
00417     // note that the contributions ym_{i,j,k} given in the tables in TkLasATModel TWiki
00418     // are defined as R*phi, so here they are divided by the R_j factors (we minimize delta phi)
00419 
00420     // unfortunately, minuit keeps parameters in a 1-dim array,
00421     // so we need to address the correct par[] for the 4 cases TIB+/TIB-/TOB+/TOB-
00422     int indexBase = 0;
00423     if( det == 2 ) { // TIB
00424       if( theSide == 0 ) indexBase = 0; // TIB+
00425       if( theSide == 1 ) indexBase = 6; // TIB-
00426     }
00427     if( det == 3 ) { // TOB
00428       if( theSide == 0 ) indexBase = 12; // TOB+
00429       if( theSide == 1 ) indexBase = 18; // TOB-
00430     }
00431 
00432     // par[0] ("subRot1"): rotation around z of first end face
00433     calculatedResidual += detReducedZ[1] * par[indexBase+0];
00434     
00435     // par[1] ("subRot2"): rotation around z of second end face
00436     calculatedResidual += detReducedZ[0] * par[indexBase+1];
00437     
00438     // par[2] ("subTransX1"): translation along x of first end face
00439     calculatedResidual += detReducedZ[1] * sin( currentPhi ) / currentR * par[indexBase+2];
00440 
00441     // par[3] ("subTransX2"): translation along x of second end face
00442     calculatedResidual += detReducedZ[0] * sin( currentPhi ) / currentR * par[indexBase+3];
00443 
00444     // par[4] ("subTransY1"): translation along y of first end face
00445     calculatedResidual += -1. * detReducedZ[1] * cos( currentPhi ) / currentR * par[indexBase+4];
00446 
00447     // par[5] ("subTransY2"): translation along y of second end face
00448     calculatedResidual += -1. * detReducedZ[0] * cos( currentPhi ) / currentR * par[indexBase+5];
00449 
00450 
00451     // now come the 8*2 beam parameters, calculate the respective parameter index base first (-> which beam)
00452     indexBase = 36 + beam * 2;
00453 
00454     // (there's no TIB/TOB/+/- distinction here for the beams)
00455 
00456     // ("beamRot1"): rotation around z at zt1
00457     calculatedResidual += beamReducedZ[1] * par[indexBase];
00458 
00459     // ("beamRot2"): rotation around z at zt2
00460     calculatedResidual += beamReducedZ[0] * par[indexBase+1];
00461  
00462 
00463     // now calculate the chisquare
00464     chisquare += pow( measuredResidual - calculatedResidual, 2 ) / pow( aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhiError(), 2 );
00465 
00466   } while( moduleLoop.TIBTOBLoop( det, beam, pos ) );
00467 
00468  
00469 
00470 
00471 
00472   // calculate residual for TEC AT
00473   det = 0; beam = 0; disk = 0;
00474   do {
00475     
00476     // define the side: TECs sides already disentangled by the "det" index, so fix this to zero
00477     const int theSide = 0;
00478     
00479     // reduced module's z position with respect to the subdetector endfaces
00480     detReducedZ[0] = aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ() - endFaceZPositions.at( det ).at( theSide ).at( 0 );
00481     detReducedZ[0] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00482     detReducedZ[1] = endFaceZPositions.at( det ).at( theSide ).at( 1 ) - aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ();
00483     detReducedZ[1] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00484 
00485     // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
00486     beamReducedZ[0] = aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ() - disk9EndFaceZPositions.at( 0 );
00487     beamReducedZ[0] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00488     beamReducedZ[1] = disk9EndFaceZPositions.at( 1 ) - aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ();
00489     beamReducedZ[1] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00490 
00491     // phi residual for this module as measured
00492     const double measuredResidual = aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhi() - //&
00493       aNominalCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhi();
00494 
00495     // shortcuts for speed
00496     const double currentPhi = aNominalCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhi();
00497     const double currentR   = aNominalCoordinates->GetTEC2TECEntry( det, beam, disk ).GetR();
00498 
00499     // phi residual for this module calculated from the parameters and nominal coordinates:
00500     // this is the sum over the contributions from all parameters
00501     double calculatedResidual = 0.;
00502 
00503     // note that the contributions ym_{i,j,k} given in the tables in TkLasATModel TWiki
00504     // are defined as R*phi, so here they are divided by the R_j factors (we minimize delta phi)
00505 
00506     // there's also a distinction between TEC+/- parameters in situ (det==0 ? <TEC+> : <TEC->)
00507 
00508     // par[0] ("subRot1"): rotation around z of first end face
00509     calculatedResidual += detReducedZ[1] * (det==0 ? par[24] : par[30]);
00510     
00511     // par[1] ("subRot2"): rotation around z of second end face
00512     calculatedResidual += detReducedZ[0] * (det==0 ? par[25] : par[31]);
00513     
00514     // par[2] ("subTransX1"): translation along x of first end face
00515     calculatedResidual += detReducedZ[1] * sin( currentPhi ) * (det==0 ? par[26] : par[32]) / currentR;
00516 
00517     // par[3] ("subTransX2"): translation along x of second end face
00518     calculatedResidual += detReducedZ[0] * sin( currentPhi ) * (det==0 ? par[27] : par[33]) / currentR;
00519 
00520     // par[4] ("subTransY1"): translation along y of first end face
00521     calculatedResidual += -1. * detReducedZ[1] * cos( currentPhi ) * (det==0 ? par[28] : par[34]) / currentR;
00522 
00523     // par[5] ("subTransY2"): translation along y of second end face
00524     calculatedResidual += -1. * detReducedZ[0] * cos( currentPhi ) * (det==0 ? par[29] : par[35]) / currentR;
00525 
00526     // now come the 8*2 beam parameters; calculate the respective parameter index base first (-> which beam)
00527     const unsigned int indexBase = 36 + beam * 2;
00528 
00529     // there's no TEC+/- distinction here
00530 
00531     // par[6] ("beamRot1"): rotation around z at zt1
00532     calculatedResidual += beamReducedZ[1] * par[indexBase];
00533 
00534     // par[7] ("beamRot2"): rotation around z at zt2
00535     calculatedResidual += beamReducedZ[0] * par[indexBase+1];
00536  
00537 
00538     // now calculate the chisquare 
00539     // TODO: check for phi != 0 !!!
00540     chisquare += pow( measuredResidual - calculatedResidual, 2 ) / pow( aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhiError(), 2 );
00541 
00542   } while( moduleLoop.TEC2TECLoop( det, beam, disk ) );
00543 
00544 
00545 
00546   // return the chisquare by ref
00547   f = chisquare;
00548   
00549 }
00550 
00551 
00552 
00553 
00558 void LASBarrelAlgorithm::Dump( void ) {
00559   
00560   if( !minuit ) {
00561     std::cerr << " [LASBarrelAlgorithm::Dump] ** WARNING: minimizer object uninitialized." << std::endl;
00562     return;
00563   }
00564 
00565   std::cout << std::endl << " [LASBarrelAlgorithm::Dump] -- Parameter dump: " << std::endl;
00566 
00567   const int subdetParMap[6] = { 24, 30, 0, 6, 12, 18 }; // map to one-dim array
00568   const std::string subdetNames[6] = { " TEC+  ", " TEC-  ", " TIB+  ", " TIB-  ", " TOB+  ", " TOB-  " };
00569   double value, error;
00570 
00571   std::cout << " Detector parameters: " << std::endl;
00572   std::cout << " -------------------" << std::endl;
00573   std::cout << " Values:     PHI1         X1          Y1         PHI2         X2          Y2   " << std::endl;
00574   for( int subdet = 0; subdet < 6; ++subdet ) {
00575     std::cout <<subdetNames[subdet];
00576     for( int par = subdetParMap[subdet]; par <= subdetParMap[subdet] + 4; par += 2 ) {
00577       minuit->GetParameter( par, value, error );
00578       std::cout << std::setw( 12 ) << std::setprecision( 6 ) << std::fixed << value;
00579     }
00580     for( int par = subdetParMap[subdet] + 1; par <= subdetParMap[subdet] + 5; par += 2 ) {
00581       minuit->GetParameter( par, value, error );
00582       std::cout << std::setw( 12 ) << std::setprecision( 6 ) << std::fixed << value;
00583     }
00584     std::cout << std::endl;
00585   }
00586 
00587   std::cout << " Errors:     PHI1         X1          Y1         PHI2         X2          Y2   " << std::endl;
00588   for( int subdet = 0; subdet < 6; ++subdet ) {
00589     std::cout <<subdetNames[subdet];
00590     for( int par = subdetParMap[subdet]; par <= subdetParMap[subdet] + 4; par += 2 ) {
00591       minuit->GetParameter( par, value, error );
00592       std::cout << std::setw( 12 ) << std::setprecision( 6 ) << std::fixed << error;
00593     }
00594     for( int par = subdetParMap[subdet] + 1; par <= subdetParMap[subdet] + 5; par += 2 ) {
00595       minuit->GetParameter( par, value, error );
00596       std::cout << std::setw( 12 ) << std::setprecision( 6 ) << std::fixed << error;
00597     }
00598     std::cout << std::endl;
00599   }
00600 
00601   std::cout << std::endl;
00602   std::cout << " Beam parameters: " << std::endl;
00603   std::cout << " ---------------" << std::endl;
00604   std::cout << " Values:  PHI1        PHI2" << std::endl;
00605   for( int beam = 0; beam < 8; ++beam ) {
00606     std::cout << " " << beam << "  ";
00607     for( int z = 0; z < 2; ++z ) {
00608       minuit->GetParameter( 36 + 2 * beam + z, value, error );
00609       std::cout << std::setw( 12 ) << std::setprecision( 6 ) << std::fixed << value;
00610     }
00611     std::cout << std::endl;
00612   }  
00613 
00614   std::cout << " Errors:  PHI1        PHI2" << std::endl;
00615   for( int beam = 0; beam < 8; ++beam ) {
00616     std::cout << " " << beam << "  ";
00617     for( int z = 0; z < 2; ++z ) {
00618       minuit->GetParameter( 36 + 2 * beam + z, value, error );
00619       std::cout << std::setw( 12 ) << std::setprecision( 6 ) << std::fixed << error;
00620     }
00621     std::cout << std::endl;
00622   }
00623 
00624 
00625 
00626   // det parameters once again without leading column (for easy read-in), into a file
00627   ofstream file( "/afs/cern.ch/user/o/olzem/public/parameters_det.txt" );
00628   for( int subdet = 0; subdet < 6; ++subdet ) {
00629     for( int par = subdetParMap[subdet]; par <= subdetParMap[subdet] + 4; par += 2 ) {
00630       minuit->GetParameter( par, value, error );
00631       file << std::setw( 12 ) << std::setprecision( 6 ) << std::fixed << value;
00632     }
00633     for( int par = subdetParMap[subdet] + 1; par <= subdetParMap[subdet] + 5; par += 2 ) {
00634       minuit->GetParameter( par, value, error );
00635       file << std::setw( 12 ) << std::setprecision( 6 ) << std::fixed << value;
00636     }
00637     file << std::endl;
00638   }
00639   file.close();
00640 
00641   // same for beam parameters
00642   file.open( "/afs/cern.ch/user/o/olzem/public/parameters_beam.txt" );
00643   for( int beam = 0; beam < 8; ++beam ) {
00644     for( int z = 0; z < 2; ++z ) {
00645       minuit->GetParameter( 36 + 2 * beam + z, value, error );
00646       file << std::setw( 12 ) << std::setprecision( 6 ) << std::fixed << value;
00647     }
00648     file << std::endl;
00649   }  
00650   file.close();
00651 
00652 
00653   std:: cout << " [LASBarrelAlgorithm::Dump] -- End parameter dump." << std::endl;
00654   std::cout << std::endl;
00655 
00656 }
00657 
00658 
00659 
00660 
00661 
00671 void LASBarrelAlgorithm::ReadMisalignmentFromFile( const char* filename, 
00672                                                    LASGlobalData<LASCoordinateSet>& measuredCoordinates,
00673                                                    LASGlobalData<LASCoordinateSet>& nominalCoordinates  ) {
00674 
00675   ifstream file( filename );
00676   if( file.bad() ) {
00677     std::cerr << " [LASBarrelAlgorithm::ReadMisalignmentFromFile] ** ERROR: cannot open file \"" << filename << "\"." << std::endl;
00678     return;
00679   }
00680 
00681   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00682   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00683   std::cerr << " [LASBarrelAlgorithm::ReadMisalignmentFromFile] ** WARNING: you are reading a fake measurement from a file!" << std::endl;
00684   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00685   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00686 
00687 
00688   // the measured coordinates will finally be overwritten;
00689   // first, set them to the nominal values
00690   measuredCoordinates = nominalCoordinates;
00691 
00692   // and put large errors on all values;
00693   {
00694     LASGlobalLoop moduleLoop;
00695     int det, ring, beam, disk, pos;
00696     
00697     det = 0; ring = 0; beam = 0; disk = 0;
00698     do {
00699       measuredCoordinates.GetTECEntry( det, ring, beam, disk ).SetPhiError( 1000. );
00700     } while( moduleLoop.TECLoop( det, ring, beam, disk ) );
00701     
00702     det = 2; beam = 0; pos = 0;
00703     do {
00704       measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).SetPhiError( 1000. );
00705     } while( moduleLoop.TIBTOBLoop( det, beam, pos  ) );
00706     
00707     det = 0; beam = 0; disk = 0;
00708     do {
00709       measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).SetPhiError( 1000. );
00710     } while( moduleLoop.TEC2TECLoop( det, beam, disk ) );
00711   }
00712 
00713 
00714   // buffers for read-in
00715   int det, beam, z, ring;
00716   double phi, phiError;
00717 
00718   while( !file.eof() ) {
00719 
00720     file >> det;
00721     if( file.eof() ) break; // do not read the last line twice, do not fill trash if file empty
00722 
00723     file >> beam; file >> z; file >> ring;
00724     file >> phi; file >> phiError;
00725 
00726     if( det > 1 ) { // TIB/TOB
00727       measuredCoordinates.GetTIBTOBEntry( det, beam, z ).SetPhi( phi );
00728       measuredCoordinates.GetTIBTOBEntry( det, beam, z ).SetPhiError( phiError );
00729     } else { // TEC or TEC(at)
00730       if( ring > -1 ) { // TEC
00731         measuredCoordinates.GetTECEntry( det, ring, beam, z ).SetPhi( phi );
00732         measuredCoordinates.GetTECEntry( det, ring, beam, z ).SetPhiError( phiError );
00733       }
00734       else { // TEC(at)
00735         measuredCoordinates.GetTEC2TECEntry( det, beam, z ).SetPhi( phi );
00736         measuredCoordinates.GetTEC2TECEntry( det, beam, z ).SetPhiError( phiError );
00737       }
00738     }
00739 
00740   }
00741   
00742   file.close();
00743 
00744 }
00745 
00746 
00747 
00748 
00749 
00758 void LASBarrelAlgorithm::ReadStartParametersFromFile( const char* filename, float values[52] ) {
00759   
00760   ifstream file( filename );
00761   if( file.bad() ) {
00762     std::cerr << " [LASBarrelAlgorithm::ReadStartParametersFromFile] ** ERROR: cannot open file \"" << filename << "\"." << std::endl;
00763     return;
00764   }
00765 
00766   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00767   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00768   std::cerr << " [LASBarrelAlgorithm::ReadStartParametersFrom File] ** WARNING: you are reading parameter start values from a file!" << std::endl;
00769   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00770   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00771 
00772   // map to the minuit par array
00773   const int subdetParMap[6] = { 24, 30, 0, 6, 12, 18 };
00774   
00775   for( int det = 0; det < 6; ++det ) {
00776     file >> values[subdetParMap[det]];   // phi1
00777     file >> values[subdetParMap[det]+2]; // x1
00778     file >> values[subdetParMap[det]+4]; // y1
00779     file >> values[subdetParMap[det]+1]; // phi2
00780     file >> values[subdetParMap[det]+3]; // x2
00781     file >> values[subdetParMap[det]+5]; // y2
00782   }
00783   
00784 }