Reference Manual

Go to the documentation of this file.
00001 
00002 #include "Alignment/LaserAlignment/interface/LASAlignmentTubeAlgorithm.h"
00003 
00004 
00008 LASAlignmentTubeAlgorithm::LASAlignmentTubeAlgorithm() {
00009 }
00010 
00011 
00012 
00013 
00014 
00018 LASBarrelAlignmentParameterSet LASAlignmentTubeAlgorithm::CalculateParameters( LASGlobalData<LASCoordinateSet>& measuredCoordinates, LASGlobalData<LASCoordinateSet>& nominalCoordinates ) {
00019 
00020   std::cout << " [LASAlignmentTubeAlgorithm::CalculateParameters] -- Starting." << std::endl;
00021 
00022 
00023   // for debugging only
00024   //######################################################################################
00025   //ReadMisalignmentFromFile( "misalign-var.txt", measuredCoordinates, nominalCoordinates );
00026   //######################################################################################
00027 
00028 
00029 
00030   // loop object
00031   LASGlobalLoop globalLoop;
00032   int det, beam, disk, pos;
00033 
00034   // phi positions of the AT beams in rad
00035   const double phiPositions[8] = { 0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784 };
00036   std::vector<double> beamPhiPositions( 8, 0. );
00037   for( beam = 0; beam < 8; ++beam ) beamPhiPositions.at( beam ) = phiPositions[beam];
00038 
00039   // the radii of the alignment tube beams for each halfbarrel.
00040   // the halfbarrels 1-6 are (see TkLasATModel TWiki): TEC+, TEC-, TIB+, TIB-. TOB+, TOB-
00041   // in TIB/TOB modules these radii differ from the beam radius..
00042   // ..due to the radial offsets (after the semitransparent mirrors)
00043   const double radii[6] = { 564., 564., 514., 514., 600., 600. };
00044   std::vector<double> beamRadii( 6, 0. );
00045   for( int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel ) beamRadii.at( aHalfbarrel ) =  radii[aHalfbarrel];
00046   
00047   // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
00048   // parameters are: det, side(0=+/1=-), z(0=lowerZ/1=higherZ). TECs have no side (use side = 0)
00049   std::vector<std::vector<std::vector<double> > > endFaceZPositions( 4, std::vector<std::vector<double> >( 2, std::vector<double>( 2, 0. ) ) );
00050   endFaceZPositions.at( 0 ).at( 0 ).at( 0 ) = 1322.5;  // TEC+, *, disk1 ///
00051   endFaceZPositions.at( 0 ).at( 0 ).at( 1 ) = 2667.5;  // TEC+, *, disk9 /// SIDE INFORMATION
00052   endFaceZPositions.at( 1 ).at( 0 ).at( 0 ) = -2667.5; // TEC-, *, disk9 /// MEANINGLESS FOR TEC -> USE .at(0)!
00053   endFaceZPositions.at( 1 ).at( 0 ).at( 1 ) = -1322.5; // TEC-, *, disk1 ///
00054   endFaceZPositions.at( 2 ).at( 1 ).at( 0 ) = -700.;   // TIB,  -, outer
00055   endFaceZPositions.at( 2 ).at( 1 ).at( 1 ) = -300.;   // TIB,  -, inner
00056   endFaceZPositions.at( 2 ).at( 0 ).at( 0 ) = 300.;    // TIB,  +, inner
00057   endFaceZPositions.at( 2 ).at( 0 ).at( 1 ) = 700.;    // TIB,  +, outer
00058   endFaceZPositions.at( 3 ).at( 1 ).at( 0 ) = -1090.;  // TOB,  -, outer
00059   endFaceZPositions.at( 3 ).at( 1 ).at( 1 ) = -300.;   // TOB,  -, inner
00060   endFaceZPositions.at( 3 ).at( 0 ).at( 0 ) = 300.;    // TOB,  +, inner
00061   endFaceZPositions.at( 3 ).at( 0 ).at( 1 ) = 1090.;   // TOB,  +, outer
00062 
00063   // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
00064   double detReducedZ[2] = { 0., 0. };
00065   // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
00066   double beamReducedZ[2] = { 0., 0. };
00067 
00068   // the z positions of the virtual planes at which the beam parameters are measured
00069   std::vector<double> disk9EndFaceZPositions( 2, 0. );
00070   disk9EndFaceZPositions.at( 0 ) = -2667.5; // TEC- disk9
00071   disk9EndFaceZPositions.at( 1 ) =  2667.5; // TEC+ disk9
00072 
00073 
00074   // define sums over measured values to "simplify" the beam parameter formulas
00075 
00076   // all these have 6 entries, one for each halfbarrel (TEC+,TEC-,TIB+,TIB-,TOB+,TOB-)
00077   std::vector<double> sumOverPhiZPrime( 6, 0. ); 
00078   std::vector<double> sumOverPhiZPrimePrime( 6, 0. ); 
00079   std::vector<double> sumOverPhiZPrimeSinTheta( 6, 0. ); 
00080   std::vector<double> sumOverPhiZPrimePrimeSinTheta( 6, 0. ); 
00081   std::vector<double> sumOverPhiZPrimeCosTheta( 6, 0. ); 
00082   std::vector<double> sumOverPhiZPrimePrimeCosTheta( 6, 0. ); 
00083 
00084   // these have 8 entries, one for each beam
00085   std::vector<double> sumOverPhiZTPrime( 8, 0. ); 
00086   std::vector<double> sumOverPhiZTPrimePrime( 8, 0. ); 
00087 
00088 
00089   // define sums over nominal values
00090 
00091   // all these have 6 entries, one for each halfbarrel (TEC+,TEC-,TIB+,TIB-,TOB+,TOB-)
00092   std::vector<double> sumOverZPrimeSquared( 6, 0. );  
00093   std::vector<double> sumOverZPrimePrimeSquared( 6, 0. ); 
00094   std::vector<double> sumOverZPrimeZPrimePrime( 6, 0. ); 
00095   std::vector<double> sumOverZPrimeZTPrime( 6, 0. ); 
00096   std::vector<double> sumOverZPrimeZTPrimePrime( 6, 0. );
00097   std::vector<double> sumOverZPrimePrimeZTPrime( 6, 0. );
00098   std::vector<double> sumOverZPrimePrimeZTPrimePrime( 6, 0. );
00099 
00100   // all these are scalars
00101   double sumOverZTPrimeSquared = 0.;
00102   double sumOverZTPrimePrimeSquared = 0.;
00103   double sumOverZTPrimeZTPrimePrime = 0.;
00104 
00105 
00106 
00107 
00108 
00109   // now calculate them for TIBTOB
00110   det = 2; beam = 0; pos = 0;
00111   do {
00112     
00113     // define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB-
00114     const int theSide = pos<3 ? 0 : 1;
00115 
00116     // define the halfbarrel number from det/side
00117     const int halfbarrel = det==2 ? det+theSide : det+1+theSide; // TIB:TOB
00118 
00119     // this is the path the beam has to travel radially after being reflected 
00120     // by the AT mirrors (TIB:50mm, TOB:36mm) -> used for beam parameters
00121     const double radialOffset = det==2 ? 50. : 36.;
00122 
00123     // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
00124     detReducedZ[0] = measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).GetZ() - endFaceZPositions.at( det ).at( theSide ).at( 0 ); // = zPrime
00125     detReducedZ[0] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00126     detReducedZ[1] = endFaceZPositions.at( det ).at( theSide ).at( 1 ) - measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).GetZ(); // = zPrimePrime
00127     detReducedZ[1] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00128 
00129     // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
00130     beamReducedZ[0] = ( measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).GetZ() - radialOffset ) - disk9EndFaceZPositions.at( 0 ); // = ZTPrime
00131     beamReducedZ[0] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00132     beamReducedZ[1] = disk9EndFaceZPositions.at( 1 ) - ( measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).GetZ() - radialOffset ); // ZTPrimePrime
00133     beamReducedZ[1] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00134 
00135     // residual in phi (in the formulas this corresponds to y_ik/R)
00136     const double phiResidual = measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).GetPhi() - nominalCoordinates.GetTIBTOBEntry( det, beam, pos ).GetPhi();
00137 
00138     // sums over measured values
00139     sumOverPhiZPrime.at( halfbarrel )               += phiResidual * detReducedZ[0];
00140     sumOverPhiZPrimePrime.at( halfbarrel )          += phiResidual * detReducedZ[1];
00141     sumOverPhiZPrimeSinTheta.at( halfbarrel )       += phiResidual * detReducedZ[0] * sin( beamPhiPositions.at( beam ) );
00142     sumOverPhiZPrimePrimeSinTheta.at( halfbarrel )  += phiResidual * detReducedZ[1] * sin( beamPhiPositions.at( beam ) );
00143     sumOverPhiZPrimeCosTheta.at( halfbarrel )       += phiResidual * detReducedZ[0] * cos( beamPhiPositions.at( beam ) );
00144     sumOverPhiZPrimePrimeCosTheta.at( halfbarrel )  += phiResidual * detReducedZ[1] * cos( beamPhiPositions.at( beam ) );
00145     
00146     sumOverPhiZTPrime.at( beam )                    += phiResidual * beamReducedZ[0]; // note the index change here..
00147     sumOverPhiZTPrimePrime.at( beam )               += phiResidual * beamReducedZ[1];
00148 
00149     // sums over nominal values
00150     sumOverZPrimeSquared.at( halfbarrel )           += pow( detReducedZ[0], 2 ) / 8.; // these are defined beam-wise, so: / 8.
00151     sumOverZPrimePrimeSquared.at( halfbarrel )      += pow( detReducedZ[1], 2 ) / 8.;
00152     sumOverZPrimeZPrimePrime.at( halfbarrel )       += detReducedZ[0] * detReducedZ[1] / 8.;
00153     sumOverZPrimeZTPrime.at( halfbarrel )           += detReducedZ[0] * beamReducedZ[0] / 8.;
00154     sumOverZPrimeZTPrimePrime.at( halfbarrel )      += detReducedZ[0] * beamReducedZ[1] / 8.;
00155     sumOverZPrimePrimeZTPrime.at( halfbarrel )      += detReducedZ[1] * beamReducedZ[0] / 8.;
00156     sumOverZPrimePrimeZTPrimePrime.at( halfbarrel ) += detReducedZ[1] * beamReducedZ[1] / 8.;
00157     
00158     sumOverZTPrimeSquared                           += pow( beamReducedZ[0], 2 ) / 8.;
00159     sumOverZTPrimePrimeSquared                      += pow( beamReducedZ[1], 2 ) / 8.;
00160     sumOverZTPrimeZTPrimePrime                      += beamReducedZ[0] * beamReducedZ[1] / 8.;
00161 
00162   } while( globalLoop.TIBTOBLoop( det, beam, pos ) );
00163   
00164 
00165 
00166 
00167 
00168 
00169 
00170 
00171   // now for TEC2TEC
00172   det = 0; beam = 0; disk = 0;
00173   do {
00174     
00175     // for the tec, the halfbarrel numbers are equal to the det numbers...
00176     const int halfbarrel = det;
00177 
00178     // ...so there's no side distinction for the TEC
00179     const int theSide = 0;
00180 
00181     // also, there's no radial offset for the TEC
00182     const double radialOffset = 0.;
00183 
00184     // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
00185     detReducedZ[0] = measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).GetZ() - endFaceZPositions.at( det ).at( theSide ).at( 0 ); // = zPrime
00186     detReducedZ[0] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00187     detReducedZ[1] = endFaceZPositions.at( det ).at( theSide ).at( 1 ) - measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).GetZ(); // = zPrimePrime
00188     detReducedZ[1] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00189 
00190     // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
00191     beamReducedZ[0] = ( measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).GetZ() - radialOffset ) - disk9EndFaceZPositions.at( 0 ); // = ZTPrime
00192     beamReducedZ[0] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00193     beamReducedZ[1] = disk9EndFaceZPositions.at( 1 ) - ( measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).GetZ() - radialOffset ); // ZTPrimePrime
00194     beamReducedZ[1] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00195 
00196     // residual in phi (in the formulas this corresponds to y_ik/R)
00197     const double phiResidual = measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).GetPhi() - nominalCoordinates.GetTEC2TECEntry( det, beam, disk ).GetPhi();
00198 
00199     // sums over measured values
00200     sumOverPhiZPrime.at( halfbarrel )               += phiResidual * detReducedZ[0];
00201     sumOverPhiZPrimePrime.at( halfbarrel )          += phiResidual * detReducedZ[1];
00202     sumOverPhiZPrimeSinTheta.at( halfbarrel )       += phiResidual * detReducedZ[0] * sin( beamPhiPositions.at( beam ) );
00203     sumOverPhiZPrimePrimeSinTheta.at( halfbarrel )  += phiResidual * detReducedZ[1] * sin( beamPhiPositions.at( beam ) );
00204     sumOverPhiZPrimeCosTheta.at( halfbarrel )       += phiResidual * detReducedZ[0] * cos( beamPhiPositions.at( beam ) );
00205     sumOverPhiZPrimePrimeCosTheta.at( halfbarrel )  += phiResidual * detReducedZ[1] * cos( beamPhiPositions.at( beam ) );
00206     
00207     sumOverPhiZTPrime.at( beam )                    += phiResidual * beamReducedZ[0]; // note the index change here..
00208     sumOverPhiZTPrimePrime.at( beam )               += phiResidual * beamReducedZ[1];
00209 
00210     // sums over nominal values
00211     sumOverZPrimeSquared.at( halfbarrel )           += pow( detReducedZ[0], 2 ) / 8.; // these are defined beam-wise, so: / 8.
00212     sumOverZPrimePrimeSquared.at( halfbarrel )      += pow( detReducedZ[1], 2 ) / 8.;
00213     sumOverZPrimeZPrimePrime.at( halfbarrel )       += detReducedZ[0] * detReducedZ[1] / 8.;
00214     sumOverZPrimeZTPrime.at( halfbarrel )           += detReducedZ[0] * beamReducedZ[0] / 8.;
00215     sumOverZPrimeZTPrimePrime.at( halfbarrel )      += detReducedZ[0] * beamReducedZ[1] / 8.;
00216     sumOverZPrimePrimeZTPrime.at( halfbarrel )      += detReducedZ[1] * beamReducedZ[0] / 8.;
00217     sumOverZPrimePrimeZTPrimePrime.at( halfbarrel ) += detReducedZ[1] * beamReducedZ[1] / 8.;
00218     
00219     sumOverZTPrimeSquared                           += pow( beamReducedZ[0], 2 ) / 8.;
00220     sumOverZTPrimePrimeSquared                      += pow( beamReducedZ[1], 2 ) / 8.;
00221     sumOverZTPrimeZTPrimePrime                      += beamReducedZ[0] * beamReducedZ[1] / 8.;
00222 
00223   } while( globalLoop.TEC2TECLoop( det, beam, disk ) );
00224 
00225 
00226 
00227   // more "simplification" terms...
00228   // these here are functions of theta and can be calculated directly
00229   double sumOverSinTheta = 0.;
00230   double sumOverCosTheta = 0.;
00231   double sumOverSinThetaSquared = 0.;
00232   double sumOverCosThetaSquared = 0.;
00233   double sumOverCosThetaSinTheta = 0.;
00234   double mixedTrigonometricTerm = 0.;
00235 
00236   for( beam = 0; beam < 8; ++beam ) {
00237     sumOverSinTheta += sin( beamPhiPositions.at( beam ) );
00238     sumOverCosTheta += cos( beamPhiPositions.at( beam ) );
00239     sumOverSinThetaSquared += pow( sin( beamPhiPositions.at( beam ) ), 2 );
00240     sumOverCosThetaSquared += pow( cos( beamPhiPositions.at( beam ) ), 2 );
00241     sumOverCosThetaSinTheta += cos( beamPhiPositions.at( beam ) ) * sin( beamPhiPositions.at( beam ) );
00242   }
00243 
00244   mixedTrigonometricTerm = 8. * ( sumOverCosThetaSquared * sumOverSinThetaSquared - pow( sumOverCosThetaSinTheta, 2 ) )
00245     - pow( sumOverCosTheta, 2 ) * sumOverSinThetaSquared
00246     - pow( sumOverSinTheta, 2 ) * sumOverCosThetaSquared
00247     + 2. * sumOverCosTheta * sumOverSinTheta * sumOverCosThetaSinTheta;
00248   
00249   
00250 
00251   // even more shortcuts before we can calculate the parameters
00252   double beamDenominator = ( pow( sumOverZTPrimeZTPrimePrime, 2 ) - sumOverZTPrimeSquared * sumOverZTPrimePrimeSquared ) * beamRadii.at( 0 );
00253   std::vector<double> alignmentDenominator( 6, 0. );
00254   std::vector<double> termA( 6, 0. ), termB( 6, 0. ), termC( 6, 0. ), termD( 6, 0. );
00255   for( unsigned int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel ) {
00256     alignmentDenominator.at ( aHalfbarrel ) = ( pow( sumOverZPrimeZPrimePrime.at( aHalfbarrel ), 2 ) - sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) ) * mixedTrigonometricTerm;
00257     termA.at( aHalfbarrel ) = sumOverZPrimeZTPrime.at( aHalfbarrel ) * sumOverZTPrimeZTPrimePrime - sumOverZPrimeZTPrimePrime.at( aHalfbarrel ) * sumOverZTPrimeSquared;
00258     termB.at( aHalfbarrel ) = sumOverZPrimePrimeZTPrime.at( aHalfbarrel ) * sumOverZTPrimeZTPrimePrime - sumOverZPrimePrimeZTPrimePrime.at( aHalfbarrel ) * sumOverZTPrimeSquared;
00259     termC.at( aHalfbarrel ) = sumOverZPrimeZTPrimePrime.at( aHalfbarrel ) * sumOverZTPrimeZTPrimePrime - sumOverZPrimeZTPrime.at( aHalfbarrel ) * sumOverZTPrimePrimeSquared;
00260     termD.at( aHalfbarrel ) = sumOverZPrimePrimeZTPrimePrime.at( aHalfbarrel ) * sumOverZTPrimeZTPrimePrime - sumOverZPrimePrimeZTPrime.at( aHalfbarrel ) * sumOverZTPrimePrimeSquared;
00261   }
00262 
00263 
00264   // have eight alignment tube beams..
00265   const int numberOfBeams = 8;
00266   
00267 
00268   // that's all for preparation, now it gets ugly:
00269   // calculate the alignment parameters
00270   LASBarrelAlignmentParameterSet theResult;
00271 
00272   // can do this in one go for all halfbarrels
00273   for( int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel ) {
00274 
00275 
00276     // no errors yet
00277 
00278     // rotation angles of the lower z endfaces (in rad)
00279     theResult.GetParameter( aHalfbarrel, 0, 0 ).first = ( sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinThetaSquared
00280                                                           - sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinThetaSquared
00281                                                           - sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinThetaSquared
00282                                                           + sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinThetaSquared
00283                                                           + sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSinTheta * sumOverSinTheta
00284                                                           - sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSinTheta * sumOverSinTheta
00285                                                           - sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinTheta
00286                                                           + sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinTheta
00287                                                           - sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * pow( sumOverCosThetaSinTheta, 2 )
00288                                                           + sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * pow( sumOverCosThetaSinTheta, 2 )
00289                                                           + sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosThetaSinTheta
00290                                                           - sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosThetaSinTheta )
00291                                                         / alignmentDenominator.at ( aHalfbarrel );
00292     
00293     // rotation angles of the upper z endfaces (in rad)
00294     theResult.GetParameter( aHalfbarrel, 1, 0 ).first = ( sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinThetaSquared 
00295                                                           -  sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinThetaSquared 
00296                                                           - sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinThetaSquared
00297                                                           + sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinThetaSquared 
00298                                                           + sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSinTheta * sumOverSinTheta 
00299                                                           - sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSinTheta * sumOverSinTheta
00300                                                           - sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinTheta 
00301                                                           + sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinTheta 
00302                                                           - sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSinTheta * sumOverCosThetaSinTheta
00303                                                           + sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSinTheta * sumOverCosThetaSinTheta 
00304                                                           + sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosThetaSinTheta 
00305                                                           - sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosThetaSinTheta )
00306                                                         / alignmentDenominator.at ( aHalfbarrel );
00307 
00308     // x deviations of the lower z endfaces (in mm)
00309     theResult.GetParameter( aHalfbarrel, 0, 1 ).first = -1. * ( sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinTheta 
00310                                                                 - sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinTheta 
00311                                                                 - sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinTheta
00312                                                                 + sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinTheta 
00313                                                                 - sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosThetaSinTheta 
00314                                                                 + sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosThetaSinTheta
00315                                                                 + sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSinTheta 
00316                                                                 - sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSinTheta 
00317                                                                 - sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSquared
00318                                                                 + sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSquared 
00319                                                                 + sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosTheta 
00320                                                                 - sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosTheta )
00321                                                               / alignmentDenominator.at ( aHalfbarrel ) * beamRadii.at( aHalfbarrel );
00322 
00323     // x deviations of the upper z endfaces (in mm)
00324     theResult.GetParameter( aHalfbarrel, 1, 1 ).first = -1. * ( sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinTheta 
00325                                                                 - sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSquared * sumOverSinTheta 
00326                                                                 - sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinTheta
00327                                                                 + sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinTheta 
00328                                                                 - sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosThetaSinTheta 
00329                                                                 + sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosThetaSinTheta
00330                                                                 + sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSinTheta 
00331                                                                 - sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSinTheta 
00332                                                                 - sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSquared
00333                                                                 + sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSquared 
00334                                                                 + sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosTheta 
00335                                                                 - sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverCosTheta )
00336                                                             / alignmentDenominator.at ( aHalfbarrel ) * beamRadii.at( aHalfbarrel );
00337 
00338     // y deviations of the lower z endfaces (in mm)
00339     theResult.GetParameter( aHalfbarrel, 0, 2 ).first = ( sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinThetaSquared 
00340                                                           - sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinThetaSquared 
00341                                                           - sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * numberOfBeams * sumOverSinThetaSquared
00342                                                           + sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * numberOfBeams * sumOverSinThetaSquared 
00343                                                           + sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverSinTheta * sumOverSinTheta 
00344                                                           - sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverSinTheta * sumOverSinTheta
00345                                                           - sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSinTheta * sumOverSinTheta 
00346                                                           + sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSinTheta * sumOverSinTheta 
00347                                                           - sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinTheta
00348                                                           + sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinTheta 
00349                                                           + sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSinTheta 
00350                                                           - sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeSquared.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSinTheta )
00351                                                         / alignmentDenominator.at ( aHalfbarrel ) * beamRadii.at( aHalfbarrel );
00352 
00353     // y deviations of the upper z endfaces (in mm)
00354     theResult.GetParameter( aHalfbarrel, 1, 2 ).first = ( sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinThetaSquared 
00355                                                           - sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinThetaSquared 
00356                                                           - sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * numberOfBeams * sumOverSinThetaSquared
00357                                                           + sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * numberOfBeams * sumOverSinThetaSquared 
00358                                                           + sumOverPhiZPrimePrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverSinTheta * sumOverSinTheta 
00359                                                           - sumOverPhiZPrimeCosTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverSinTheta * sumOverSinTheta
00360                                                           - sumOverPhiZPrimePrime.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosThetaSinTheta * sumOverSinTheta 
00361                                                           + sumOverPhiZPrime.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosThetaSinTheta * sumOverSinTheta 
00362                                                           - sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinTheta
00363                                                           + sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * sumOverCosTheta * sumOverSinTheta 
00364                                                           + sumOverPhiZPrimePrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimeZPrimePrime.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSinTheta 
00365                                                           - sumOverPhiZPrimeSinTheta.at( aHalfbarrel ) * sumOverZPrimePrimeSquared.at( aHalfbarrel ) * numberOfBeams * sumOverCosThetaSinTheta )
00366                                                         / alignmentDenominator.at ( aHalfbarrel ) * beamRadii.at( aHalfbarrel );
00367 
00368   }
00369 
00370 
00371 
00372 
00373   // another loop is needed here to calculate some terms for the beam parameters
00374   double vsumA = 0., vsumB = 0., vsumC = 0., vsumD = 0., vsumE = 0., vsumF = 0.;
00375   for( unsigned int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel ) {
00376     vsumA += theResult.GetParameter( aHalfbarrel, 1, 2 ).first * termA.at( aHalfbarrel ) + theResult.GetParameter( aHalfbarrel, 0, 2 ).first * termB.at( aHalfbarrel );
00377     vsumB += theResult.GetParameter( aHalfbarrel, 1, 1 ).first * termA.at( aHalfbarrel ) + theResult.GetParameter( aHalfbarrel, 0, 1 ).first * termB.at( aHalfbarrel );
00378     vsumC += beamRadii.at( aHalfbarrel ) * ( theResult.GetParameter( aHalfbarrel, 1, 0 ).first * termA.at( aHalfbarrel ) + theResult.GetParameter( aHalfbarrel, 0, 0 ).first * termB.at( aHalfbarrel ) );
00379     vsumD += theResult.GetParameter( aHalfbarrel, 1, 2 ).first * termC.at( aHalfbarrel ) + theResult.GetParameter( aHalfbarrel, 0, 2 ).first * termD.at( aHalfbarrel );
00380     vsumE += theResult.GetParameter( aHalfbarrel, 1, 1 ).first * termC.at( aHalfbarrel ) + theResult.GetParameter( aHalfbarrel, 0, 1 ).first * termD.at( aHalfbarrel );
00381     vsumF += beamRadii.at( aHalfbarrel ) * ( theResult.GetParameter( aHalfbarrel, 1, 0 ).first * termC.at( aHalfbarrel ) + theResult.GetParameter( aHalfbarrel, 0, 0 ).first * termD.at( aHalfbarrel ) );
00382   }
00383 
00384 
00385 
00386   // calculate the beam parameters
00387   for( unsigned int beam = 0; beam < 8; ++beam ) {
00388 
00389     // parameter A, defined at lower z
00390     theResult.GetBeamParameter( beam, 0 ).first = ( cos( beamPhiPositions.at( beam ) ) * vsumA
00391                                                     - sin( beamPhiPositions.at( beam ) ) * vsumB
00392                                                     - vsumC
00393                                                     + sumOverPhiZTPrime.at( beam ) * sumOverZTPrimeZTPrimePrime - sumOverPhiZTPrimePrime.at( beam ) * sumOverZTPrimeSquared )
00394                                                   / beamDenominator;
00395 
00397     std::cout << "BBBBBBBB: " << cos( beamPhiPositions.at( beam ) ) * vsumA << "  "
00398               << -1. * sin( beamPhiPositions.at( beam ) ) * vsumB << "  "
00399               << -1. * vsumC << "  "
00400               << sumOverPhiZTPrime.at( beam ) * sumOverZTPrimeZTPrimePrime - sumOverPhiZTPrimePrime.at( beam ) * sumOverZTPrimeSquared << "  "
00401               << beamDenominator
00402               << std::endl;
00404 
00405 
00406     // parameter B, defined at upper z
00407     theResult.GetBeamParameter( beam, 1 ).first = ( cos( beamPhiPositions.at( beam ) ) * vsumD
00408                                                     - sin( beamPhiPositions.at( beam ) ) * vsumE
00409                                                     - vsumF
00410                                                     + sumOverPhiZTPrimePrime.at( beam ) * sumOverZTPrimeZTPrimePrime - sumOverPhiZTPrime.at( beam ) * sumOverZTPrimePrimeSquared )
00411                                                   / beamDenominator;
00412 
00413   }
00414 
00415 
00416   return theResult;
00417 
00418 }
00419 
00420 
00421 
00422 
00423 
00428 double LASAlignmentTubeAlgorithm::GetTIBTOBAlignmentParameterCorrection( int det, int beam, int pos, LASGlobalData<LASCoordinateSet>& nominalCoordinates, LASBarrelAlignmentParameterSet& alignmentParameters ) {
00429 
00430   // INITIALIZATION;
00431   // ALL THIS IS DUPLICATED FOR THE MOMENT, SHOULD FINALLY BE CALCULATED ONLY ONCE
00432   // AND HARD CODED NUMBERS SHOULD CENTRALLY BE IMPORTED FROM src/LASConstants.h
00433 
00434 
00435   // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
00436   // parameters are: det, side(0=+/1=-), z(0=lowerZ/1=higherZ). TECs have no side (use side = 0)
00437   std::vector<std::vector<std::vector<double> > > endFaceZPositions( 4, std::vector<std::vector<double> >( 2, std::vector<double>( 2, 0. ) ) );
00438   endFaceZPositions.at( 0 ).at( 0 ).at( 0 ) = 1322.5;  // TEC+, *, disk1 ///
00439   endFaceZPositions.at( 0 ).at( 0 ).at( 1 ) = 2667.5;  // TEC+, *, disk9 /// SIDE INFORMATION
00440   endFaceZPositions.at( 1 ).at( 0 ).at( 0 ) = -2667.5; // TEC-, *, disk9 /// MEANINGLESS FOR TEC -> USE .at(0)!
00441   endFaceZPositions.at( 1 ).at( 0 ).at( 1 ) = -1322.5; // TEC-, *, disk1 ///
00442   endFaceZPositions.at( 2 ).at( 1 ).at( 0 ) = -700.;   // TIB,  -, outer
00443   endFaceZPositions.at( 2 ).at( 1 ).at( 1 ) = -300.;   // TIB,  -, inner
00444   endFaceZPositions.at( 2 ).at( 0 ).at( 0 ) = 300.;    // TIB,  +, inner
00445   endFaceZPositions.at( 2 ).at( 0 ).at( 1 ) = 700.;    // TIB,  +, outer
00446   endFaceZPositions.at( 3 ).at( 1 ).at( 0 ) = -1090.;  // TOB,  -, outer
00447   endFaceZPositions.at( 3 ).at( 1 ).at( 1 ) = -300.;   // TOB,  -, inner
00448   endFaceZPositions.at( 3 ).at( 0 ).at( 0 ) = 300.;    // TOB,  +, inner
00449   endFaceZPositions.at( 3 ).at( 0 ).at( 1 ) = 1090.;   // TOB,  +, outer
00450 
00451   // the z positions of the virtual planes at which the beam parameters are measured
00452   std::vector<double> disk9EndFaceZPositions( 2, 0. );
00453   disk9EndFaceZPositions.at( 0 ) = -2667.5; // TEC- disk9
00454   disk9EndFaceZPositions.at( 1 ) =  2667.5; // TEC+ disk9
00455 
00456   // define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB-
00457   const int theSide = pos<3 ? 0 : 1;
00458   
00459   // define the halfbarrel number from det/side
00460   const int halfbarrel = det==2 ? det+theSide : det+1+theSide; // TIB:TOB
00461   
00462   // this is the path the beam has to travel radially after being reflected 
00463   // by the AT mirrors (TIB:50mm, TOB:36mm) -> used for beam parameters
00464   const double radialOffset = det==2 ? 50. : 36.;
00465 
00466   // phi positions of the AT beams in rad
00467   const double phiPositions[8] = { 0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784 };
00468   std::vector<double> beamPhiPositions( 8, 0. );
00469   for( unsigned int aBeam = 0; aBeam < 8; ++aBeam ) beamPhiPositions.at( aBeam ) = phiPositions[aBeam];
00470 
00471   // the radii of the alignment tube beams for each halfbarrel.
00472   // the halfbarrels 1-6 are (see TkLasATModel TWiki): TEC+, TEC-, TIB+, TIB-. TOB+, TOB-
00473   // in TIB/TOB modules these radii differ from the beam radius..
00474   // ..due to the radial offsets (after the semitransparent mirrors)
00475   const double radii[6] = { 564., 564., 514., 514., 600., 600. };
00476   std::vector<double> beamRadii( 6, 0. );
00477   for( int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel ) beamRadii.at( aHalfbarrel ) =  radii[aHalfbarrel];
00478 
00479   // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
00480   double detReducedZ[2] = { 0., 0. };
00481   // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
00482   double beamReducedZ[2] = { 0., 0. };
00483 
00484   // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
00485   detReducedZ[0] = nominalCoordinates.GetTIBTOBEntry( det, beam, pos ).GetZ() - endFaceZPositions.at( det ).at( theSide ).at( 0 ); // = zPrime
00486   detReducedZ[0] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00487   detReducedZ[1] = endFaceZPositions.at( det ).at( theSide ).at( 1 ) - nominalCoordinates.GetTIBTOBEntry( det, beam, pos ).GetZ(); // = zPrimePrime
00488   detReducedZ[1] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00489   
00490   // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
00491   beamReducedZ[0] = ( nominalCoordinates.GetTIBTOBEntry( det, beam, pos ).GetZ() - radialOffset ) - disk9EndFaceZPositions.at( 0 ); // = ZTPrime
00492   beamReducedZ[0] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00493   beamReducedZ[1] = disk9EndFaceZPositions.at( 1 ) - ( nominalCoordinates.GetTIBTOBEntry( det, beam, pos ).GetZ() - radialOffset ); // ZTPrimePrime
00494   beamReducedZ[1] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00495 
00496   // the correction to phi from the endcap algorithm;
00497   // it is defined such that the correction is to be subtracted ///////////////////////////////// ???
00498   double phiCorrection = 0.;
00499   
00500   // contribution from phi rotation of first end face
00501   phiCorrection += detReducedZ[1] * alignmentParameters.GetParameter( halfbarrel, 0, 0 ).first;
00502 
00503   // contribution from phi rotation of second end face
00504   phiCorrection += detReducedZ[0] * alignmentParameters.GetParameter( halfbarrel, 1, 0 ).first;
00505   
00506   // contribution from translation along x of first endface
00507   phiCorrection += detReducedZ[1] * sin( beamPhiPositions.at( beam ) ) * alignmentParameters.GetParameter( halfbarrel, 0, 1 ).first / beamRadii.at( halfbarrel );
00508 
00509   // contribution from translation along x of second endface
00510   phiCorrection += detReducedZ[0] * sin( beamPhiPositions.at( beam ) ) * alignmentParameters.GetParameter( halfbarrel, 1, 1 ).first / beamRadii.at( halfbarrel );
00511 
00512   // contribution from translation along y of first endface
00513   phiCorrection -= detReducedZ[1] * cos( beamPhiPositions.at( beam ) ) * alignmentParameters.GetParameter( halfbarrel, 0, 2 ).first / beamRadii.at( halfbarrel );
00514 
00515   // contribution from translation along y of second endface
00516   phiCorrection -= detReducedZ[0] * cos( beamPhiPositions.at( beam ) ) * alignmentParameters.GetParameter( halfbarrel, 1, 2 ).first / beamRadii.at( halfbarrel );
00517   
00518   // contribution from beam parameters;
00519   // originally, the contribution in meter is proportional to the radius of the beams: beamRadii.at( 0 )
00520   // the additional factor: beamRadii.at( halfbarrel ) converts from meter to radian on the module
00521   phiCorrection += beamReducedZ[1] * alignmentParameters.GetBeamParameter( beam, 0 ).first * beamRadii.at( 0 ) / beamRadii.at( halfbarrel );
00522   phiCorrection += beamReducedZ[0] * alignmentParameters.GetBeamParameter( beam, 1 ).first * beamRadii.at( 0 ) / beamRadii.at( halfbarrel );
00523 
00524 
00525   return phiCorrection;
00526 
00527 }
00528 
00529 
00530 
00531 
00532 
00537 double LASAlignmentTubeAlgorithm::GetTEC2TECAlignmentParameterCorrection( int det, int beam, int disk, LASGlobalData<LASCoordinateSet>& nominalCoordinates, LASBarrelAlignmentParameterSet& alignmentParameters ) {
00538 
00539   // INITIALIZATION;
00540   // ALL THIS IS DUPLICATED FOR THE MOMENT, SHOULD FINALLY BE CALCULATED ONLY ONCE
00541   // AND HARD CODED NUMBERS SHOULD CENTRALLY BE IMPORTED FROM src/LASConstants.h
00542 
00543 
00544   // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
00545   // parameters are: det, side(0=+/1=-), z(0=lowerZ/1=higherZ). TECs have no side (use side = 0)
00546   std::vector<std::vector<std::vector<double> > > endFaceZPositions( 4, std::vector<std::vector<double> >( 2, std::vector<double>( 2, 0. ) ) );
00547   endFaceZPositions.at( 0 ).at( 0 ).at( 0 ) = 1322.5;  // TEC+, *, disk1 ///
00548   endFaceZPositions.at( 0 ).at( 0 ).at( 1 ) = 2667.5;  // TEC+, *, disk9 /// SIDE INFORMATION
00549   endFaceZPositions.at( 1 ).at( 0 ).at( 0 ) = -2667.5; // TEC-, *, disk9 /// MEANINGLESS FOR TEC -> USE .at(0)!
00550   endFaceZPositions.at( 1 ).at( 0 ).at( 1 ) = -1322.5; // TEC-, *, disk1 ///
00551   endFaceZPositions.at( 2 ).at( 1 ).at( 0 ) = -700.;   // TIB,  -, outer
00552   endFaceZPositions.at( 2 ).at( 1 ).at( 1 ) = -300.;   // TIB,  -, inner
00553   endFaceZPositions.at( 2 ).at( 0 ).at( 0 ) = 300.;    // TIB,  +, inner
00554   endFaceZPositions.at( 2 ).at( 0 ).at( 1 ) = 700.;    // TIB,  +, outer
00555   endFaceZPositions.at( 3 ).at( 1 ).at( 0 ) = -1090.;  // TOB,  -, outer
00556   endFaceZPositions.at( 3 ).at( 1 ).at( 1 ) = -300.;   // TOB,  -, inner
00557   endFaceZPositions.at( 3 ).at( 0 ).at( 0 ) = 300.;    // TOB,  +, inner
00558   endFaceZPositions.at( 3 ).at( 0 ).at( 1 ) = 1090.;   // TOB,  +, outer
00559 
00560   // the z positions of the virtual planes at which the beam parameters are measured
00561   std::vector<double> disk9EndFaceZPositions( 2, 0. );
00562   disk9EndFaceZPositions.at( 0 ) = -2667.5; // TEC- disk9
00563   disk9EndFaceZPositions.at( 1 ) =  2667.5; // TEC+ disk9
00564 
00565   // for the tec, the halfbarrel numbers are equal to the det numbers...
00566   const int halfbarrel = det;
00567   
00568   // ...so there's no side distinction for the TEC
00569   const int theSide = 0;
00570   
00571   // also, there's no radial offset for the TEC
00572   const double radialOffset = 0.;
00573   
00574   // phi positions of the AT beams in rad
00575   const double phiPositions[8] = { 0.392699, 1.289799, 1.851794, 2.748894, 3.645995, 4.319690, 5.216791, 5.778784 };
00576   std::vector<double> beamPhiPositions( 8, 0. );
00577   for( unsigned int aBeam = 0; aBeam < 8; ++aBeam ) beamPhiPositions.at( aBeam ) = phiPositions[aBeam];
00578 
00579   // the radii of the alignment tube beams for each halfbarrel.
00580   // the halfbarrels 1-6 are (see TkLasATModel TWiki): TEC+, TEC-, TIB+, TIB-. TOB+, TOB-
00581   // in TIB/TOB modules these radii differ from the beam radius..
00582   // ..due to the radial offsets (after the semitransparent mirrors)
00583   const double radii[6] = { 564., 564., 514., 514., 600., 600. };
00584   std::vector<double> beamRadii( 6, 0. );
00585   for( int aHalfbarrel = 0; aHalfbarrel < 6; ++aHalfbarrel ) beamRadii.at( aHalfbarrel ) =  radii[aHalfbarrel];
00586 
00587   // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
00588   double detReducedZ[2] = { 0., 0. };
00589   // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
00590   double beamReducedZ[2] = { 0., 0. };
00591 
00592   // reduced module's z position with respect to the subdetector endfaces (zPrime, zPrimePrime)
00593   detReducedZ[0] = nominalCoordinates.GetTEC2TECEntry( det, beam, disk ).GetZ() - endFaceZPositions.at( det ).at( theSide ).at( 0 ); // = zPrime
00594   detReducedZ[0] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00595   detReducedZ[1] = endFaceZPositions.at( det ).at( theSide ).at( 1 ) - nominalCoordinates.GetTEC2TECEntry( det, beam, disk ).GetZ(); // = zPrimePrime
00596   detReducedZ[1] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
00597   
00598   // reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
00599   beamReducedZ[0] = ( nominalCoordinates.GetTEC2TECEntry( det, beam, disk ).GetZ() - radialOffset ) - disk9EndFaceZPositions.at( 0 ); // = ZTPrime
00600   beamReducedZ[0] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00601   beamReducedZ[1] = disk9EndFaceZPositions.at( 1 ) - ( nominalCoordinates.GetTEC2TECEntry( det, beam, disk ).GetZ() - radialOffset ); // ZTPrimePrime
00602   beamReducedZ[1] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
00603 
00604 
00605   // the correction to phi from the endcap algorithm;
00606   // it is defined such that the correction is to be subtracted ///////////////////////////////// ???
00607   double phiCorrection = 0.;
00608   
00609   // contribution from phi rotation of first end face
00610   phiCorrection += detReducedZ[1] * alignmentParameters.GetParameter( halfbarrel, 0, 0 ).first;
00611 
00612   // contribution from phi rotation of second end face
00613   phiCorrection += detReducedZ[0] * alignmentParameters.GetParameter( halfbarrel, 1, 0 ).first;
00614   
00615   // contribution from translation along x of first endface
00616   phiCorrection += detReducedZ[1] * sin( beamPhiPositions.at( beam ) ) * alignmentParameters.GetParameter( halfbarrel, 0, 1 ).first / beamRadii.at( halfbarrel );
00617 
00618   // contribution from translation along x of second endface
00619   phiCorrection += detReducedZ[0] * sin( beamPhiPositions.at( beam ) ) * alignmentParameters.GetParameter( halfbarrel, 1, 1 ).first / beamRadii.at( halfbarrel );
00620 
00621   // contribution from translation along y of first endface
00622   phiCorrection -= detReducedZ[1] * cos( beamPhiPositions.at( beam ) ) * alignmentParameters.GetParameter( halfbarrel, 0, 2 ).first / beamRadii.at( halfbarrel );
00623 
00624   // contribution from translation along y of second endface
00625   phiCorrection -= detReducedZ[0] * cos( beamPhiPositions.at( beam ) ) * alignmentParameters.GetParameter( halfbarrel, 1, 2 ).first / beamRadii.at( halfbarrel );
00626   
00627   // contribution from beam parameters;
00628   // originally, the contribution in meter is proportional to the radius of the beams: beamRadii.at( 0 )
00629   // the additional factor: beamRadii.at( halfbarrel ) converts from meter to radian on the module
00630   phiCorrection += beamReducedZ[1] * alignmentParameters.GetBeamParameter( beam, 0 ).first * beamRadii.at( 0 ) / beamRadii.at( halfbarrel );
00631   phiCorrection += beamReducedZ[0] * alignmentParameters.GetBeamParameter( beam, 1 ).first * beamRadii.at( 0 ) / beamRadii.at( halfbarrel );
00632 
00633 
00634   return phiCorrection;
00635 
00636 }
00637 
00638 
00639 
00640 
00641 
00651 void LASAlignmentTubeAlgorithm::ReadMisalignmentFromFile( const char* filename, 
00652                                                           LASGlobalData<LASCoordinateSet>& measuredCoordinates,
00653                                                           LASGlobalData<LASCoordinateSet>& nominalCoordinates  ) {
00654 
00655   std::ifstream file( filename );
00656   if( file.bad() ) {
00657     std::cerr << " [LASAlignmentTubeAlgorithm::ReadMisalignmentFromFile] ** ERROR: cannot open file \"" << filename << "\"." << std::endl;
00658     return;
00659   }
00660 
00661   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00662   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00663   std::cerr << " [LASAlignmentTubeAlgorithm::ReadMisalignmentFromFile] ** WARNING: you are reading a fake measurement from a file!" << std::endl;
00664   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00665   std::cerr << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
00666 
00667 
00668   // the measured coordinates will finally be overwritten;
00669   // first, set them to the nominal values
00670   measuredCoordinates = nominalCoordinates;
00671 
00672   // and put large errors on all values;
00673   {
00674     LASGlobalLoop moduleLoop;
00675     int det, ring, beam, disk, pos;
00676     
00677     det = 0; ring = 0; beam = 0; disk = 0;
00678     do {
00679       measuredCoordinates.GetTECEntry( det, ring, beam, disk ).SetPhiError( 1000. );
00680     } while( moduleLoop.TECLoop( det, ring, beam, disk ) );
00681     
00682     det = 2; beam = 0; pos = 0;
00683     do {
00684       measuredCoordinates.GetTIBTOBEntry( det, beam, pos ).SetPhiError( 1000. );
00685     } while( moduleLoop.TIBTOBLoop( det, beam, pos  ) );
00686     
00687     det = 0; beam = 0; disk = 0;
00688     do {
00689       measuredCoordinates.GetTEC2TECEntry( det, beam, disk ).SetPhiError( 1000. );
00690     } while( moduleLoop.TEC2TECLoop( det, beam, disk ) );
00691   }
00692 
00693 
00694   // buffers for read-in
00695   int det, beam, z, ring;
00696   double phi, phiError;
00697 
00698   while( !file.eof() ) {
00699 
00700     file >> det;
00701     if( file.eof() ) break; // do not read the last line twice, do not fill trash if file empty
00702 
00703     file >> beam; file >> z; file >> ring;
00704     file >> phi; file >> phiError;
00705 
00706     if( det > 1 ) { // TIB/TOB
00707       measuredCoordinates.GetTIBTOBEntry( det, beam, z ).SetPhi( phi );
00708       measuredCoordinates.GetTIBTOBEntry( det, beam, z ).SetPhiError( phiError );
00709     } else { // TEC or TEC(at)
00710       if( ring > -1 ) { // TEC
00711         measuredCoordinates.GetTECEntry( det, ring, beam, z ).SetPhi( phi );
00712         measuredCoordinates.GetTECEntry( det, ring, beam, z ).SetPhiError( phiError );
00713       }
00714       else { // TEC(at)
00715         measuredCoordinates.GetTEC2TECEntry( det, beam, z ).SetPhi( phi );
00716         measuredCoordinates.GetTEC2TECEntry( det, beam, z ).SetPhiError( phiError );
00717       }
00718     }
00719 
00720   }
00721   
00722   file.close();
00723 
00724 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/Alignment/LaserAlignment/src/LASAlignmentTubeAlgorithm.cc
