MassConstraint.C

Go to the documentation of this file.

// Class MassConstraint
//
// Author: Jenny Boehme
// Last update: $Date: 2005/01/12 10:11:45 $
//          by: $Author: blist $
// 
// Description: M(p1,p2)=M_W constraint
//               

#include "jbltools/kinfit/MassConstraint.h"
#include "jbltools/kinfit/ParticleFitObject.h"

#include<iostream>
#include<cmath>
#include<cassert>

using std::cerr;
using std::cout;
using std::endl;

// constructor
MassConstraint::MassConstraint (double mass_) : mass(mass_) {}

// destructor
MassConstraint::~MassConstraint () {}

// calulate current value of constraint function
double MassConstraint::getValue() const {
  double totE[2] = {0,0};
  double totpx[2] = {0,0}; 
  double totpy[2] = {0,0}; 
  double totpz[2] = {0,0}; 
  for (unsigned int i = 0; i < fitobjects.size(); i++) {
    int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
    totE[index] += fitobjects[i]->getE(); 
    totpx[index] += fitobjects[i]->getPx(); 
    totpy[index] += fitobjects[i]->getPy(); 
    totpz[index] += fitobjects[i]->getPz(); 
  }
  double result = -mass;
  result += std::sqrt(std::abs(totE[0]*totE[0]-totpx[0]*totpx[0]-totpy[0]*totpy[0]-totpz[0]*totpz[0]));
  result -= std::sqrt(std::abs(totE[1]*totE[1]-totpx[1]*totpx[1]-totpy[1]*totpy[1]-totpz[1]*totpz[1]));
  return result;
}

// calculate vector/array of derivatives of this contraint 
// w.r.t. to ALL parameters of all fitobjects
// here: d M /d par(j) 
//          = d M /d p(i) * d p(i) /d par(j)
//          =  +-1/M * p(i) * d p(i) /d par(j)
void MassConstraint::getDerivatives(int idim, double der[]) const {
  double totE[2] = {0,0};
  double totpx[2] = {0,0}; 
  double totpy[2] = {0,0}; 
  double totpz[2] = {0,0}; 
  for (unsigned int i = 0; i < fitobjects.size(); i++) {
    int index = (flags[i]==1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
    totE[index] += fitobjects[i]->getE(); 
    totpx[index] += fitobjects[i]->getPx(); 
    totpy[index] += fitobjects[i]->getPy(); 
    totpz[index] += fitobjects[i]->getPz(); 
  }
  for (unsigned int i = 0; i < fitobjects.size(); i++) {
    int index = (flags[i]==1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
    for (int ilocal = 0; ilocal < fitobjects[i]->getNPar(); ilocal++) {
      if (!fitobjects[i]->isParamFixed(ilocal)) {
        int iglobal = fitobjects[i]->getGlobalParNum (ilocal);
        assert (iglobal >= 0 && iglobal < idim);
        der[iglobal] = 0;
        der[iglobal] += totE[index] * fitobjects[i]->getDE (ilocal);
        der[iglobal] -= totpx[index] * fitobjects[i]->getDPx (ilocal);
        der[iglobal] -= totpy[index] * fitobjects[i]->getDPy (ilocal);
        der[iglobal] -= totpz[index] * fitobjects[i]->getDPz (ilocal);
        double m2 = totE[index]*totE[index] - totpx[index]*totpx[index]
                    - totpy[index]*totpy[index] - totpz[index]*totpz[index];
        if (m2 < 0) cerr << "MassConstraint::getDerivatives: m2<0!" << endl;
        if (m2 == 0) cerr << "MassConstraint::getDerivatives: m2==0!" << endl;
        else der[iglobal] /= std::sqrt (std::abs(m2));
        if (index == 1) der[iglobal] *= -1.;
      }
    }
  }
}
  
double MassConstraint::getMass (int flag) {
  double totE = 0;
  double totpx = 0; 
  double totpy = 0; 
  double totpz = 0; 
  for (unsigned int i = 0; i < fitobjects.size(); i++) {
    if (flags[i] == flag) {
      totE += fitobjects[i]->getE(); 
      totpx += fitobjects[i]->getPx(); 
      totpy += fitobjects[i]->getPy(); 
      totpz += fitobjects[i]->getPz(); 
    }
  }
  return std::sqrt(std::abs(totE*totE-totpx*totpx-totpy*totpy-totpz*totpz));
}

---