|
|
|
#include "Alignment/LaserAlignment/src/LASBarrelAlgorithm.h"Go to the source code of this file.
Functions | |
| void | fcn (int &npar, double *gin, double &f, double *par, int iflag) |
| minuit chisquare func | |
Variables | |
| LASGlobalData< LASCoordinateSet > * | aMeasuredCoordinates |
| LASGlobalData< LASCoordinateSet > * | aNominalCoordinates |
minuit chisquare func
Definition at line 334 of file LASBarrelAlgorithm.cc.
References funct::cos(), muonGeometry::disk, LASCoordinateSet::GetPhi(), LASCoordinateSet::GetPhiError(), LASCoordinateSet::GetR(), LASGlobalData< T >::GetTEC2TECEntry(), LASGlobalData< T >::GetTIBTOBEntry(), LASCoordinateSet::GetZ(), funct::pow(), funct::sin(), LASGlobalLoop::TEC2TECLoop(), and LASGlobalLoop::TIBTOBLoop().
Referenced by LASBarrelAlgorithm::CalculateParameters().
00334 { 00335 00336 double chisquare = 0.; 00337 00338 // the loop object and its variables 00339 LASGlobalLoop moduleLoop; 00340 int det, beam, pos, disk; 00341 00343 // ADJUST THIS ALSO IN LASGeometryUpdater 00345 00346 // the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm); 00347 // parameters are: det, side(0=+/1=-), z(lower/upper). TECs have no side (use side = 0) 00348 std::vector<std::vector<std::vector<double> > > endFaceZPositions( 4, std::vector<std::vector<double> >( 2, std::vector<double>( 2, 0. ) ) ); 00349 endFaceZPositions.at( 0 ).at( 0 ).at( 0 ) = 1250.; // TEC+, *, disk1 /// 00350 endFaceZPositions.at( 0 ).at( 0 ).at( 1 ) = 2595.; // TEC+, *, disk9 /// SIDE INFORMATION 00351 endFaceZPositions.at( 1 ).at( 0 ).at( 0 ) = -2595.; // TEC-, *, disk1 /// MEANINGLESS FOR TEC -> USE .at(0)! 00352 endFaceZPositions.at( 1 ).at( 0 ).at( 1 ) = -1250.; // TEC-, *, disk9 /// 00353 endFaceZPositions.at( 2 ).at( 1 ).at( 0 ) = -700.; // TIB, -, small z 00354 endFaceZPositions.at( 2 ).at( 1 ).at( 1 ) = -300.; // TIB, -, large z 00355 endFaceZPositions.at( 2 ).at( 0 ).at( 0 ) = 300.; // TIB, +, small z 00356 endFaceZPositions.at( 2 ).at( 0 ).at( 1 ) = 700.; // TIB, +, large z 00357 endFaceZPositions.at( 3 ).at( 1 ).at( 0 ) = -1090.; // TOB, -, small z 00358 endFaceZPositions.at( 3 ).at( 1 ).at( 1 ) = -300.; // TOB, -, large z 00359 endFaceZPositions.at( 3 ).at( 0 ).at( 0 ) = 300.; // TOB, +, small z 00360 endFaceZPositions.at( 3 ).at( 0 ).at( 1 ) = 1090.; // TOB, +, large z 00361 00362 // the z positions of the TEC outer disks (9) in mm 00363 // (in priciple one could also use the above vector set here, but it's more compact) 00364 std::vector<double> disk9EndFaceZPositions( 2, 0. ); 00365 disk9EndFaceZPositions.at( 0 ) = -2595.; // TEC- disk9 00366 disk9EndFaceZPositions.at( 1 ) = 2595.; // TEC+ disk9 00367 00368 // reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} ) 00369 double detReducedZ[2] = { 0., 0. }; 00370 // reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} ) 00371 double beamReducedZ[2] = { 0., 0. }; 00372 00373 // calculate residual for TIBTOB 00374 det = 2; beam = 0; pos = 0; 00375 do { 00376 00377 // define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB- 00378 const int theSide = pos<3 ? 0 : 1; 00379 00380 // reduced module's z position with respect to the subdetector endfaces 00381 detReducedZ[0] = aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ() - endFaceZPositions.at( det ).at( theSide ).at( 0 ); 00382 detReducedZ[0] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) ); 00383 detReducedZ[1] = endFaceZPositions.at( det ).at( theSide ).at( 1 ) - aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ(); 00384 detReducedZ[1] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) ); 00385 00386 // reduced module's z position with respect to the tec disks +-9 (for the beam parameters) 00387 beamReducedZ[0] = aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ() - disk9EndFaceZPositions.at( 0 ); 00388 beamReducedZ[0] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) ); 00389 beamReducedZ[1] = disk9EndFaceZPositions.at( 1 ) - aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ(); 00390 beamReducedZ[1] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) ); 00391 00392 // phi residual for this module as measured 00393 const double measuredResidual = aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhi() - //& 00394 aNominalCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhi(); 00395 00396 // shortcuts for speed 00397 const double currentPhi = aNominalCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhi(); 00398 const double currentR = aNominalCoordinates->GetTIBTOBEntry( det, beam, pos ).GetR(); 00399 00400 // phi residual for this module calculated from the parameters and nominal coordinates: 00401 // this is the sum over the contributions from all parameters 00402 double calculatedResidual = 0.; 00403 00404 // note that the contributions ym_{i,j,k} given in the tables in TkLasATModel TWiki 00405 // are defined as R*phi, so here they are divided by the R_j factors (we minimize delta phi) 00406 00407 // unfortunately, minuit keeps parameters in a 1-dim array, 00408 // so we need to address the correct par[] for the 4 cases TIB+/TIB-/TOB+/TOB- 00409 int indexBase = 0; 00410 if( det == 2 ) { // TIB 00411 if( theSide == 0 ) indexBase = 0; // TIB+ 00412 if( theSide == 1 ) indexBase = 6; // TIB- 00413 } 00414 if( det == 3 ) { // TOB 00415 if( theSide == 0 ) indexBase = 12; // TOB+ 00416 if( theSide == 1 ) indexBase = 18; // TOB- 00417 } 00418 00419 // par[0] ("subRot1"): rotation around z of first end face 00420 calculatedResidual += detReducedZ[1] * par[indexBase+0]; 00421 00422 // par[1] ("subRot2"): rotation around z of second end face 00423 calculatedResidual += detReducedZ[0] * par[indexBase+1]; 00424 00425 // par[2] ("subTransX1"): translation along x of first end face 00426 calculatedResidual += detReducedZ[1] * sin( currentPhi ) / currentR * par[indexBase+2]; 00427 00428 // par[3] ("subTransX2"): translation along x of second end face 00429 calculatedResidual += detReducedZ[0] * sin( currentPhi ) / currentR * par[indexBase+3]; 00430 00431 // par[4] ("subTransY1"): translation along y of first end face 00432 calculatedResidual += -1. * detReducedZ[1] * cos( currentPhi ) / currentR * par[indexBase+4]; 00433 00434 // par[5] ("subTransY2"): translation along y of second end face 00435 calculatedResidual += -1. * detReducedZ[0] * cos( currentPhi ) / currentR * par[indexBase+5]; 00436 00437 00438 // now come the 8*2 beam parameters, calculate the respective parameter index base first (-> which beam) 00439 indexBase = 36 + beam * 2; 00440 00441 // (there's no TIB/TOB/+/- distinction here for the beams) 00442 00443 // ("beamRot1"): rotation around z at zt1 00444 calculatedResidual += beamReducedZ[1] * par[indexBase]; 00445 00446 // ("beamRot2"): rotation around z at zt2 00447 calculatedResidual += beamReducedZ[0] * par[indexBase+1]; 00448 00449 00450 // now calculate the chisquare 00451 chisquare += pow( measuredResidual - calculatedResidual, 2 ) / pow( aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhiError(), 2 ); 00452 00453 } while( moduleLoop.TIBTOBLoop( det, beam, pos ) ); 00454 00455 00456 00457 00458 00459 // calculate residual for TEC AT 00460 det = 0; beam = 0; disk = 0; 00461 do { 00462 00463 // define the side: TECs sides already disentangled by the "det" index, so fix this to zero 00464 const int theSide = 0; 00465 00466 // reduced module's z position with respect to the subdetector endfaces 00467 detReducedZ[0] = aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ() - endFaceZPositions.at( det ).at( theSide ).at( 0 ); 00468 detReducedZ[0] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) ); 00469 detReducedZ[1] = endFaceZPositions.at( det ).at( theSide ).at( 1 ) - aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ(); 00470 detReducedZ[1] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) ); 00471 00472 // reduced module's z position with respect to the tec disks +-9 (for the beam parameters) 00473 beamReducedZ[0] = aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ() - disk9EndFaceZPositions.at( 0 ); 00474 beamReducedZ[0] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) ); 00475 beamReducedZ[1] = disk9EndFaceZPositions.at( 1 ) - aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ(); 00476 beamReducedZ[1] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) ); 00477 00478 // phi residual for this module as measured 00479 const double measuredResidual = aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhi() - //& 00480 aNominalCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhi(); 00481 00482 // shortcuts for speed 00483 const double currentPhi = aNominalCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhi(); 00484 const double currentR = aNominalCoordinates->GetTEC2TECEntry( det, beam, disk ).GetR(); 00485 00486 // phi residual for this module calculated from the parameters and nominal coordinates: 00487 // this is the sum over the contributions from all parameters 00488 double calculatedResidual = 0.; 00489 00490 // note that the contributions ym_{i,j,k} given in the tables in TkLasATModel TWiki 00491 // are defined as R*phi, so here they are divided by the R_j factors (we minimize delta phi) 00492 00493 // there's also a distinction between TEC+/- parameters in situ (det==0 ? <TEC+> : <TEC->) 00494 00495 // par[0] ("subRot1"): rotation around z of first end face 00496 calculatedResidual += detReducedZ[1] * (det==0 ? par[24] : par[30]); 00497 00498 // par[1] ("subRot2"): rotation around z of second end face 00499 calculatedResidual += detReducedZ[0] * (det==0 ? par[25] : par[31]); 00500 00501 // par[2] ("subTransX1"): translation along x of first end face 00502 calculatedResidual += detReducedZ[1] * sin( currentPhi ) * (det==0 ? par[26] : par[32]) / currentR; 00503 00504 // par[3] ("subTransX2"): translation along x of second end face 00505 calculatedResidual += detReducedZ[0] * sin( currentPhi ) * (det==0 ? par[27] : par[33]) / currentR; 00506 00507 // par[4] ("subTransY1"): translation along y of first end face 00508 calculatedResidual += -1. * detReducedZ[1] * cos( currentPhi ) * (det==0 ? par[28] : par[34]) / currentR; 00509 00510 // par[5] ("subTransY2"): translation along y of second end face 00511 calculatedResidual += -1. * detReducedZ[0] * cos( currentPhi ) * (det==0 ? par[29] : par[35]) / currentR; 00512 00513 // now come the 8*2 beam parameters; calculate the respective parameter index base first (-> which beam) 00514 const unsigned int indexBase = 36 + beam * 2; 00515 00516 // there's no TEC+/- distinction here 00517 00518 // par[6] ("beamRot1"): rotation around z at zt1 00519 calculatedResidual += beamReducedZ[1] * par[indexBase]; 00520 00521 // par[7] ("beamRot2"): rotation around z at zt2 00522 calculatedResidual += beamReducedZ[0] * par[indexBase+1]; 00523 00524 00525 // now calculate the chisquare 00526 // TODO: check for phi != 0 !!! 00527 chisquare += pow( measuredResidual - calculatedResidual, 2 ) / pow( aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhiError(), 2 ); 00528 00529 } while( moduleLoop.TEC2TECLoop( det, beam, disk ) ); 00530 00531 00532 00533 // return the chisquare by ref 00534 f = chisquare; 00535 00536 }
Definition at line 8 of file LASBarrelAlgorithm.cc.
Definition at line 9 of file LASBarrelAlgorithm.cc.
1.5.4