Reference Manual

Go to the documentation of this file.
00001 
00002 #include "Alignment/LaserAlignment/interface/LASEndcapAlgorithm.h"
00003 
00004 
00005 
00009 LASEndcapAlgorithm::LASEndcapAlgorithm() {
00010 }
00011 
00012 
00013 
00014 
00015 
00022 LASEndcapAlignmentParameterSet LASEndcapAlgorithm::CalculateParameters( LASGlobalData<LASCoordinateSet>& measuredCoordinates, 
00023                                                                         LASGlobalData<LASCoordinateSet>& nominalCoordinates ) {
00024   
00025   std::cout << " [LASEndcapAlgorithm::CalculateParameters] -- Starting." << std::endl;
00026 
00027   // loop object
00028   LASGlobalLoop globalLoop;
00029   int det, ring, beam, disk;
00030 
00031   // vector containing the z positions of the disks in mm;
00032   // outer vector: TEC+/-, inner vector: 9 disks
00033   const double zPositions[9] = { 1322.5,  1462.5,  1602.5,  1742.5,  1882.5,  2057.5,  2247.5,  2452.5,  2667.5 };
00034   std::vector<std::vector<double> > diskZPositions( 2, std::vector<double>( 9, 0. ) );
00035   for( det = 0; det < 2; ++det ) {
00036     for( disk = 0; disk < 9; ++disk ) {
00037       // sign depends on side of course
00038       diskZPositions.at( det ).at( disk ) = ( det==0 ? zPositions[disk] : -1. * zPositions[disk] ); 
00039     }
00040   }
00041 
00042   // vector containing the phi positions of the beam in rad;
00043   // outer vector: TEC+/-, inner vector: 8 beams
00044   const double phiPositions[8] = { 0.392699, 1.178097, 1.963495, 2.748894, 3.534292, 4.319690, 5.105088, 5.890486 };
00045   std::vector<std::vector<double> > beamPhiPositions( 2, std::vector<double>( 8, 0. ) );
00046   for( det = 0; det < 2; ++ det ) {
00047     for( beam = 0; beam < 8; ++beam ) {
00048       beamPhiPositions.at( det ).at( beam ) = phiPositions[beam];
00049     }
00050   }
00051   
00052   // vector containing the radius of ring4,ring6 = (0,1)
00053   std::vector<double> radius( 2, 0. );
00054   radius.at( 0 ) = 564.; radius.at( 1 ) = 840.;
00055 
00056   // constants
00057   const double endcapLength = 1345.; // mm
00058 
00059   // now come some sums which are later used in the formulas for the parameters.
00060   // the rings are implicitly summed over, however, this brings a little complication:
00061   // to make the calculation of the parameters independent of the ring (=radius),
00062   // we define some of the sums twice, once for phi and once for y=r*phi
00063 
00064   // sum over r*phi or phi for each endcap and for each disk (both rings)
00065   // outer vector: TEC+/-, inner vector: 9 disks
00066   std::vector<std::vector<double> > sumOverY( 2, std::vector<double>( 9, 0. ) );
00067   std::vector<std::vector<double> > sumOverPhi( 2, std::vector<double>( 9, 0. ) );
00068 
00069   // sum over phi for each endcap and for each beam (both rings)
00070   // outer vector: TEC+/-, inner vector: 8 beams
00071   std::vector<std::vector<double> > kSumOverPhi( 2, std::vector<double>( 8, 0. ) );
00072 
00073   // double sum over r*phi or phi, for each endcap (both rings)
00074   // outer vector: TEC+/-
00075   std::vector<double> doubleSumOverY( 2, 0. );
00076   std::vector<double> doubleSumOverPhi( 2, 0. );
00077 
00078   // sum over r*phi*z or phi*z, for each endcap and for each beam (both rings)
00079   // outer vector: TEC+/-, inner vector: 8 beams
00080   std::vector<std::vector<double> > kSumOverPhiZ( 2, std::vector<double>( 8, 0. ) );
00081 
00082   // sum over r*phi*z or phi*z, for each endcap (both rings)
00083   // outer vector: TEC+/-
00084   std::vector<double> doubleSumOverYZ( 2, 0. );
00085   std::vector<double> doubleSumOverPhiZ( 2, 0. );
00086 
00087   // sum over sin(phi_nominal)*R*phi for each endcap and for each disk (both rings)
00088   std::vector<std::vector<double> > sumOverSinThetaY( 2, std::vector<double>( 9, 0. ) );
00089 
00090   // sum over cos(phi_nominal)*R*phi for each endcap and for each disk (both rings)
00091   std::vector<std::vector<double> > sumOverCosThetaY( 2, std::vector<double>( 9, 0. ) );
00092 
00093   // double sum over sin or cos(phi_nominal)*phi, for each endcap
00094   std::vector<double> doubleSumOverSinThetaY( 2, 0. );
00095   std::vector<double> doubleSumOverCosThetaY( 2, 0. );
00096 
00097   // double sum over sin or cos(phi_nominal)*phi*z, for each endcap
00098   std::vector<double> doubleSumOverSinThetaYZ( 2, 0. );
00099   std::vector<double> doubleSumOverCosThetaYZ( 2, 0. );
00100 
00101   // sum over z values / sum over z^2 values
00102   std::vector<double> sumOverZ( 2, 0. );
00103   std::vector<double> sumOverZSquared( 2, 0. );
00104 
00105 
00106   // now calculate the sums
00107   det = 0; ring = 0; beam = 0; disk = 0;
00108   do {
00109     if( ring == 1 ) continue; //################################################################################################### BOUND TO RING6
00110     // current radius, depends on the ring
00111     const double radius = nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetR();
00112 
00113     // residual in r*phi (in the formulas this corresponds to y_ik)
00114     const double residual = ( measuredCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() - nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * radius;
00115     
00116     sumOverY.at( det ).at( disk ) += residual;
00117     sumOverPhi.at( det ).at( disk ) += residual / radius;
00118     kSumOverPhi.at( det ).at( beam ) += residual / radius;
00119 
00120     doubleSumOverY.at( det ) += residual;
00121     doubleSumOverPhi.at( det ) += residual / radius;
00122 
00123     kSumOverPhiZ.at( det ).at( beam ) += diskZPositions.at( det ).at( disk ) * residual / radius;
00124 
00125     doubleSumOverYZ.at( det ) += diskZPositions.at( det ).at( disk ) * residual;
00126     doubleSumOverPhiZ.at( det ) +=  diskZPositions.at( det ).at( disk ) * residual / radius;
00127 
00128     sumOverSinThetaY.at( det ).at( disk ) += sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * residual;
00129     sumOverCosThetaY.at( det ).at( disk ) += cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * residual;
00130 
00131     doubleSumOverSinThetaY.at( det ) += sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * residual;
00132     doubleSumOverCosThetaY.at( det ) += cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * residual;
00133 
00134     doubleSumOverSinThetaYZ.at( det ) += sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * diskZPositions.at( det ).at( disk ) * residual;
00135     doubleSumOverCosThetaYZ.at( det ) += cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * diskZPositions.at( det ).at( disk ) * residual;
00136 
00137   } while( globalLoop.TECLoop( det, ring, beam, disk ) );
00138 
00139 
00140   // disk-wise sums
00141   for( disk = 0; disk < 9; ++disk ) {
00142     sumOverZ.at( 0 ) += diskZPositions.at( 0 ).at( disk ); sumOverZ.at( 1 ) += diskZPositions.at( 1 ).at( disk );
00143     sumOverZSquared.at( 0 ) += pow( diskZPositions.at( 0 ).at( disk ), 2 ); sumOverZSquared.at( 1 ) += pow( diskZPositions.at( 1 ).at( disk ), 2 );
00144   }
00145   
00146 
00147   // now we can calculate the parameters for both TECs simultaneously,
00148   // so they're all vectors( 2 ) for TEC+/- (global parameters), or dim 2*9 (disk parameters),
00149   // or dim 2*8 (beam parameters)
00150 
00151   // define them..
00152 
00153   // deltaPhi_0
00154   std::vector<double> deltaPhi0( 2, 0. );
00155   
00156   // deltaPhi_t
00157   std::vector<double> deltaPhiT( 2, 0. );
00158   
00159   // deltaPhi_k (k=0..8)
00160   std::vector<std::vector<double> > deltaPhiK( 2, std::vector<double>( 9, 0. ) );
00161 
00162   // deltaX_0
00163   std::vector<double> deltaX0( 2, 0. );
00164   
00165   // deltaX_t
00166   std::vector<double> deltaXT( 2, 0. );
00167   
00168   // deltaX_k (k=0..8)
00169   std::vector<std::vector<double> > deltaXK( 2, std::vector<double>( 9, 0. ) );
00170 
00171   // deltaY_0
00172   std::vector<double> deltaY0( 2, 0. );
00173   
00174   // deltaY_t
00175   std::vector<double> deltaYT( 2, 0. );
00176   
00177   // deltaY_k (k=0..8)
00178   std::vector<std::vector<double> > deltaYK( 2, std::vector<double>( 9, 0. ) );
00179 
00180   // beam parameters: deltaTheta_A, deltaTheta_B (i=0..7)
00181   std::vector<std::vector<double> > deltaTA( 2, std::vector<double>( 8, 0. ) );
00182   std::vector<std::vector<double> > deltaTB( 2, std::vector<double>( 8, 0. ) );
00183   
00184 
00185 
00186   // ..and fill them;
00187   // the additional divisors "/ 2." come from the fact that we average over both rings
00188   for( det = 0; det < 2; ++det ) { // TEC+/- loop
00189     
00190     // deltaPhi_0
00191     // here we use the phi-sums to eliminate the radius
00192     deltaPhi0.at( det ) = ( sumOverZSquared.at( det ) * doubleSumOverPhi.at( det ) - sumOverZ.at( det ) * doubleSumOverPhiZ.at( det ) )
00193       / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) ); // / 2.; // @@@@@@@
00194 
00195     // deltaPhi_t
00196     // again use the phi-sums
00197     deltaPhiT.at( det ) = endcapLength * ( 9. * doubleSumOverPhiZ.at( det ) - sumOverZ.at( det ) * doubleSumOverPhi.at( det ) )
00198       / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) ); // / 2.; // @@@@@@@
00199 
00200     // deltaPhi_k (k=0..8)
00201     // again use the phi-sums
00202     for( disk = 0; disk < 9; ++disk ) {
00203       deltaPhiK.at( det ).at( disk ) = ( -1. * diskZPositions.at( det ).at( disk ) * deltaPhiT.at( det ) / endcapLength )
00204         -  ( deltaPhi0.at( det ) )  -  sumOverPhi.at( det ).at( disk ) / 8.;// / 2.; // @@@@@@@
00205     }
00206     
00207     // deltaX_0
00208     deltaX0.at( det ) = 2. * ( sumOverZ.at( det ) * doubleSumOverSinThetaYZ.at( det ) - sumOverZSquared.at( det ) * doubleSumOverSinThetaY.at( det ) )
00209       / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) );// / 2.; // @@@@@@@
00210 
00211     // deltaX_t
00212     deltaXT.at( det ) = 2. * endcapLength * ( sumOverZ.at( det ) * doubleSumOverSinThetaY.at( det ) - 9. * doubleSumOverSinThetaYZ.at( det ) )
00213       / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) );// / 2.; // @@@@@@@
00214 
00215     // deltaX_k (k=0..8)
00216     for( disk = 0; disk < 9; ++disk ) {
00217       deltaXK.at( det ).at( disk ) =  ( -1. * diskZPositions.at( det ).at( disk ) * deltaXT.at( det ) / endcapLength )
00218         -  ( deltaX0.at( det ) )  +  2. * sumOverSinThetaY.at( det ).at( disk ) / 8.;// / 2.; // @@@@@@@
00219     }
00220 
00221     // deltaY_0
00222     deltaY0.at( det ) = 2. * ( sumOverZSquared.at( det ) * doubleSumOverCosThetaY.at( det ) - sumOverZ.at( det ) * doubleSumOverCosThetaYZ.at( det ) )
00223       / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) );// / 2.; // @@@@@@@
00224 
00225     // deltaY_t
00226     deltaYT.at( det ) = 2. * endcapLength * ( 9. * doubleSumOverCosThetaYZ.at( det ) - sumOverZ.at( det ) * doubleSumOverCosThetaY.at( det ) )
00227       / ( 8. * ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) ) );// / 2.; // @@@@@@@
00228 
00229     // deltaY_k (k=0..8)
00230     for( disk = 0; disk < 9; ++disk ) {
00231       deltaYK.at( det ).at( disk ) =  ( -1. * diskZPositions.at( det ).at( disk ) * deltaYT.at( det ) / endcapLength )
00232         -  ( deltaY0.at( det ) )  -  2. * sumOverCosThetaY.at( det ).at( disk ) / 8.;// / 2.; // @@@@@@@
00233     }
00234 
00235 
00236     // the beam parameters deltaTA & deltaTB
00237     for( beam = 0; beam < 8; ++beam ) {
00238 
00239       deltaTA.at( det ).at( beam ) = deltaPhi0.at( det ) 
00240         - ( kSumOverPhi.at( det ).at( beam ) * sumOverZSquared.at( det ) - kSumOverPhiZ.at( det ).at( beam ) * sumOverZ.at( det ) )
00241         / ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) )
00242         + ( cos( beamPhiPositions.at( det ).at( beam ) ) * deltaY0.at( det ) - sin( beamPhiPositions.at( det ).at( beam ) ) * deltaX0.at( det ) ) / radius.at( 0 ); // for ring 4..
00243         // + ( cos( beamPhiPositions.at( det ).at( beam ) ) * deltaY0.at( det ) - sin( beamPhiPositions.at( det ).at( beam ) ) * deltaX0.at( det ) ) / radius.at( 1 ); // ...and ring 6
00244 
00245 
00246       deltaTB.at( det ).at( beam ) = -1. * deltaPhiT.at( det ) - deltaPhi0.at( det )
00247         - ( kSumOverPhi.at( det ).at( beam ) * sumOverZ.at( det ) - 9. * kSumOverPhiZ.at( det ).at( beam ) )
00248         * endcapLength / ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) )
00249         - ( kSumOverPhiZ.at( det ).at( beam ) * sumOverZ.at( det ) -  kSumOverPhi.at( det ).at( beam ) * sumOverZSquared.at( det ) )
00250         / ( pow( sumOverZ.at( det ), 2 ) - 9. * sumOverZSquared.at( det ) )
00251         + ( ( deltaXT.at( det ) + deltaX0.at( det ) ) * sin( beamPhiPositions.at( det ).at( beam ) ) - ( deltaYT.at( det ) + deltaY0.at( det ) ) * cos( beamPhiPositions.at( det ).at( beam ) ) )
00252         / radius.at( 0 ); // for ring4..
00253       // + ( ( deltaXT.at( det ) + deltaX0.at( det ) ) * sin( beamPhiPositions.at( det ).at( beam ) ) - ( deltaYT.at( det ) + deltaY0.at( det ) ) * cos( beamPhiPositions.at( det ).at( beam ) ) )
00254       // / radius.at( 1 ); // ..and ring6
00255 
00256     }
00257 
00258 
00259   }
00260 
00261 
00262   // fill the result
00263   LASEndcapAlignmentParameterSet theResult;
00264 
00265   // for the moment we fill only the values, not the errors
00266 
00267   // disk parameters
00268   for( det = 0; det < 2; ++det ) {
00269     for( disk = 0; disk < 9; ++disk ) {
00270       // the rotation parameters: deltaPhi_k
00271       theResult.GetDiskParameter( det, disk, 0 ).first = deltaPhiK.at( det ).at( disk );
00272       // the x offsets: deltaX_k
00273       theResult.GetDiskParameter( det, disk, 1 ).first = deltaXK.at( det ).at( disk );
00274       // the y offsets: deltaY_k
00275       theResult.GetDiskParameter( det, disk, 2 ).first = deltaYK.at( det ).at( disk );
00276     }
00277   }
00278 
00279   // global parameters
00280   for( int det = 0; det < 2; ++det ) {
00281     theResult.GetGlobalParameter( det, 0 ).first = deltaPhi0.at( det );
00282     theResult.GetGlobalParameter( det, 1 ).first = deltaPhiT.at( det );
00283     theResult.GetGlobalParameter( det, 2 ).first = deltaX0.at( det );
00284     theResult.GetGlobalParameter( det, 3 ).first = deltaXT.at( det );
00285     theResult.GetGlobalParameter( det, 4 ).first = deltaY0.at( det );
00286     theResult.GetGlobalParameter( det, 5 ).first = deltaYT.at( det );
00287   }
00288 
00289   // beam parameters
00290   for( int det = 0; det < 2; ++det ) {
00291     for( int beam = 0; beam < 8; ++beam ) {
00292       theResult.GetBeamParameter( det, 1 /*R6*/, beam, 0 ).first = deltaTA.at( det ).at( beam ); 
00293       theResult.GetBeamParameter( det, 1 /*R6*/, beam, 1 ).first = deltaTB.at( det ).at( beam );
00294     }
00295   }
00296 
00297   std::cout << " [LASEndcapAlgorithm::CalculateParameters] -- Done." << std::endl;
00298 
00299   return( theResult );
00300 
00301 }
00302 
00303 
00304 
00305 
00306 
00312 double LASEndcapAlgorithm::GetAlignmentParameterCorrection( int det, int ring, int beam, int disk, LASGlobalData<LASCoordinateSet>& nominalCoordinates, LASEndcapAlignmentParameterSet& endcapParameters ) {
00313   
00314   // ring dependent radius, to be softcoded...
00315   const double radius = ring==0 ? 564. : 840.;
00316   const double endcapLength = 1345.; // mm
00317   
00318   // the correction to phi from the endcap algorithm;
00319   // it is defined such that the correction is to be subtracted
00320   double phiCorrection = 0.;
00321   
00322   // plain disk phi
00323   phiCorrection += endcapParameters.GetDiskParameter( det, disk, 0 ).first;
00324   
00325   // phi component from x deviation
00326   phiCorrection -= sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) / radius * endcapParameters.GetDiskParameter( det, disk, 1 ).first;
00327   
00328   // phi component from y deviation
00329   phiCorrection += cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) / radius * endcapParameters.GetDiskParameter( det, disk, 2 ).first;
00330   
00331   // phi correction from global phi
00332   phiCorrection += endcapParameters.GetGlobalParameter( det, 0 ).first;
00333   
00334   // correction from global x deviation
00335   phiCorrection -= sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) / radius * endcapParameters.GetGlobalParameter( det, 2 ).first;
00336   
00337   // correction from global y deviation
00338   phiCorrection += cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) / radius * endcapParameters.GetGlobalParameter( det, 4 ).first;
00339   
00340   // correction from global torsion
00341   phiCorrection += nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetZ() / endcapLength * endcapParameters.GetGlobalParameter( det, 1 ).first;
00342   
00343   // correction from global x shear
00344   phiCorrection -= nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetZ() / endcapLength / radius *
00345     sin( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * endcapParameters.GetGlobalParameter( det, 3 ).first;
00346   
00347   // correction from global y shear
00348   phiCorrection += nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetZ() / endcapLength / radius *
00349     cos( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetPhi() ) * endcapParameters.GetGlobalParameter( det, 5 ).first;
00350   
00351   // correction from beam parameters
00352   phiCorrection += ( nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetZ() / endcapLength - 1. ) * endcapParameters.GetBeamParameter( det, 1, beam, 0 ).first;
00353   phiCorrection += nominalCoordinates.GetTECEntry( det, ring, beam, disk ).GetZ() / endcapLength * endcapParameters.GetBeamParameter( det, 1, beam, 1 ).first;
00354   
00355   return phiCorrection;
00356 
00357 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/Alignment/LaserAlignment/src/LASEndcapAlgorithm.cc
