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 388 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(), CTPPSOpticalFunctionsESSource::produce(), npstat::ArrayND< Numeric, StackLen, StackDim >::projectLoop(), npstat::ArrayND< Numeric, StackLen, StackDim >::projectLoop2(), PVFitter::runBXFitter(), PVFitter::runFitter(), and MuonResidualsFitter::selectedResidualsFlags().

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