MuonErrorMatrix Class Reference

#include <MuonErrorMatrix.h>

Public Types
enum	action { use, constructor }
	enum type to define if the class is used as a tool or to be created More...

Public Member Functions
void	adjust (FreeTrajectoryState &state)
	adjust the error matrix on the state More...
void	adjust (TrajectoryStateOnSurface &state)
	adjust the error matrix on the state More...
void	close ()
	close the root file attached to the class More...
void	complicatedTerm (const AlgebraicSymMatrix55 &input, AlgebraicSymMatrix55 &output)
	convert sigma/rho -> sigma2/COV More...
int	findBin (TAxis *axis, double value)
	method to get the bin index, taking care of under/overlow: first(1)/last(GetNbins())returned More...
CurvilinearTrajectoryError	get (GlobalVector momentum, bool convolute=true)
	main method to be used. Retrieve a 5x5 symetrical matrix according to parametrization of error or scale factor More...
TProfile3D *	get (int i, int j)
	actually get access to the TProfile3D used for the parametrization More...
CurvilinearTrajectoryError	getFast (GlobalVector momentum)
unsigned int	index (int i, int j)
	MuonErrorMatrix (const edm::ParameterSet &pset)
	constructor from a parameter set More...
void	simpleTerm (const AlgebraicSymMatrix55 &input, AlgebraicSymMatrix55 &output)
	convert sigma2/COV -> sigma/rho More...
	~MuonErrorMatrix ()
	destructor More...

Static Public Member Functions
static bool	divide (CurvilinearTrajectoryError &num_error, const CurvilinearTrajectoryError &denom_error)
	divide term by term the two matrix More...
static void	multiply (CurvilinearTrajectoryError &initial_error, const CurvilinearTrajectoryError &scale_error)
	multiply term by term the two matrix More...
static double	Term (const AlgebraicSymMatrix55 &curv, int i, int j)
	provide the numerical value used. sigma or correlation factor More...

Static Public Attributes
static const TString	vars [5] = {"#frac{q}{|p|}", "#lambda", "#varphi_{0}", "X_{T}", "Y_{T}"}
	names of the variables of the 5x5 error matrix More...

Private Types
enum	TermAction { error, scale, assign }
	decide whether to scale of to assigne terms More...

Private Member Functions
TProfile3D *	Index (int i, int j)
	internal method to get access to the profiles More...
int	Pindex (int i, int j)
	internal methods to get the index of a matrix term. More...
double	Rms (GlobalVector &momentum, int i, int j)
	internal method that retreives the error on the value of the parametrization for term i,j More...
double	Value (GlobalVector &momentum, int i, int j, bool convolute=true)
	internal method that retreives the value of the parametrization for term i,j More...

Private Attributes
std::string	theCategory
	log category: "MuonErrorMatrix" More...
TDirectory *	theD
	the attached root file, where the parametrization is saved More...
TProfile3D *	theData [15]
	15 TProfile, each holding he parametrization of each term of the 5x5 More...
TProfile3D *	theData_fast [5][5]
TermAction	theTermAction [15]

Detailed Description

This class holds an error matrix paramertization. Either the error matrix value or scale factors. The implementation uses TProfile3D with pt, eta, phi axis. Error/Scale factor matrix is obtained using get(GlobalVector momentum)

Author

Jean-Roch Vlimant UCSB

Finn Rebassoo UCSB

Definition at line 24 of file MuonErrorMatrix.h.

Member Enumeration Documentation

enum MuonErrorMatrix::action

enum type to define if the class is used as a tool or to be created

Enumerator
use
constructor

Definition at line 27 of file MuonErrorMatrix.h.

{ use, constructor };

enum MuonErrorMatrix::TermAction

private

decide whether to scale of to assigne terms

Enumerator
error
scale
assign

Definition at line 83 of file MuonErrorMatrix.h.

{ error, scale, assign };

Constructor & Destructor Documentation

MuonErrorMatrix::MuonErrorMatrix

(

const edm::ParameterSet &

pset

)

constructor from a parameter set

Definition at line 17 of file MuonErrorMatrix.cc.

