LASBarrelAlgorithm.h File Reference

#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)

Function Documentation

void fcn	(	int &	npar,
		double *	gin,
		double &	f,
		double *	par,
		int	iflag
	)

minuit chisquare func

Definition at line 343 of file LASBarrelAlgorithm.cc.

References EcalCondDBWriter_cfi::beam, funct::cos(), 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(), 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(), PVFitter::runFitter(), and HBHENegativeFlagSetter::setPulseShapeFlags().

                                                                       {

  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;
  
}