test/BSMTriggerCheck/EffPullCalculator.cc

0001 using namespace std;
0002 #include <TH1D.h>
0003 #include <TAxis.h>
0004 #include <TMath.h>
0005 #include "EffPullCalculator.hh"
0006
0007 EffPullcalculator::EffPullcalculator(TH1D* pathHisto1_v, TH1D* pathHisto2_v, vector<TH1D*> sortedHisto1_v, vector<TH1D*> sortedHisto2_v, string error_v) {
0008   pathHisto.push_back(pathHisto1_v);
0009   pathHisto.push_back(pathHisto2_v);
0010   sortedHisto1 = sortedHisto1_v;
0011   sortedHisto2 = sortedHisto2_v;
0012   error = error_v;
0013 }
0014
0015 void EffPullcalculator::CalculatePulls() {
0016
0017   if(pathHisto[0]->GetNbinsX() != pathHisto[1]->GetNbinsX()) {
0018     lines.push_back("--------------------------------------------------------------------------------");
0019     lines.push_back("WARNING in EffPullcalculator::CalculatePulls: the vectors to compare ");
0020     lines.push_back("                                              have different size."   );
0021     lines.push_back("                                              Some HLT path was added or removed");
0022     lines.push_back(" ");
0023   }
0024
0025   TAxis* axis1 = pathHisto[0]->GetXaxis();
0026   TAxis* axis2 = pathHisto[1]->GetXaxis();
0027   double nEv1 = 0;
0028   double nEv2 = 0;
0029
0030   for(int i1=1; i1 <= pathHisto[0]->GetNbinsX(); i1++) {
0031     string label1 = axis1->GetBinLabel(i1);
0032     if(label1 == "Total") nEv1 = pathHisto[0]->GetBinContent(i1);
0033     for(int i2=1; i2 <= pathHisto[1]->GetNbinsX(); i2++) {
0034       string label2 = axis2->GetBinLabel(i2);
0035       if(label2 == "Total") nEv2 = pathHisto[1]->GetBinContent(i2);
0036       if(label1 == label2) { // we are comparing the same path
0037     eff1.push_back(pathHisto[0]->GetBinContent(i1));
0038     eff2.push_back(pathHisto[1]->GetBinContent(i2));
0039     name.push_back(label1);
0040       }
0041     }
0042   }
0043
0044   for(int i = 0; i< int(eff1.size()); i++) {
0045     eff1[i] *= (nEv1 == 0 ? 0. : 1./nEv1);
0046     eff2[i] *= (nEv2 == 0 ? 0. : 1./nEv2);
0047      err_eff1.push_back(sqrt(eff1[i]*(1.-eff1[i])/nEv1));
0048      err_eff2.push_back(sqrt(eff2[i]*(1.-eff2[i])/nEv2));
0049 //     err_eff1.push_back(sqrt(eff1[i]/nEv1));
0050 //     err_eff2.push_back(sqrt(eff2[i]/nEv2));
0051   }
0052
0053   // create 10x2 histograms with efficiencies
0054   int nbin = 10;//eff1.size()/8+1;
0055   vector<int> nbins; nbins.resize(20);
0056   nbins[0]=nbin;
0057   nbins[1]=nbin;
0058   nbins[10]=nbin;
0059   nbins[11]=nbin;
0060   for(int i=2; i<10; ++i)
0061     {
0062       if(i-2<int(sortedHisto1.size()))
0063     {
0064       nbins[i] = sortedHisto1[i-2]->GetNbinsX();
0065       nbins[i+10] = sortedHisto1[i-2]->GetNbinsX();
0066     }
0067       else
0068     {
0069       nbins[i]=0;
0070       nbins[i+10]=0;
0071     }
0072     }
0073   for(int iplot=0; iplot<20;iplot++) {
0074     char name[256];
0075     sprintf(name,"%s_eff_%i",pathHisto[0]->GetName(),iplot);
0076     effhisto.push_back(new TH1D(name, "Trigger Efficiency", nbins[iplot], 0., double(nbins[iplot])));
0077   }
0078   for(int i=2; i<10 && i-2<int(sortedHisto1.size()); ++i)
0079     {
0080       effhisto[i]->SetTitle(sortedHisto1[i-2]->GetTitle());
0081       effhisto[i+10]->SetTitle(sortedHisto1[i-2]->GetTitle());
0082     }
0083
0084   for(int k=0; k<int(sortedHisto1.size()); ++k)
0085     {
0086       int i=0, j=0;
0087       cout<<"111\n";
0088       for(i=2; i<=sortedHisto1[k]->GetNbinsX(); ++i)
0089     {
0090       if(std::abs(sortedHisto1[k]->GetBinContent(i) - sortedHisto2[k]->GetBinContent(i)) > std::abs(sortedHisto1[k]->GetBinContent(i-1)-sortedHisto2[k]->GetBinContent(i-1)))
0091         {
0092           double content1=sortedHisto1[k]->GetBinContent(i), content2=sortedHisto2[k]->GetBinContent(i), error1=sortedHisto1[k]->GetBinError(i), error2=sortedHisto2[k]->GetBinError(i);
0093           char name1[256], name2[256], char_name[256];
0094           sprintf(name1,"%s",sortedHisto1[k]->GetXaxis()->GetBinLabel(i));
0095           sprintf(name2,"%s",sortedHisto2[k]->GetXaxis()->GetBinLabel(i));
0096           for(j=i-1; j>=1; --j)
0097         {
0098           if(std::abs(content1-content2) > std::abs(sortedHisto1[k]->GetBinContent(j) - sortedHisto2[k]->GetBinContent(j)))
0099             {
0100               sortedHisto1[k]->SetBinContent(j+1, sortedHisto1[k]->GetBinContent(j));
0101               sortedHisto1[k]->SetBinError(j+1, sortedHisto1[k]->GetBinError(j));
0102               sprintf(char_name,"%s",sortedHisto1[k]->GetXaxis()->GetBinLabel(j));
0103               sortedHisto1[k]->GetXaxis()->SetBinLabel(j+1, char_name);
0104               sortedHisto2[k]->SetBinContent(j+1, sortedHisto2[k]->GetBinContent(j));
0105               sortedHisto2[k]->SetBinError(j+1, sortedHisto2[k]->GetBinError(j));
0106               sprintf(char_name,"%s",sortedHisto2[k]->GetXaxis()->GetBinLabel(j));
0107               sortedHisto2[k]->GetXaxis()->SetBinLabel(j+1, char_name);
0108             }
0109           else
0110             break;
0111         }
0112           sortedHisto1[k]->SetBinContent(j+1, content1);
0113           sortedHisto1[k]->SetBinError(j+1, error1);
0114           sortedHisto1[k]->GetXaxis()->SetBinLabel(j+1, name1);
0115           sortedHisto2[k]->SetBinContent(j+1, content2);
0116           sortedHisto2[k]->SetBinError(j+1, error2);
0117           sortedHisto2[k]->GetXaxis()->SetBinLabel(j+1, name2);
0118         }
0119     }
0120     }
0121
0122   // set plot style
0123   for(int iplot=0; iplot<10;iplot++) {
0124     effhisto[iplot]->GetXaxis()->SetLabelSize(0.03);
0125     effhisto[iplot+10]->GetXaxis()->SetLabelSize(0.03);
0126     effhisto[iplot]->SetStats(0);
0127     effhisto[iplot+10]->SetStats(0);
0128     effhisto[iplot]->SetYTitle("Efficiency");
0129     effhisto[iplot+10]->SetYTitle("Efficiency");
0130     effhisto[iplot]->SetFillColor(38);
0131     effhisto[iplot+10]->SetFillColor(28);
0132     effhisto[iplot]->SetBarWidth(0.4*effhisto[iplot]->GetXaxis()->GetBinWidth(1));
0133     effhisto[iplot]->SetBarOffset(0.1*effhisto[iplot]->GetXaxis()->GetBinWidth(1));
0134     effhisto[iplot+10]->SetBarWidth(0.4*effhisto[iplot+10]->GetXaxis()->GetBinWidth(1));
0135     effhisto[iplot+10]->SetBarOffset(0.5*effhisto[iplot+10]->GetXaxis()->GetBinWidth(1));
0136   }
0137
0138   for(int i=2; i<10;++i)
0139     {
0140       if(i-2<int(sortedHisto1.size()))
0141     for(int j=0; j<nbins[i]; ++j)
0142       {
0143         effhisto[i]->SetBinContent(j+1,sortedHisto1[i-2]->GetBinContent(j+1));
0144         effhisto[i]->SetBinError(j+1,sortedHisto1[i-2]->GetBinError(j+1));
0145         effhisto[i]->GetXaxis()->SetBinLabel(j+1,sortedHisto1[i-2]->GetXaxis()->GetBinLabel(j+1));
0146         effhisto[i+10]->SetBinContent(j+1,sortedHisto2[i-2]->GetBinContent(j+1));
0147         effhisto[i+10]->SetBinError(j+1,sortedHisto2[i-2]->GetBinError(j+1));
0148         effhisto[i+10]->GetXaxis()->SetBinLabel(j+1,sortedHisto2[i-2]->GetXaxis()->GetBinLabel(j+1));
0149       }
0150     }
0151
0152   // Histograms with largest efficiencies
0153   vector<int> SortedEff = SortVec(eff2);
0154
0155   int nb=1;
0156   for(int i=0; i< int(SortedEff.size()); i++) {
0157     if(eff1[SortedEff[i]] != 0. && name[SortedEff[i]] != "Total") {
0158       effhisto[0]->SetBinContent(nb,eff1[SortedEff[i]]);
0159       effhisto[0]->SetBinError(nb,err_eff1[SortedEff[i]]);
0160       effhisto[10]->SetBinContent(nb,eff2[SortedEff[i]]);
0161       effhisto[10]->SetBinError(nb,err_eff2[SortedEff[i]]);
0162       effhisto[0]->GetXaxis()->SetBinLabel(nb,name[SortedEff[i]].c_str());
0163       effhisto[10]->GetXaxis()->SetBinLabel(nb,name[SortedEff[i]].c_str());
0164       nb++;
0165     }
0166     if(nb == effhisto[0]->GetXaxis()->GetNbins()+1) break;
0167   }
0168
0169   // Histograms with largest differences in efficiencies
0170   vector<double> deltaeff;
0171   for(int i=0; i< int(eff1.size()); i++)
0172     deltaeff.push_back(fabs(eff1[i]-eff2[i]));
0173   SortedEff = SortVec(deltaeff);
0174   nb=1;
0175   for(int i=0; i< int(SortedEff.size()); i++) {
0176     if(eff2[SortedEff[i]] != 0. && eff1[SortedEff[i]] != 0.) {
0177       effhisto[1]->SetBinContent(nb,eff1[SortedEff[i]]);
0178       effhisto[1]->SetBinError(nb,err_eff1[SortedEff[i]]);
0179       effhisto[11]->SetBinContent(nb,eff2[SortedEff[i]]);
0180       effhisto[11]->SetBinError(nb,err_eff2[SortedEff[i]]);
0181       effhisto[1]->GetXaxis()->SetBinLabel(nb,name[SortedEff[i]].c_str());
0182       effhisto[11]->GetXaxis()->SetBinLabel(nb,name[SortedEff[i]].c_str());
0183       nb++;
0184     }
0185     if(nb == effhisto[1]->GetXaxis()->GetNbins()+1) break;
0186   }
0187
0188   // pull histogram
0189   pullDist = new TH1D(string(string(pathHisto[0]->GetName())+"_pullDist").c_str(),
0190               "Efficiency pulls", 20, -10., 10.);
0191   TAxis* pullXaxis = pullDist->GetXaxis();
0192   pullXaxis->SetTitle("Pull of trigger efficiency");
0193
0194   //  residual histogram
0195   resHisto = new TH1D(string(string(pathHisto[0]->GetName())+"_resHisto").c_str(),
0196               "Efficiency residuals", eff1.size(), 0., eff1.size()*1.);
0197   TAxis* resXaxis = resHisto->GetXaxis();
0198   resXaxis->SetTitle("");
0199   TAxis* resYaxis = resHisto->GetYaxis();
0200   resYaxis->SetTitle("Efficiency Residual");
0201
0202   for(int i = 0; i< int(eff1.size()); i++) {
0203     if(eff1[i] == 0. && eff2[i] != 0) {
0204       lines.push_back("---------------------------------------------------------------------------------------------");
0205       lines.push_back("WARNING: " + name[i] + " Trigger Path was masked and now it is used");
0206       lines.push_back(" ");
0207     } else if(eff2[i] == 0. && eff1[i] != 0) {
0208       lines.push_back("---------------------------------------------------------------------------------------------");
0209       lines.push_back("WARNING: " + name[i] + " Trigger Path was used and now it is masked");
0210       lines.push_back(" ");
0211     } else {
0212       double err = (error == "uncorrelated" ?
0213             sqrt(std::pow(err_eff1[i],2.)+ std::pow(err_eff2[i],2.)) :
0214             err_eff1[i]-err_eff2[i]);
0215       double pull = (eff1[i]-eff2[i])/err;
0216       if(fabs(pull) > 3.) {
0217     lines.push_back("---------------------------------------------------------------------------------------------");
0218     lines.push_back("WARNING: A discrepancy bigger than 3 sigmas was found for " + name[i] + " Trigger Path");
0219     lines.push_back(" ");
0220       }
0221
0222       resXaxis->SetBinLabel(i+1,name[i].c_str());
0223       pullDist->Fill(pull);
0224       resHisto->SetBinContent(i+1,eff1[i]-eff2[i]);
0225       resHisto->SetBinError(i+1,err);
0226     }
0227   }
0228 }
0229
0230 double EffPullcalculator::GetEff(string label, int ind) {
0231   double eff = 0.;
0232   if(ind<2) {
0233     TAxis* axis = pathHisto[ind]->GetXaxis();
0234     for(int i=1; i<=pathHisto[ind]->GetNbinsX(); i++) {
0235       if(axis->GetBinLabel(i) == label) {
0236     // last bin is total
0237     eff = pathHisto[ind]->GetBinContent(i)/pathHisto[ind]->GetBinContent(pathHisto[ind]->GetNbinsX());
0238       }
0239     }
0240   }
0241   return eff;
0242 }
0243
0244 void EffPullcalculator::WriteLogFile(string namefile) {
0245   FILE*  f=fopen(namefile.c_str(),"w");
0246   for(int j = 0; j< int(lines.size()); j++)
0247     fprintf(f,"%s\n",lines[j].c_str());
0248   fclose(f);
0249 }
0250
0251 vector<int> EffPullcalculator::SortVec(vector<double> eff) {
0252   vector<double> eff_tmp = eff;
0253   vector<double> eff_tmp2 = eff;
0254   vector<int> indexes;
0255
0256   std::sort(eff_tmp.begin(), eff_tmp.end());
0257
0258   for(int i=0; i< int(eff_tmp.size()); i++)
0259     for(int j=0; j< int(eff.size()); j++)
0260       if(eff_tmp2[j] == eff_tmp[eff_tmp.size()-1-i]) {
0261     indexes.push_back(j);
0262     eff_tmp2[j] = 1.1; // to ignore it the next round
0263     //  eff_tmp2.remove(i);
0264     j=0;
0265     i++;
0266       }
0267
0268   return indexes;
0269 }
0270
0271 void EffPullcalculator::AddGoldenPath(string name) {
0272   goldenpaths.push_back(name);
0273 }
0274
0275 void EffPullcalculator::PrintTwikiTable(string filename) {
0276
0277   FILE* f=fopen(filename.c_str(),"w");
0278
0279   for(int j=0; j< int(goldenpaths.size()); j++) {
0280     for(int i=0; i < int(name.size()); i++) {
0281       if(name[i] == goldenpaths[j]) {
0282     fprintf(f,"%s|%i +/- %i|%i +/- %i |\n",name[i].c_str(),eff1[i],err_eff1[i],eff2[i],err_eff2[i]);
0283       }
0284     }
0285   }
0286   fclose(f);
0287 }
0288
0289 double EffPullcalculator::abs(double value)
0290 {
0291   if(value<0)
0292     return -1.0*value;
0293   else
0294     return value;
0295 }