References a, assign, constructor, data, TrackerOfflineValidation_Dqm_cff::dirName, error, edm::ParameterSet::exists(), beamvalidation::exit(), MillePedeFileConverter_cfg::fileName, edm::FileInPath::fullPath(), das_client::fullpath(), edm::ParameterSet::getParameter(), mps_fire::i, dqmiolumiharvest::j, LogDebug, LogTrace, maxEta, HLT_FULL_cff::maxPhi, HLT_FULL_cff::minPhi, HLT_FULL_cff::minPt, Pi, Pindex(), scale, str, AlCaHLTBitMon_QueryRunRegistry::string, theCategory, theD, theData, theData_fast, theTermAction, use, and makeHLTPrescaleTable::values.

                                                               : theD(nullptr) {
  theCategory = "MuonErrorMatrix";
  std::string action = iConfig.getParameter<std::string>("action");

  bool madeFromCff = iConfig.exists("errorMatrixValuesPSet");
  edm::ParameterSet errorMatrixValuesPSet;

  std::string fileName;
  if (!madeFromCff) {
    fileName = iConfig.getParameter<std::string>("rootFileName");
  } else {
    errorMatrixValuesPSet = iConfig.getParameter<edm::ParameterSet>("errorMatrixValuesPSet");
  }
  MuonErrorMatrix::action a = use;

  int NPt = 5;
  std::vector<double> xBins;
  double *xBinsArray = nullptr;
  double minPt = 1;
  double maxPt = 200;
  int NEta = 5;
  std::vector<double> yBins;
  double *yBinsArray = nullptr;
  double minEta = 0;
  double maxEta = 2.5;
  int NPhi = 1;
  double minPhi = -TMath::Pi();
  double maxPhi = TMath::Pi();

  if (action != "use") {
    a = constructor;

    NPt = iConfig.getParameter<int>("NPt");
    if (NPt != 0) {
      minPt = iConfig.getParameter<double>("minPt");
      maxPt = iConfig.getParameter<double>("maxPt");
    } else {
      xBins = iConfig.getParameter<std::vector<double> >("PtBins");
      if (xBins.empty()) {
        edm::LogError(theCategory) << "Npt=0 and no entries in the vector. I will do aseg fault soon.";
      }
      NPt = xBins.size() - 1;
      xBinsArray = &(xBins.front());
      minPt = xBins.front();
      maxPt = xBins.back();
    }

    NEta = iConfig.getParameter<int>("NEta");
    if (NEta != 0) {
      minEta = iConfig.getParameter<double>("minEta");
      maxEta = iConfig.getParameter<double>("maxEta");
    } else {
      yBins = iConfig.getParameter<std::vector<double> >("EtaBins");
      if (yBins.empty()) {
        edm::LogError(theCategory) << "NEta=0 and no entries in the vector. I will do aseg fault soon.";
      }
      NEta = yBins.size() - 1;
      yBinsArray = &(yBins.front());
      minEta = yBins.front();
      maxEta = yBins.back();
    }

    NPhi = iConfig.getParameter<int>("NPhi");
    std::stringstream get(iConfig.getParameter<std::string>("minPhi"));
    if (get.str() == "-Pi") {
      minPhi = -TMath::Pi();
    } else if (get.str() == "Pi") {
      minPhi = TMath::Pi();
    } else {
      get >> minPhi;
    }
    get.str(iConfig.getParameter<std::string>("maxPhi"));
    if (get.str() == "-Pi") {
      maxPhi = -TMath::Pi();
    } else if (get.str() == "Pi") {
      maxPhi = TMath::Pi();
    } else {
      get >> maxPhi;
    }
  }  //action!=use
  else if (madeFromCff) {
    xBins = errorMatrixValuesPSet.getParameter<std::vector<double> >("xAxis");
    NPt = xBins.size() - 1;
    xBinsArray = &(xBins.front());
    minPt = xBins.front();
    maxPt = xBins.back();

    yBins = errorMatrixValuesPSet.getParameter<std::vector<double> >("yAxis");
    NEta = yBins.size() - 1;
    yBinsArray = &(yBins.front());
    minEta = yBins.front();
    maxEta = yBins.back();

    std::vector<double> zBins = errorMatrixValuesPSet.getParameter<std::vector<double> >("zAxis");
    NPhi = 1;
    if (zBins.size() != 2) {
      edm::LogError(theCategory) << "please specify a zAxis with 2 entries only. more bins not implemented yet.";
    }
    minPhi = zBins.front();
    maxPhi = zBins.back();
  }

  if (a == use) {
    if (!madeFromCff) {
      edm::LogInfo(theCategory) << "using an error matrix object from: " << fileName;
      edm::FileInPath data(fileName);
      std::string fullpath = data.fullPath();
      gROOT->cd();
      theD = new TFile(fullpath.c_str());
      theD->SetWritable(false);
    } else {
      static std::atomic<unsigned int> neverTheSame{0};
      std::stringstream dirName("MuonErrorMatrixDirectory");
      dirName << neverTheSame++;
      edm::LogInfo(theCategory)
          << "using an error matrix object from configuration file. putting memory histograms to: " << dirName.str();
      gROOT->cd();
      theD = new TDirectory(dirName.str().c_str(), "transient directory to host MuonErrorMatrix TProfile3D");
      theD->SetWritable(false);
    }
  } else {
    edm::LogInfo(theCategory) << "creating  an error matrix object: " << fileName;
    theD = new TFile(fileName.c_str(), "recreate");
  }

  if (a == use && !madeFromCff) {
    gROOT->cd();
  } else {
    theD->cd();
  }

  for (int i = 0; i != 5; i++) {
    for (int j = i; j != 5; j++) {
      TString pfname(Form("pf3_V%1d%1d", i + 1, j + 1));
      TProfile3D *pf = nullptr;
      if (a == use && !madeFromCff) {
        //read from the rootfile
        edm::LogVerbatim(theCategory) << "getting " << pfname << " from " << fileName;
        pf = (TProfile3D *)theD->Get(pfname);
        theData[Pindex(i, j)] = pf;
        theData_fast[i][j] = pf;
        theData_fast[j][i] = pf;
      } else {
        //    curvilinear coordinate system
        //need to make some input parameter to be to change the number of bins

        TString pftitle;
        if (i == j) {
          pftitle = "#sigma_{" + vars[i] + "}";
        } else {
          pftitle = "#rho(" + vars[i] + "," + vars[j] + ")";
        }
        edm::LogVerbatim(theCategory) << "booking " << pfname << " into " << fileName;
        pf = new TProfile3D(pfname, pftitle, NPt, minPt, maxPt, NEta, minEta, maxEta, NPhi, minPhi, maxPhi, "S");
        pf->SetXTitle("muon p_{T} [GeV]");
        pf->SetYTitle("muon |#eta|");
        pf->SetZTitle("muon #varphi");

        //set variable size binning
        if (xBinsArray) {
          pf->GetXaxis()->Set(NPt, xBinsArray);
        }
        if (yBinsArray) {
          pf->GetYaxis()->Set(NEta, yBinsArray);
        }

        if (madeFromCff) {
          edm::ParameterSet pSet = errorMatrixValuesPSet.getParameter<edm::ParameterSet>(pfname.Data());
          //set the values from the configuration file itself
          std::vector<double> values = pSet.getParameter<std::vector<double> >("values");
          unsigned int iX = pf->GetNbinsX();
          unsigned int iY = pf->GetNbinsY();
          unsigned int iZ = pf->GetNbinsZ();
          unsigned int continuous_i = 0;
          for (unsigned int ix = 1; ix <= iX; ++ix) {
            for (unsigned int iy = 1; iy <= iY; ++iy) {
              for (unsigned int iz = 1; iz <= iZ; ++iz) {
                LogTrace(theCategory) << "filling profile:"
                                      << "\n pt (x)= " << pf->GetXaxis()->GetBinCenter(ix)
                                      << "\n eta (y)= " << pf->GetYaxis()->GetBinCenter(iy)
                                      << "\n phi (z)= " << pf->GetZaxis()->GetBinCenter(iz)
                                      << "\n value= " << values[continuous_i];
                pf->Fill(pf->GetXaxis()->GetBinCenter(ix),
                         pf->GetYaxis()->GetBinCenter(iy),
                         pf->GetZaxis()->GetBinCenter(iz),
                         values[continuous_i++]);
              }
            }
          }
          //term action
          std::string tAction = pSet.getParameter<std::string>("action");
          if (tAction == "scale")
            theTermAction[Pindex(i, j)] = scale;
          else if (tAction == "assign")
            theTermAction[Pindex(i, j)] = assign;
          else {
            edm::LogError(theCategory) << " wrong configuration: term action: " << tAction << " is not recognized.";
            theTermAction[Pindex(i, j)] = error;
          }
        }
      }

      LogDebug(theCategory) << " index " << i << ":" << j << " -> " << Pindex(i, j);
      theData[Pindex(i, j)] = pf;
      theData_fast[i][j] = pf;
      theData_fast[j][i] = pf;
      if (!pf) {
        edm::LogError(theCategory) << " profile " << pfname << " in file " << fileName << " is not valid. exiting.";
        exit(1);
      }
    }
  }

  //verify it
  for (int i = 0; i != 15; i++) {
    if (theData[i]) {
      edm::LogVerbatim(theCategory) << i << " :" << theData[i]->GetName() << " " << theData[i]->GetTitle() << std::endl;
    }
  }
}

