EgammaTools/src/ConversionLikelihoodCalculator.cc

0001 #include "RecoEgamma/EgammaTools/interface/ConversionLikelihoodCalculator.h"
0002
0003 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0004
0005 ConversionLikelihoodCalculator::ConversionLikelihoodCalculator() {
0006   reader_ = std::make_unique<TMVA::Reader>("!Color:Silent");
0007
0008   //   std::cout << "Init Reader()" << std::endl;
0009
0010   reader_->AddVariable("log(e_over_p)", &log_e_over_p_);
0011   reader_->AddVariable("log(abs(cot_theta))", &log_abs_cot_theta_);
0012   reader_->AddVariable("log(abs(delta_phi))", &log_abs_delta_phi_);
0013   reader_->AddVariable("log(chi2_max_pt)", &log_chi2_max_pt_);
0014   reader_->AddVariable("log(chi2_min_pt)", &log_chi2_min_pt_);
0015 }
0016
0017 void ConversionLikelihoodCalculator::setWeightsFile(const char* weightsFile) {
0018   //   std::cout << "Before BookMVA " << weightsFile << std::endl;
0019   reader_->BookMVA("Likelihood", weightsFile);
0020   //   std::cout << "After  BookMVA" << std::endl;
0021 }
0022
0023 double ConversionLikelihoodCalculator::calculateLikelihood(reco::ConversionRef conversion) {
0024   if (conversion->nTracks() != 2)
0025     return -1.;
0026
0027   log_e_over_p_ = log(conversion->EoverP());
0028
0029   log_abs_cot_theta_ = log(fabs(conversion->pairCotThetaSeparation()));
0030
0031   double delta_phi = conversion->tracks()[0]->innerMomentum().phi() - conversion->tracks()[1]->innerMomentum().phi();
0032   double pi = 3.14159265;
0033   // phi normalization
0034   while (delta_phi > pi)
0035     delta_phi -= 2 * pi;
0036   while (delta_phi < -pi)
0037     delta_phi += 2 * pi;
0038   log_abs_delta_phi_ = log(fabs(delta_phi));
0039
0040   double chi2_1 = conversion->tracks()[0]->normalizedChi2();
0041   double pt_1 = conversion->tracks()[0]->pt();
0042
0043   double chi2_2 = conversion->tracks()[1]->normalizedChi2();
0044   double pt_2 = conversion->tracks()[1]->pt();
0045
0046   double chi2_max_pt = chi2_1;
0047   double chi2_min_pt = chi2_2;
0048
0049   if (pt_2 > pt_1) {
0050     chi2_max_pt = chi2_2;
0051     chi2_min_pt = chi2_1;
0052   }
0053
0054   log_chi2_max_pt_ = log(chi2_max_pt);
0055   log_chi2_min_pt_ = log(chi2_min_pt);
0056
0057   //   std::cout << "log_e_over_p_ " << log_e_over_p_ << std::endl;
0058   //   std::cout << "log_abs_cot_theta_ " << log_abs_cot_theta_ << std::endl;
0059   //   std::cout << "log_abs_delta_phi_ " << log_abs_delta_phi_ << std::endl;
0060   //   std::cout << "log_chi2_max_pt_ " << log_chi2_max_pt_ << std::endl;
0061   //   std::cout << "log_chi2_min_pt_ " << log_chi2_min_pt_ << std::endl;
0062   std::vector<Float_t> inputVec;
0063   inputVec.push_back(log_e_over_p_);
0064   inputVec.push_back(log_abs_cot_theta_);
0065   inputVec.push_back(log_abs_delta_phi_);
0066   inputVec.push_back(log_chi2_max_pt_);
0067   inputVec.push_back(log_chi2_min_pt_);
0068   float like = reader_->EvaluateMVA(inputVec, "Likelihood");
0069   //   std::cout << "reader_->EvaluateMVA(\"Likelihood\") " << reader_->EvaluateMVA(inputVec,"Likelihood") << std::endl;
0070
0071   return like;
0072 }
0073
0074 double ConversionLikelihoodCalculator::calculateLikelihood(reco::Conversion& conversion) {
0075   if (conversion.nTracks() != 2)
0076     return -1.;
0077
0078   log_e_over_p_ = log(conversion.EoverP());
0079
0080   log_abs_cot_theta_ = log(fabs(conversion.pairCotThetaSeparation()));
0081
0082   double delta_phi = conversion.tracksPin()[0].phi() - conversion.tracksPin()[1].phi();
0083   double pi = 3.14159265;
0084   // phi normalization
0085   while (delta_phi > pi)
0086     delta_phi -= 2 * pi;
0087   while (delta_phi < -pi)
0088     delta_phi += 2 * pi;
0089   log_abs_delta_phi_ = log(fabs(delta_phi));
0090
0091   double chi2_1 = conversion.tracks()[0]->normalizedChi2();
0092   double pt_1 = conversion.tracks()[0]->pt();
0093
0094   double chi2_2 = conversion.tracks()[1]->normalizedChi2();
0095   double pt_2 = conversion.tracks()[1]->pt();
0096
0097   double chi2_max_pt = chi2_1;
0098   double chi2_min_pt = chi2_2;
0099
0100   if (pt_2 > pt_1) {
0101     chi2_max_pt = chi2_2;
0102     chi2_min_pt = chi2_1;
0103   }
0104
0105   log_chi2_max_pt_ = log(chi2_max_pt);
0106   log_chi2_min_pt_ = log(chi2_min_pt);
0107
0108   //   std::cout << "log_e_over_p_ " << log_e_over_p_ << std::endl;
0109   //   std::cout << "log_abs_cot_theta_ " << log_abs_cot_theta_ << std::endl;
0110   //   std::cout << "log_abs_delta_phi_ " << log_abs_delta_phi_ << std::endl;
0111   //   std::cout << "log_chi2_max_pt_ " << log_chi2_max_pt_ << std::endl;
0112   //   std::cout << "log_chi2_min_pt_ " << log_chi2_min_pt_ << std::endl;
0113   std::vector<Float_t> inputVec;
0114   inputVec.push_back(log_e_over_p_);
0115   inputVec.push_back(log_abs_cot_theta_);
0116   inputVec.push_back(log_abs_delta_phi_);
0117   inputVec.push_back(log_chi2_max_pt_);
0118   inputVec.push_back(log_chi2_min_pt_);
0119   float like = reader_->EvaluateMVA(inputVec, "Likelihood");
0120   //   std::cout << "reader_->EvaluateMVA(\"Likelihood\") " << reader_->EvaluateMVA(inputVec,"Likelihood") << std::endl;
0121
0122   return like;
0123 }