MuonRPCGeometry/test/RPCHistos.cc

0001 #include <iostream>
0002 #include <sstream>
0003 #include <fstream>
0004 #include <algorithm>
0005 #include <map>
0006 #include <cmath>
0007 #include "TH2D.h"
0008 #include "TCanvas.h"
0009 #include "TStyle.h"
0010 #include "TLegend.h"
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "L1Trigger/RPCTrigger/interface/RPCConst.h"
0013 #include "FWCore/Utilities/interface/Exception.h"
0014
0015 const double ptRanges[RPCConst::m_PT_CODE_MAX + 2] = {0.0, 0.01, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0,  4.5,  5.,   6.,
0016                                                       7.,  8.,   10., 12., 14., 16., 18., 20.,  25.,  30.,  35.,
0017                                                       40., 45.,  50., 60., 70., 80., 90., 100., 120., 140., 160.};
0018
0019 class RPCHistos {
0020 public:
0021   RPCHistos();
0022   ~RPCHistos();
0023   void go();
0024   TH2D *eff();
0025   void drawActCurv(int towerFrom, int towerTo);
0026   TH1D *drawActCurv(int ptCut, int towerFrom, int towerTo);
0027   void rateVsPt();
0028   double frate(double p);
0029   double rate(int pt, int tower);
0030   double rate(int pt);
0031   void phiEff();
0032   void etaCorr();
0033
0034   void setFormat(std::string fm) { m_ext = fm; };
0035
0036 private:
0037   std::ifstream m_ifile;
0038   std::string m_ext;
0039   std::map<int, TH2D *> m_tower2pt;
0040
0041   TH1D *m_phiEffRec;  // efficienct vs eta
0042   TH1D *m_phiEffAll;  // efficienct vs eta
0043   TH2F *m_etaCorr;
0044
0045   RPCConst m_const;
0046 };
0047
0048 RPCHistos::RPCHistos() {
0049   m_ext = "png";
0050   m_ifile.open("muons.txt");
0051   if (!m_ifile) {
0052     throw cms::Exception("IOError") << "Bkah\n";
0053   }
0054
0055   for (int i = -RPCConst::m_TOWER_COUNT + 1; i < RPCConst::m_TOWER_COUNT; ++i) {
0056     std::stringstream ss;
0057     ss << i;
0058     m_tower2pt[i] = new TH2D(ss.str().c_str(),
0059                              "",
0060                              RPCConst::m_PT_CODE_MAX + 1,
0061                              -0.5,
0062                              RPCConst::m_PT_CODE_MAX + 0.5,  // gen
0063                              RPCConst::m_PT_CODE_MAX + 1,
0064                              -0.5,
0065                              RPCConst::m_PT_CODE_MAX + 0.5);  // rec
0066   }
0067
0068   m_etaCorr = new TH2F("etaCorr", "", 35, -17, 17, 35, -17, 17);
0069
0070   //int bins = 440;
0071   int bins = 110;
0072   m_phiEffRec = new TH1D("phiEffRec", "efficiency vs eta", bins, -2.2, 2.2);
0073   m_phiEffAll = new TH1D("phiEffAll", "eff vs eta", bins, -2.2, 2.2);
0074 }
0075
0076 /****************************************************************
0077 *
0078 * dtor
0079 *
0080 *****************************************************************/
0081 RPCHistos::~RPCHistos() {
0082   // std::for_each(m_tower2pt.begin(), m_tower2pt.end(), std::with_data(delete));
0083
0084   std::map<int, TH2D *>::iterator it = m_tower2pt.begin();
0085   for (; it != m_tower2pt.end(); ++it) {
0086     delete it->second;
0087   }
0088
0089   delete m_phiEffRec;
0090   delete m_phiEffAll;
0091 }
0092 /****************************************************************
0093 *
0094 * `Do the work` function
0095 *
0096 *****************************************************************/
0097 void RPCHistos::go() {
0098   float etaGen = 0, ptGen = 0, phiGen = 0;
0099   int ptCodeRec = 0, towerRec = 0, phiRec = 0, ghost = 0, qual = 0;
0100
0101   std::cout << "Readout start" << std::endl;
0102   int count = 0;
0103   while (1) {
0104     if ((++count) % 20000 == 0)
0105       std::cout << " Read: " << count << std::endl;
0106
0107     if (count == 1000000)
0108       break;
0109     // -1.01966 -1.67928 99.9912 7 34 31 2
0110
0111     /*
0112                 m_outfileR << etaGen << " " << phiGen << " " << ptGen << " "
0113                             << towerRec << " " << phiRec << " " << ptCodeRec << " " << qual<< " "
0114                             << ghost
0115
0116    */
0117     m_ifile >> etaGen >> phiGen >> ptGen >> towerRec >> phiRec >> ptCodeRec >> qual >> ghost;
0118
0119     if (!m_ifile.good())
0120       break;
0121     int ptCodeGen = m_const.iptFromPt(ptGen);
0122     int towerGen = m_const.towerNumFromEta(etaGen);
0123     m_tower2pt[towerGen]->Fill(ptCodeGen, ptCodeRec);
0124
0125     // PhiEff plot, rate plot
0126     //    float etaGenAbs = etaGen > 0 ? etaGen : -etaGen;
0127     //m_phiEffAll->Fill(etaGenAbs);
0128     m_phiEffAll->Fill(etaGen);
0129     if (ghost == 0 && ptCodeRec != 0) {
0130       // m_phiEffRec->Fill(etaGenAbs);
0131       m_phiEffRec->Fill(etaGen);
0132       m_etaCorr->Fill(towerGen, towerRec);
0133     }
0134   }
0135
0136   std::cout << "Readout finished. Read " << count << std::endl;
0137 }
0138 // integrated rate for pt > p, whole detector (|eta|<2.1)
0139 double RPCHistos::frate(double p) {
0140   double a = -0.235802;  //        +/- 0.01494      (6.335%)
0141   double b = -2.82345;   //         +/- 0.09143      (3.238%)
0142   double c = 17.162;     //           +/- 0.1342       (0.7822%)
0143
0144   float ret = pow(p, (a * log(p))) * pow(p, b) * exp(c);
0145
0146   if (ret <= 0)
0147     throw cms::Exception("float") << "Lower then zero \n";
0148
0149   return ret;
0150 }
0151
0152 double RPCHistos::rate(int ptCode, int tower) {
0153   if (tower > 16 || tower < -16)
0154     throw cms::Exception("data") << " Bad tower\n ";
0155
0156   double etaNorm = 2.1 * 2;
0157
0158   std::vector<float> towerEtaSize;
0159   towerEtaSize.push_back((0.07 - 0.00) * 2);  //T0
0160   towerEtaSize.push_back(0.27 - 0.07);        //1
0161   towerEtaSize.push_back(0.44 - 0.27);        //2
0162   towerEtaSize.push_back(0.58 - 0.44);        //3
0163   towerEtaSize.push_back(0.72 - 0.58);        //4
0164   towerEtaSize.push_back(0.83 - 0.72);        //5
0165   towerEtaSize.push_back(0.93 - 0.83);        //6
0166   towerEtaSize.push_back(1.04 - 0.93);        //7
0167   towerEtaSize.push_back(1.14 - 1.04);        //8
0168   towerEtaSize.push_back(1.24 - 1.14);        //9
0169   towerEtaSize.push_back(1.36 - 1.24);        //10
0170   towerEtaSize.push_back(1.48 - 1.36);        //11
0171   towerEtaSize.push_back(1.61 - 1.48);        //12
0172   towerEtaSize.push_back(1.73 - 1.61);        //13
0173   towerEtaSize.push_back(1.85 - 1.73);        //14
0174   towerEtaSize.push_back(1.97 - 1.85);        //15
0175   towerEtaSize.push_back(2.10 - 1.97);        //16
0176
0177   return rate(ptCode) * towerEtaSize.at(std::abs(tower)) / etaNorm;
0178 }
0179 double RPCHistos::rate(int ptCode) {
0180   // Is it ok??
0181   if (ptCode == 0)
0182     return 0;
0183
0184   if (ptCode > 31 || ptCode < 0)
0185     throw cms::Exception("data") << " Bad ptCode " << ptCode << "\n";
0186
0187   double pts[] = {
0188       0.0, 0.01, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.,  6.,  7.,  8.,   10.,  12.,  14.,  16.,
0189       18., 20.,  25., 30., 35., 40., 45., 50., 60., 70., 80., 90., 100., 120., 140., 1000.};  // Make 1 TeV highest (inf ;)
0190
0191   //std::cout<< pts[32] << std::endl;
0192   double r2 = frate(pts[ptCode + 1]);
0193   double r1 = frate(pts[ptCode]);
0194
0195   return r1 - r2;
0196 }
0197 /****************************************************************
0198 *
0199 * rate vs ptCut histo
0200 *
0201 ****************************************************************/
0202 void RPCHistos::rateVsPt() {
0203   //  m_rateVsptCut->Divide(m_rateVsptCutAll);
0204
0205   TH2D *th2eff = eff();
0206
0207   for (int x = 1; x <= 16 + 1; ++x)  // tower
0208   {
0209     // 1 means ptCode 0
0210     for (int y = 1; y <= 31 + 1; ++y) {  // ptCode
0211       double bincont = th2eff->GetBinContent(x, y);
0212       //bincont*=m_const.vxIntegMuRate(y-1,x-1);
0213       bincont *= rate(y - 1, x - 1);
0214       th2eff->SetBinContent(x, y, bincont);
0215     }
0216   }
0217
0218   // Sum rate in all towers;
0219   TH1D *th1rate = th2eff->ProjectionY("ptRate", 1, 18);
0220   TH1D *genRate = new TH1D("genRate", "", RPCConst::m_PT_CODE_MAX + 1, -0.5, RPCConst::m_PT_CODE_MAX + 0.5);
0221
0222   double gratesum = 0;
0223   for (int y = 1; y <= 31 + 1; ++y) {
0224     gratesum += rate(y - 1);
0225   }
0226
0227   for (int y = 1; y <= 31 + 1; ++y) {
0228     //     genRate->SetBinContent(y, m_const.vxMuRate(y-1));
0229     gratesum -= rate(y - 1);
0230     genRate->SetBinContent(y, gratesum);
0231   }
0232
0233   //th1rate->Multiply(genRate);
0234   // Make integrated rate
0235   double records = 0, sum = 0;
0236
0237   //*
0238   for (int i = 1; i <= 31 + 1; ++i) {
0239     records += th1rate->GetBinContent(i);
0240   }
0241
0242   for (int i = 1; i <= 31 + 1; ++i) {
0243     double sumDelta = th1rate->GetBinContent(i);
0244     double cont = (records - sum);
0245     sum += sumDelta;
0246     th1rate->SetBinContent(i, cont);
0247   }  //*/
0248      //////////////////
0249      ///////////////
0250
0251   for (int ptCut = 0; ptCut <= 31; ++ptCut) {
0252     double arate = 0;
0253     //for (int tower = 0; tower < 17; ++tower){
0254     for (int tower = 0; tower < 13; ++tower) {
0255       TH1D *cur = drawActCurv(ptCut, tower, tower);
0256       for (int i = 2; i <= RPCConst::m_PT_CODE_MAX + 1; ++i) {
0257         double cont = cur->GetBinContent(i);
0258         arate += cont * rate(i - 1, tower);
0259       }
0260       delete cur;
0261     }
0262     std::cout << "TT " << ptCut << " " << arate << std::endl;
0263     th1rate->SetBinContent(ptCut + 1, arate);
0264   }
0265
0266   TCanvas c1;
0267   c1.SetLogy();
0268   c1.SetLogx();
0269
0270   th1rate->SetBins(RPCConst::m_PT_CODE_MAX + 1, ptRanges);
0271   genRate->SetBins(RPCConst::m_PT_CODE_MAX + 1, ptRanges);
0272   th1rate->SetLineColor(1);
0273   genRate->SetLineColor(2);
0274
0275   th1rate->SetStats(false);
0276   genRate->SetStats(false);
0277
0278   //th1rate->GetXaxis()->SetRange(1,100);
0279   //genRate->GetXaxis()->SetRange(1,100);
0280   th1rate->Draw("");
0281   genRate->Draw("SAME");
0282   th1rate->GetXaxis()->SetRangeUser(2, 100);
0283
0284   th1rate->SetMaximum(6000000);
0285   th1rate->SetMinimum(0.1);
0286   //genRate->GetXaxis()->SetRangeUser(1,100);
0287
0288   th1rate->GetYaxis()->SetTitle("rate [Hz]");
0289   th1rate->GetXaxis()->SetTitle("ptCodeCut");
0290
0291   TLegend *legend = new TLegend(0.66, 0.15 + 0.5, 0.86, 0.37 + 0.5);
0292   legend->AddEntry(th1rate, "RPC", "l");
0293   legend->AddEntry(genRate, "Generated rate", "l");
0294   legend->Draw();
0295   std::string fname = "rate." + m_ext;
0296   c1.Print(fname.c_str());
0297 }
0298 /****************************************************************
0299 *
0300 * Reco vs eta histo
0301 *
0302 ****************************************************************/
0303 void RPCHistos::phiEff() {
0304   m_phiEffRec->Divide(m_phiEffAll);
0305   TCanvas c1;
0306   m_phiEffRec->SetStats(false);
0307   m_phiEffRec->Draw("");
0308   m_phiEffRec->SetMaximum(1);
0309   m_phiEffRec->SetMinimum(0);
0310   m_phiEffRec->GetYaxis()->SetTitle("rpc trg efficiency");
0311   m_phiEffRec->GetXaxis()->SetTitle("eta generated");
0312   std::string fname = "phiEff." + m_ext;
0313   c1.Print(fname.c_str());
0314 }
0315 /****************************************************************
0316 *
0317 * Eta corelation
0318 *
0319 ****************************************************************/
0320 void RPCHistos::etaCorr() {
0321   TCanvas c1;
0322   m_etaCorr->SetFillColor(2);
0323   m_etaCorr->GetXaxis()->SetTitle("tower generated");
0324   m_etaCorr->GetYaxis()->SetTitle("tower reconstructed");
0325   m_etaCorr->SetStats(false);
0326   m_etaCorr->Draw("BOX");
0327
0328   std::string filename = "etacorr." + m_ext;
0329   c1.Print(filename.c_str());
0330 }
0331 /****************************************************************
0332 *
0333 * Draw efficiency histo (tower vs ptCode). Split draw and return.
0334 *
0335 ****************************************************************/
0336 TH2D *RPCHistos::eff() {
0337   TH2D *reco2d = new TH2D("reco2d",
0338                           "",
0339                           //RPCConst::m_TOWER_COUNT, -0.5, RPCConst::m_TOWER_COUNT-0.5, //Towers
0340                           RPCConst::m_TOWER_COUNT * 2 - 1,
0341                           -RPCConst::m_TOWER_COUNT + 0.5,
0342                           RPCConst::m_TOWER_COUNT - 0.5,  //Towers
0343                           RPCConst::m_PT_CODE_MAX + 1,
0344                           -0.5,
0345                           RPCConst::m_PT_CODE_MAX + 0.5);  // Pt's
0346   TH2D *lost2d = new TH2D("lost2d",
0347                           "",
0348                           //RPCConst::m_TOWER_COUNT, -0.5, RPCConst::m_TOWER_COUNT-0.5, //Towers
0349                           RPCConst::m_TOWER_COUNT * 2 - 1,
0350                           -RPCConst::m_TOWER_COUNT + 0.5,
0351                           RPCConst::m_TOWER_COUNT - 0.5,  //Towers
0352                           RPCConst::m_PT_CODE_MAX + 1,
0353                           -0.5,
0354                           RPCConst::m_PT_CODE_MAX + 0.5);  // Pt's
0355
0356   std::map<int, TH2D *>::iterator it = m_tower2pt.begin();
0357   for (; it != m_tower2pt.end(); ++it) {
0358     TH1D *reco = it->second->ProjectionX("dd", 2, RPCConst::m_PT_CODE_MAX + 1);
0359     TH1D *lost = it->second->ProjectionX("lost", 1, 1);  // ptCode 0 is in first bin
0360     int tower = it->first;
0361     for (int i = 1; i <= RPCConst::m_PT_CODE_MAX + 1; ++i) {
0362       //double rcur=reco2d->GetBinContent(std::abs(tower)+1,i);
0363       double rcur = reco2d->GetBinContent(tower + RPCConst::m_TOWER_COUNT, i);
0364       double radd = reco->GetBinContent(i);
0365       //reco2d->SetBinContent(std::abs(tower)+1,i,rcur+radd);
0366       reco2d->SetBinContent(tower + RPCConst::m_TOWER_COUNT, i, rcur + radd);
0367
0368       //double lcur=lost2d->GetBinContent(std::abs(tower)+1,i);
0369       double lcur = lost2d->GetBinContent(tower + RPCConst::m_TOWER_COUNT, i);
0370       double ladd = lost->GetBinContent(i);
0371       //lost2d->SetBinContent(std::abs(tower)+1,i,lcur+ladd);
0372       lost2d->SetBinContent(tower + RPCConst::m_TOWER_COUNT, i, lcur + ladd);
0373     }
0374
0375     delete reco;
0376     delete lost;
0377   }
0378
0379   TH2D *all2d = (TH2D *)reco2d->Clone();
0380   all2d->Add(lost2d);
0381
0382   reco2d->Divide(all2d);
0383   TCanvas c1;
0384   gStyle->SetPalette(1, 0);
0385   reco2d->SetStats(false);
0386   reco2d->Draw("COLZ");
0387   reco2d->SetMaximum(1);
0388   reco2d->SetMinimum(0);
0389   reco2d->GetXaxis()->SetTitle("tower generated");
0390   reco2d->GetYaxis()->SetTitle("ptCode generated");
0391   std::string fname = "effkar." + m_ext;
0392   c1.Print(fname.c_str());
0393
0394   return reco2d;
0395 }
0396
0397 /****************************************************************
0398 *
0399 * Returns act. curve for single ptCut for given towers range
0400 *  (inspiered by Karols implementation)
0401 *
0402 ****************************************************************/
0403 TH1D *RPCHistos::drawActCurv(int ptCut, int towerFrom, int towerTo) {
0404   TH2D *sum = new TH2D("sums",
0405                        "",
0406                        RPCConst::m_PT_CODE_MAX + 1,
0407                        -0.5,
0408                        RPCConst::m_PT_CODE_MAX + 0.5,  // gen
0409                        RPCConst::m_PT_CODE_MAX + 1,
0410                        -0.5,
0411                        RPCConst::m_PT_CODE_MAX + 0.5);  // rec
0412
0413   std::map<int, TH2D *>::iterator it = m_tower2pt.begin();
0414   for (; it != m_tower2pt.end(); ++it) {
0415     int absTower = std::abs(it->first);
0416     if (absTower < towerFrom || absTower > towerTo)
0417       continue;
0418     sum->Add(it->second);
0419   }
0420
0421   std::stringstream name;
0422   name << ptCut << "_" << towerFrom << "_" << towerTo;
0423   TH1D *passed = sum->ProjectionX(name.str().c_str(), ptCut + 1, RPCConst::m_PT_CODE_MAX + 1);
0424   TH1D *all = sum->ProjectionX("lost", 1, RPCConst::m_PT_CODE_MAX + 1);
0425
0426   std::stringstream title;
0427   title << "towers " << towerFrom << " - " << towerTo;
0428   passed->SetTitle(title.str().c_str());
0429
0430   passed->Divide(all);
0431
0432   passed->SetBins(RPCConst::m_PT_CODE_MAX + 1, ptRanges);
0433   delete sum;
0434   return passed;
0435 }
0436 /****************************************************************
0437 *
0438 * Draws act. curves for given tower range
0439 *
0440 ****************************************************************/
0441 void RPCHistos::drawActCurv(int twFrom, int twTo) {
0442   std::vector<int> ptTresh;
0443   ptTresh.push_back(9);
0444   ptTresh.push_back(15);
0445   ptTresh.push_back(18);
0446   ptTresh.push_back(22);
0447
0448   TCanvas c1;
0449   c1.SetTicks(1, 1);
0450   c1.SetLogx();
0451   c1.SetGrid();
0452
0453   gStyle->SetPalette(1, 0);
0454
0455   Color_t color = kRed;
0456   TLegend *legend = new TLegend(0.66, 0.15, 0.86, 0.37);
0457
0458   for (unsigned int i = 0; i < ptTresh.size(); ++i) {
0459     ++color;
0460     TH1D *a1 = drawActCurv(ptTresh.at(i), twFrom, twTo);
0461     a1->SetLineColor(color);
0462     if (color == kYellow)
0463       a1->SetLineColor(46);
0464
0465     a1->SetStats(false);
0466     a1->GetXaxis()->SetRangeUser(2, 200);
0467     a1->Draw("SAME");
0468     a1->SetMaximum(1);
0469     a1->SetMinimum(0);
0470     a1->GetXaxis()->SetTitle("pt generated [GeV]");
0471     a1->GetYaxis()->SetTitle("efficiency");
0472
0473     // Legend
0474     std::stringstream leg;
0475     leg << "pt >= " << m_const.ptFromIpt(ptTresh.at(i)) << " GeV";
0476     //leg << "ptCode >= " <<ptTresh.at(i) << " GeV";
0477     legend->AddEntry(a1, leg.str().c_str(), "l");
0478   }
0479
0480   legend->Draw();
0481   c1.Update();
0482
0483   std::stringstream fname;
0484   fname << "akt" << twFrom << "_" << twTo << "." << m_ext;
0485   c1.Print(fname.str().c_str());
0486 }
0487 /****************************************************************
0488 *
0489 *
0490 *
0491 ****************************************************************/
0492 int main() {
0493   gStyle->SetCanvasBorderMode(0);
0494   gStyle->SetPadBorderMode(0);
0495   gStyle->SetPadColor(0);
0496   gStyle->SetCanvasColor(0);
0497   gStyle->SetTitleBorderSize(0);
0498   gStyle->SetStatColor(0);
0499   gStyle->SetOptDate(0);
0500   gStyle->SetPalette(1);
0501   gStyle->SetOptFit(0111);
0502   gStyle->SetOptStat(1000000);
0503
0504   RPCHistos rh;
0505   rh.setFormat("png");
0506   rh.go();
0507   rh.eff();
0508   rh.phiEff();
0509   rh.etaCorr();
0510   //rh.rateVsPt();
0511   rh.drawActCurv(0, 7);
0512   rh.drawActCurv(8, 12);
0513   return 0;
0514 }