MuonErrorMatrix::~MuonErrorMatrix

(

)

destructor

Definition at line 256 of file MuonErrorMatrix.cc.

References close().

{ close(); }

Member Function Documentation

void MuonErrorMatrix::adjust

(

FreeTrajectoryState &

state

)

adjust the error matrix on the state

Definition at line 498 of file MuonErrorMatrix.cc.

References assign, complicatedTerm(), FreeTrajectoryState::curvilinearError(), error, mps_fire::i, dqmiolumiharvest::j, LogDebug, makeMuonMisalignmentScenario::matrix, FreeTrajectoryState::momentum(), FreeTrajectoryState::parameters(), Pindex(), scale, simpleTerm(), theCategory, and theTermAction.

Referenced by TSGForRoadSearch::IPfts(), and TSGForRoadSearch::notAtIPtsos().

                                                       {
  LogDebug(theCategory + "|Adjust") << "state: \n" << state;
  AlgebraicSymMatrix55 simpleTerms;
  simpleTerm(state.curvilinearError(), simpleTerms);
  //the above contains sigma(i), rho(i,j)
  LogDebug(theCategory + "|Adjust") << "state sigma(i), rho(i,j): \n" << simpleTerms;

  AlgebraicSymMatrix55 simpleValues = get(state.momentum(), false).matrix();
  LogDebug(theCategory + "|Adjust") << "config: (i,i), (i,j): \n" << simpleValues;

  for (int i = 0; i != 5; i++) {
    for (int j = i; j != 5; j++) {
      //check on each term for desired action
      switch (theTermAction[Pindex(i, j)]) {
        case scale: {
          simpleTerms(i, j) *= simpleValues(i, j);
          break;
        }
        case assign: {
          simpleTerms(i, j) = simpleValues(i, j);
          break;
        }
        case error: {
          edm::LogError(theCategory + "|Adjust") << " cannot properly adjust for term: " << i << "," << j;
        }
      }
    }
  }
  LogDebug(theCategory + "|Adjust") << "updated state sigma(i), rho(i,j): \n" << simpleTerms;

  AlgebraicSymMatrix55 finalTerms;
  complicatedTerm(simpleTerms, finalTerms);
  LogDebug(theCategory + "|Adjust") << "updated state COV(i,j): \n" << finalTerms;

  CurvilinearTrajectoryError oMat(finalTerms);
  state = FreeTrajectoryState(state.parameters(), oMat);
  LogDebug(theCategory + "|Adjust") << "updated state:\n" << state;
}

