flavour
Classes | Functions
draw.cpp File Reference
#include "MCMC.h"
#include "BGL.h"
#include "TF2.h"
#include "TProfile3D.h"
#include "THStack.h"
#include "TColor.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TGraph.h"
#include "TLatex.h"
#include "TFile.h"
#include "TH2F.h"
#include "TVector.h"
#include "TCanvas.h"
#include "TMath.h"
#include <iostream>
#include <fstream>
#include <cln/cln.h>
#include <cln/float.h>
Include dependency graph for draw.cpp:

Go to the source code of this file.

Classes

class  BGL2
 A second implementation of the BGL model, for testing purposes. More...
 

Functions

int main (int argc, char *argv[])
 the main function takes the arguments inputfile gL gQ lup qup which specify the file containing the simulation results for a BGL model and draws the plots for that model More...
 

Function Documentation

int main ( int  argc,
char *  argv[] 
)

the main function takes the arguments inputfile gL gQ lup qup which specify the file containing the simulation results for a BGL model and draws the plots for that model

Definition at line 357 of file draw.cpp.

References BGLmodels::Vud().

357  {
358  // Check the number of parameters
359 
360  if(argc<6){
361  std::cerr<<"Usage: "<<argv[0]<<" inputfile gL gQ lup qup"<<std::endl;
362  return 1;}
363  CD cmu=conj(Vud[1][0])*Vud[1][1];
364  CD umu=conj(Vud[0][0])*Vud[0][1];
365 
366 
367  cout<<"RATIO "<<cmu*cmu<<endl;
368  cout<<umu*umu<<endl;
369 
370  int gL=atoi(argv[2]);
371  int gQ=atoi(argv[3]);
372  int lup=atoi(argv[4]);
373  int qup=atoi(argv[5]);
374  char name[5]="0000";
375 
376  name[0]+=gL;
377  name[1]+=gQ;
378  name[2]+=lup;
379  name[3]+=qup;
380  string ll[2][3]={{"#nu_{1}","#nu_{2}","#nu_{3}"},{"e","#mu","#tau"}};
381  string qq[2][3]={{"u","c","t"},{"d","s","b"}};
382  //Int_t MyPalette[100];
383  Double_t r[] = {1, 0.3};
384  Double_t g[] = {1, 0.3};
385  Double_t b[] = {1, 0.3};
386  Double_t stop[] = {0., 1.0};
387  TColor::CreateGradientColorTable(2, stop, r, g,b, 100);
388  //TH1F * pdf1=new TH1F("pdf1","pdf1",npoints,10,500);
389  //TGraph * chi2=new TGraph(npoints);
390 
391  uint npoints=200;
392  double init1=-3, final1=3;
393  double init2=10, final2=1000;
394  double initBmumu=0, finalBmumu=3;
395  double initBsmumu=0, finalBsmumu=6;
396 
397  double llmax=-1000,McHmax=1000,MRmax=1000,MImax=1000,tbmax=1;
398 
399  TFile *f=new TFile(argv[1],"read");
400  if(!f->IsOpen()) cout<<"NOFILE"<<endl;
401  //f->ShowStreamerInfo();
402 
403  TH2F *limits4,*Bmumu_Bsmumu,*limits_tb_MR,*limits_tb_MI;
404  TH2F *limits_MR_MI,*limits_MR_McH,*limits_MI_McH;
405 
406  f->GetObject("limits4;1",limits4);
407  f->GetObject("Bmumu_Bsmumu;1",Bmumu_Bsmumu);
408  f->GetObject("limits_tb_MR;1",limits_tb_MR);
409  f->GetObject("limits_tb_MI;1",limits_tb_MI);
410  f->GetObject("limits_MR_MI;1",limits_MR_MI);
411  f->GetObject("limits_MR_McH;1",limits_MR_McH);
412  f->GetObject("limits_MI_McH;1",limits_MI_McH);
413 
414  TVectorD* vllmax=NULL;
415 
416  f->GetObject("vllmax;1",vllmax);
417  if(!vllmax) cout<<"ERROR"<<endl;
418  llmax=(*vllmax)[0];
419  //tbmax=(*vllmax)[1];
420  //McHmax=(*vllmax)[2];
421  //MRmax=McHmax+(*vllmax)[3];
422  //MImax=MRmax+(*vllmax)[4];
423  cout<<llmax<<" "<<tbmax<<" "<<McHmax<<" "<<MRmax<<" "<<MImax<<" "<<endl;
424 
425  /*BGL2* m=new BGL2(gL,gQ,lup,qup);
426  double sm_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,5));
427  double charged_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,1));
428  double neutral_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,2));
429  double neutralR_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,3));
430  double neutralI_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,4));
431  double eK_(m->epsK(tbmax,McHmax,MRmax,MImax));
432 
433  double all_(m->bsgammawidth(tbmax,McHmax,MRmax,MImax,0));
434  */
435  //for(int gL=2;gL>=0;gL--)
436  //for(int gQ=2;gQ>=0;gQ--)
437  //for(uint lup=0;lup<2;lup++)
438  //for(uint qup=0;qup<2;qup++)
439  uint min1=npoints, min2=npoints, min3=npoints;
440  uint min11=npoints, min21=npoints, min31=npoints;
441  uint min12=npoints, min22=npoints, min32=npoints;
442 
443  for(uint i=0;i<npoints;i++)
444  for(uint j=0;j<npoints;j++){
445  int binmax=limits4->GetBin(i+1,j+1);
446  double rest=limits4->GetBinContent(binmax);
447  if(rest>=llmax) rest=1;
448  else rest=TMath::Prob(-2*(rest-llmax),2);
449  if(rest>=0.05 && j<min1){min1=j;}
450  if(rest>=0.05 && j<min11 && j>uint(180*npoints/990)){min11=j;}
451  if(rest>=0.05 && j<min12 && j>uint(400*npoints/990)){min12=j;}
452  limits4->SetBinContent(i+1,j+1,rest);
453 
454  rest=Bmumu_Bsmumu->GetBinContent(binmax);
455  if(rest>=llmax) rest=1;
456  else rest=TMath::Prob(-2*(rest-llmax),2);
457  //int nn=4;
458  //int ii=(i/nn)*nn, jj=(j/nn)*nn;
459  //for(int iii=ii;iii<ii+n;++iii)
460  //for(int iii=ii;iii<ii+n;++iii)
461  Bmumu_Bsmumu->SetBinContent(i+1,j+1,rest);
462 
463  rest=limits_MR_MI->GetBinContent(binmax);
464  if(rest>=llmax) rest=1;
465  else rest=TMath::Prob(-2*(rest-llmax),2);
466  limits_MR_MI->SetBinContent(i+1,j+1,rest);
467 
468  rest=limits_MR_McH->GetBinContent(binmax);
469  if(rest>=llmax) rest=1;
470  else rest=TMath::Prob(-2*(rest-llmax),2);
471  if(rest>=0.05 && i<min2){min2=i;}
472  if(rest>=0.05 && i<min21 && j>uint(180*npoints/990)){min21=i;}
473  if(rest>=0.05 && i<min22 && j>uint(400*npoints/990)){min22=i;}
474  limits_MR_McH->SetBinContent(i+1,j+1,rest);
475 
476  rest=limits_MI_McH->GetBinContent(binmax);
477  if(rest>=llmax) rest=1;
478  else rest=TMath::Prob(-2*(rest-llmax),2);
479  if(rest>=0.05 && i<min3){min3=i;}
480  if(rest>=0.05 && i<min31 && j>uint(180*npoints/990)){min31=i;}
481  if(rest>=0.05 && i<min32 && j>uint(400*npoints/990)){min32=i;}
482  limits_MI_McH->SetBinContent(i+1,j+1,rest);
483 
484 
485  rest=limits_tb_MR->GetBinContent(binmax);
486  if(rest>=llmax) rest=1;
487  else rest=TMath::Prob(-2*(rest-llmax),2);
488  limits_tb_MR->SetBinContent(i+1,j+1,rest);
489  }
490  double mmin1=init2+((min1+0.5)*(final2-init2))/npoints;
491  double mmin2=init2+((min2+0.5)*(final2-init2))/npoints;
492  double mmin3=init2+((min3+0.5)*(final2-init2))/npoints;
493 
494  double mmin11=init2+((min11+0.5)*(final2-init2))/npoints;
495  double mmin21=init2+((min21+0.5)*(final2-init2))/npoints;
496  double mmin31=init2+((min31+0.5)*(final2-init2))/npoints;
497 
498  double mmin12=init2+((min12+0.5)*(final2-init2))/npoints;
499  double mmin22=init2+((min22+0.5)*(final2-init2))/npoints;
500  double mmin32=init2+((min32+0.5)*(final2-init2))/npoints;
501 
502  //mass<<name<<" "<<mmin1<<" "<<mmin2<<" "<<mmin3<<endl;
503 
504  ofstream maxs((string("maxs_")+string(name)+string(".out")).c_str());
505  maxs<<mmin1<<" "<<mmin2<<" "<<mmin3<<endl;
506  maxs<<mmin11<<" "<<mmin21<<" "<<mmin31<<endl;
507  maxs<<mmin12<<" "<<mmin22<<" "<<mmin32<<endl;
508  maxs<<llmax<<" "<<tbmax<<" "<<McHmax<<" "<<MRmax<<" "<<MImax<<" "<<endl;
509  //maxs<<eK_<<endl;
510  //maxs<<sm_<<" "<<charged_<<" "<<neutral_<<" "<<neutralR_<<" "<<neutralI_<<" "<<all_<<endl;
511 
512  //for(uint j=0;j<npoints;j++)
513  //for(uint i=0;i<npoints;i++){
514  // int binmax=limits4->GetBin(i+1,j+1);
515  // maxs<<"("<<i<<","<<j<<"):"<<limits4->GetBinContent(binmax)<<endl;
516  // }
517 
518  maxs.close();
519 
520  double ma=0,me=.2, x0=1,y0=120;
521  gStyle->SetOptTitle(0);
522  gStyle->SetPaperSize(10.,10.);
523  TCanvas * c21=new TCanvas("c21","",int(800*(1+ma+me)),int(600*(1+ma+me)));
524  c21->SetMargin(me,ma,me,ma);
525  c21->SetGrid();
526 
527  limits4->SetStats(0);
528  limits4->GetXaxis()->SetTitle("log_{10}(tan\\beta)");
529  limits4->GetYaxis()->SetTitle("M_{H+} (GeV)");
530 
531  Bmumu_Bsmumu->SetStats(0);
532  Bmumu_Bsmumu->GetXaxis()->SetTitle("Br(B\\to\\mu\\mu)/10^{-10}");
533  Bmumu_Bsmumu->GetYaxis()->SetTitle("Br(B_{s}\\to\\mu\\mu)/10^{-9}");
534 
535  limits_MR_MI->SetStats(0);
536  limits_MR_MI->GetYaxis()->SetTitle("M_{I} (GeV)");
537  limits_MR_MI->GetXaxis()->SetTitle("M_{R} (GeV)");
538 
539 
540  limits_MR_McH->SetStats(0);
541  limits_MR_McH->GetYaxis()->SetTitle("M_{H+} (GeV)");
542  limits_MR_McH->GetXaxis()->SetTitle("M_{R} (GeV)");
543 
544  limits_MI_McH->SetStats(0);
545  limits_MI_McH->GetYaxis()->SetTitle("M_{H+} (GeV)");
546  limits_MI_McH->GetXaxis()->SetTitle("M_{I} (GeV)");
547 
548  limits_tb_MR->SetStats(0);
549  limits_tb_MR->GetXaxis()->SetTitle("log_{10}(tan\\beta)");
550  limits_tb_MR->GetYaxis()->SetTitle("M_{R} (GeV)");
551 
552 
553 
554  Double_t contours[3];
555  contours[0] = 0.003;
556  contours[1] = 0.05;
557  contours[2] = 0.32;
558 
559 
560 
561 
562  limits4->SetContour(3, contours);
563  //limits4->GetYaxis()->SetLabelOffset(0.02);
564  limits4->GetYaxis()->SetLabelSize(0.08);
565  limits4->GetYaxis()->SetTitleSize(0.08);
566  limits4->GetYaxis()->SetTitleOffset(1.2);
567  limits4->GetYaxis()->SetLimits(1,999);
568 
569 
570 
571  //limits4->GetXaxis()->SetLabelOffset(0.02);
572  limits4->GetXaxis()->SetLabelSize(0.08);
573  limits4->GetXaxis()->SetTitleSize(0.08);
574  limits4->GetXaxis()->SetTitleOffset(1.2);
575  limits4->GetXaxis()->SetLimits(-2.99,2.99);
576 
577 
578 
579  TLatex l;
580  l.SetTextSize(0.08);
581  string ss=qq[qup][gQ]+","+ll[lup][gL];
582 
583 
584  limits4->Draw("CONT Z LIST");
585  //limits4->Draw("CONT LIST");
586  //limits4->Draw("colz");
587 
588  l.DrawLatex(x0,y0,ss.c_str());
589 
590  c21->SaveAs((string("pdf_")+string(name)+string(".png")).c_str());
591 
592  delete c21;
593  //Bmumu_Bsmumu->SetBit(TH1::kCanRebin);
594  Bmumu_Bsmumu->Rebin2D(2,2);
595 
596  //Bmumu_Bsmumu->SetContour(3, contours);
597  //limits4->GetYaxis()->SetLabelOffset(0.02);
598  Bmumu_Bsmumu->GetYaxis()->SetLabelSize(0.08);
599  Bmumu_Bsmumu->GetYaxis()->SetTitleSize(0.08);
600  Bmumu_Bsmumu->GetYaxis()->SetTitleOffset(0.8);
601  Bmumu_Bsmumu->GetYaxis()->SetLimits(0.01,4.99);
602  //Bmumu_Bsmumu->GetYaxis()->SetRangeUser(0.01, 3.49);
603  // Bmumu_Bsmumu->GetYaxis()->SetLimits(0.01,3.49);
604  //limits4->GetXaxis()->SetLabelOffset(0.02);
605  Bmumu_Bsmumu->GetYaxis()->SetNdivisions(5, kTRUE);
606  Bmumu_Bsmumu->GetXaxis()->SetNdivisions(5, kTRUE);
607 
608  Bmumu_Bsmumu->GetXaxis()->SetLabelSize(0.08);
609  Bmumu_Bsmumu->GetXaxis()->SetTitleSize(0.08);
610  Bmumu_Bsmumu->GetXaxis()->SetTitleOffset(1.2);
611  Bmumu_Bsmumu->GetXaxis()->SetLimits(0.01,1.99);
612  //Bmumu_Bsmumu->GetXaxis()->SetRangeUser(0., 2);
613 
614  TCanvas * cB=new TCanvas("cB","",int(800*(1+ma+me)),int(600*(1+ma+me)));
615  cB->SetMargin(.14,ma,me,ma);
616  cB->SetGrid();
617  //limits4->Draw("CONT Z LIST");
618  Bmumu_Bsmumu->Draw("COLZ");
619 
620  l.DrawLatex(1.5,1,ss.c_str());
621 
622  cB->SaveAs((string("Bmumu_Bsmumu_")+string(name)+string(".png")).c_str());
623 
624  delete cB;
625 
626  limits_MR_MI->SetContour(3, contours);
627  //limits4->GetYaxis()->SetLabelOffset(0.02);
628  limits_MR_MI->GetYaxis()->SetLabelSize(0.06);
629  limits_MR_MI->GetYaxis()->SetTitleSize(0.06);
630  limits_MR_MI->GetYaxis()->SetTitleOffset(1.1);
631  limits_MR_MI->GetYaxis()->SetLimits(1,999);
632  //limits4->GetXaxis()->SetLabelOffset(0.02);
633  limits_MR_MI->GetXaxis()->SetLabelSize(0.06);
634  limits_MR_MI->GetXaxis()->SetTitleSize(0.06);
635  limits_MR_MI->GetXaxis()->SetTitleOffset(1.1);
636  limits_MR_MI->GetXaxis()->SetLimits(1,999);
637 
638  TCanvas * c3=new TCanvas("c3","",800,600);
639  c3->SetMargin(me,ma,me,ma);
640  c3->SetGrid();
641 
642  limits_MR_MI->Draw("CONT LIST");
643 
644  c3->SaveAs((string("pdf_")+string(name)+string("_MRMI.png")).c_str());
645 
646  limits_MR_McH->SetContour(3, contours);
647  //limits4->GetYaxis()->SetLabelOffset(0.02);
648  limits_MR_McH->GetYaxis()->SetLabelSize(0.06);
649  limits_MR_McH->GetYaxis()->SetTitleSize(0.06);
650  limits_MR_McH->GetYaxis()->SetTitleOffset(1.1);
651  limits_MR_McH->GetYaxis()->SetLimits(1,999);
652  //limits4->GetXaxis()->SetLabelOffset(0.02);
653  limits_MR_McH->GetXaxis()->SetLabelSize(0.06);
654  limits_MR_McH->GetXaxis()->SetTitleSize(0.06);
655  limits_MR_McH->GetXaxis()->SetTitleOffset(1.1);
656  limits_MR_McH->GetXaxis()->SetLimits(1,999);
657 
658  TCanvas * c4=new TCanvas("c4","",800,600);
659  limits_MR_McH->Draw("CONT LIST");
660  c4->SetMargin(me,ma,me,ma);
661  c4->SetGrid();
662  c4->SaveAs((string("pdf_")+string(name)+string("_MRMcH.png")).c_str());
663 
664  limits_MI_McH->SetContour(3, contours);
665  //limits4->GetYaxis()->SetLabelOffset(0.02);
666  limits_MI_McH->GetYaxis()->SetLabelSize(0.06);
667  limits_MI_McH->GetYaxis()->SetTitleSize(0.06);
668  limits_MI_McH->GetYaxis()->SetTitleOffset(1.1);
669  limits_MI_McH->GetYaxis()->SetLimits(1,999);
670  //limits4->GetXaxis()->SetLabelOffset(0.02);
671  limits_MI_McH->GetXaxis()->SetLabelSize(0.06);
672  limits_MI_McH->GetXaxis()->SetTitleSize(0.06);
673  limits_MI_McH->GetXaxis()->SetTitleOffset(1.1);
674  limits_MI_McH->GetXaxis()->SetLimits(1,999);
675 
676  TCanvas * c6=new TCanvas("c6","",800,600);
677  limits_MI_McH->Draw("CONT LIST");
678  c6->SetMargin(me,ma,me,ma);
679  c6->SetGrid();
680  c6->SaveAs((string("pdf_")+string(name)+string("_MIMcH.png")).c_str());
681 
682  TCanvas * c5=new TCanvas("c5","",800,600);
683  limits_tb_MR->Draw("colz");
684 
685  c5->SaveAs((string("pdf_")+string(name)+string("_tbMR.png")).c_str());
686 
687  //delete m;
688  //mass.close();
689  f->Close();
690  delete f;
691  return 0;
692 
693 }
uint
unsigned int uint
Definition: script.cpp:4
BGLmodels::CD
std::complex< double > CD
Definition: Formulas.h:65
BGLmodels::Vud
const Matrixx Vud(13.04 *M_PI/180, 0.201 *M_PI/180, 2.38 *M_PI/180, 1.2)

Here is the call graph for this function: