File indexing completed on 2024-04-06 12:32:55
0001
0002
0003
0004
0005
0006 #include <iostream>
0007 #include <string>
0008 #include <sstream>
0009 #include <fstream>
0010
0011
0012 #include "Validation/MuonRPCGeometry/interface/PTStatistics.h"
0013 #include "Validation/MuonRPCGeometry/test/style.h"
0014
0015 #include "TH2D.h"
0016 #include "THStack.h"
0017 #include "TCanvas.h"
0018 #include "TStyle.h"
0019 #include "TLegend.h"
0020 #include "TFile.h"
0021
0022
0023
0024
0025
0026
0027
0028 double etas[] = {-0.07, 0.07, 0.27, 0.44, 0.58, 0.72,
0029
0030 0.83, 0.93, 1.04, 1.14, 1.24, 1.36,
0031
0032 1.48, 1.61, 1.73, 1.85, 1.97, 2.10};
0033
0034 double pts[] = {
0035 0.0, 0.01,
0036
0037 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5,
0038
0039 5., 6., 7., 8.,
0040
0041 10., 12., 14., 16., 18.,
0042 20., 25., 30., 35., 40., 45.,
0043 50., 60., 70., 80., 90., 100., 120., 140.};
0044
0045
0046
0047
0048
0049 struct SMuonBinPos{
0050 SMuonBinPos(): m_binNo(-1), m_tower(0) {};
0051 SMuonBinPos(int tower, int bin): m_binNo(bin), m_tower(tower) {};
0052 int m_binNo;
0053 int m_tower;
0054
0055 bool operator<(const SMuonBinPos & mp1) const{
0056
0057 if (this->m_binNo != mp1.m_binNo){
0058 return (this->m_binNo < mp1.m_binNo);
0059 }
0060
0061 return this->m_tower < mp1.m_tower;
0062 };
0063 };
0064
0065 struct SMuonBin {
0066 SMuonBin(): m_rate(-1) {};
0067 PTStatistics m_stats;
0068 std::map<int, PTStatistics> m_qual2stats;
0069 long double m_rate;
0070 };
0071
0072 typedef std::map<SMuonBinPos,SMuonBin> TEffMap;
0073
0074
0075
0076
0077
0078 double rate(double x);
0079 float getPt(unsigned int ptCode);
0080 double binRate(double pt1, double pt2);
0081 int getBinNo(double pt);
0082 int getPtCode(double pt);
0083 void getBinBounds(int binNo, double & pt1, double & pt2);
0084 void initRates(TEffMap & effMap);
0085 void ratesVsPtcut(TEffMap & effMap);
0086 void rateVsTower(TEffMap & effMap,int ptCut);
0087 void draw(TEffMap & effMap, int tower1, int tower2, TCanvas & c1, int padNo=0);
0088 void actCurves(TEffMap & effMap);
0089 void drawActCur(TEffMap & effMap, int towerMin, int towerMax, TCanvas & c1, int ptCut, int padNo=0);
0090 int eta2tower(double eta);
0091 void drawKar(TEffMap & effMap, int ptCut, int option=0);
0092 void effAndPurity(TEffMap & effMap, int option, int towerBeg, int towerEnd);
0093
0094
0095
0096
0097
0098 int main(int argc, char *argv[]){
0099
0100 std::string fname;
0101 if (argc!=2){
0102 fname = "muons.txt";
0103 std::cout<< " Trying to use defualt: " << fname << std::endl;
0104 } else {
0105 fname = argv[1];
0106 std::cout<< " File: " << fname << std::endl;
0107 }
0108
0109
0110 std::ifstream fIn(fname.c_str());
0111 if (!fIn.is_open()) {
0112 std::cout << "could not open the file" << std::endl;
0113 return 1;
0114 }
0115
0116 double etaGen;
0117 double phiGen;
0118 double ptGen;
0119 int towerRec, phiRec, ptCodeRec, qual;
0120 TEffMap effMap;
0121 TEffMap effMapMC;
0122 int read = 0;
0123
0124 int toRead = -1;
0125 int ghost = 0;
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138 while( fIn >> etaGen >> phiGen >> ptGen
0139 >> towerRec >> phiRec >> ptCodeRec >> qual
0140 >> ghost)
0141
0142
0143 {
0144 if (read>toRead && toRead!=-1)
0145 break;
0146
0147 if (ghost!=0) continue;
0148 if (std::abs(etaGen)>2.1) continue;
0149 if (ptGen < 1.7) continue ;
0150
0151 ++read;
0152 if (read%100000==0){
0153 std::cout << " Read: " << read << std::endl;
0154 }
0155 int t = towerRec;
0156 if (ptCodeRec == 0) {
0157 t = eta2tower(etaGen);
0158 }
0159
0160
0161
0162
0163
0164
0165
0166 SMuonBinPos p(std::abs(t), getBinNo(ptGen) );
0167 SMuonBinPos pMC( eta2tower(etaGen) , getBinNo(ptGen) );
0168
0169 ++(effMap[p].m_qual2stats[qual].at(ptCodeRec));
0170 ++(effMapMC[pMC].m_qual2stats[qual].at(ptCodeRec));
0171 }
0172 std::cout << " Finished. Read: " << read << std::endl;
0173
0174 TFile * hfile = new TFile("plots.root","RECREATE","Demo ROOT file with histograms");
0175
0176
0177 initRates(effMap);
0178 initRates(effMapMC);
0179
0180
0181 setTDRStyle();
0182
0183
0184 ratesVsPtcut(effMap);
0185
0186
0187 rateVsTower(effMap,9);
0188 rateVsTower(effMap,13);
0189 rateVsTower(effMap,18);
0190 rateVsTower(effMap,24);
0191 rateVsTower(effMap,31);
0192
0193 actCurves(effMap);
0194
0195
0196
0197 drawKar(effMapMC, 1);
0198 drawKar(effMapMC, 5);
0199 drawKar(effMapMC, 10);
0200 drawKar(effMapMC, 15);
0201 drawKar(effMapMC, 20);
0202 drawKar(effMapMC, 25);
0203 drawKar(effMapMC, 30);
0204 drawKar(effMapMC, 31);
0205
0206
0207 drawKar(effMapMC, 1,1);
0208 drawKar(effMapMC, 5,1);
0209 drawKar(effMapMC, 10,1);
0210 drawKar(effMapMC, 15,1);
0211 drawKar(effMapMC, 20,1);
0212 drawKar(effMapMC, 25,1);
0213 drawKar(effMapMC, 30,1);
0214 drawKar(effMapMC, 31,1);
0215
0216
0217 effAndPurity(effMapMC,1,0,8);
0218 effAndPurity(effMapMC,1,9,16);
0219
0220
0221 hfile->Write();
0222
0223 return 0;
0224 }
0225
0226
0227
0228
0229
0230 void actCurves(TEffMap & effMap){
0231
0232 std::vector<int> cuts;
0233 cuts.push_back(5);
0234 cuts.push_back(10);
0235 cuts.push_back(15);
0236 cuts.push_back(20);
0237 cuts.push_back(25);
0238 cuts.push_back(29);
0239 cuts.push_back(30);
0240 cuts.push_back(31);
0241
0242 TCanvas c1("actCanvas","",600,500);
0243 std::stringstream ss;
0244 std::vector<int>::iterator itC = cuts.begin();
0245
0246 for(;itC!=cuts.end();++itC){
0247
0248
0249 c1.Clear();
0250 drawActCur(effMap, 0, 7, c1, *itC);
0251 ss.clear();
0252 ss.str("");
0253 ss << "akt" << *itC << "t0_7.pdf";
0254 c1.Print(ss.str().c_str());
0255
0256
0257
0258 c1.Clear();
0259 drawActCur(effMap, 8, 16, c1, *itC);
0260 ss.clear();
0261 ss.str("");
0262 ss << "akt" << *itC << "t8_16.pdf";
0263 c1.Print(ss.str().c_str());
0264
0265
0266 for (int tower = 0; tower < 17; ++tower){
0267 c1.Clear();
0268 drawActCur(effMap, tower, tower, c1, *itC);
0269 ss.clear();
0270 ss.str("");
0271 ss << "akt" << *itC << "t" << tower <<".pdf";
0272 c1.Print(ss.str().c_str());
0273 }
0274
0275 }
0276
0277
0278 }
0279
0280
0281
0282
0283
0284
0285 void ratesVsPtcut(TEffMap & effMap){
0286
0287
0288 TCanvas c1;
0289 draw(effMap, 0, 7, c1);
0290 c1.Print("t0_7.pdf");
0291
0292 c1.Clear();
0293 draw(effMap, 8, 16, c1);
0294 c1.Print("t8_16.pdf");
0295
0296 for (int i = 0; i<17; ++i){
0297 c1.Clear();
0298 draw(effMap, i, i, c1);
0299 std::stringstream ss;
0300 ss << "t" << i << ".pdf";
0301 c1.Print(ss.str().c_str());
0302 }
0303
0304 c1.Clear();
0305 c1.Divide(2,3);
0306 for (int i = 0; i<6; ++i){
0307 draw(effMap, i, i, c1, i+1);
0308 }
0309 c1.Print("t0t5.pdf");
0310
0311 c1.Clear();
0312 c1.Divide(2,3);
0313 for (int i = 6; i<12; ++i){
0314 draw(effMap, i, i, c1, i+1-6);
0315 }
0316 c1.Print("t6t11.pdf");
0317
0318 c1.Clear();
0319 c1.Divide(2,3);
0320 for (int i = 12; i<17; ++i){
0321 draw(effMap, i, i, c1, i+1-12);
0322 }
0323 c1.Print("t12t16.pdf");
0324
0325 return;
0326 }
0327
0328
0329
0330
0331
0332 void rateVsTower(TEffMap & effMap,int ptCodeCut){
0333
0334
0335 double pt = pts[ptCodeCut];
0336
0337
0338
0339 std::vector<double> rateVsTower(17,0);
0340
0341
0342 TEffMap::iterator it=effMap.begin();
0343 TEffMap::iterator itend=effMap.end();
0344
0345
0346
0347
0348
0349
0350
0351
0352 typedef std::vector<int> TQinTowerCnt;
0353 TQinTowerCnt qInTCnt(17,0);
0354 typedef std::vector<TQinTowerCnt> TQl2Tow;
0355 TQl2Tow tow2qlCnt(4,qInTCnt);
0356
0357 for(;it!=itend;++it){
0358 int tower = it->first.m_tower;
0359
0360 rateVsTower.at(tower)+=it->second.m_rate*it->second.m_stats.eff(ptCodeCut);
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372 std::map<int, PTStatistics>::iterator it1 = it->second.m_qual2stats.begin();
0373 std::map<int, PTStatistics>::iterator it1end = it->second.m_qual2stats.end();
0374 for (;it1!=it1end;++it1){
0375
0376 tow2qlCnt.at(it1->first).at(tower)+=it1->second.sum(ptCodeCut);
0377 }
0378 }
0379
0380 std::stringstream name;
0381 name << "aRate_ptCodeCut_" << ptCodeCut;
0382
0383 TH1D *hist = new TH1D(name.str().c_str(), "a", 17, -0.5, 16.5);
0384
0385 TCanvas c1;
0386
0387
0388 c1.SetTicks(1,1);
0389 c1.SetLogy();
0390 c1.SetGrid();
0391
0392
0393 hist->SetStats(false);
0394 hist->GetYaxis()->SetTitle("Rate [Hz]");
0395 hist->GetXaxis()->SetTitle("Tower");
0396
0397 for(unsigned int i=0; i<rateVsTower.size();++i){
0398 hist->SetBinContent(i+1,rateVsTower.at(i));
0399 }
0400 TLegend * leg = new TLegend(0.2,0.87,0.6,0.94);
0401 std::stringstream stitle;
0402 stitle << "Rate in towers with ptCode cut " << ptCodeCut << "(" << pt << "GeV/c)";
0403 leg->SetHeader(stitle.str().c_str());
0404
0405 hist->Draw();
0406 leg->Draw();
0407
0408
0409 std::stringstream ss;
0410 ss << "rateVsTow_cut" << ptCodeCut << ".pdf";
0411 c1.Print(ss.str().c_str());
0412
0413 delete leg;
0414
0415
0416 TQinTowerCnt sum(17,0);
0417 for (unsigned int qual = 0; qual< tow2qlCnt.size();++qual){
0418 for (unsigned int tower = 0; tower< tow2qlCnt.at(qual).size();++tower){
0419 sum.at(tower)+=tow2qlCnt.at(qual).at(tower);
0420 }
0421 }
0422
0423 THStack *hs = new THStack(ss.str().c_str(),"blabbla");
0424
0425
0426 leg = new TLegend(0.2,0.6,0.4,0.94);
0427
0428 leg->SetHeader("Jakosc:");
0429
0430
0431 for (unsigned int qual = 0; qual< tow2qlCnt.size();++qual){
0432
0433 std::stringstream name;
0434 name << "Rate_ptCodeCut_" << ptCodeCut << "_qual_" << qual;
0435 TH1D *hist = new TH1D(name.str().c_str(), "a", 17, -0.5, 16.5);
0436 bool empty = true;
0437 for (unsigned int tower = 0; tower< tow2qlCnt.at(qual).size();++tower){
0438 if (sum.at(tower)!=0){
0439 double bc = double(tow2qlCnt.at(qual).at(tower))/sum.at(tower)*rateVsTower.at(tower);
0440 hist->SetBinContent(tower+1,bc);
0441 empty=false;
0442 }
0443 }
0444 if (!empty){
0445 hist->SetFillColor(qual+1);
0446 hs->Add(hist);
0447 std::stringstream lname;
0448 lname << qual;
0449 leg->AddEntry(hist,lname.str().c_str(),"f");
0450 }
0451
0452 }
0453
0454
0455 c1.Clear();
0456 hs->Draw();
0457 leg->Draw();
0458 hs->GetYaxis()->SetTitle("Czestosc [Hz]");
0459 hs->GetXaxis()->SetTitle("Wieza");
0460 std::stringstream ssa;
0461 ssa << "crateVsTow_cut" << ptCodeCut << ".pdf";
0462 c1.Print(ssa.str().c_str());
0463
0464 return;
0465 }
0466
0467
0468
0469
0470
0471 #include <cmath>
0472 double rate(double x){
0473
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488 const double lum = 2.0e33;
0489 const double dabseta = 1.0;
0490 const double dpt = 1.0;
0491 const double afactor = 1.0e-34*lum*dabseta*dpt;
0492 const double a = 2*1.3084E6;
0493 const double mu=-0.725;
0494 const double sigma=0.4333;
0495 const double s2=2*sigma*sigma;
0496 double ptlog10;
0497 ptlog10 = std::log10(x);
0498 double ex = (ptlog10-mu)*(ptlog10-mu)/s2;
0499 double rate = (a * exp(-ex) * afactor);
0500
0501 return rate;
0502
0503 }
0504
0505
0506
0507
0508
0509 double binRate(double pt1, double pt2, int tower){
0510
0511
0512
0513 tower = std::abs(tower);
0514 if (tower > 16){
0515 std::cout << "Bad tower " << tower << std::endl;
0516 throw "Bad tower";
0517 }
0518
0519
0520 double etaTowerSize = etas[tower+1]-etas[tower];
0521 double detEta = 2.1;
0522
0523
0524 double ret =( rate(pt2)-rate(pt1) ) /2/(pt2-pt1) ;
0525
0526 if (ret<0){
0527 ret = -ret;
0528 }
0529
0530 return ret*etaTowerSize/detEta;
0531
0532 }
0533
0534
0535
0536
0537
0538
0539
0540
0541 int getBinNo(double pt){
0542
0543
0544
0545 int ret = 0;
0546 int nbins = 1000;
0547
0548 double ptLow = 1.49;
0549
0550 double ptHigh = 200.;
0551 double delta = ptHigh-ptLow;
0552 double binSize = delta/nbins;
0553
0554 while (ret < nbins){
0555
0556 double lowEdge = ptLow+ret*binSize;
0557 if (lowEdge>pt)
0558 break;
0559
0560 ++ret;
0561 }
0562
0563 return ret;
0564 }
0565
0566
0567
0568
0569
0570
0571
0572
0573 void getBinBounds(int binNo, double & pt1, double & pt2){
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585 int nbins = 1000;
0586
0587 double ptLow = 1.49;
0588 double ptHigh = 200.;
0589 double delta = ptHigh-ptLow;
0590 double binSize = delta/nbins;
0591 pt1 = ptLow+binNo*binSize;
0592 pt2 = ptLow+(binNo+1)*binSize;
0593
0594 }
0595
0596
0597
0598
0599
0600 void initRates(TEffMap & effMap){
0601
0602 TEffMap::iterator it=effMap.begin();
0603 TEffMap::iterator itend=effMap.end();
0604
0605 for(;it!=itend;++it){
0606 double p1 = -1;
0607 double p2 = -1;
0608 getBinBounds(it->first.m_binNo,p1,p2);
0609 it->second.m_rate = binRate(p1,p2,it->first.m_tower);
0610
0611
0612
0613
0614
0615
0616
0617 std::map<int, PTStatistics>::iterator it1 = it->second.m_qual2stats.begin();
0618 std::map<int, PTStatistics>::iterator it1end = it->second.m_qual2stats.end();
0619 for (;it1!=it1end;++it1){
0620 it->second.m_stats.update(it1->second);
0621 }
0622
0623
0624
0625 }
0626
0627
0628 }
0629
0630
0631
0632
0633
0634 void draw(TEffMap & effMap, int towerMin, int towerMax, TCanvas & c1, int padNo){
0635
0636
0637
0638 std::vector<double> rateVsPtCut(32,0);
0639
0640 TEffMap::iterator it=effMap.begin();
0641 TEffMap::iterator itend=effMap.end();
0642
0643 for(;it!=itend;++it){
0644 int tower = it->first.m_tower;
0645 if ( tower < towerMin || tower > towerMax){
0646 continue;
0647 }
0648 for(int i=0;i<32;++i){
0649 rateVsPtCut.at(i)+=it->second.m_rate*it->second.m_stats.eff(i);
0650 }
0651
0652 }
0653
0654
0655 std::stringstream ss,fname;
0656
0657 if (towerMin==towerMax){
0658 ss << "Tower: " << towerMin;
0659 } else {
0660 ss << "Towers: " << towerMin << "..." << towerMax;
0661 }
0662
0663
0664
0665 c1.cd(padNo);
0666 c1.SetTicks(1,1);
0667 c1.SetLogy();
0668 c1.SetGrid();
0669
0670 fname << "RateInT" << towerMin << "_" << towerMax << ".pdf";
0671 TH1D *hist = new TH1D(fname.str().c_str(), ss.str().c_str(), 31, 0.5, 31.5);
0672 hist->SetStats(false);
0673 hist->GetYaxis()->SetTitle("czestosc [Hz]");
0674 hist->GetXaxis()->SetTitle("ptCodeCut");
0675
0676 for(unsigned int i=1; i<rateVsPtCut.size();++i){
0677
0678 hist->SetBinContent(i,rateVsPtCut.at(i));
0679 }
0680
0681 TLegend * leg = new TLegend(0.8,0.88,0.97,0.93);
0682 leg->SetHeader(ss.str().c_str());
0683
0684
0685 hist->Draw();
0686 leg->Draw();
0687
0688
0689
0690 }
0691
0692
0693
0694
0695
0696 #include "TAxis.h"
0697 void drawActCur(TEffMap & effMap, int towerMin, int towerMax, TCanvas & c1, int ptCut, int padNo){
0698
0699
0700
0701
0702 int bins = getBinNo(10000000)+1;
0703
0704
0705
0706 std::vector<int> rec(bins,0);
0707 std::vector<int> all(bins,0);
0708
0709
0710
0711 TEffMap::iterator it=effMap.begin();
0712 TEffMap::iterator itend=effMap.end();
0713
0714 for(;it!=itend;++it){
0715 int tower = it->first.m_tower;
0716 if ( tower < towerMin || tower > towerMax){
0717 continue;
0718 }
0719 int bin = it->first.m_binNo;
0720
0721
0722 rec.at(bin)+=it->second.m_stats.sum(ptCut);
0723 all.at(bin)+=it->second.m_stats.sum();
0724 }
0725
0726
0727 std::stringstream ss,fname;
0728
0729 if (towerMin==towerMax){
0730 ss << "Tower: " << towerMin;
0731 } else {
0732 ss << "Towers: " << towerMin << "..." << towerMax;
0733 }
0734
0735 c1.cd(padNo);
0736 c1.SetTicks(1,1);
0737
0738 c1.SetLogx();
0739 c1.SetGrid();
0740
0741 fname << "akt" << ptCut << "inT" << towerMin << "_" << towerMax << ".pdf";
0742
0743
0744
0745
0746
0747
0748 double xbins[1000];
0749 double ptLow=0, ptHigh=0;
0750 for(int i = 0; i < bins; ++i){
0751 getBinBounds(i,ptLow,ptHigh);
0752 xbins[i]=ptLow;
0753 }
0754 xbins[bins]=ptHigh;
0755
0756 TH1D *hist = new TH1D(fname.str().c_str(), ss.str().c_str(), bins, xbins);
0757 hist->SetStats(false);
0758 hist->GetYaxis()->SetTitle("efektywnosc");
0759 hist->GetXaxis()->SetTitle("pT [GeV/c]");
0760 hist->SetAxisRange(1.5,200.);
0761
0762
0763 for(unsigned int i=0; i<all.size();++i){
0764
0765
0766
0767
0768 int div = all.at(i);
0769 if (div > 0){
0770 hist->SetBinContent(i+1,double(rec.at(i))/div);
0771 }
0772 }
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782
0783
0784
0785
0786 std::stringstream fnew;
0787 fnew << "aktR" << ptCut << "inT" << towerMin << "_" << towerMax << ".pdf";
0788
0789
0790 hist->Draw();
0791 hist ->GetYaxis()->SetRangeUser(0.,1.);
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802
0803
0804
0805
0806
0807 }
0808
0809
0810
0811
0812
0813 void drawKar(TEffMap & effMap, int ptCut, int option){
0814
0815 static TEffMap newmap;
0816 static bool firstRun = true;
0817
0818 if (firstRun){
0819
0820 firstRun = false;
0821 TEffMap::iterator it=effMap.begin();
0822 TEffMap::iterator itend=effMap.end();
0823 int tower = 0, bin = 0;
0824 double r1 = 0, r2 = 0;
0825 for(;it!=itend;++it){
0826 tower = it->first.m_tower;
0827 bin = it->first.m_binNo;
0828 getBinBounds(bin, r1, r2);
0829 int ptCode = getPtCode( (r1+r2)/2 );
0830
0831
0832
0833
0834
0835 SMuonBinPos p(tower, ptCode );
0836
0837 newmap[p].m_stats.update(it->second.m_stats);
0838 }
0839 }
0840
0841
0842 std::stringstream sname;
0843 if (option == 0){
0844 sname << "effkar" << ptCut << ".pdf";
0845 } else {
0846 sname << "effkarTT" << ptCut << ".pdf";
0847 }
0848
0849 int m_TOWER_COUNT = 17, m_PT_CODE_MAX = 31;
0850 TH2D *eff = new TH2D( sname.str().c_str(), "",
0851 m_TOWER_COUNT, -0.5, m_TOWER_COUNT-0.5,
0852 m_PT_CODE_MAX+1, -0.5, m_PT_CODE_MAX+0.5);
0853
0854 TEffMap::iterator it=newmap.begin();
0855 TEffMap::iterator itend=newmap.end();
0856 for(;it!=itend;++it){
0857
0858
0859
0860
0861
0862 eff->SetBinContent(it->first.m_tower+1, it->first.m_binNo+1,it->second.m_stats.eff(ptCut));
0863 }
0864
0865
0866
0867 TCanvas c1;
0868 gStyle->SetPalette(1,0);
0869 eff->SetStats(false);
0870 if (option == 0){
0871 eff->Draw("COLZ");
0872 } else {
0873 eff->Draw("TEXT");
0874 }
0875 eff-> SetMaximum(1);
0876 eff-> SetMinimum(0);
0877 eff->GetXaxis()->SetTitle("tower generated");
0878 eff->GetYaxis()->SetTitle("ptCode generated");
0879 c1.Print(sname.str().c_str());
0880
0881
0882 }
0883
0884
0885
0886
0887
0888 int eta2tower(double eta){
0889
0890 bool neg = (eta < 0 ? true : false);
0891 if (neg) {
0892 eta = - eta;
0893 }
0894 if (eta < 0){
0895 std::cout << "Eta neg " << eta << std::endl;
0896 throw "Eta negative";
0897 }
0898
0899 if (eta > 2.101) {
0900 std::cout << "Eta big " << eta << std::endl;
0901 throw "Eta to big";
0902 }
0903
0904
0905 int tower = 0;
0906 while (tower < 16 && eta > etas[tower+1] ) ++tower;
0907
0908 return (neg ? -tower : tower);
0909
0910 }
0911
0912
0913
0914
0915
0916 int getPtCode(double pt){
0917
0918 int ret = 0;
0919 while ( pt>pts[ret+1] && ret < 31) ++ret;
0920
0921 return ret;
0922 }
0923
0924 float getPt(unsigned int ptCode){
0925
0926 if (ptCode > 31){
0927 std::cout << "Ptcode big " << ptCode << std::endl;
0928 throw "Bad ptCode";
0929 }
0930
0931 return pts[ptCode];
0932 }
0933
0934
0935
0936
0937
0938
0939
0940 void effAndPurity(TEffMap & effMap, int option, int towerBeg, int towerEnd){
0941
0942 TEffMap::iterator it=effMap.begin();
0943 TEffMap::iterator itend=effMap.end();
0944
0945
0946
0947 typedef std::vector<double> TPtCode2Rate;
0948 typedef std::vector<TPtCode2Rate> TTower2PtCode;
0949
0950 TPtCode2Rate ini(32,0);
0951 TTower2PtCode accepted(17,ini);
0952 TTower2PtCode shouldBeAccepted(17,ini);
0953 TTower2PtCode wasAndShouldBeAccepted(17,ini);
0954
0955 float ptCut = 0;
0956 double pt1 = 0, pt2 = 0;
0957 int tower = -1;
0958 for(;it!=itend;++it){
0959
0960 for (int ptCodeCut = 1; ptCodeCut < 32; ++ptCodeCut){
0961 ptCut = getPt(ptCodeCut);
0962 tower = std::abs(it->first.m_tower);
0963 getBinBounds(it->first.m_binNo,pt1,pt2);
0964
0965 double aRate = 0;
0966 int sum = it->second.m_stats.sum();
0967 if (sum != 0){
0968 aRate = (double)it->second.m_stats.sum(ptCodeCut);
0969 aRate/=sum;
0970 aRate*=it->second.m_rate;
0971 accepted[tower][ptCodeCut]+=aRate;
0972
0973 }
0974
0975
0976 if ( (pt1+pt2)/2 > ptCut){
0977 shouldBeAccepted[tower][ptCodeCut]+=it->second.m_rate;
0978 wasAndShouldBeAccepted[tower][ptCodeCut]+=aRate;
0979
0980
0981 }
0982 }
0983
0984
0985 }
0986
0987
0988
0989 std::stringstream pname;
0990 if (option == 0){
0991 pname << "pur"<< towerBeg << "_" << towerEnd << ".pdf";
0992 } else {
0993 pname << "purTT"<< towerBeg << "_" << towerEnd << ".pdf";;
0994 }
0995
0996 std::stringstream ename;
0997 if (option == 0){
0998 ename << "eff"<< towerBeg << "_" << towerEnd << ".pdf";;
0999 } else {
1000 ename << "effTT"<< towerBeg << "_" << towerEnd << ".pdf";;
1001 }
1002
1003
1004
1005 int m_TOWER_COUNT = towerEnd-towerBeg+1;
1006 int m_PT_CODE_MAX = 31;
1007 int ptCodeMin = 5;
1008 TH2D *eff = new TH2D( ename.str().c_str(), "",
1009 m_TOWER_COUNT, towerBeg-0.5, towerEnd+0.5,
1010 m_PT_CODE_MAX+1-ptCodeMin, ptCodeMin-0.5, m_PT_CODE_MAX+0.5);
1011
1012 TH2D *pur = new TH2D( pname.str().c_str(), "",
1013 m_TOWER_COUNT, towerBeg-0.5, towerEnd+0.5,
1014 m_PT_CODE_MAX+1-ptCodeMin, ptCodeMin-0.5, m_PT_CODE_MAX+0.5);
1015
1016
1017
1018 double div = 0;
1019 for (int ptCodeCut = ptCodeMin; ptCodeCut < 32; ++ptCodeCut){
1020 for (int tower = towerBeg; tower <= towerEnd; ++tower){
1021
1022 double dEff = -1;
1023
1024 div = shouldBeAccepted[tower][ptCodeCut];
1025 if (div > 0) {
1026 dEff = wasAndShouldBeAccepted[tower][ptCodeCut];
1027 dEff /= div;
1028 }
1029
1030 double dPur = -1;
1031
1032 div = accepted[tower][ptCodeCut];
1033 if (div > 0){
1034 dPur = wasAndShouldBeAccepted[tower][ptCodeCut];
1035 dPur /= div;
1036 }
1037
1038 eff->SetBinContent(tower-towerBeg+1, ptCodeCut+1-ptCodeMin, dEff);
1039 pur->SetBinContent(tower-towerBeg+1, ptCodeCut+1-ptCodeMin, dPur);
1040
1041 }
1042
1043 }
1044
1045
1046 TCanvas c1;
1047 gStyle->SetPalette(1,0);
1048 eff->SetStats(false);
1049 if (option == 0){
1050 eff->Draw("COLZ");
1051 } else {
1052 eff->Draw("TEXT");
1053 }
1054 eff-> SetMaximum(1);
1055 eff-> SetMinimum(0);
1056 eff-> SetMarkerSize(2);
1057 eff->GetXaxis()->SetTitle("Wieza");
1058 eff->GetXaxis()->SetTickLength(0);
1059 eff->GetYaxis()->SetTitle("ptCodeCut");
1060 eff->GetYaxis()->SetTickLength(0.01);
1061 c1.Print(ename.str().c_str());
1062
1063 TCanvas c2("purCanv", "", 800, 500);
1064
1065 c2.SetGridy();
1066 gStyle->SetPalette(1,0);
1067 pur->SetStats(false);
1068 if (option == 0){
1069 pur->Draw("COLZ");
1070 } else {
1071 pur->Draw("TEXT");
1072 }
1073 pur-> SetMaximum(1);
1074 pur-> SetMinimum(0);
1075 pur-> SetMarkerSize(1.7);
1076
1077 pur->GetXaxis()->SetTitle("Wieza");
1078 pur->GetXaxis()->SetTickLength(0);
1079 pur->GetYaxis()->SetTitle("ptCodeCut");
1080 pur->GetYaxis()->SetTickLength(0.01);
1081 c2.Print(pname.str().c_str());
1082
1083
1084 }