void MuonErrorMatrix::adjust

(

TrajectoryStateOnSurface &

state

)

adjust the error matrix on the state

Definition at line 537 of file MuonErrorMatrix.cc.

References assign, complicatedTerm(), TrajectoryStateOnSurface::curvilinearError(), error, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalParameters(), mps_fire::i, dqmiolumiharvest::j, LogDebug, makeMuonMisalignmentScenario::matrix, Pindex(), scale, simpleTerm(), TrajectoryStateOnSurface::surface(), TrajectoryStateOnSurface::surfaceSide(), theCategory, theTermAction, and TrajectoryStateOnSurface::weight().

                                                            {
  AlgebraicSymMatrix55 simpleTerms;
  simpleTerm(state.curvilinearError(), simpleTerms);
  LogDebug(theCategory + "|Adjust") << "state sigma(i), rho(i,j): \n" << simpleTerms;

  AlgebraicSymMatrix55 simpleValues = get(state.globalMomentum(), false).matrix();
  LogDebug(theCategory + "|Adjust") << "config: (i,i), (i,j):\n" << simpleValues;

  for (int i = 0; i != 5; i++) {
    for (int j = i; j != 5; j++) {
      //check on each term for desired action
      switch (theTermAction[Pindex(i, j)]) {
        case scale: {
          simpleTerms(i, j) *= simpleValues(i, j);
          break;
        }
        case assign: {
          simpleTerms(i, j) = simpleValues(i, j);
          break;
        }
        case error: {
          edm::LogError(theCategory + "|Adjust") << " cannot properly adjust for term: " << i << "," << j;
        }
      }
    }
  }
  LogDebug(theCategory + "|Adjust") << "updated state sigma(i), rho(i,j): \n" << simpleTerms;

  AlgebraicSymMatrix55 finalTerms;
  complicatedTerm(simpleTerms, finalTerms);
  LogDebug(theCategory + "|Adjust") << "updated state COV(i,j): \n" << finalTerms;

  CurvilinearTrajectoryError oMat(finalTerms);
  state =
      TrajectoryStateOnSurface(state.weight(), state.globalParameters(), oMat, state.surface(), state.surfaceSide());
  LogDebug(theCategory + "|Adjust") << "updated state:\n" << state;
}

void MuonErrorMatrix::close

(

void

)

close the root file attached to the class

Definition at line 240 of file MuonErrorMatrix.cc.

References mps_fire::i, theD, and theData.

Referenced by esMonitoring.AsyncLineReaderMixin::handle_close(), esMonitoring.FDJsonServer::handle_close(), and ~MuonErrorMatrix().

                            {
  //close the file
  if (theD) {
    theD->cd();
    //write to it first if constructor
    if (theD->IsWritable()) {
      for (int i = 0; i != 15; i++) {
        if (theData[i]) {
          theData[i]->Write();
        }
      }
    }
    theD->Close();
  }
}

void MuonErrorMatrix::complicatedTerm	(	const AlgebraicSymMatrix55 &	input,
		AlgebraicSymMatrix55 &	output
	)

convert sigma/rho -> sigma2/COV

Definition at line 487 of file MuonErrorMatrix.cc.

References mps_fire::i, input, dqmiolumiharvest::j, and convertSQLitetoXML_cfg::output.

Referenced by adjust().

                                                                                                     {
  for (int i = 0; i != 5; i++) {
    for (int j = i; j != 5; j++) {
      if (i == j)
        output(i, j) = input(i, j) * input(i, j);  //sigma squared
      else
        output(i, j) = input(i, j) * input(i, i) * input(j, j);  //rho*sigma*sigma
    }
  }
}

bool MuonErrorMatrix::divide ( CurvilinearTrajectoryError &  num_error, const CurvilinearTrajectoryError &  denom_error  )

static

divide term by term the two matrix

Definition at line 459 of file MuonErrorMatrix.cc.

References mps_fire::i, dqmiolumiharvest::j, and CurvilinearTrajectoryError::matrix().

                                                                                                                 {
  //divide term by term the matrix
  const AlgebraicSymMatrix55 &denom_matrix = denom_error.matrix();
  AlgebraicSymMatrix55 num_matrix = num_error.matrix();
  // the true type of the matrix is such that [i][j] is the same memory object as [j][i]: looping i:0-5, j:0-5 double multiply the terms
  // need to loop only on i:0-5, j:i-5
  for (int i = 0; i != 5; i++) {
    for (int j = i; j != 5; j++) {
      if (denom_matrix(i, j) == 0)
        return false;
      num_matrix(i, j) /= denom_matrix(i, j);
    }
  }
  num_error = CurvilinearTrajectoryError(num_matrix);
  return true;
}

int MuonErrorMatrix::findBin	(	TAxis *	axis,
		double	value
	)

method to get the bin index, taking care of under/overlow: first(1)/last(GetNbins())returned

Definition at line 322 of file MuonErrorMatrix.cc.

References mps_fire::result.

Referenced by getFast(), and Value().

                                                      {
  //find the proper bin, protecting against under/over flow
  int result = axis->FindBin(value);
  if (result <= 0)
    result = 1;  //protect against under flow
  else if (result > axis->GetNbins())
    result = axis->GetNbins();
  return result;
}

CurvilinearTrajectoryError MuonErrorMatrix::get	(	GlobalVector	momentum,
		bool	convolute = true
	)

main method to be used. Retrieve a 5x5 symetrical matrix according to parametrization of error or scale factor

Definition at line 258 of file MuonErrorMatrix.cc.

References mps_fire::i, dqmiolumiharvest::j, cms::cuda::V, and Value().

Referenced by Options.Options::__getitem__(), betterConfigParser.BetterConfigParser::__updateDict(), submitPVValidationJobs.BetterConfigParser::__updateDict(), MuonErrorMatrixAdjuster::fix_cov_matrix(), betterConfigParser.BetterConfigParser::getCompares(), betterConfigParser.BetterConfigParser::getGeneral(), betterConfigParser.BetterConfigParser::getResultingSection(), and submitPVValidationJobs.BetterConfigParser::getResultingSection().

                                                                                     {
  AlgebraicSymMatrix55 V;
  //retrieves a 55 matrix containing (i,i)^2 and (i,j)*(i,i)*(j,j)
  for (int i = 0; i != 5; i++) {
    for (int j = i; j != 5; j++) {
      V(i, j) = Value(momentum, i, j, convolute);
    }
  }
  return CurvilinearTrajectoryError(V);
}

TProfile3D* MuonErrorMatrix::get ( int  i, int  j  )

inline

actually get access to the TProfile3D used for the parametrization

Definition at line 48 of file MuonErrorMatrix.h.

References mps_fire::i, Index(), and dqmiolumiharvest::j.

Referenced by Options.Options::__getitem__(), betterConfigParser.BetterConfigParser::__updateDict(), submitPVValidationJobs.BetterConfigParser::__updateDict(), betterConfigParser.BetterConfigParser::getCompares(), betterConfigParser.BetterConfigParser::getGeneral(), betterConfigParser.BetterConfigParser::getResultingSection(), and submitPVValidationJobs.BetterConfigParser::getResultingSection().

{ return Index(i, j); }

CurvilinearTrajectoryError MuonErrorMatrix::getFast

(

GlobalVector

momentum

)

Definition at line 269 of file MuonErrorMatrix.cc.

References PVValHelper::eta, PV3DBase< T, PVType, FrameType >::eta(), findBin(), mps_fire::i, dqmiolumiharvest::j, PV3DBase< T, PVType, FrameType >::perp(), phi, PV3DBase< T, PVType, FrameType >::phi(), PVValHelper::pT, theData_fast, cms::cuda::V, and makeHLTPrescaleTable::values.

                                                                         {
  //will be faster but make assumptions that could be broken at some point
  //  same bining for all TProfile
  AlgebraicSymMatrix55 V;

  double pT = momentum.perp();
  double eta = fabs(momentum.eta());
  double phi = momentum.phi();

  //assume all the same axis in X,Y,Z
  int iBin_x = findBin(theData_fast[0][0]->GetXaxis(), pT);
  int iBin_y = findBin(theData_fast[0][0]->GetYaxis(), eta);
  int iBin_z = findBin(theData_fast[0][0]->GetZaxis(), phi);

  //retreive values
  double values[5][5];  //sigma_i and rho_ij
  for (int i = 0; i != 5; i++) {
    for (int j = i; j != 5; j++) {
      values[i][j] = theData_fast[i][j]->GetBinContent(iBin_x, iBin_y, iBin_z);
    }
  }

  for (int i = 0; i != 5; i++) {
    for (int j = i; j != 5; j++) {
      if (i == j) {
        //sigma_i * sigma_i
        V(i, j) = values[i][j];
        V(i, j) *= V(i, j);
      } else {
        //sigma_i * sigma_j * rho_ij
        V(i, j) = values[i][i] * values[j][j] * values[i][j];
      }
    }
  }

  return CurvilinearTrajectoryError(V);
}

unsigned int MuonErrorMatrix::index ( int  i, int  j  )

inline

Definition at line 49 of file MuonErrorMatrix.h.

References Pindex().

Referenced by BeautifulSoup.PageElement::_invert().

{ return Pindex(i, j); }

TProfile3D* MuonErrorMatrix::Index ( int  i, int  j  )

inlineprivate

internal method to get access to the profiles

Definition at line 92 of file MuonErrorMatrix.h.

References Pindex(), and theData.

Referenced by get(), Rms(), and Value().

{ return theData[Pindex(i, j)]; }

void MuonErrorMatrix::multiply ( CurvilinearTrajectoryError &  initial_error, const CurvilinearTrajectoryError &  scale_error  )

static

multiply term by term the two matrix

Definition at line 445 of file MuonErrorMatrix.cc.

References mps_fire::i, dqmiolumiharvest::j, and CurvilinearTrajectoryError::matrix().

Referenced by TSGFromPropagation::adjust(), and FastTSGFromPropagation::adjust().

                                                                              {
  //scale term by term the matrix
  const AlgebraicSymMatrix55 &scale_matrix = scale_error.matrix();
  AlgebraicSymMatrix55 revised_matrix = initial_error.matrix();
  // the true type of the matrix is such that [i][j] is the same memory object as [j][i]: looping i:0-5, j:0-5 double multiply the terms
  // need to loop only on i:0-5, j:i-5
  for (int i = 0; i != 5; i++) {
    for (int j = i; j != 5; j++) {
      revised_matrix(i, j) *= scale_matrix(i, j);
    }
  }
  initial_error = CurvilinearTrajectoryError(revised_matrix);
}

