25 #include "TDecompSVD.h"
26 #include "TMatrixDSymEigen.h"
42 using std::setprecision;
47 void validate_reduction(
int matrixsize, TMatrixDSym& corr, std::vector<std::pair<TH1*, TH1*>>& eigenvars, TH1*
result,
const std::string&
scheme,
const std::string&
tagger,
const std::string&
wp,
const std::string& jetauthor){
55 TMatrixDSym approx_cov(matrixsize);
58 for (
const std::pair<TH1*,TH1*>&
var : eigenvars){
61 TH1* pure_var =
static_cast<TH1*
>(
var.first->Clone()); pure_var->SetDirectory(0);
62 pure_var->Add(
result, -1.0);
64 std::vector<double> pure_var_vec;
65 for (
int binz = 1 ; binz <= pure_var->GetNbinsZ() ; ++binz){
66 for (
int biny = 1 ; biny <= pure_var->GetNbinsY() ; ++biny){
67 for (
int binx = 1 ; binx <= pure_var->GetNbinsX() ; ++binx){
68 pure_var_vec.push_back(pure_var->GetBinContent(pure_var->GetBin(binx,biny,binz)));
73 for(
int i = 0 ;
i < corr.GetNrows() ;
i++){
74 for(
int j = 0 ; j < corr.GetNcols() ; j++){
75 approx_cov(
i,j) += (pure_var_vec[
i])*(pure_var_vec[j]);
80 }
else if (
scheme ==
"SFGlobalEigen"){
83 for (
const std::pair<TH1*,TH1*>&
var : eigenvars){
86 TH1* pure_var =
static_cast<TH1*
>(
var.first->Clone()); pure_var->SetDirectory(0);
87 pure_var->Add(comb_result, -1.0);
88 for(
int i = 0 ;
i < approx_cov.GetNrows() ;
i++){
89 for(
int j = 0 ; j < approx_cov.GetNcols() ; j++){
90 approx_cov(
i,j) += (pure_var->GetBinContent(
i))*(pure_var->GetBinContent(j));
96 std::cout <<
" Error: The scheme " <<
scheme <<
" didn't match SFEigen or SFGlobalEigen." << std::endl;
101 TMatrixDSym approx_corr(approx_cov);
102 for(
int row = 0 ;
row < approx_cov.GetNrows() ;
row++){
103 for(
int col = 0 ;
col < approx_cov.GetNcols() ;
col++){
104 double rowvar = sqrt(approx_cov(
row,
row));
105 double colvar = sqrt(approx_cov(
col,
col));
106 approx_corr(
row,
col) = approx_cov(
row,
col)/(rowvar * colvar);
110 TMatrixDSym numerator = (approx_corr - corr);
112 TMatrixDSym correlation_comparison(matrixsize);
113 for(
int row = 0 ;
row < approx_corr.GetNrows() ;
row++){
114 for(
int col = 0 ;
col < approx_corr.GetNcols() ;
col++){
116 correlation_comparison(
row,
col) = 0.0;
123 string outputname =
"Validate_" +
scheme +
"_" +
tagger +
"_" +
wp +
"_" + jetauthor +
"_" +
std::to_string(eigenvars.size()) +
"vars" ;
124 string eospath =
"/your/path/here/plotting/";
125 string outputfilename =
eospath + outputname +
".root";
126 std::cout <<
"ATTEMPTING TO WRITE " << outputfilename <<
" TO FILE..." << std::endl;
127 TFile
fout(outputfilename.c_str(),
"recreate");
129 std::unique_ptr<TH2D> hout(
new TH2D(correlation_comparison));
132 hout->SetTitle(
"Correlation Matrix Comparison");
148 TMatrixDSym getStatCovarianceMatrix(
const TH1*
hist,
bool zero) {
149 Int_t nbinx =
hist->GetNbinsX()+2, nbiny =
hist->GetNbinsY()+2, nbinz =
hist->GetNbinsZ()+2;
151 if (
hist->GetDimension() > 1)
rows *= nbiny;
152 if (
hist->GetDimension() > 2)
rows *= nbinz;
160 for (Int_t binx = 1; binx < nbinx-1; ++binx)
161 for (Int_t biny = 1; biny < nbiny-1; ++biny)
162 for (Int_t binz = 1; binz < nbinz-1; ++binz) {
163 Int_t
bin =
hist->GetBin(binx, biny, binz);
175 TMatrixDSym getSystCovarianceMatrix(
const TH1*
ref,
const TH1* unc,
bool doCorrelated,
int tagWeightAxis) {
176 Int_t nbinx =
ref->GetNbinsX()+2, nbiny =
ref->GetNbinsY()+2, nbinz =
ref->GetNbinsZ()+2;
178 if (
ref->GetDimension() > 1)
rows *= nbiny;
179 if (
ref->GetDimension() > 2)
rows *= nbinz;
184 if (unc->GetNbinsX()+2 != nbinx || unc->GetNbinsY()+2 != nbiny || unc->GetNbinsZ()+2 != nbinz || unc->GetDimension() !=
ref->GetDimension()) {
185 std::cout <<
"getSystCovarianceMatrix: inconsistency found in histograms " <<
ref->GetName() <<
" and " << unc->GetName() << std::endl;
196 if (! doCorrelated) {
197 if (tagWeightAxis < 0) {
199 for (Int_t binx = 1; binx < nbinx-1; ++binx)
200 for (Int_t biny = 1; biny < nbiny-1; ++biny)
201 for (Int_t binz = 1; binz < nbinz-1; ++binz) {
202 Int_t
bin =
ref->GetBin(binx, biny, binz);
203 double err = unc->GetBinContent(
bin);
207 }
else if (tagWeightAxis == 0) {
209 for (Int_t biny = 1; biny < nbiny-1; ++biny)
210 for (Int_t binz = 1; binz < nbinz-1; ++binz)
211 for (Int_t binx = 1; binx < nbinx-1; ++binx) {
212 Int_t
bin =
ref->GetBin(binx, biny, binz);
213 double err = unc->GetBinContent(
bin);
214 for (Int_t binx2 = 1; binx2 < nbinx-1; ++binx2) {
215 Int_t
bin2 =
ref->GetBin(binx2, biny, binz);
216 double err2 = unc->GetBinContent(
bin2);
221 }
else if (tagWeightAxis == 1) {
223 for (Int_t binx = 1; binx < nbinx-1; ++binx)
224 for (Int_t binz = 1; binz < nbinz-1; ++binz)
225 for (Int_t biny = 1; biny < nbiny-1; ++biny) {
226 Int_t
bin =
ref->GetBin(binx, biny, binz);
227 double err = unc->GetBinContent(
bin);
228 for (Int_t biny2 = 1; biny2 < nbiny-1; ++biny2) {
229 Int_t
bin2 =
ref->GetBin(binx, biny2, binz);
230 double err2 = unc->GetBinContent(
bin2);
235 }
else if (tagWeightAxis == 2) {
237 for (Int_t binx = 1; binx < nbinx-1; ++binx)
238 for (Int_t biny = 1; biny < nbiny-1; ++biny)
239 for (Int_t binz = 1; binz < nbinz-1; ++binz) {
240 Int_t
bin =
ref->GetBin(binx, biny, binz);
241 double err = unc->GetBinContent(
bin);
242 for (Int_t binz2 = 1; binz2 < nbinz-1; ++binz2) {
243 Int_t
bin2 =
ref->GetBin(binx, biny, binz2);
244 double err2 = unc->GetBinContent(
bin2);
252 for (Int_t binx = 1; binx < nbinx-1; ++binx)
253 for (Int_t biny = 1; biny < nbiny-1; ++biny)
254 for (Int_t binz = 1; binz < nbinz-1; ++binz) {
255 Int_t
bin =
ref->GetBin(binx, biny, binz);
256 double err = unc->GetBinContent(
bin);
257 for (Int_t binx2 = 1; binx2 < nbinx-1; ++binx2)
258 for (Int_t biny2 = 1; biny2 < nbiny-1; ++biny2)
259 for (Int_t binz2 = 1; binz2 < nbinz-1; ++binz2) {
260 Int_t
bin2 =
ref->GetBin(binx2, biny2, binz2);
261 double err2 = unc->GetBinContent(
bin2);
271 int rows = unc.size();
274 for(
int i = 0 ;
i <
rows ;
i++){
275 for(
int j = 0 ; j <
rows ; j++){
276 double err1 = unc[
i];
277 double err2 = unc[j];
278 cov(
i,j) = err1*err2;
285 TMatrixDSym getStatCovarianceMatrix(
const std::vector<double>& unc,
bool zero) {
287 int rows = unc.size();
289 for (
int i = 0;
i <
rows;
i++) {
331 CalibrationDataEigenVariations::CalibrationDataEigenVariations(
const std::string& cdipath,
const std::string&
tagger,
const std::string&
wp,
const std::string& jetcollection, CalibrationDataHistogramContainer*
cnt,
bool excludeRecommendedUncertaintySet,
bool base) :
332 m_cnt(
cnt), m_initialized(false), m_validate(false), m_namedExtrapolation(-1), m_statVariations(false), m_cdipath(cdipath), m_taggername(
tagger), m_wp(
wp), m_jetauthor(jetcollection), m_totalvariance(0), m_capturedvariance(0), m_verbose(false)
335 if (m_verbose) std::cout <<
" CDEV Constructor : info " << cdipath <<
" " <<
tagger <<
" " <<
wp <<
" " << jetcollection << std::endl;
337 if (excludeRecommendedUncertaintySet &&
base) {
338 std::vector<std::string> to_exclude =
split(m_cnt->getExcludedUncertainties());
339 for (
const auto&
name : to_exclude) {
340 excludeNamedUncertainty(
name,
const_cast<CalibrationDataHistogramContainer*
>(
cnt));
342 if (to_exclude.size() == 0) {
343 std::cerr <<
"CalibrationDataEigenVariations warning: exclusion of pre-set uncertainties list requested but no (or empty) list found" <<std::endl;
348 for (
const auto&
name : uncs) {
349 if (
name.find(
"stat_np") != string::npos) m_statVariations =
true;
355 CalibrationDataEigenVariations::~CalibrationDataEigenVariations()
382 std::cerr <<
"CalibrationDataEigenVariations::excludeNamedUncertainty error:" <<
" initialization already done" << std::endl;
383 }
else if (
name ==
"comment" ||
name ==
"result" ||
name ==
"systematics" ||
384 name ==
"statistics" ||
name ==
"combined" ||
name ==
"extrapolation" ||
385 name ==
"MCreference" ||
name ==
"MChadronisation" ||
name ==
"ReducedSets" ||
388 std::cerr <<
"CalibrationDataEigenVariations::excludeNamedUncertainty error:" <<
" name " <<
name <<
" not allowed" << std::endl;
391 else if (
name.back() ==
'*'){
392 std::string temp_name =
name.substr(0,
name.size()-1);
394 std::vector<std::string> unc_subgroup;
395 std::copy_if(uncs.begin(), uncs.end(), back_inserter(unc_subgroup),
396 [&temp_name](
const std::string&
el) {
397 return el.compare(0, temp_name.size(), temp_name) == 0;
399 if (
m_verbose) std::cout <<
"Found a group of uncertainties to exclude: " <<
name <<
" found " <<unc_subgroup.size() <<
" uncertainties corresponding to the query" <<std::endl;
400 for (
const auto& single_name : unc_subgroup){
403 if (
m_verbose) std::cout <<
"Name : " << single_name << std::endl;
404 m_named.push_back(std::pair<TH1*, TH1*>(0, 0));
409 else if (!
cnt->GetValue(
name.c_str())){
410 std::cerr <<
"CalibrationDataEigenVariations::excludeNamedUncertainty error:" <<
" uncertainty named " <<
name <<
" not found" << std::endl;
433 TH1*
result =
dynamic_cast<TH1*
>(
m_cnt->GetValue(
"result"));
435 std::cerr<<
"CalibrationDataEigenVariations::getEigenCovarianceMatrix(): failed dynamic cast to TH1*\n";
436 return TMatrixDSym();
444 TMatrixDSym* sCov =
dynamic_cast<TMatrixDSym*
>(
m_cnt->GetValue(
"statistics"));
452 for (
unsigned int t = 0;
t < uncs.size(); ++
t) {
454 if (uncs[
t] ==
"comment" || uncs[
t] ==
"result" || uncs[
t] ==
"combined" ||
455 uncs[
t] ==
"statistics" || uncs[
t] ==
"systematics" || uncs[
t] ==
"MCreference" ||
456 uncs[
t] ==
"MChadronisation" || uncs[
t] ==
"extrapolation" || uncs[
t] ==
"ReducedSets" ||
457 uncs[
t] ==
"excluded_set")
continue;
460 TH1* hunc =
dynamic_cast<TH1*
>(
m_cnt->GetValue(uncs[
t].c_str()));
462 std::cerr<<
"CalibrationDataEigenVariations::getEigenCovarianceMatrix(): dynamic cast failed\n";
479 TH1 *
result =
dynamic_cast<TH1*
>(
m_cnt->GetValue(
"result"));
481 int nbins = jac.GetNcols();
483 auto variation = std::make_unique<double[]>(
nbins);
485 for (
const std::pair<TH1*, TH1*>&
it :
m_eigen){
486 TH1* resultVariedUp =
it.first;
487 for (
unsigned int u = 0;
u < (
unsigned int)
nbins; ++
u) variation[
u] = resultVariedUp->GetBinContent(
u) -
result->GetBinContent(
u);
490 cov(
u,
v) += variation[
u]*variation[
v];
510 TH1*
result =
dynamic_cast<TH1*
>(
m_cnt->GetValue(
"result"));
512 std::cerr<<
"CalibrationDataEigenVariations::getJacobianReductionMatrix(): dynamic cast failed\n";
525 int ndim =
result->GetDimension();
531 std::vector<int> zeroComponents;
533 std::cerr <<
" error: covariance matrix size (" <<
cov.GetNrows() <<
") doesn't match histogram size (" <<
nbins <<
")" << std::endl;
534 return TMatrixDSym();
540 Int_t binx, biny, binz;
541 result->GetBinXYZ(
i, binx, biny, binz);
542 if ((binx == 0 || binx ==
result->GetNbinsX()+1) ||
543 (ndim > 1 && (biny == 0 || biny ==
result->GetNbinsY()+1)) ||
544 (ndim > 2 && (binz == 0 || binz ==
result->GetNbinsZ()+1)))
547 zeroComponents.push_back(
i);
551 else if (fabs(
cov(
i,0)) < 1
e-10) {
552 bool isThereANonZero(
false);
553 for (
int j = 0; j <
nbins; ++j) {
554 if (fabs(
cov(
i,j)) > 1
e-10) {
555 isThereANonZero=
true;
break;
558 if (! isThereANonZero) {
560 zeroComponents.push_back(
i);
608 if (nZeros >=
nbins) {
609 std::cerr <<
" Problem. Found n. " << nZeros <<
" while size of matrix is " <<
nbins << std::endl;
610 return TMatrixDSym();
615 TMatrixT<double> matrixVariationJacobian(
size,
nbins);
619 for (
unsigned int s = 0;
s < zeroComponents.size(); ++
s) {
620 if (zeroComponents.at(
s) ==
i) {
629 matrixVariationJacobian(
i-nMissed,
i)=1;
632 return matrixVariationJacobian;
646 TH1*
result =
dynamic_cast<TH1*
>(
m_cnt->GetValue(
"result"));
648 std::cerr<<
"CalibrationDataEigenVariations::initialize: dynamic cast failed\n";
664 corr(rowc,colc) =
cov(
row,
col)/(rowvar * colvar);
675 TH1* hunc = (TH1*)
m_cnt->GetValue(
it->first.c_str());
676 TString namedvar(
"namedVar");
677 namedvar +=
it->first.c_str();
678 TString namedvarUp(namedvar); namedvarUp +=
"_up";
679 TString namedvarDown(namedvar); namedvarDown +=
"_down";
680 TH1* resultVariedUp = (TH1*)
result->Clone(namedvarUp); resultVariedUp->Add(hunc, 1.0); resultVariedUp->SetDirectory(0);
681 TH1* resultVariedDown = (TH1*)
result->Clone(namedvarDown); resultVariedDown->Add(hunc, -1.0); resultVariedDown->SetDirectory(0);
682 p.first = resultVariedUp;
683 p.second = resultVariedDown;
687 if (TH1* hunc = (TH1*)
m_cnt->GetValue(
"extrapolation")) {
688 TH1* resultVariedUp = (TH1*) hunc->Clone(
"extrapolation_up"); resultVariedUp->SetDirectory(0);
689 TH1* resultVariedDown = (TH1*) hunc->Clone(
"extrapolation_down"); resultVariedDown->SetDirectory(0);
690 Int_t nbinx = hunc->GetNbinsX()+2, nbiny = hunc->GetNbinsY()+2, nbinz = hunc->GetNbinsZ()+2;
691 Int_t firstbinx = hunc->GetXaxis()->FindFixBin(
result->GetXaxis()->GetBinCenter(1));
692 Int_t firstbiny =
result->GetDimension() > 1 ? hunc->GetYaxis()->FindFixBin(
result->GetYaxis()->GetBinCenter(1)) : 1;
693 Int_t firstbinz =
result->GetDimension() > 2 ? hunc->GetZaxis()->FindFixBin(
result->GetZaxis()->GetBinCenter(1)) : 1;
694 for (Int_t binx = 1; binx < nbinx-1; ++binx) {
695 Int_t binxResult = binx - firstbinx + 1;
696 if (binxResult < 1) binxResult = 1;
697 if (binxResult >
result->GetNbinsX()) binxResult =
result->GetNbinsX();
698 for (Int_t biny = 1; biny < nbiny-1; ++biny) {
699 Int_t binyResult = biny - firstbiny + 1;
700 if (binyResult >
result->GetNbinsY()) binyResult =
result->GetNbinsY();
701 for (Int_t binz = 1; binz < nbinz-1; ++binz) {
702 Int_t binzResult = binz - firstbinz + 1;
703 if (binzResult >
result->GetNbinsZ()) binzResult =
result->GetNbinsZ();
704 Int_t
bin = hunc->GetBin(binx, biny, binz);
705 double contResult =
result->GetBinContent(binxResult, binyResult, binzResult);
706 resultVariedUp->SetBinContent(
bin, contResult + hunc->GetBinContent(
bin));
707 resultVariedDown->SetBinContent(
bin, contResult - hunc->GetBinError(
bin));
711 m_named.push_back(std::make_pair(resultVariedUp, resultVariedDown));
717 int ndim =
result->GetDimension();
747 int size = matrixVariationJacobian.GetNrows();
751 const TMatrixDSym matrixCovariance =
cov.Similarity(matrixVariationJacobian);
754 TMatrixDSymEigen eigenValueMaker (matrixCovariance);
755 TVectorT<double> eigenValues = eigenValueMaker.GetEigenValues();
756 TMatrixT<double> eigenVectors = eigenValueMaker.GetEigenVectors();
757 TMatrixT<double> matrixVariations (
size,
size);
762 for (
int i = 0;
i <
size; ++
i) {
766 matrixVariations(
i,
r) = -1.0*eigenVectors[
r][
i]*sqrt(fabs(eigenValues[
i]));
769 TMatrixT<double> matrixVariationsWithZeros = matrixVariations * matrixVariationJacobian;
772 for (
int i = 0;
i < matrixVariationsWithZeros.GetNrows(); ++
i) {
775 TString superstring(
"eigenVar");
778 TString nameUp(superstring); nameUp +=
"_up";
779 TString nameDown(superstring); nameDown +=
"_down";
782 TH1* resultVariedUp = (TH1*)
result->Clone(nameUp); resultVariedUp->SetDirectory(0);
783 TH1* resultVariedDown = (TH1*)
result->Clone(nameDown); resultVariedDown->SetDirectory(0);
786 resultVariedUp->SetBinContent(
u,
result->GetBinContent(
u) + matrixVariationsWithZeros(
i,
u));
787 resultVariedDown->SetBinContent(
u,
result->GetBinContent(
u) - matrixVariationsWithZeros(
i,
u));
790 m_eigen.push_back(std::make_pair(resultVariedUp, resultVariedDown));
796 size_t current_set = 0;
798 bool keep_variation =
false;
799 TH1 *
up =
static_cast<TH1*
>(
m_eigen[
index].first->Clone());
up->SetDirectory(0);
804 if (fabs(
up->GetBinContent(
bin)) > min_variance) {
805 keep_variation =
true;
810 if (!keep_variation){
811 final_set.insert(current_set);
817 if (final_set.size() > 0)
818 if (
m_verbose) std::cout <<
"CalibrationDataEigenVariations: Removing " << final_set.size() <<
" eigenvector variations leading to sub-tolerance effects, retaining " <<
m_eigen.size()-final_set.size() <<
" variations" << std::endl;
837 if (
set.size() == 0)
return;
839 std::vector<std::pair<TH1*, TH1*> > new_eigen;
854 for (
IndexSet::iterator subset_it = set_it->begin(); subset_it != set_it->end(); ++subset_it)
855 simple_set.insert(*subset_it);
871 simple_set.insert(
it);
891 for (
IndexSet::iterator subset_it = set_it->begin(); subset_it != set_it->end(); ++subset_it){
892 if (
checker.count(*subset_it) == 0 && *subset_it <=
m_eigen.size())
895 std::cerr <<
"Error in CalibrationDataEigenVariations::mergeVariations: \
896 IndexSets must not overlap and must lie between 1 and " <<
m_eigen.size() <<
". Aborting!" << std::endl;
903 TH1 *
result =
static_cast<TH1*
>(
m_cnt->GetValue(
"result"));
906 int ndim =
result->GetDimension();
918 if (set_it->empty())
continue;
920 double sum_H_up = 0.0, sum_H_down = 0.0;
921 size_t lowest_index = *set_it->lower_bound(0);
922 TH1 *total_var_up =
static_cast<TH1*
>(
m_eigen[lowest_index].first->Clone()),
923 *total_var_down =
static_cast<TH1*
>(
m_eigen[lowest_index].
second->Clone());
924 total_var_up->SetDirectory(0);
925 total_var_down->SetDirectory(0);
927 total_var_up->Reset();
928 total_var_down->Reset();
932 subset_it != set_it->end(); ++subset_it) {
933 size_t actual_index = *subset_it;
935 if (actual_index != lowest_index) toDelete.insert(*subset_it);
937 TH1 *partial_var_up =
static_cast<TH1*
>(
m_eigen[actual_index].first->Clone()),
938 *partial_var_down =
static_cast<TH1*
>(
m_eigen[actual_index].
second->Clone());
939 partial_var_up->SetDirectory(0);
940 partial_var_down->SetDirectory(0);
942 partial_var_up->Add(
result, -1.0);
943 partial_var_down->Add(
result, -1.0);
945 partial_var_down->SetBinContent(
i, -1.0*partial_var_down->GetBinContent(
i));
949 double sum_up = total_var_up->GetBinContent(
u),
950 sum_down = total_var_down->GetBinContent(
u);
952 sum_up += partial_var_up->GetBinContent(
u)*partial_var_up->GetBinContent(
v);
953 sum_H_up += partial_var_up->GetBinContent(
u)*partial_var_up->GetBinContent(
v);
954 sum_down += partial_var_down->GetBinContent(
u)*partial_var_down->GetBinContent(
v);
955 sum_H_down += partial_var_down->GetBinContent(
u)*partial_var_down->GetBinContent(
v);
957 total_var_up->SetBinContent(
u, sum_up);
958 total_var_down->SetBinContent(
u, sum_down);
960 delete partial_var_up;
961 delete partial_var_down;
967 total_var_up->SetBinContent(
i, total_var_up->GetBinContent(
i)/sqrt(sum_H_up));
969 total_var_up->SetBinContent(
i, 0.0);
970 if (sum_H_down != 0.0)
971 total_var_down->SetBinContent(
i, -1.0*total_var_down->GetBinContent(
i)/sqrt(sum_H_down));
973 total_var_down->SetBinContent(
i, 0.0);
976 total_var_up->Add(
result);
977 total_var_down->Add(
result);
979 m_eigen[lowest_index].first = total_var_up;
980 m_eigen[lowest_index].second = total_var_down;
1019 TH1*&
up, TH1*& down)
1031 if (variation <
m_eigen.size()) {
1033 down =
m_eigen[variation].second;
1044 TH1*&
up, TH1*& down)
1064 TH1*&
up, TH1*& down)
1076 if (nameIndex <
m_named.size()) {
1078 down =
m_named[nameIndex].second;
1110 std::map<std::string, std::map<std::string, float>> &coefficientMap)
1121 std::vector<TH1*> originSF_hvec;
1122 std::vector<TH1*> eigenSF_hvec;
1126 std::vector<string> uncList;
1127 for (
unsigned int t = 0;
t < fullUncList.size(); ++
t) {
1129 if (fullUncList[
t] ==
"comment" || fullUncList[
t] ==
"result" || fullUncList[
t] ==
"combined" ||
1130 fullUncList[
t] ==
"statistics" || fullUncList[
t] ==
"systematics" || fullUncList[
t] ==
"MCreference" ||
1131 fullUncList[
t] ==
"MChadronisation" || fullUncList[
t] ==
"extrapolation" || fullUncList[
t] ==
"ReducedSets" ||
1132 fullUncList[
t] ==
"excluded_set")
continue;
1137 TH1* hunc =
dynamic_cast<TH1*
>(
m_cnt->GetValue(fullUncList[
t].c_str()));
1139 std::cerr<<
"CalibrationDataEigenVariations::EigenVectorRecomposition: dynamic cast failed\n";
1143 Int_t nx = hunc->GetNbinsX();
1144 Int_t ny = hunc->GetNbinsY();
1145 Int_t nz = hunc->GetNbinsZ();
1147 bool retain =
false;
1151 for(Int_t binx = 1; binx <= nx; binx++)
1152 for(Int_t biny = 1; biny <= ny; biny++)
1153 for(Int_t binz = 1; binz <= nz; binz++){
1154 bin = hunc->GetBin(binx, biny, binz);
1155 if (fabs(hunc->GetBinContent(
bin)) > 1
E-20){
1161 if (
m_verbose) std::cout<<
"Eigenvector Recomposition: Empty uncertainty "<<fullUncList.at(
t)<<
" is discarded."<<std::endl;
1165 uncList.push_back(fullUncList.at(
t));
1166 originSF_hvec.push_back(hunc);
1169 TH1*
nom =
dynamic_cast<TH1*
>(
m_cnt->GetValue(
"result"));
1171 if (
m_verbose) std::cout<<
"Eigenvector Recomposition: dynamic cast failed\n";
1174 int dim =
nom->GetDimension();
1175 Int_t nx =
nom->GetNbinsX();
1176 Int_t ny =
nom->GetNbinsY();
1177 Int_t nz =
nom->GetNbinsZ();
1181 TMatrixD matSF(uncList.size(),
nbins);
1184 for(
unsigned int i = 0;
i < uncList.size();
i++){
1187 for(
int binz = 1; binz <= nz; binz++)
1188 for(
int biny = 1; biny <= ny; biny++)
1189 for(
int binx = 1; binx <= nx; binx++){
1190 Int_t
bin = originSF_hvec.at(
i)->GetBin(binx, biny, binz);
1191 TMatrixDRow(matSF,
i)[
col] = originSF_hvec[
i]->GetBinContent(
bin);
1199 TH1* down =
nullptr;
1200 for (
int i = 0;
i < nEigen;
i++){
1202 std::cerr<<
"EigenVectorRecomposition: Error on retrieving eigenvector "<<
i<<std::endl;
1207 eigenSF_hvec.push_back(
up);
1209 TMatrixD matEigen(nEigen,
nbins);
1212 for(
int i = 0;
i < nEigen;
i++){
1215 for(
int binz = 1; binz <= nz; binz++)
1216 for(
int biny = 1; biny <= ny; biny++)
1217 for(
int binx = 1; binx <= nx; binx++){
1218 Int_t
bin = eigenSF_hvec.at(
i)->GetBin(binx, biny, binz);
1219 TMatrixDRow(matEigen,
i)[
col] = eigenSF_hvec[
i]->GetBinContent(
bin);
1225 TMatrixD matEigenOriginal = matEigen;
1226 TMatrixD matEigenTranspose = matEigen;
1227 matEigenTranspose = matEigenTranspose.T();
1228 TMatrixD matOriginalTimesTranspose = matEigenOriginal*matEigenTranspose;
1229 TMatrixD matEigenInvert = matEigenTranspose*matOriginalTimesTranspose.Invert();
1232 TMatrixD matCoeff = matSF*matEigenInvert;
1233 int nRows = matCoeff.GetNrows();
1234 int nCols = matCoeff.GetNcols();
1235 std::map<std::string, float> temp_map;
1239 temp_map[uncList[
row]] = TMatrixDRow(matCoeff,
row)[
col];
1295 TFile*
f = TFile::Open(
fileName.Data(),
"READ");
1302 TString contName(
dir);
1304 f->GetObject(contName.Data(),
c);
1306 if (
m_verbose) std::cout <<
"Found " << contName.Data() << std::endl;
1308 std::vector<std::string> uncs =
c->listUncertainties();
1309 TH1*
result =
dynamic_cast<TH1*
>(
c->GetValue(
"result"));
1311 if (
m_verbose) std::cout <<
"Dynamic cast failed at "<<__LINE__<<
"\n";
1316 for (
const std::string& unc : uncs){
1318 if ((unc==
"result")||(unc==
"comment")||(unc==
"ReducedSets")||(unc==
"systematics")||(unc==
"statistics")||(unc==
"extrapolation")||(unc==
"MChadronisation")||(unc==
"combined")||(unc==
"extrapolation from charm")) {
1326 if (excludeRecommendedUncertaintySet) {
1327 std::vector<std::string> to_exclude =
split(
c->getExcludedUncertainties());
1328 if (to_exclude.size() == 0) {
1329 std::cerr <<
"CalibrationDataEigenVariations warning: exclusion of pre-set uncertainties list requested but no (or empty) list found" << std::endl;
1331 for (
const auto&
name : to_exclude) {
1332 if (
name ==
"")
continue;
1343 if (
m_verbose) std::cout <<
"Printing out all shared uncertainties for " <<
tagger <<
"/" << jetcollection <<
"/" <<
wp << std::endl;
1344 if (
m_verbose) std::cout <<
"| " << std::endl;
1346 if (
m_verbose) std::cout <<
"|-- " << (*it) << std::endl;
1350 if (
m_verbose) std::cout <<
"| no shared systematics between ";
1353 }
if (
m_verbose) std::cout << std::endl;
1393 std::map<std::string, TMatrixDSym> cov_matrices;
1397 std::vector<int> flavs_in_common;
1399 std::vector<double> comb_syst;
1404 int flavour_size = 0;
1407 TH1* hunc =
dynamic_cast<TH1*
>(
c->GetValue(tunc.Data()));
1408 TH1*
ref =
dynamic_cast<TH1*
>(
c->GetValue(
"result"));
1410 if (
m_verbose) std::cout <<
" There was no uncertainty OR SF/EFF results... Are you sure you have the right CDIContainer path?" << std::endl;
1413 int tagweightax =
c->getTagWeightAxis();
1415 flavs_in_common.push_back(1);
1416 Int_t nbinx = hunc->GetNbinsX(), nbiny = hunc->GetNbinsY(), nbinz = hunc->GetNbinsZ();
1418 if (hunc->GetDimension() > 1)
rows *= nbiny;
1419 if (hunc->GetDimension() > 2)
rows *= nbinz;
1420 flavour_size =
rows;
1422 flavs_in_common.push_back(0);
1424 Int_t nbinx =
ref->GetNbinsX(), nbiny =
ref->GetNbinsY(), nbinz =
ref->GetNbinsZ();
1426 if (
ref->GetDimension() > 1)
rows *= nbiny;
1427 if (
ref->GetDimension() > 2)
rows *= nbinz;
1428 flavour_size =
rows;
1432 if (tagweightax == -1){
1433 for (
int i = 1 ;
i <= flavour_size ;
i++){
1435 Int_t
bin = hunc->GetBin(1,
i,1);
1436 double unc_val = hunc->GetBinContent(
bin);
1437 comb_syst.push_back(unc_val);
1439 comb_syst.push_back(0.0);
1442 }
else if (tagweightax == 0){
1445 for (Int_t zbins = 1 ; zbins <=
ref->GetNbinsZ() ; zbins++){
1449 double unc_val = hunc->GetBinContent(
bin);
1450 comb_syst.push_back(unc_val);
1452 comb_syst.push_back(0.0);
1458 }
else if (tagweightax == 1){
1461 for (Int_t zbins = 1 ; zbins <=
ref->GetNbinsZ() ; zbins++){
1465 double unc_val = hunc->GetBinContent(
bin);
1466 comb_syst.push_back(unc_val);
1468 comb_syst.push_back(0.0);
1474 }
else if (tagweightax == 2){
1476 for (Int_t zbins = 1 ; zbins <=
ref->GetNbinsZ() ; zbins++){
1481 double unc_val = hunc->GetBinContent(
bin);
1482 comb_syst.push_back(unc_val);
1484 comb_syst.push_back(0.0);
1494 TMatrixDSym unc_cov(comb_syst.size());
1495 if (unc ==
"statistics"){
1498 unc_cov = getSystCovarianceMatrix(comb_syst);
1500 cov_matrices.insert({unc, unc_cov});
1503 if (flavs_in_common[0]+flavs_in_common[1]+flavs_in_common[2]+flavs_in_common[3] > 1){
1507 global_covariance += unc_cov;
1509 return global_covariance;
1525 if (
m_initialized) std::cerr <<
"CalibrationDataGlobalEigenVariations::excludeNamedUncertainty error:" <<
" initialization already done" << std::endl;
1527 else if (
name ==
"comment" ||
name ==
"result" ||
name ==
"systematics" ||
1528 name ==
"statistics" ||
name ==
"combined" ||
name ==
"extrapolation" ||
1529 name ==
"MCreference" ||
name ==
"MChadronisation" ||
name ==
"ReducedSets" ||
1531 std::cerr <<
"CalibrationDataGlobalEigenVariations::excludeNamedUncertainty error:" <<
" name [" <<
name <<
"] not allowed" << std::endl;
1536 m_flav_named[flavour].push_back(std::pair<TH1*, TH1*>(0, 0));
1543 std::vector<std::string>
1548 std::vector<std::string>
names;
1550 names.push_back(namedar.first);
1570 TMatrixDSym corr(
cov);
1584 std::vector<double> combined_result;
1585 std::map<std::string, int> flav_bins;
1592 TH1*
result =
dynamic_cast<TH1*
>(
c->GetValue(
"result"));
1594 std::cerr<<
"CalibrationDataGlobalEigenVariations::initialize: dynamic cast failed\n ";
1598 if (
c->getTagWeightAxis() == -1){
1599 flav_bins[flavour] =
result->GetNbinsY();
1600 if (
m_verbose) std::cout <<
"flav_bins["<<flavour<<
"] = " << flav_bins[flavour] << std::endl;
1601 for(
int i = 0 ;
i < flav_bins[flavour] ;
i++){
1604 double res_value =
result->GetBinContent(
bin);
1605 combined_result.push_back(res_value);
1607 }
else if (
c->getTagWeightAxis() == 0) {
1608 flav_bins[flavour] =
result->GetNbinsX()*
result->GetNbinsY();
1609 int tagbins =
result->GetNbinsX();
1610 int ptbins =
result->GetNbinsY();
1611 if (
m_verbose) std::cout <<
"flav_bins["<<flavour<<
"] = " << flav_bins[flavour] << std::endl;
1612 for(
int i = 0 ;
i < tagbins ;
i++){
1613 for(
int j = 0 ; j < ptbins ; j++){
1616 double res_value =
result->GetBinContent(
bin);
1617 combined_result.push_back(res_value);
1620 }
else if (
c->getTagWeightAxis() == 1) {
1621 flav_bins[flavour] =
result->GetNbinsX()*
result->GetNbinsY();
1622 int tagbins =
result->GetNbinsY();
1623 int ptbins =
result->GetNbinsX();
1624 if (
m_verbose) std::cout <<
"flav_bins["<<flavour<<
"] = " << flav_bins[flavour] << std::endl;
1625 for(
int i = 0 ;
i < tagbins ;
i++){
1626 for(
int j = 0 ; j < ptbins ; j++){
1629 double res_value =
result->GetBinContent(
bin);
1630 combined_result.push_back(res_value);
1633 }
else if (
c->getTagWeightAxis() == 2) {
1634 flav_bins[flavour] =
result->GetNbinsX()*
result->GetNbinsZ();
1635 int tagbins =
result->GetNbinsZ();
1636 int ptbins =
result->GetNbinsX();
1637 if (
m_verbose) std::cout <<
"flav_bins["<<flavour<<
"] = " << flav_bins[flavour] << std::endl;
1638 for(
int i = 0 ;
i < tagbins ;
i++){
1639 for(
int j = 0 ; j < ptbins ; j++){
1642 double res_value =
result->GetBinContent(
bin);
1643 combined_result.push_back(res_value);
1652 TH1* hunc = (TH1*)
c->GetValue(
it->first.c_str());
1656 TString namedvar(
"namedVar");
1657 namedvar +=
it->first.c_str();
1658 TString namedvarUp(namedvar); namedvarUp +=
"_up";
1659 TString namedvarDown(namedvar); namedvarDown +=
"_down";
1660 TH1* resultVariedUp = (TH1*)
result->Clone(namedvarUp);
1661 TH1* resultVariedDown = (TH1*)
result->Clone(namedvarDown);
1663 resultVariedUp->Add(hunc, 1.0); resultVariedUp->SetDirectory(0);
1664 resultVariedDown->Add(hunc, -1.0); resultVariedDown->SetDirectory(0);
1666 resultVariedUp->SetDirectory(0);
1667 resultVariedDown->SetDirectory(0);
1669 p.first = resultVariedUp;
1670 p.second = resultVariedDown;
1676 if (TH1* hunc = (TH1*)
c->GetValue(
"extrapolation")) {
1677 TH1* resultVariedUp = (TH1*) hunc->Clone(
"extrapolation_up"); resultVariedUp->SetDirectory(0);
1678 TH1* resultVariedDown = (TH1*) hunc->Clone(
"extrapolation_down"); resultVariedDown->SetDirectory(0);
1679 Int_t nbinx = hunc->GetNbinsX()+2, nbiny = hunc->GetNbinsY()+2, nbinz = hunc->GetNbinsZ()+2;
1680 Int_t firstbinx = hunc->GetXaxis()->FindFixBin(
result->GetXaxis()->GetBinCenter(1));
1681 Int_t firstbiny =
result->GetDimension() > 1 ? hunc->GetYaxis()->FindFixBin(
result->GetYaxis()->GetBinCenter(1)) : 1;
1682 Int_t firstbinz =
result->GetDimension() > 2 ? hunc->GetZaxis()->FindFixBin(
result->GetZaxis()->GetBinCenter(1)) : 1;
1683 for (Int_t binx = 1; binx < nbinx-1; ++binx) {
1684 Int_t binxResult = binx - firstbinx + 1;
1685 if (binxResult < 1) binxResult = 1;
1686 if (binxResult >
result->GetNbinsX()) binxResult =
result->GetNbinsX();
1687 for (Int_t biny = 1; biny < nbiny-1; ++biny) {
1688 Int_t binyResult = biny - firstbiny + 1;
1689 if (binyResult >
result->GetNbinsY()) binyResult =
result->GetNbinsY();
1690 for (Int_t binz = 1; binz < nbinz-1; ++binz) {
1691 Int_t binzResult = binz - firstbinz + 1;
1692 if (binzResult >
result->GetNbinsZ()) binzResult =
result->GetNbinsZ();
1693 Int_t
bin = hunc->GetBin(binx, biny, binz);
1694 double contResult =
result->GetBinContent(binxResult, binyResult, binzResult);
1695 resultVariedUp->SetBinContent(
bin, contResult + hunc->GetBinContent(
bin));
1696 resultVariedDown->SetBinContent(
bin, contResult - hunc->GetBinError(
bin));
1700 m_flav_named[flavour].push_back(std::make_pair(resultVariedUp, resultVariedDown));
1711 std::unique_ptr<TH1> comb_result(
new TH1D(
"combined_result",
"", combined_result.size(), 0., 1.));
1712 int nbins = comb_result->GetNbinsX()+2;
1713 int ndim = comb_result->GetDimension();
1714 if (ndim > 1)
nbins*= (comb_result->GetNbinsY()+2);
1715 if (ndim > 2)
nbins*= (comb_result->GetNbinsZ()+2);
1718 for (
unsigned int i=0 ;
i<combined_result.size() ;
i++){
1720 comb_result->SetBinContent(
i+1, combined_result[
i]);
1726 int size = matrixVariationJacobian.GetNrows();
1729 const TMatrixDSym matrixCovariance =
cov.Similarity(matrixVariationJacobian);
1732 TMatrixDSymEigen eigenValueMaker (matrixCovariance);
1733 TVectorT<double> eigenValues = eigenValueMaker.GetEigenValues();
1734 TMatrixT<double> eigenVectors = eigenValueMaker.GetEigenVectors();
1735 TMatrixT<double> matrixVariations (
size,
size);
1740 for (
int i = 0;
i <
size; ++
i) {
1741 for (
int r = 0;
r <
size; ++
r) {
1743 matrixVariations(
i,
r) = -1.0*eigenVectors[
r][
i]*sqrt(fabs(eigenValues[
i]));
1747 TMatrixT<double> matrixVariationsWithZeros = matrixVariations * matrixVariationJacobian;
1750 for (
int i = 0;
i < matrixVariationsWithZeros.GetNrows(); ++
i) {
1751 TString superstring(
"eigenVar");
1754 TString nameUp(superstring); nameUp +=
"_up";
1755 TString nameDown(superstring); nameDown +=
"_down";
1758 TH1* resultVariedUp = (TH1*)comb_result->Clone(nameUp); resultVariedUp->SetDirectory(0);
1759 TH1* resultVariedDown = (TH1*)comb_result->Clone(nameDown); resultVariedDown->SetDirectory(0);
1761 for (
int u = 0;
u < comb_result->GetNbinsX(); ++
u) {
1762 resultVariedUp->SetBinContent(
u,(comb_result->GetBinContent(
u) + matrixVariationsWithZeros(
i,
u)));
1763 resultVariedDown->SetBinContent(
u,(comb_result->GetBinContent(
u) - matrixVariationsWithZeros(
i,
u)));
1766 m_eigen.push_back(std::make_pair(resultVariedUp, resultVariedDown));
1784 size_t current_set = 0;
1787 min_variance = 1.0E-6;
1789 bool keep_variation =
false;
1790 TH1*
up =
static_cast<TH1*
>(
m_eigen[
index].first->Clone());
up->SetDirectory(0);
1791 up->Add(comb_result.get(), -1.0);
1794 if (fabs(
up->GetBinContent(
bin)) > min_variance) {
1795 keep_variation =
true;
1799 if (!keep_variation){
1800 final_set.insert(current_set);
1807 if (final_set.size() > 0){
1808 if (
m_verbose) std::cout <<
"CalibrationDataEigenVariations: Removing " << final_set.size() <<
" eigenvector variations leading to sub-tolerance effects, retaining " <<
m_eigen.size()-final_set.size() <<
" variations" << std::endl;
1816 std::streamsize
ss = std::cout.precision();
1818 std::cout.precision(
ss);
1837 for(
const std::pair<TH1*,TH1*>&
var :
m_eigen){
1839 TString eigenvarup =
var.first->GetName();
1840 TString eigenvardown =
var.second->GetName();
1841 int bin_baseline = 0;
1842 for (
const std::string& flavour :
m_flavours){
1844 TH1*
result =
dynamic_cast<TH1*
>(
c->GetValue(
"result"));
1846 std::cerr<<
"CalibrationDataGlobalEigenVariations::initialize: dynamic cast failed\n";
1849 TH1* resultVariedUp = (TH1*)
result->Clone(eigenvarup); resultVariedUp->SetDirectory(0);
1850 TH1* resultVariedDown = (TH1*)
result->Clone(eigenvardown); resultVariedDown->SetDirectory(0);
1851 int up_to_bin = flav_bins[flavour];
1852 int current_bin = 1;
1853 for(
int flav_bin = bin_baseline+1 ; flav_bin < up_to_bin+1 ; flav_bin++){
1854 Int_t
bin =
result->GetBin(1,current_bin);
1855 resultVariedUp->SetBinContent(
bin,
var.first->GetBinContent(flav_bin));
1856 resultVariedDown->SetBinContent(
bin,
var.second->GetBinContent(flav_bin));
1859 bin_baseline+=up_to_bin;
1860 m_flav_eigen[flavour].push_back(std::make_pair(resultVariedUp, resultVariedDown));
1885 std::vector<int> zeroComponents;
1889 if (fabs(
cov(
i,0)) < 1
e-10){
1890 bool isThereANonZero =
false;
1894 if (fabs(
cov(
i,j)) > 1
e-10){
1895 isThereANonZero =
true;
1899 if (! isThereANonZero){
1901 zeroComponents.push_back(
i) ;
1907 std::cerr <<
" Problem. Found n. " << nZeros <<
" while size of matrix is " <<
m_blockmatrixsize << std::endl;
1908 return TMatrixDSym();
1916 bool missed =
false;
1917 for (
unsigned int s = 0 ;
s < zeroComponents.size(); ++
s) {
1918 if (zeroComponents.at(
s) ==
i) {
1928 matrixVariationJacobian(
i-nMissed,
i)=1;
1931 return matrixVariationJacobian;
1945 if (
m_verbose) std::cout <<
"No m_flav_eigen for flavour " << flavour << std::endl;
1964 std::vector<std::pair<TH1*,TH1*>> flav_variations =
m_flav_eigen.find(flavour)->second;
1965 if (variation < flav_variations.size()){
1966 up = flav_variations[variation].first;
1967 down = flav_variations[variation].second;
2038 return (extrapIndex ==
int(nameIndex));
2054 simple_set.insert(
it);
2079 for (
IndexSet::iterator subset_it = set_it->begin(); subset_it != set_it->end(); ++subset_it){
2083 std::cerr <<
"Error in CalibrationDataGlobalEigenVariations::mergeVariations: \
2084 IndexSets must not overlap and must lie between 1 and " <<
m_eigen.size() <<
". Aborting!" << std::endl;
2090 std::string flavour = flav;
2094 TH1*
result =
dynamic_cast<TH1*
>(
c->GetValue(
"result"));
2096 std::cerr <<
"Error in CalibrationDataGlobalEigenVariations::mergeVariations: failed dynamic cast\n";
2103 int ndim =
result->GetDimension();
2109 if (set_it->empty())
continue;
2111 double sum_H_up = 0.0, sum_H_down = 0.0;
2112 size_t lowest_index = *set_it->lower_bound(0);
2113 TH1 *total_var_up =
static_cast<TH1*
>(
m_flav_eigen[flavour][lowest_index].first->Clone());
2114 TH1 *total_var_down =
static_cast<TH1*
>(
m_flav_eigen[flavour][lowest_index].second->Clone());
2115 total_var_up->SetDirectory(0);
2116 total_var_down->SetDirectory(0);
2118 total_var_up->Reset();
2119 total_var_down->Reset();
2122 for (
IndexSet::iterator subset_it = set_it->begin(); subset_it != set_it->end(); ++subset_it) {
2123 size_t actual_index = *subset_it;
2125 if (actual_index != lowest_index) toDelete.insert(*subset_it);
2127 TH1 *partial_var_up =
static_cast<TH1*
>(
m_flav_eigen[flavour][actual_index].first->Clone());
2128 TH1 *partial_var_down =
static_cast<TH1*
>(
m_flav_eigen[flavour][actual_index].second->Clone());
2129 partial_var_up->SetDirectory(0);
2130 partial_var_down->SetDirectory(0);
2132 partial_var_up->Add(
result, -1.0);
2133 partial_var_down->Add(
result, -1.0);
2135 partial_var_down->SetBinContent(
i, -1.0*partial_var_down->GetBinContent(
i));
2139 double sum_up = total_var_up->GetBinContent(
u);
2140 double sum_down = total_var_down->GetBinContent(
u);
2142 sum_up += partial_var_up->GetBinContent(
u)*partial_var_up->GetBinContent(
v);
2143 sum_H_up += partial_var_up->GetBinContent(
u)*partial_var_up->GetBinContent(
v);
2144 sum_down += partial_var_down->GetBinContent(
u)*partial_var_down->GetBinContent(
v);
2145 sum_H_down += partial_var_down->GetBinContent(
u)*partial_var_down->GetBinContent(
v);
2147 total_var_up->SetBinContent(
u, sum_up);
2148 total_var_down->SetBinContent(
u, sum_down);
2150 delete partial_var_up;
2151 delete partial_var_down;
2156 if (sum_H_up != 0.0)
2157 total_var_up->SetBinContent(
i, total_var_up->GetBinContent(
i)/sqrt(sum_H_up));
2159 total_var_up->SetBinContent(
i, 0.0);
2160 if (sum_H_down != 0.0)
2161 total_var_down->SetBinContent(
i, -1.0*total_var_down->GetBinContent(
i)/sqrt(sum_H_down));
2163 total_var_down->SetBinContent(
i, 0.0);
2166 total_var_up->Add(
result);
2167 total_var_down->Add(
result);
2169 m_flav_eigen[flavour][lowest_index].first = total_var_up;
2170 m_flav_eigen[flavour][lowest_index].second = total_var_down;
2181 if (
set.size() == 0)
return;
2182 std::vector<std::pair<TH1*, TH1*> > new_eigen;
2183 std::vector<std::pair<TH1*, TH1*>> eigen =
m_flav_eigen[flav];
2186 new_eigen.push_back(eigen[
index]);
2188 delete eigen[
index].first;
delete eigen[
index].second;
2203 for (
IndexSet::iterator subset_it = set_it->begin(); subset_it != set_it->end(); ++subset_it)
2204 simple_set.insert(*subset_it);