Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:18 

0001 #include "RecoMuon/TrackingTools/interface/MuonErrorMatrix.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/ParameterSet/interface/FileInPath.h"
0004 
0005 #include "TROOT.h"
0006 #include "TString.h"
0007 #include "TRandom2.h"
0008 #include "TMath.h"
0009 
0010 #include <sstream>
0011 #include <atomic>
0012 
0013 using namespace std;
0014 
0015 const TString MuonErrorMatrix::vars[5] = {"#frac{q}{|p|}", "#lambda", "#varphi_{0}", "X_{T}", "Y_{T}"};
0016 
0017 MuonErrorMatrix::MuonErrorMatrix(const edm::ParameterSet &iConfig) : theD(nullptr) {
0018   theCategory = "MuonErrorMatrix";
0019   std::string action = iConfig.getParameter<std::string>("action");
0020 
0021   bool madeFromCff = iConfig.exists("errorMatrixValuesPSet");
0022   edm::ParameterSet errorMatrixValuesPSet;
0023 
0024   std::string fileName;
0025   if (!madeFromCff) {
0026     fileName = iConfig.getParameter<std::string>("rootFileName");
0027   } else {
0028     errorMatrixValuesPSet = iConfig.getParameter<edm::ParameterSet>("errorMatrixValuesPSet");
0029   }
0030   MuonErrorMatrix::action a = use;
0031 
0032   int NPt = 5;
0033   std::vector<double> xBins;
0034   double *xBinsArray = nullptr;
0035   double minPt = 1;
0036   double maxPt = 200;
0037   int NEta = 5;
0038   std::vector<double> yBins;
0039   double *yBinsArray = nullptr;
0040   double minEta = 0;
0041   double maxEta = 2.5;
0042   int NPhi = 1;
0043   double minPhi = -TMath::Pi();
0044   double maxPhi = TMath::Pi();
0045 
0046   if (action != "use") {
0047     a = constructor;
0048 
0049     NPt = iConfig.getParameter<int>("NPt");
0050     if (NPt != 0) {
0051       minPt = iConfig.getParameter<double>("minPt");
0052       maxPt = iConfig.getParameter<double>("maxPt");
0053     } else {
0054       xBins = iConfig.getParameter<std::vector<double> >("PtBins");
0055       if (xBins.empty()) {
0056         edm::LogError(theCategory) << "Npt=0 and no entries in the vector. I will do aseg fault soon.";
0057       }
0058       NPt = xBins.size() - 1;
0059       xBinsArray = &(xBins.front());
0060       minPt = xBins.front();
0061       maxPt = xBins.back();
0062     }
0063 
0064     NEta = iConfig.getParameter<int>("NEta");
0065     if (NEta != 0) {
0066       minEta = iConfig.getParameter<double>("minEta");
0067       maxEta = iConfig.getParameter<double>("maxEta");
0068     } else {
0069       yBins = iConfig.getParameter<std::vector<double> >("EtaBins");
0070       if (yBins.empty()) {
0071         edm::LogError(theCategory) << "NEta=0 and no entries in the vector. I will do aseg fault soon.";
0072       }
0073       NEta = yBins.size() - 1;
0074       yBinsArray = &(yBins.front());
0075       minEta = yBins.front();
0076       maxEta = yBins.back();
0077     }
0078 
0079     NPhi = iConfig.getParameter<int>("NPhi");
0080     std::stringstream get(iConfig.getParameter<std::string>("minPhi"));
0081     if (get.str() == "-Pi") {
0082       minPhi = -TMath::Pi();
0083     } else if (get.str() == "Pi") {
0084       minPhi = TMath::Pi();
0085     } else {
0086       get >> minPhi;
0087     }
0088     get.str(iConfig.getParameter<std::string>("maxPhi"));
0089     if (get.str() == "-Pi") {
0090       maxPhi = -TMath::Pi();
0091     } else if (get.str() == "Pi") {
0092       maxPhi = TMath::Pi();
0093     } else {
0094       get >> maxPhi;
0095     }
0096   }  //action!=use
0097   else if (madeFromCff) {
0098     xBins = errorMatrixValuesPSet.getParameter<std::vector<double> >("xAxis");
0099     NPt = xBins.size() - 1;
0100     xBinsArray = &(xBins.front());
0101     minPt = xBins.front();
0102     maxPt = xBins.back();
0103 
0104     yBins = errorMatrixValuesPSet.getParameter<std::vector<double> >("yAxis");
0105     NEta = yBins.size() - 1;
0106     yBinsArray = &(yBins.front());
0107     minEta = yBins.front();
0108     maxEta = yBins.back();
0109 
0110     std::vector<double> zBins = errorMatrixValuesPSet.getParameter<std::vector<double> >("zAxis");
0111     NPhi = 1;
0112     if (zBins.size() != 2) {
0113       edm::LogError(theCategory) << "please specify a zAxis with 2 entries only. more bins not implemented yet.";
0114     }
0115     minPhi = zBins.front();
0116     maxPhi = zBins.back();
0117   }
0118 
0119   if (a == use) {
0120     if (!madeFromCff) {
0121       edm::LogInfo(theCategory) << "using an error matrix object from: " << fileName;
0122       edm::FileInPath data(fileName);
0123       std::string fullpath = data.fullPath();
0124       gROOT->cd();
0125       theD = new TFile(fullpath.c_str());
0126       theD->SetWritable(false);
0127     } else {
0128       static std::atomic<unsigned int> neverTheSame{0};
0129       std::stringstream dirName("MuonErrorMatrixDirectory");
0130       dirName << neverTheSame++;
0131       edm::LogInfo(theCategory)
0132           << "using an error matrix object from configuration file. putting memory histograms to: " << dirName.str();
0133       gROOT->cd();
0134       theD = new TDirectory(dirName.str().c_str(), "transient directory to host MuonErrorMatrix TProfile3D");
0135       theD->SetWritable(false);
0136     }
0137   } else {
0138     edm::LogInfo(theCategory) << "creating  an error matrix object: " << fileName;
0139     theD = new TFile(fileName.c_str(), "recreate");
0140   }
0141 
0142   if (a == use && !madeFromCff) {
0143     gROOT->cd();
0144   } else {
0145     theD->cd();
0146   }
0147 
0148   for (int i = 0; i != 5; i++) {
0149     for (int j = i; j != 5; j++) {
0150       TString pfname(Form("pf3_V%1d%1d", i + 1, j + 1));
0151       TProfile3D *pf = nullptr;
0152       if (a == use && !madeFromCff) {
0153         //read from the rootfile
0154         edm::LogVerbatim(theCategory) << "getting " << pfname << " from " << fileName;
0155         pf = (TProfile3D *)theD->Get(pfname);
0156         theData[Pindex(i, j)] = pf;
0157         theData_fast[i][j] = pf;
0158         theData_fast[j][i] = pf;
0159       } else {
0160         //    curvilinear coordinate system
0161         //need to make some input parameter to be to change the number of bins
0162 
0163         TString pftitle;
0164         if (i == j) {
0165           pftitle = "#sigma_{" + vars[i] + "}";
0166         } else {
0167           pftitle = "#rho(" + vars[i] + "," + vars[j] + ")";
0168         }
0169         edm::LogVerbatim(theCategory) << "booking " << pfname << " into " << fileName;
0170         pf = new TProfile3D(pfname, pftitle, NPt, minPt, maxPt, NEta, minEta, maxEta, NPhi, minPhi, maxPhi, "S");
0171         pf->SetXTitle("muon p_{T} [GeV]");
0172         pf->SetYTitle("muon |#eta|");
0173         pf->SetZTitle("muon #varphi");
0174 
0175         //set variable size binning
0176         if (xBinsArray) {
0177           pf->GetXaxis()->Set(NPt, xBinsArray);
0178         }
0179         if (yBinsArray) {
0180           pf->GetYaxis()->Set(NEta, yBinsArray);
0181         }
0182 
0183         if (madeFromCff) {
0184           edm::ParameterSet pSet = errorMatrixValuesPSet.getParameter<edm::ParameterSet>(pfname.Data());
0185           //set the values from the configuration file itself
0186           std::vector<double> values = pSet.getParameter<std::vector<double> >("values");
0187           unsigned int iX = pf->GetNbinsX();
0188           unsigned int iY = pf->GetNbinsY();
0189           unsigned int iZ = pf->GetNbinsZ();
0190           unsigned int continuous_i = 0;
0191           for (unsigned int ix = 1; ix <= iX; ++ix) {
0192             for (unsigned int iy = 1; iy <= iY; ++iy) {
0193               for (unsigned int iz = 1; iz <= iZ; ++iz) {
0194                 LogTrace(theCategory) << "filling profile:"
0195                                       << "\n pt (x)= " << pf->GetXaxis()->GetBinCenter(ix)
0196                                       << "\n eta (y)= " << pf->GetYaxis()->GetBinCenter(iy)
0197                                       << "\n phi (z)= " << pf->GetZaxis()->GetBinCenter(iz)
0198                                       << "\n value= " << values[continuous_i];
0199                 pf->Fill(pf->GetXaxis()->GetBinCenter(ix),
0200                          pf->GetYaxis()->GetBinCenter(iy),
0201                          pf->GetZaxis()->GetBinCenter(iz),
0202                          values[continuous_i++]);
0203               }
0204             }
0205           }
0206           //term action
0207           std::string tAction = pSet.getParameter<std::string>("action");
0208           if (tAction == "scale")
0209             theTermAction[Pindex(i, j)] = scale;
0210           else if (tAction == "assign")
0211             theTermAction[Pindex(i, j)] = assign;
0212           else {
0213             edm::LogError(theCategory) << " wrong configuration: term action: " << tAction << " is not recognized.";
0214             theTermAction[Pindex(i, j)] = error;
0215           }
0216         }
0217       }
0218 
0219       LogDebug(theCategory) << " index " << i << ":" << j << " -> " << Pindex(i, j);
0220       theData[Pindex(i, j)] = pf;
0221       theData_fast[i][j] = pf;
0222       theData_fast[j][i] = pf;
0223       if (!pf) {
0224         edm::LogError(theCategory) << " profile " << pfname << " in file " << fileName << " is not valid. exiting.";
0225         exit(1);
0226       }
0227     }
0228   }
0229 
0230   //verify it
0231   for (int i = 0; i != 15; i++) {
0232     if (theData[i]) {
0233       edm::LogVerbatim(theCategory) << i << " :" << theData[i]->GetName() << " " << theData[i]->GetTitle() << std::endl;
0234     }
0235   }
0236 }
0237 
0238 //void MuonErrorMatrix::writeIntoCFF(){}
0239 
0240 void MuonErrorMatrix::close() {
0241   //close the file
0242   if (theD) {
0243     theD->cd();
0244     //write to it first if constructor
0245     if (theD->IsWritable()) {
0246       for (int i = 0; i != 15; i++) {
0247         if (theData[i]) {
0248           theData[i]->Write();
0249         }
0250       }
0251     }
0252     theD->Close();
0253   }
0254 }
0255 
0256 MuonErrorMatrix::~MuonErrorMatrix() { close(); }
0257 
0258 CurvilinearTrajectoryError MuonErrorMatrix::get(GlobalVector momentum, bool convolute) {
0259   AlgebraicSymMatrix55 V;
0260   //retrieves a 55 matrix containing (i,i)^2 and (i,j)*(i,i)*(j,j)
0261   for (int i = 0; i != 5; i++) {
0262     for (int j = i; j != 5; j++) {
0263       V(i, j) = Value(momentum, i, j, convolute);
0264     }
0265   }
0266   return CurvilinearTrajectoryError(V);
0267 }
0268 
0269 CurvilinearTrajectoryError MuonErrorMatrix::getFast(GlobalVector momentum) {
0270   //will be faster but make assumptions that could be broken at some point
0271   //  same bining for all TProfile
0272   AlgebraicSymMatrix55 V;
0273 
0274   double pT = momentum.perp();
0275   double eta = fabs(momentum.eta());
0276   double phi = momentum.phi();
0277 
0278   //assume all the same axis in X,Y,Z
0279   int iBin_x = findBin(theData_fast[0][0]->GetXaxis(), pT);
0280   int iBin_y = findBin(theData_fast[0][0]->GetYaxis(), eta);
0281   int iBin_z = findBin(theData_fast[0][0]->GetZaxis(), phi);
0282 
0283   //retreive values
0284   double values[5][5];  //sigma_i and rho_ij
0285   for (int i = 0; i != 5; i++) {
0286     for (int j = i; j != 5; j++) {
0287       values[i][j] = theData_fast[i][j]->GetBinContent(iBin_x, iBin_y, iBin_z);
0288     }
0289   }
0290 
0291   for (int i = 0; i != 5; i++) {
0292     for (int j = i; j != 5; j++) {
0293       if (i == j) {
0294         //sigma_i * sigma_i
0295         V(i, j) = values[i][j];
0296         V(i, j) *= V(i, j);
0297       } else {
0298         //sigma_i * sigma_j * rho_ij
0299         V(i, j) = values[i][i] * values[j][j] * values[i][j];
0300       }
0301     }
0302   }
0303 
0304   return CurvilinearTrajectoryError(V);
0305 }
0306 
0307 /*CurvilinearTrajectoryError MuonErrorMatrix::get_random(GlobalVector momentum)  { 
0308   static TRandom2 rand;
0309   AlgebraicSymMatrix55 V;//result
0310   //first proceed with diagonal elements
0311   for (int i=0;i!=5;i++){
0312   V(i,i)=rand.Gaus( Value(momentum,i,i), Rms(momentum,i,i));}
0313 
0314   //now proceed with the correlations
0315   for (int i=0;i!=5;i++){for (int j=i+1;j<5;j++){
0316       double corr = rand.Gaus( Value(momentum,i,j), Rms(momentum,i,j));
0317       //assign the covariance from correlation and sigmas
0318       V(i,j)= corr * sqrt( V[i][i] * V[j][j]);}}
0319       return CurvilinearTrajectoryError(V);  }
0320 */
0321 
0322 int MuonErrorMatrix::findBin(TAxis *axis, double value) {
0323   //find the proper bin, protecting against under/over flow
0324   int result = axis->FindBin(value);
0325   if (result <= 0)
0326     result = 1;  //protect against under flow
0327   else if (result > axis->GetNbins())
0328     result = axis->GetNbins();
0329   return result;
0330 }
0331 
0332 double MuonErrorMatrix::Value(GlobalVector &momentum, int i, int j, bool convolute) {
0333   double result = 0;
0334   TProfile3D *ij = Index(i, j);
0335   if (!ij) {
0336     edm::LogError(theCategory) << "cannot get the profile (" << i << ":" << j << ")";
0337     return result;
0338   }
0339 
0340   double pT = momentum.perp();
0341   double eta = fabs(momentum.eta());
0342   double phi = momentum.phi();
0343 
0344   int iBin_x = findBin(ij->GetXaxis(), pT);
0345   int iBin_y = findBin(ij->GetYaxis(), eta);
0346   int iBin_z = findBin(ij->GetZaxis(), phi);
0347 
0348   if (convolute) {
0349     if (i != j) {
0350       //return the covariance = correlation*sigma_1 *sigma_2;
0351       TProfile3D *ii = Index(i, i);
0352       TProfile3D *jj = Index(j, j);
0353       if (!ii) {
0354         edm::LogError(theCategory) << "cannot get the profile (" << i << ":" << i << ")";
0355         return result;
0356       }
0357       if (!jj) {
0358         edm::LogError(theCategory) << "cannot get the profile (" << j << ":" << j << ")";
0359         return result;
0360       }
0361 
0362       int iBin_i_x = findBin(ii->GetXaxis(), pT);
0363       int iBin_i_y = findBin(ii->GetYaxis(), eta);
0364       int iBin_i_z = findBin(ii->GetZaxis(), phi);
0365       int iBin_j_x = findBin(jj->GetXaxis(), pT);
0366       int iBin_j_y = findBin(jj->GetYaxis(), eta);
0367       int iBin_j_z = findBin(jj->GetZaxis(), phi);
0368 
0369       double corr = ij->GetBinContent(iBin_x, iBin_y, iBin_z);
0370       double sigma_1 = (ii->GetBinContent(iBin_i_x, iBin_i_y, iBin_i_z));
0371       double sigma_2 = (jj->GetBinContent(iBin_j_x, iBin_j_y, iBin_j_z));
0372 
0373       result = corr * sigma_1 * sigma_2;
0374       LogDebug(theCategory) << "for: (pT,eta,phi)=(" << pT << ", " << eta << ", " << phi << ") nterms are"
0375                             << "\nrho[" << i << "," << j << "]: " << corr << " [" << iBin_x << ", " << iBin_y << ", "
0376                             << iBin_z << "]"
0377                             << "\nsigma[" << i << "," << i << "]: " << sigma_1 << "\nsigma[" << j << "," << j
0378                             << "]: " << sigma_2 << "Covariance[" << i << "," << j << "] is: " << result;
0379       return result;
0380     } else {
0381       //return the variance = sigma_1 **2
0382       //    result=ij->GetBinContent(iBin);
0383       result = ij->GetBinContent(iBin_x, iBin_y, iBin_z);
0384       result *= result;
0385       LogDebug(theCategory) << "for: (pT,eta,phi)=(" << pT << ", " << eta << ", " << phi << ") sigma^2[" << i << ","
0386                             << j << "] is: " << result;
0387       return result;
0388     }
0389   } else {
0390     //do not convolute
0391     result = ij->GetBinContent(iBin_x, iBin_y, iBin_z);
0392     return result;
0393   }
0394 }
0395 
0396 double MuonErrorMatrix::Rms(GlobalVector &momentum, int i, int j) {
0397   double result = 0;
0398   TProfile3D *ij = Index(i, j);
0399   if (!ij) {
0400     edm::LogError(theCategory) << "cannot get the profile (" << i << ":" << i << ")";
0401     return result;
0402   }
0403   double pT = momentum.perp();
0404   double eta = fabs(momentum.eta());
0405   double phi = momentum.phi();
0406 
0407   int iBin_x = ij->GetXaxis()->FindBin(pT);
0408   int iBin_y = ij->GetYaxis()->FindBin(eta);
0409   int iBin_z = ij->GetZaxis()->FindBin(phi);
0410   result = ij->GetBinError(iBin_x, iBin_y, iBin_z);
0411 
0412   LogDebug(theCategory) << "for: (pT,eta,phi)=(" << pT << ", " << eta << ", " << phi << ") error[" << i << "," << j
0413                         << "] is: " << result;
0414   return result;
0415 }
0416 
0417 double MuonErrorMatrix::Term(const AlgebraicSymMatrix55 &curv, int i, int j) {
0418   //return sigma or correlation factor
0419   double result = 0;
0420   if (i == j) {
0421     result = curv(i, j);
0422     if (result < 0) {
0423       //check validity of this guy
0424       edm::LogError("MuonErrorMatrix") << "invalid term in the error matrix.\n sii: " << result;
0425       return 0;
0426     }
0427     return sqrt(result);
0428   } else {
0429     double si = curv(i, i);
0430     double sj = curv(j, j);
0431     if (si <= 0 || sj <= 0) {
0432       //check validity
0433       edm::LogError("MuonErrorMatrix") << "invalid term in the error matrix.\n si: " << si << " sj: " << sj
0434                                        << ". result will be corrupted\n"
0435                                        << curv;
0436       return 0;
0437     }
0438     result = curv(i, j) / sqrt(si * sj);
0439     return result;
0440   }
0441   //by default
0442   return 0;
0443 }
0444 
0445 void MuonErrorMatrix::multiply(CurvilinearTrajectoryError &initial_error,
0446                                const CurvilinearTrajectoryError &scale_error) {
0447   //scale term by term the matrix
0448   const AlgebraicSymMatrix55 &scale_matrix = scale_error.matrix();
0449   AlgebraicSymMatrix55 revised_matrix = initial_error.matrix();
0450   // 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
0451   // need to loop only on i:0-5, j:i-5
0452   for (int i = 0; i != 5; i++) {
0453     for (int j = i; j != 5; j++) {
0454       revised_matrix(i, j) *= scale_matrix(i, j);
0455     }
0456   }
0457   initial_error = CurvilinearTrajectoryError(revised_matrix);
0458 }
0459 bool MuonErrorMatrix::divide(CurvilinearTrajectoryError &num_error, const CurvilinearTrajectoryError &denom_error) {
0460   //divide term by term the matrix
0461   const AlgebraicSymMatrix55 &denom_matrix = denom_error.matrix();
0462   AlgebraicSymMatrix55 num_matrix = num_error.matrix();
0463   // 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
0464   // need to loop only on i:0-5, j:i-5
0465   for (int i = 0; i != 5; i++) {
0466     for (int j = i; j != 5; j++) {
0467       if (denom_matrix(i, j) == 0)
0468         return false;
0469       num_matrix(i, j) /= denom_matrix(i, j);
0470     }
0471   }
0472   num_error = CurvilinearTrajectoryError(num_matrix);
0473   return true;
0474 }
0475 
0476 void MuonErrorMatrix::simpleTerm(const AlgebraicSymMatrix55 &input, AlgebraicSymMatrix55 &output) {
0477   for (int i = 0; i != 5; i++) {
0478     for (int j = i; j != 5; j++) {
0479       if (i == j)
0480         output(i, j) = sqrt(input(i, j));  //sigma
0481       else
0482         output(i, j) = input(i, j) / sqrt(input(i, i) * input(j, j));  //rho
0483     }
0484   }
0485 }
0486 
0487 void MuonErrorMatrix::complicatedTerm(const AlgebraicSymMatrix55 &input, AlgebraicSymMatrix55 &output) {
0488   for (int i = 0; i != 5; i++) {
0489     for (int j = i; j != 5; j++) {
0490       if (i == j)
0491         output(i, j) = input(i, j) * input(i, j);  //sigma squared
0492       else
0493         output(i, j) = input(i, j) * input(i, i) * input(j, j);  //rho*sigma*sigma
0494     }
0495   }
0496 }
0497 
0498 void MuonErrorMatrix::adjust(FreeTrajectoryState &state) {
0499   LogDebug(theCategory + "|Adjust") << "state: \n" << state;
0500   AlgebraicSymMatrix55 simpleTerms;
0501   simpleTerm(state.curvilinearError(), simpleTerms);
0502   //the above contains sigma(i), rho(i,j)
0503   LogDebug(theCategory + "|Adjust") << "state sigma(i), rho(i,j): \n" << simpleTerms;
0504 
0505   AlgebraicSymMatrix55 simpleValues = get(state.momentum(), false).matrix();
0506   LogDebug(theCategory + "|Adjust") << "config: (i,i), (i,j): \n" << simpleValues;
0507 
0508   for (int i = 0; i != 5; i++) {
0509     for (int j = i; j != 5; j++) {
0510       //check on each term for desired action
0511       switch (theTermAction[Pindex(i, j)]) {
0512         case scale: {
0513           simpleTerms(i, j) *= simpleValues(i, j);
0514           break;
0515         }
0516         case assign: {
0517           simpleTerms(i, j) = simpleValues(i, j);
0518           break;
0519         }
0520         case error: {
0521           edm::LogError(theCategory + "|Adjust") << " cannot properly adjust for term: " << i << "," << j;
0522         }
0523       }
0524     }
0525   }
0526   LogDebug(theCategory + "|Adjust") << "updated state sigma(i), rho(i,j): \n" << simpleTerms;
0527 
0528   AlgebraicSymMatrix55 finalTerms;
0529   complicatedTerm(simpleTerms, finalTerms);
0530   LogDebug(theCategory + "|Adjust") << "updated state COV(i,j): \n" << finalTerms;
0531 
0532   CurvilinearTrajectoryError oMat(finalTerms);
0533   state = FreeTrajectoryState(state.parameters(), oMat);
0534   LogDebug(theCategory + "|Adjust") << "updated state:\n" << state;
0535 }
0536 
0537 void MuonErrorMatrix::adjust(TrajectoryStateOnSurface &state) {
0538   AlgebraicSymMatrix55 simpleTerms;
0539   simpleTerm(state.curvilinearError(), simpleTerms);
0540   LogDebug(theCategory + "|Adjust") << "state sigma(i), rho(i,j): \n" << simpleTerms;
0541 
0542   AlgebraicSymMatrix55 simpleValues = get(state.globalMomentum(), false).matrix();
0543   LogDebug(theCategory + "|Adjust") << "config: (i,i), (i,j):\n" << simpleValues;
0544 
0545   for (int i = 0; i != 5; i++) {
0546     for (int j = i; j != 5; j++) {
0547       //check on each term for desired action
0548       switch (theTermAction[Pindex(i, j)]) {
0549         case scale: {
0550           simpleTerms(i, j) *= simpleValues(i, j);
0551           break;
0552         }
0553         case assign: {
0554           simpleTerms(i, j) = simpleValues(i, j);
0555           break;
0556         }
0557         case error: {
0558           edm::LogError(theCategory + "|Adjust") << " cannot properly adjust for term: " << i << "," << j;
0559         }
0560       }
0561     }
0562   }
0563   LogDebug(theCategory + "|Adjust") << "updated state sigma(i), rho(i,j): \n" << simpleTerms;
0564 
0565   AlgebraicSymMatrix55 finalTerms;
0566   complicatedTerm(simpleTerms, finalTerms);
0567   LogDebug(theCategory + "|Adjust") << "updated state COV(i,j): \n" << finalTerms;
0568 
0569   CurvilinearTrajectoryError oMat(finalTerms);
0570   state =
0571       TrajectoryStateOnSurface(state.weight(), state.globalParameters(), oMat, state.surface(), state.surfaceSide());
0572   LogDebug(theCategory + "|Adjust") << "updated state:\n" << state;
0573 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!