int MuonErrorMatrix::Pindex ( int  i, int  j  )

inlineprivate

internal methods to get the index of a matrix term.

Definition at line 87 of file MuonErrorMatrix.h.

References funct::abs(), mps_fire::i, and hltrates_dqm_sourceclient-live_cfg::offset.

Referenced by adjust(), index(), Index(), and MuonErrorMatrix().

                                  {
    static const int offset[5] = {0, 5, 5 + 4, 5 + 4 + 3, 5 + 4 + 3 + 2};
    return offset[i] + abs(j - i);
  }

double MuonErrorMatrix::Rms ( GlobalVector &  momentum, int  i, int  j  )

private

internal method that retreives the error on the value of the parametrization for term i,j

Definition at line 396 of file MuonErrorMatrix.cc.

References PVValHelper::eta, PV3DBase< T, PVType, FrameType >::eta(), Index(), LogDebug, PV3DBase< T, PVType, FrameType >::perp(), phi, PV3DBase< T, PVType, FrameType >::phi(), PVValHelper::pT, mps_fire::result, and theCategory.

                                                                {
  double result = 0;
  TProfile3D *ij = Index(i, j);
  if (!ij) {
    edm::LogError(theCategory) << "cannot get the profile (" << i << ":" << i << ")";
    return result;
  }
  double pT = momentum.perp();
  double eta = fabs(momentum.eta());
  double phi = momentum.phi();

  int iBin_x = ij->GetXaxis()->FindBin(pT);
  int iBin_y = ij->GetYaxis()->FindBin(eta);
  int iBin_z = ij->GetZaxis()->FindBin(phi);
  result = ij->GetBinError(iBin_x, iBin_y, iBin_z);

  LogDebug(theCategory) << "for: (pT,eta,phi)=(" << pT << ", " << eta << ", " << phi << ") error[" << i << "," << j
                        << "] is: " << result;
  return result;
}

void MuonErrorMatrix::simpleTerm	(	const AlgebraicSymMatrix55 &	input,
		AlgebraicSymMatrix55 &	output
	)

convert sigma2/COV -> sigma/rho

Definition at line 476 of file MuonErrorMatrix.cc.

References mps_fire::i, input, dqmiolumiharvest::j, convertSQLitetoXML_cfg::output, and mathSSE::sqrt().

Referenced by adjust().

                                                                                                {
  for (int i = 0; i != 5; i++) {
    for (int j = i; j != 5; j++) {
      if (i == j)
        output(i, j) = sqrt(input(i, j));  //sigma
      else
        output(i, j) = input(i, j) / sqrt(input(i, i) * input(j, j));  //rho
    }
  }
}

double MuonErrorMatrix::Term ( const AlgebraicSymMatrix55 &  curv, int  i, int  j  )

static

provide the numerical value used. sigma or correlation factor

Definition at line 417 of file MuonErrorMatrix.cc.

References mps_fire::result, and mathSSE::sqrt().

                                                                           {
  //return sigma or correlation factor
  double result = 0;
  if (i == j) {
    result = curv(i, j);
    if (result < 0) {
      //check validity of this guy
      edm::LogError("MuonErrorMatrix") << "invalid term in the error matrix.\n sii: " << result;
      return 0;
    }
    return sqrt(result);
  } else {
    double si = curv(i, i);
    double sj = curv(j, j);
    if (si <= 0 || sj <= 0) {
      //check validity
      edm::LogError("MuonErrorMatrix") << "invalid term in the error matrix.\n si: " << si << " sj: " << sj
                                       << ". result will be corrupted\n"
                                       << curv;
      return 0;
    }
    result = curv(i, j) / sqrt(si * sj);
    return result;
  }
  //by default
  return 0;
}

double MuonErrorMatrix::Value ( GlobalVector &  momentum, int  i, int  j, bool  convolute = true  )

private

internal method that retreives the value of the parametrization for term i,j

Definition at line 332 of file MuonErrorMatrix.cc.

References alignCSCRings::corr, PVValHelper::eta, PV3DBase< T, PVType, FrameType >::eta(), findBin(), cuy::ii, Index(), findQualityFiles::jj, LogDebug, PV3DBase< T, PVType, FrameType >::perp(), phi, PV3DBase< T, PVType, FrameType >::phi(), PVValHelper::pT, mps_fire::result, and theCategory.

