Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:55 

0001 /*

0002 *

0003 *    Tomasz Fruboes,  IPJ Warsaw 2007

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 * Ugly global data definitions

0025 *

0026 **************************************************/
0027 //                     0     1     2     3     4      5

0028 double etas[] = {-0.07, 0.07, 0.27, 0.44, 0.58, 0.72,
0029    //                  6     7     8     9     10     11 

0030                   0.83, 0.93, 1.04, 1.14, 1.24, 1.36,
0031                   //   12    13    14    15    16                

0032                  1.48,  1.61, 1.73, 1.85, 1.97, 2.10};
0033 
0034 double pts[] = {
0035    0.0,  0.01,    
0036 //1    2     3    4     5      6     7    8

0037    1.5,  2.0, 2.5,  3.0,  3.5,  4.0,  4.5, 
0038 //    9     10    11    12

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 * Data structures filled when reading muons.txt

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; // m_stats = sum of m_qual2stats->seconds 

0069    long double m_rate;
0070 };
0071 
0072 typedef std::map<SMuonBinPos,SMuonBin> TEffMap;
0073 /*************************************************

0074 *

0075 * Fwd definitions

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 * Read the data, call functions for drawing

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;//, muonsFound;

0120    TEffMap effMap;
0121    TEffMap effMapMC;
0122    int read = 0;
0123    //int toRead = 1000000;

0124    int toRead = -1;
0125    int ghost = 0;
0126 //   while (fIn >> etaGen >> phiGen >> ptGen

0127 //              >> towerRec >> phiRec >> ptCodeRec >> qual 

0128 //              >> muonsFound)

0129 
0130   // 0.453722 -2.06621 53.1866 3 95 31 3 0

0131 
0132  /*                        m_outfileR << etaGen << " " << phiGen << " " << ptGen

0133                                     << " "    << towerRec << " " << phiRec << " " << ptCodeRec << " " << qual

0134                                     << " "    << ghost

0135                                     << std::endl;

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; // TODO: add ghost histo

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) { // avoid putting events without muons in tower = 0

0157          t = eta2tower(etaGen);
0158       }
0159       /* 

0160       if (read < 1000) {

0161        std::cout << "p "<<  t << " " << std::abs(t) << " " << " " << ptGen << " " << getBinNo(ptGen) << std::endl;

0162        std::cout << "pMC "<< etaGen << " " << eta2tower(etaGen) << " " << ptGen << " " << getBinNo(ptGen) << std::endl;

0163        std::cout << " qual "  << qual << " ptCodeRec " << ptCodeRec << std::endl;

0164       }*/
0165 
0166       SMuonBinPos p(std::abs(t), getBinNo(ptGen) );
0167       SMuonBinPos pMC( eta2tower(etaGen) , getBinNo(ptGen) );
0168    //   ++(effMap[p].m_stats.at(ptCodeRec)); // ptCodeRec==0 - no muon found

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    // initilize some data after readout 

0177    initRates(effMap);
0178    initRates(effMapMC);
0179 
0180    // set style for plots

0181    setTDRStyle();   
0182 
0183    // Draw rate Vs pt cut plots. Plenty of plots 

0184    ratesVsPtcut(effMap);
0185    
0186    // draw rate vs tower for variuos cuts (9,13,...) 

0187    rateVsTower(effMap,9);
0188    rateVsTower(effMap,13);
0189    rateVsTower(effMap,18);
0190    rateVsTower(effMap,24);
0191    rateVsTower(effMap,31);
0192    // draw efficiency for various cuts of towers. Plenty of plots

0193    actCurves(effMap);
0194  
0195    // draw efficiency for different cuts 

0196    //  (note 1 means no cut) 

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    // Same as above, draw with "text" option

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    // effciency and purity

0217    effAndPurity(effMapMC,1,0,8); // barrel

0218    effAndPurity(effMapMC,1,9,16);// endacp

0219 
0220    // save to root file 

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       // Barrel

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       // Endcap

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       // All towers    

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 * Rates vs ptCut  

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 * Rates vs tower

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    typedef std::vector<int> TQlCnt; 

0347    TQlCnt qlCnt(8,0); // 8 qualities max

