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 29 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 32 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 88 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, MuonErrorMatrixValues_cff::errorMatrixValuesPSet, edm::ParameterSet::exists(), cmsRelvalreport::exit, MillePedeFileConverter_cfg::fileName, edm::FileInPath::fullPath(), edm::ParameterSet::getParameter(), mps_fire::i, LogDebug, LogTrace, maxEta, trackingParticleSelector_cfi::maxPhi, MuonErrorMatrixAnalyzer_cfi::maxPt, cutBasedElectronID_CSA14_50ns_V0_cff::minEta, trackingParticleSelector_cfi::minPhi, lostTracks_cfi::minPt, MuonErrorMatrixAnalyzer_cfi::NEta, MuonErrorMatrixAnalyzer_cfi::NPhi, MuonErrorMatrixAnalyzer_cfi::NPt, packedPFCandidateRefMixer_cfi::pf, Pi, Pindex(), scale, harvestTrackValidationPlots::str, AlCaHLTBitMon_QueryRunRegistry::string, theCategory, theD, theData, theData_fast, theTermAction, use, MuonErrorMatrixValues_cff::values, and vars.

                                                               :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());
      minPt = yBins.front();
      maxPt = 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());
    minPt = yBins.front();
    maxPt = 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 225 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 446 of file MuonErrorMatrix.cc.

References assign, complicatedTerm(), FreeTrajectoryState::curvilinearError(), error, mps_fire::i, 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 489 of file MuonErrorMatrix.cc.

References assign, complicatedTerm(), TrajectoryStateOnSurface::curvilinearError(), error, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalParameters(), mps_fire::i, 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 215 of file MuonErrorMatrix.cc.

References mps_fire::i, theD, and theData.

Referenced by lumiQTWidget.ApplicationWindow::fileQuit(), esMonitoring.AsyncLineReaderMixin::handle_close(), esMonitoring.FDJsonServer::handle_close(), Vispa.Gui.BoxContentDialog.BoxContentDialog::keyPressEvent(), Vispa.Gui.FindDialog.FindDialog::keyPressEvent(), 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 434 of file MuonErrorMatrix.cc.

References mps_fire::i, input, 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 411 of file MuonErrorMatrix.cc.

References mps_fire::i, 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 289 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 231 of file MuonErrorMatrix.cc.

References mps_fire::i, and Value().

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

                                                                                     {
  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 53 of file MuonErrorMatrix.h.

References Index().

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

{return Index(i,j);}

CurvilinearTrajectoryError MuonErrorMatrix::getFast

(

GlobalVector

momentum

)

Definition at line 238 of file MuonErrorMatrix.cc.

References PVValHelper::eta, PV3DBase< T, PVType, FrameType >::eta(), findBin(), mps_fire::i, PV3DBase< T, PVType, FrameType >::perp(), phi, PV3DBase< T, PVType, FrameType >::phi(), PVValHelper::pT, theData_fast, and MuonErrorMatrixValues_cff::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 54 of file MuonErrorMatrix.h.

References Pindex().

Referenced by BeautifulSoup.PageElement::insert().

{return Pindex(i,j);}

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

inlineprivate

internal method to get access to the profiles

Definition at line 96 of file MuonErrorMatrix.h.

References Pindex(), Rms(), and Value().

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 400 of file MuonErrorMatrix.cc.

References mps_fire::i, 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 92 of file MuonErrorMatrix.h.

References funct::abs(), mps_fire::i, and PFRecoTauDiscriminationByIsolation_cfi::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 354 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.

Referenced by Index().

                                                                  {
  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 425 of file MuonErrorMatrix.cc.

References mps_fire::i, input, 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 371 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){
      edm::LogError("MuonErrorMatrix")<<"invalid term in the error matrix.\n si: "
                <<si<<" sj: "<<sj<<". result will be corrupted\n"<<curv;
      return 0;}

    si=sqrt(si);
    sj=sqrt(sj);
    //check validity

    return result = curv(i,j)/(si*sj);}
  //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 297 of file MuonErrorMatrix.cc.

References 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(), and Index().

                                                                                   {
  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 79 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 82 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 84 of file MuonErrorMatrix.h.

Referenced by close(), and MuonErrorMatrix().

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

private

Definition at line 85 of file MuonErrorMatrix.h.

Referenced by getFast(), and MuonErrorMatrix().

TermAction MuonErrorMatrix::theTermAction[15]

private

Definition at line 89 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 57 of file MuonErrorMatrix.h.

Referenced by MuonErrorMatrix().