Referenced by get().

                                                                                  {
  double result = 0;
  TProfile3D *ij = Index(i, j);
  if (!ij) {
    edm::LogError(theCategory) << "cannot get the profile (" << i << ":" << j << ")";
    return result;
  }

  double pT = momentum.perp();
  double eta = fabs(momentum.eta());
  double phi = momentum.phi();

  int iBin_x = findBin(ij->GetXaxis(), pT);
  int iBin_y = findBin(ij->GetYaxis(), eta);
  int iBin_z = findBin(ij->GetZaxis(), phi);

  if (convolute) {
    if (i != j) {
      //return the covariance = correlation*sigma_1 *sigma_2;
      TProfile3D *ii = Index(i, i);
      TProfile3D *jj = Index(j, j);
      if (!ii) {
        edm::LogError(theCategory) << "cannot get the profile (" << i << ":" << i << ")";
        return result;
      }
      if (!jj) {
        edm::LogError(theCategory) << "cannot get the profile (" << j << ":" << j << ")";
        return result;
      }

      int iBin_i_x = findBin(ii->GetXaxis(), pT);
      int iBin_i_y = findBin(ii->GetYaxis(), eta);
      int iBin_i_z = findBin(ii->GetZaxis(), phi);
      int iBin_j_x = findBin(jj->GetXaxis(), pT);
      int iBin_j_y = findBin(jj->GetYaxis(), eta);
      int iBin_j_z = findBin(jj->GetZaxis(), phi);

      double corr = ij->GetBinContent(iBin_x, iBin_y, iBin_z);
      double sigma_1 = (ii->GetBinContent(iBin_i_x, iBin_i_y, iBin_i_z));
      double sigma_2 = (jj->GetBinContent(iBin_j_x, iBin_j_y, iBin_j_z));

      result = corr * sigma_1 * sigma_2;
      LogDebug(theCategory) << "for: (pT,eta,phi)=(" << pT << ", " << eta << ", " << phi << ") nterms are"
                            << "\nrho[" << i << "," << j << "]: " << corr << " [" << iBin_x << ", " << iBin_y << ", "
                            << iBin_z << "]"
                            << "\nsigma[" << i << "," << i << "]: " << sigma_1 << "\nsigma[" << j << "," << j
                            << "]: " << sigma_2 << "Covariance[" << i << "," << j << "] is: " << result;
      return result;
    } else {
      //return the variance = sigma_1 **2
      //    result=ij->GetBinContent(iBin);
      result = ij->GetBinContent(iBin_x, iBin_y, iBin_z);
      result *= result;
      LogDebug(theCategory) << "for: (pT,eta,phi)=(" << pT << ", " << eta << ", " << phi << ") sigma^2[" << i << ","
                            << j << "] is: " << result;
      return result;
    }
  } else {
    //do not convolute
    result = ij->GetBinContent(iBin_x, iBin_y, iBin_z);
    return result;
  }
}

Member Data Documentation

std::string MuonErrorMatrix::theCategory

private

log category: "MuonErrorMatrix"

Definition at line 74 of file MuonErrorMatrix.h.

Referenced by adjust(), MuonErrorMatrix(), Rms(), and Value().

TDirectory* MuonErrorMatrix::theD

private

the attached root file, where the parametrization is saved

Definition at line 77 of file MuonErrorMatrix.h.

Referenced by close(), and MuonErrorMatrix().

TProfile3D* MuonErrorMatrix::theData[15]

private

15 TProfile, each holding he parametrization of each term of the 5x5

Definition at line 79 of file MuonErrorMatrix.h.

Referenced by close(), Index(), and MuonErrorMatrix().

TProfile3D* MuonErrorMatrix::theData_fast[5][5]

private

Definition at line 80 of file MuonErrorMatrix.h.

Referenced by getFast(), and MuonErrorMatrix().

TermAction MuonErrorMatrix::theTermAction[15]

private

Definition at line 84 of file MuonErrorMatrix.h.

Referenced by adjust(), and MuonErrorMatrix().

const TString MuonErrorMatrix::vars = {"#frac{q}{|p|}", "#lambda", "#varphi_{0}", "X_{T}", "Y_{T}"}

static

names of the variables of the 5x5 error matrix

Definition at line 52 of file MuonErrorMatrix.h.