|
|
|
#include <vector>#include <cmath>#include <string>#include <iostream>#include <iomanip>#include <fstream>#include <TMinuit.h>#include "Alignment/LaserAlignment/interface/LASBarrelAlignmentParameterSet.h"#include "Alignment/LaserAlignment/interface/LASCoordinateSet.h"#include "Alignment/LaserAlignment/interface/LASGlobalData.h"#include "Alignment/LaserAlignment/interface/LASGlobalLoop.h"Go to the source code of this file.
Classes | |
| class | LASBarrelAlgorithm |
Functions | |
| void | fcn (int &, double *, double &, double *, int) |
| void fcn | ( | int & | npar, |
| double * | gin, | ||
| double & | f, | ||
| double * | par, | ||
| int | iflag | ||
| ) |
minuit chisquare func
Definition at line 343 of file LASBarrelAlgorithm.cc.
References funct::cos(), LASCoordinateSet::GetPhi(), LASCoordinateSet::GetPhiError(), LASCoordinateSet::GetR(), LASGlobalData< T >::GetTEC2TECEntry(), LASGlobalData< T >::GetTIBTOBEntry(), LASCoordinateSet::GetZ(), pos, funct::pow(), funct::sin(), LASGlobalLoop::TEC2TECLoop(), and LASGlobalLoop::TIBTOBLoop().
Referenced by LASBarrelAlgorithm::CalculateParameters(), MuonResidualsFitter::dofit(), fftjetcms::fftjet_JetFunctor_parser(), fftjetcms::fftjet_PeakFunctor_parser(), npstat::ArrayND< Numeric, StackLen, StackDim >::jointSliceLoop(), npstat::ArrayND< Numeric, StackLen, StackDim >::projectLoop(), npstat::ArrayND< Numeric, StackLen, StackDim >::projectLoop2(), PVFitter::runBXFitter(), and PVFitter::runFitter().
{
double chisquare = 0.;
// the loop object and its variables
LASGlobalLoop moduleLoop;
int det, beam, pos, disk;
// ADJUST THIS ALSO IN LASGeometryUpdater
// the z positions of the halfbarrel_end_faces / outer_TEC_disks (in mm);
// parameters are: det, side(0=+/1=-), z(lower/upper). TECs have no side (use side = 0)
std::vector<std::vector<std::vector<double> > > endFaceZPositions( 4, std::vector<std::vector<double> >( 2, std::vector<double>( 2, 0. ) ) );
endFaceZPositions.at( 0 ).at( 0 ).at( 0 ) = 1322.5; // TEC+, *, disk1 ///
endFaceZPositions.at( 0 ).at( 0 ).at( 1 ) = 2667.5; // TEC+, *, disk9 /// SIDE INFORMATION
endFaceZPositions.at( 1 ).at( 0 ).at( 0 ) = -2667; // TEC-, *, disk1 /// MEANINGLESS FOR TEC -> USE .at(0)!
endFaceZPositions.at( 1 ).at( 0 ).at( 1 ) = -1322.5; // TEC-, *, disk9 ///
endFaceZPositions.at( 2 ).at( 1 ).at( 0 ) = -700.; // TIB, -, small z
endFaceZPositions.at( 2 ).at( 1 ).at( 1 ) = -300.; // TIB, -, large z
endFaceZPositions.at( 2 ).at( 0 ).at( 0 ) = 300.; // TIB, +, small z
endFaceZPositions.at( 2 ).at( 0 ).at( 1 ) = 700.; // TIB, +, large z
endFaceZPositions.at( 3 ).at( 1 ).at( 0 ) = -1090.; // TOB, -, small z
endFaceZPositions.at( 3 ).at( 1 ).at( 1 ) = -300.; // TOB, -, large z
endFaceZPositions.at( 3 ).at( 0 ).at( 0 ) = 300.; // TOB, +, small z
endFaceZPositions.at( 3 ).at( 0 ).at( 1 ) = 1090.; // TOB, +, large z
// the z positions of the virtual planes at which the beam parameters are measured
std::vector<double> disk9EndFaceZPositions( 2, 0. );
disk9EndFaceZPositions.at( 0 ) = -2667.5; // TEC- disk9
disk9EndFaceZPositions.at( 1 ) = 2667.5; // TEC+ disk9
// reduced z positions of the beam spots ( z'_{k,j}, z"_{k,j} )
double detReducedZ[2] = { 0., 0. };
// reduced beam splitter positions ( zt'_{k,j}, zt"_{k,j} )
double beamReducedZ[2] = { 0., 0. };
// calculate residual for TIBTOB
det = 2; beam = 0; pos = 0;
do {
// define the side: 0 for TIB+/TOB+ and 1 for TIB-/TOB-
const int theSide = pos<3 ? 0 : 1;
// this is the path the beam has to travel radially after being reflected
// by the AT mirrors (TIB:50mm, TOB:36mm) -> used for beam parameters
const double radialOffset = det==2 ? 50. : 36.;
// reduced module's z position with respect to the subdetector endfaces
detReducedZ[0] = aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ() - endFaceZPositions.at( det ).at( theSide ).at( 0 );
detReducedZ[0] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
detReducedZ[1] = endFaceZPositions.at( det ).at( theSide ).at( 1 ) - aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ();
detReducedZ[1] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
// reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
beamReducedZ[0] = ( aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ() - radialOffset ) - disk9EndFaceZPositions.at( 0 );
beamReducedZ[0] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
beamReducedZ[1] = disk9EndFaceZPositions.at( 1 ) - ( aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetZ() - radialOffset );
beamReducedZ[1] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
// phi residual for this module as measured
const double measuredResidual = aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhi() - //&
aNominalCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhi();
// shortcuts for speed
const double currentPhi = aNominalCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhi();
const double currentR = aNominalCoordinates->GetTIBTOBEntry( det, beam, pos ).GetR();
// phi residual for this module calculated from the parameters and nominal coordinates:
// this is the sum over the contributions from all parameters
double calculatedResidual = 0.;
// note that the contributions ym_{i,j,k} given in the tables in TkLasATModel TWiki
// are defined as R*phi, so here they are divided by the R_j factors (we minimize delta phi)
// unfortunately, minuit keeps parameters in a 1-dim array,
// so we need to address the correct par[] for the 4 cases TIB+/TIB-/TOB+/TOB-
int indexBase = 0;
if( det == 2 ) { // TIB
if( theSide == 0 ) indexBase = 0; // TIB+
if( theSide == 1 ) indexBase = 6; // TIB-
}
if( det == 3 ) { // TOB
if( theSide == 0 ) indexBase = 12; // TOB+
if( theSide == 1 ) indexBase = 18; // TOB-
}
// par[0] ("subRot1"): rotation around z of first end face
calculatedResidual += detReducedZ[1] * par[indexBase+0];
// par[1] ("subRot2"): rotation around z of second end face
calculatedResidual += detReducedZ[0] * par[indexBase+1];
// par[2] ("subTransX1"): translation along x of first end face
calculatedResidual += detReducedZ[1] * sin( currentPhi ) / currentR * par[indexBase+2];
// par[3] ("subTransX2"): translation along x of second end face
calculatedResidual += detReducedZ[0] * sin( currentPhi ) / currentR * par[indexBase+3];
// par[4] ("subTransY1"): translation along y of first end face
calculatedResidual += -1. * detReducedZ[1] * cos( currentPhi ) / currentR * par[indexBase+4];
// par[5] ("subTransY2"): translation along y of second end face
calculatedResidual += -1. * detReducedZ[0] * cos( currentPhi ) / currentR * par[indexBase+5];
// now come the 8*2 beam parameters, calculate the respective parameter index base first (-> which beam)
indexBase = 36 + beam * 2;
// (there's no TIB/TOB/+/- distinction here for the beams)
// ("beamRot1"): rotation around z at zt1
calculatedResidual += beamReducedZ[1] * par[indexBase];
// ("beamRot2"): rotation around z at zt2
calculatedResidual += beamReducedZ[0] * par[indexBase+1];
// now calculate the chisquare
chisquare += pow( measuredResidual - calculatedResidual, 2 ) / pow( aMeasuredCoordinates->GetTIBTOBEntry( det, beam, pos ).GetPhiError(), 2 );
} while( moduleLoop.TIBTOBLoop( det, beam, pos ) );
// calculate residual for TEC AT
det = 0; beam = 0; disk = 0;
do {
// define the side: TECs sides already disentangled by the "det" index, so fix this to zero
const int theSide = 0;
// reduced module's z position with respect to the subdetector endfaces
detReducedZ[0] = aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ() - endFaceZPositions.at( det ).at( theSide ).at( 0 );
detReducedZ[0] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
detReducedZ[1] = endFaceZPositions.at( det ).at( theSide ).at( 1 ) - aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ();
detReducedZ[1] /= ( endFaceZPositions.at( det ).at( theSide ).at( 1 ) - endFaceZPositions.at( det ).at( theSide ).at( 0 ) );
// reduced module's z position with respect to the tec disks +-9 (for the beam parameters)
beamReducedZ[0] = aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ() - disk9EndFaceZPositions.at( 0 );
beamReducedZ[0] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
beamReducedZ[1] = disk9EndFaceZPositions.at( 1 ) - aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetZ();
beamReducedZ[1] /= ( disk9EndFaceZPositions.at( 1 ) - disk9EndFaceZPositions.at( 0 ) );
// phi residual for this module as measured
const double measuredResidual = aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhi() - //&
aNominalCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhi();
// shortcuts for speed
const double currentPhi = aNominalCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhi();
const double currentR = aNominalCoordinates->GetTEC2TECEntry( det, beam, disk ).GetR();
// phi residual for this module calculated from the parameters and nominal coordinates:
// this is the sum over the contributions from all parameters
double calculatedResidual = 0.;
// note that the contributions ym_{i,j,k} given in the tables in TkLasATModel TWiki
// are defined as R*phi, so here they are divided by the R_j factors (we minimize delta phi)
// there's also a distinction between TEC+/- parameters in situ (det==0 ? <TEC+> : <TEC->)
// par[0] ("subRot1"): rotation around z of first end face
calculatedResidual += detReducedZ[1] * (det==0 ? par[24] : par[30]);
// par[1] ("subRot2"): rotation around z of second end face
calculatedResidual += detReducedZ[0] * (det==0 ? par[25] : par[31]);
// par[2] ("subTransX1"): translation along x of first end face
calculatedResidual += detReducedZ[1] * sin( currentPhi ) * (det==0 ? par[26] : par[32]) / currentR;
// par[3] ("subTransX2"): translation along x of second end face
calculatedResidual += detReducedZ[0] * sin( currentPhi ) * (det==0 ? par[27] : par[33]) / currentR;
// par[4] ("subTransY1"): translation along y of first end face
calculatedResidual += -1. * detReducedZ[1] * cos( currentPhi ) * (det==0 ? par[28] : par[34]) / currentR;
// par[5] ("subTransY2"): translation along y of second end face
calculatedResidual += -1. * detReducedZ[0] * cos( currentPhi ) * (det==0 ? par[29] : par[35]) / currentR;
// now come the 8*2 beam parameters; calculate the respective parameter index base first (-> which beam)
const unsigned int indexBase = 36 + beam * 2;
// there's no TEC+/- distinction here
// par[6] ("beamRot1"): rotation around z at zt1
calculatedResidual += beamReducedZ[1] * par[indexBase];
// par[7] ("beamRot2"): rotation around z at zt2
calculatedResidual += beamReducedZ[0] * par[indexBase+1];
// now calculate the chisquare
// TODO: check for phi != 0 !!!
chisquare += pow( measuredResidual - calculatedResidual, 2 ) / pow( aMeasuredCoordinates->GetTEC2TECEntry( det, beam, disk ).GetPhiError(), 2 );
} while( moduleLoop.TEC2TECLoop( det, beam, disk ) );
// return the chisquare by ref
f = chisquare;
}
1.7.3