0348    typedef std::vector<TQlCnt> TTow2QlCnt; 

0349    TTow2QlCnt tow2qlCnt(17,qlCnt);

0350    */
0351    
0352    typedef std::vector<int> TQinTowerCnt; // count occurences of given quality in towers

0353    TQinTowerCnt qInTCnt(17,0); 
0354    typedef std::vector<TQinTowerCnt> TQl2Tow; 
0355    TQl2Tow tow2qlCnt(4,qInTCnt);// 4 qualities max

0356    
0357    for(;it!=itend;++it){
0358       int tower = it->first.m_tower;
0359       //      double oldR = rateVsTower.at(tower);

0360       rateVsTower.at(tower)+=it->second.m_rate*it->second.m_stats.eff(ptCodeCut);
0361       /*

0362       if (rateVsTower.at(tower)/oldR > 10 ) {

0363            std::cout << oldR << " " << rateVsTower.at(tower) 

0364                      << " m_rate " << it->second.m_rate 

0365                      << " eff " << it->second.m_stats.eff(ptCodeCut)

0366                      << " m_binNo " << it->first.m_binNo

0367                      << " m_tower " << it->first.m_tower

0368                      << std::endl;

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          //tow2qlCnt.at(tower).at(it1->first)+=it1->second.sum(ptCodeCut);

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    //c1.cd(padNo);

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    //delete hist;

0413    delete leg;
0414    
0415    // draw stacked rates with qualities shown

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    //leg = new TLegend(0.4,0.6,0.89,0.89);

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 //   for (unsigned int qual = tow2qlCnt.size()-1; qual >= 0;--qual){

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 * Muon rate vs ptCut in whole detector region

0469 *

0470 **************************************************/
0471 #include <cmath>
0472 double rate(double x){
0473 
0474 
0475    /*

0476    double ret = 0;

0477    double a = -0.235801;

0478    double b = -2.82346;

0479    double c = 17.162;

0480 

0481    //f1(x) =  x**(a*log(x)) *(x**b)*exp(c)

0482   //ret = std::pow( x,a*std::log(x) ) * std::pow(x,b)*std::exp(c)l

0483    ret = std::pow( x,a*std::log(x) ) * std::pow(x,b)*std::exp(c);

0484    return ret;

0485 */
0486 
0487 
0488  const double lum = 2.0e33; //defoult is 1.0e34;

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 * Int. rate  between pt1 and pt2 for given tower

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    //double ret = rate(pt2)-rate(pt1);

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 * Converts pt into bin number. 

0536 *

0537 *     -> Should be method of SMuonBinPos

0538 *     -> getBinBounds() has the same internal data 

0539 *

0540 **************************************************/
0541 int  getBinNo(double pt){
0542    
0543    //return getPtCode(pt);

0544    
0545    int ret = 0;
0546    int nbins = 1000;
0547    //int nbins = 200;

0548    double ptLow = 1.49;
0549    //double ptHigh = 201.;

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 * Converts bin number into pt range.

0568 *

0569 *     -> Should be method of SMuonBinPos

0570 *     -> getBinNo() has the same internal data

0571 *

0572 **************************************************/
0573 void getBinBounds(int binNo, double & pt1, double & pt2){
0574 
0575    /*

0576    pt1 = getPt(binNo);

0577    if(binNo == 31) 

0578       pt2 = 200.;

0579    else  

0580       pt2 = getPt(binNo+1);

0581          

0582    return;      

0583    */
0584 
0585    int nbins = 1000;
0586    //int nbins = 200;

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 * Initilize data after readout

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       if (it->second.m_rate > 1000000) {

0612         std::cout << " initRates "<< p1 << " " << p2 << " " << it->first.m_tower << " " << binRate(p1,p2,it->first.m_tower) << std::endl; 

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 *  Draws rate vs ptCut hist

0632 *

0633 *******************************************************************************************/
0634 void draw(TEffMap & effMap, int towerMin, int towerMax, TCanvas  & c1, int padNo){
0635 
0636    //std::cout << "D " << towerMin << " " << towerMax << std::endl;

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){ // `i` is ptCut

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       //std::cout << i << " " << rateVsPtCut.at(i) << std::endl;

0678       hist->SetBinContent(i,rateVsPtCut.at(i));
0679    }
0680    //TLegend * leg = new TLegend(0.84,0.88,0.97,0.93);

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    //delete hist;

0689    //delete leg;

0690 }
0691 /*******************************************************************************************

0692 *

0693 *  Draws effIciency vs ptGen for various ptCuts

0694 *

0695 *******************************************************************************************/
0696 #include "TAxis.h"
0697 void drawActCur(TEffMap & effMap, int towerMin, int towerMax, TCanvas  & c1, int ptCut, int padNo){
0698 
0699    
0700    //std::cout << "Akt " << towerMin << " " << towerMax << std::endl;

0701 
0702    int bins = getBinNo(10000000)+1; // returns max_bin bin no 0 is the first one

0703    //std::vector<double> effVsPt(bins,0);

0704    //std::vector<int> cnt(bins,0);

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       //effVsPt.at(bin)+=it->second.m_stats.eff(ptCut);

0721       //++cnt.at(bin);

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 //   c1.SetLogy();

0738    c1.SetLogx();
0739    c1.SetGrid();
0740 
0741    fname << "akt" <<  ptCut << "inT" << towerMin << "_" << towerMax << ".pdf";
0742    
0743    //double ptLow = -1, ptHigh = -1, trash = -1;

0744    //getBinBounds(0,ptLow,trash);

0745    //getBinBounds(bins-1,trash,ptHigh);

0746    //TH1D *hist = new TH1D(fname.str().c_str(), ss.str().c_str(), bins, ptLow, ptHigh);

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; // Fill the high edge of last bin

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    //for(unsigned int i=0; i<effVsPt.size();++i){

0763    for(unsigned int i=0; i<all.size();++i){
0764       /*

0765       if(cnt.at(i)!=0){

0766          hist->SetBinContent(i+1,effVsPt.at(i)/cnt.at(i));

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    double pts[30] = {

0776 //      0.0,  0.01,    

0777 //      1.5,  2.0, 2.5,  3.0,  3.5,  4.0,  4.5, 

0778       2.0, 2.5,  3.0,  3.5,  4.0,  4.5,  // 6

0779       5.,   6.,   7.,   8.,   // 4

0780       10.,  12., 14.,  16.,  18., // 5 

0781       20.,  25.,  30., 35.,  40.,  45., // 5

0782       50.,  60.,  70., 80.,  90.,  100., 120., 140., 200.};//9

0783 

0784    */
0785 
0786    std::stringstream fnew;
0787    fnew << "aktR" <<  ptCut << "inT" << towerMin << "_" << towerMax << ".pdf";
0788    //hist = (TH1D *)hist->Rebin(5, fnew.str().c_str());

0789    //hist->Scale(0.2);

0790    hist->Draw();
0791    hist ->GetYaxis()->SetRangeUser(0.,1.);
0792 
0793    /*

0794 //   hist ->GetXaxis()->SetRangeUser(2.0,200.0);

0795    //*/
0796 
0797    
0798    
0799    //TLegend * leg = new TLegend(0.84,0.88,0.97,0.93);

0800    //TLegend * leg = new TLegend(0.8,0.88,0.97,0.93);

0801    //leg->SetHeader(ss.str().c_str());

0802    //leg->Draw();

0803    
0804    
0805    //delete hist;

0806    //delete leg;

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       //std::cout << " FIRST run" << std::endl;

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          /*std::cout << r1 

0831                    << " " << r2

0832                    << " " << (r1+r2)/2

0833                    << " " << ptCode 

0834                    << std::endl;*/
0835          SMuonBinPos p(tower, ptCode );
0836          // copy the statistics into proper bin of newmap 

0837          newmap[p].m_stats.update(it->second.m_stats);
0838       }
0839    }
0840    
0841    // Fill the histogram

0842    std::stringstream sname;
0843    if (option == 0){
0844       sname << "effkar" << ptCut << ".pdf";
0845    } else { // (option == 1){

0846       sname << "effkarTT" << ptCut << ".pdf";
0847    }
0848    
0849    int m_TOWER_COUNT = 17, m_PT_CODE_MAX = 31; // From RPCConst

0850    TH2D *eff = new TH2D( sname.str().c_str(), "",
0851                          m_TOWER_COUNT, -0.5, m_TOWER_COUNT-0.5, //Towers

0852                          m_PT_CODE_MAX+1, -0.5, m_PT_CODE_MAX+0.5);// Pt's

0853 
0854    TEffMap::iterator it=newmap.begin();
0855    TEffMap::iterator itend=newmap.end();
0856    for(;it!=itend;++it){
0857       // XXX Tower may be negative...

0858       /*std::cout << it->first.m_tower

0859                 << " " << it->first.m_binNo

0860                 << " " <<it->second.m_stats.eff(ptCut)

0861                 << std::endl;*/
0862       eff->SetBinContent(it->first.m_tower+1, it->first.m_binNo+1,it->second.m_stats.eff(ptCut));
0863    }
0864 
0865    // draw it and save to file

0866    //TCanvas c1("canbaaas","",800,600);

0867    TCanvas c1;
0868    gStyle->SetPalette(1,0);
0869    eff->SetStats(false);
0870    if (option == 0){
0871      eff->Draw("COLZ");
0872    } else { // (option == 1){

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 * Draw efficiency and purity. Both use rate() for 

0936 *  weighting. Since rate goes exp. down with pt

0937 *  it may be wrong for efficiency plot 

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    // Rate wieghted calculations

0947    typedef std::vector<double> TPtCode2Rate;
0948    typedef std::vector<TPtCode2Rate> TTower2PtCode;
0949    
0950    TPtCode2Rate ini(32,0); // 31 ptCodes

0951    TTower2PtCode accepted(17,ini); // 17 towers  // #1

0952    TTower2PtCode shouldBeAccepted(17,ini); // 17 towers // #2

0953    TTower2PtCode wasAndShouldBeAccepted(17,ini); // 17 towers  // #3

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 { // (option == 1){

0993       pname << "purTT"<< towerBeg << "_" << towerEnd << ".pdf";;
0994    }
0995    
0996    std::stringstream ename;
0997    if (option == 0){
0998       ename << "eff"<< towerBeg << "_" << towerEnd << ".pdf";;
0999    } else { // (option == 1){

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; // From RPCConst

1007    int ptCodeMin = 5;
1008    TH2D *eff = new TH2D( ename.str().c_str(), "",
1009                          m_TOWER_COUNT, towerBeg-0.5, towerEnd+0.5, //Towers

1010                          m_PT_CODE_MAX+1-ptCodeMin, ptCodeMin-0.5, m_PT_CODE_MAX+0.5);// Pt's

1011    
1012    TH2D *pur = new TH2D( pname.str().c_str(), "",
1013                          m_TOWER_COUNT, towerBeg-0.5, towerEnd+0.5, //Towers

1014                          m_PT_CODE_MAX+1-ptCodeMin, ptCodeMin-0.5, m_PT_CODE_MAX+0.5);// Pt's

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          // use unweighted values for efficiency calculations

1024          div = shouldBeAccepted[tower][ptCodeCut];
1025          if (div > 0) {
1026            dEff = wasAndShouldBeAccepted[tower][ptCodeCut];
1027            dEff /= div;
1028          } 
1029          
1030          double dPur = -1;
1031          // use weighted values for efficiency calculations

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    //TCanvas c1("effCanv", "", 800, 500);

1046    TCanvas c1;
1047    gStyle->SetPalette(1,0);
1048    eff->SetStats(false);
1049    if (option == 0){
1050       eff->Draw("COLZ");
1051    } else { // (option == 1){

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    //TCanvas c2;

1065    c2.SetGridy();
1066    gStyle->SetPalette(1,0);
1067    pur->SetStats(false);
1068    if (option == 0){
1069       pur->Draw("COLZ");
1070    } else { // (option == 1){

1071       pur->Draw("TEXT");
1072    }
1073    pur-> SetMaximum(1);
1074    pur-> SetMinimum(0);
1075    pur-> SetMarkerSize(1.7);
1076    //pur-> Scale(1000);

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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!