SusyScan/PlotScript/TheLimits.cc

#include "TheLimits.h"
#include "SusyScan.h"
#include "GeneratorMasses.h"

#include <fstream>
#include <iostream>
#include <cmath>

void TheLimits::Fill(int argc, char** argv)
{
   for (int i = 1; i<argc; ++i)
   {
     add( new SusyScan(argv[i]) );
   }
   std::cout << "Read " << _scan.size() << " points from argument list." <<std::endl;
}

void TheLimits::Fill(const std::string filelist)
{
   std::ifstream masses_file;
   masses_file.open(filelist.c_str());
   std::string file;
   while (1) {
      GeneratorMasses * p = new GeneratorMasses;
      masses_file >> file;
      if (!masses_file.good()) break;
      add( new SusyScan(file));
   }
   std::cout << "Read " << _scan.size() << " points from file " << filelist <<std::endl;
   masses_file.close();
}

void TheLimits::match()
{
  for (std::vector<SusyScan*>::iterator it=_scan.begin();it!=_scan.end();++it){
    bool match = false;
    for (std::vector<GeneratorMasses*>::const_iterator gt=_masses.begin();gt!=_masses.end();++gt){
      if ((*it)->Mzero==(*gt)->Mzero && 
          (*it)->Mhalf==(*gt)->Mhalf && 
          (*it)->Azero==(*gt)->Azero &&
          (*it)->TanBeta==(*gt)->TanBeta &&
          (*it)->Mu==(*gt)->Mu ) {
        (*it)->M1 = (*gt)->M1;
        (*it)->M2 = (*gt)->M2;
        (*it)->M3 = (*gt)->M3;
        (*it)->MGL = (*gt)->MGL;
        (*it)->MUL = (*gt)->MUL;
        (*it)->MB1 = (*gt)->MB1;
        (*it)->MSN = (*gt)->MSN;
        (*it)->MNTAU = (*gt)->MNTAU;
        (*it)->MZ1 = (*gt)->MZ1;
        (*it)->MW1 = (*gt)->MW1;
        (*it)->MHL = (*gt)->MHL;
        (*it)->MUR = (*gt)->MUR;
        (*it)->MB2 = (*gt)->MB2;
        (*it)->MEL = (*gt)->MEL;
        (*it)->MTAU1 = (*gt)->MTAU1;
        (*it)->MZ2 = (*gt)->MZ2;
        (*it)->MW2 = (*gt)->MW2;
        (*it)->MHH = (*gt)->MHH;
        (*it)->MDL = (*gt)->MDL;
        (*it)->MT1 = (*gt)->MT1;
        (*it)->MER = (*gt)->MER;
        (*it)->MTAU2 = (*gt)->MTAU2;
        (*it)->MZ3 = (*gt)->MZ3;
        (*it)->MHA = (*gt)->MHA;
        (*it)->MDR = (*gt)->MDR;
        (*it)->MT2 = (*gt)->MT2;
        (*it)->MZ4 = (*gt)->MZ4;
        (*it)->MHp = (*gt)->MHp;
        match = true;
        //std::cout <<"m0="<<(*it)->Mzero<<", m1/2="<<(*it)->Mhalf<< ", sq="<<(*it)->MUL << ", gl="<<(*it)->MGL <<std::endl; 
      }   
    }
    //if (!match) std::cout << "No match for M0="<<(*it)->Mzero
    //                      << ", M12="<<(*it)->Mhalf<<std::endl;
  }
}

void TheLimits::FillGeneratorMasses(std::string file)
{
   std::ifstream masses_file;
   masses_file.open(file.c_str());
   while (1) {
      GeneratorMasses * p = new GeneratorMasses;
      masses_file >> p->Mzero
                  >> p->Mhalf
                  >> p->TanBeta
                  >> p->Mu  
                  >> p->Azero
                  >> p->Mtop  
                  >> p->muQ       
                  >> p->Q         
                  >> p->M1        
                  >> p->M2
                  >> p->M3        
                  >> p->MGL       
                  >> p->MUL       
                  >> p->MB1       
                  >> p->MSN       
                  >> p->MNTAU     
                  >> p->MZ1       
                  >> p->MW1       
                  >> p->MHL       
                  >> p->MUR       
                  >> p->MB2       
                  >> p->MEL       
                  >> p->MTAU1     
                  >> p->MZ2       
                  >> p->MW2       
                  >> p->MHH       
                  >> p->MDL       
                  >> p->MT1       
                  >> p->MER       
                  >> p->MTAU2     
                  >> p->MZ3       
                  >> p->MHA       
                  >> p->MDR       
                  >> p->MT2       
                  >> p->MZ4       
                  >> p->MHp;

      if (!masses_file.good()) break;
      //std::cout << p->Mzero<<", "<<p->Mhalf<<", "<<p->TanBeta<<std::endl;
      if (fabs(p->Mu)!=1.) {
         std::cerr << "check lines near m0=" << p->Mzero << ", m1/2=" << p->Mhalf << std::endl;
         break;
      } 
      //std::cout <<"m0="<< p->Mzero<<", m1/2="<< p->Mhalf<< ", sq="<<p->MUL << ", gl="<<p->MGL <<std::endl; 
      _masses.push_back( p );
   }
   std::cout << "Read " << _masses.size() << " points from file " << file <<std::endl;

}


void TheLimits::OverwriteLimits(std::string flag)
{
  std::cout<<"WARNING: OVERWRITING LIMITS!!!!"<<std::endl;
  for (std::vector<SusyScan*>::iterator it=_scan.begin();it!=_scan.end();++it){
    if (flag=="ABCD_MHT") {
      (*it)->ExpNsigLimit = 17.85;
      (*it)->PLExpNsigLimit = 12.558;
      (*it)->FCExpNsigLimit = 15.75;
      (*it)->MCMCExpNsigLimit = 15.8085;
      (*it)->ObsNsigLimit = 10.0637;
      (*it)->PLObsNsigLimit = 4.22124;
      (*it)->FCObsNsigLimit  = 0.75;
      (*it)->MCMCObsNsigLimit  = 8.75403;
    }
    else if (flag=="ABCD_HT") {
      (*it)->ExpNsigLimit = 24.0438;
      (*it)->PLExpNsigLimit = 17.8537;
      (*it)->FCExpNsigLimit = 22.05;
      (*it)->MCMCExpNsigLimit = 20.781;
      (*it)->ObsNsigLimit    = 18.24;
      (*it)->PLObsNsigLimit = 9.70375;
      (*it)->FCObsNsigLimit  = 11.85;
      (*it)->MCMCObsNsigLimit  = 14.34;
    }
    
    
      (*it)->ExpXsecLimit    = (*it)->ExpNsigLimit * (*it)->Xsection / (*it)->signal;
      (*it)->PLExpXsecLimit  = (*it)->PLExpNsigLimit  * (*it)->Xsection / (*it)->signal;
      (*it)->FCExpXsecLimit  = (*it)->FCExpNsigLimit  * (*it)->Xsection / (*it)->signal;
      (*it)->MCMCExpXsecLimit= (*it)->MCMCExpNsigLimit  * (*it)->Xsection / (*it)->signal;
      (*it)->ObsXsecLimit    = (*it)->ObsNsigLimit  * (*it)->Xsection / (*it)->signal;
      (*it)->PLObsXsecLimit  = (*it)->PLObsNsigLimit  * (*it)->Xsection / (*it)->signal;
      (*it)->FCObsXsecLimit  = (*it)->FCObsNsigLimit   * (*it)->Xsection / (*it)->signal;
      (*it)->MCMCObsXsecLimit= (*it)->MCMCObsNsigLimit  * (*it)->Xsection / (*it)->signal;
    
  
  } 


}

void TheLimits::ExpandGrid(int s)
{
  std::vector<SusyScan*> new_grid;
  for (std::vector<SusyScan*>::iterator it=_scan.begin(); it!=_scan.end(); ++it){
    double x=9999, y=9999;
    std::vector<SusyScan*>::iterator next_x=_scan.end(), next_y=_scan.end();
    for (std::vector<SusyScan*>::iterator nx=_scan.begin(); nx!=_scan.end(); ++nx)
      if ((*it)->Mzero<(*nx)->Mzero && fabs((*it)->Mzero-(*nx)->Mzero)<x) {
        x=fabs((*it)->Mzero-(*nx)->Mzero);
        next_x=nx;
      }
    for (std::vector<SusyScan*>::iterator ny=_scan.begin(); ny!=_scan.end(); ++ny)
      if ( (*it)->Mzero==(*ny)->Mzero && (*it)->Mhalf<(*ny)->Mhalf && fabs((*it)->Mhalf-(*ny)->Mhalf)<y) {
        y=fabs((*it)->Mhalf-(*ny)->Mhalf);
        next_y=ny;
      }
    if (next_y!=_scan.end()) {  
        SusyScan*  n = new SusyScan((**it + **next_y) * 0.5);
        new_grid.push_back( n );
        //std::cout 
        //<< "l:"<<(*it)->Mzero<<","<<(*it)->Mhalf
        //<< " <>"<<    n->Mzero<<","<<n->Mhalf
        //<< " <> r:"<<(*next_y)->Mzero<<","<<(*next_y)->Mhalf
        //<<std::endl;
    }  
      
  } 
  _scan.insert(_scan.end(), new_grid.begin(), new_grid.end());
}
Revision:	1.10
Committed:	Wed Jun 22 15:03:36 2011 UTC (13 years, 10 months ago) by auterman
Content type:	text/plain
Branch:	MAIN
CVS Tags:	JHEP2010, HEAD
Changes since 1.9:	+7 -1 lines
Log Message:	2010 RA2 paper
#	Content
1	#include "TheLimits.h"
2	#include "SusyScan.h"
3	#include "GeneratorMasses.h"
4
5	#include <fstream>
6	#include <iostream>
7	#include <cmath>
8
9	void TheLimits::Fill(int argc, char** argv)
10	{
11	for (int i = 1; i<argc; ++i)
12	{
13	add( new SusyScan(argv[i]) );
14	}
15	std::cout << "Read " << _scan.size() << " points from argument list." <<std::endl;
16	}
17
18	void TheLimits::Fill(const std::string filelist)
19	{
20	std::ifstream masses_file;
21	masses_file.open(filelist.c_str());
22	std::string file;
23	while (1) {
24	GeneratorMasses * p = new GeneratorMasses;
25	masses_file >> file;
26	if (!masses_file.good()) break;
27	add( new SusyScan(file));
28	}
29	std::cout << "Read " << _scan.size() << " points from file " << filelist <<std::endl;
30	masses_file.close();
31	}
32
33	void TheLimits::match()
34	{
35	for (std::vector<SusyScan*>::iterator it=_scan.begin();it!=_scan.end();++it){
36	bool match = false;
37	for (std::vector<GeneratorMasses*>::const_iterator gt=_masses.begin();gt!=_masses.end();++gt){
38	if ((it)->Mzero==(gt)->Mzero &&
39	(it)->Mhalf==(gt)->Mhalf &&
40	(it)->Azero==(gt)->Azero &&
41	(it)->TanBeta==(gt)->TanBeta &&
42	(it)->Mu==(gt)->Mu ) {
43	(it)->M1 = (gt)->M1;
44	(it)->M2 = (gt)->M2;
45	(it)->M3 = (gt)->M3;
46	(it)->MGL = (gt)->MGL;
47	(it)->MUL = (gt)->MUL;
48	(it)->MB1 = (gt)->MB1;
49	(it)->MSN = (gt)->MSN;
50	(it)->MNTAU = (gt)->MNTAU;
51	(it)->MZ1 = (gt)->MZ1;
52	(it)->MW1 = (gt)->MW1;
53	(it)->MHL = (gt)->MHL;
54	(it)->MUR = (gt)->MUR;
55	(it)->MB2 = (gt)->MB2;
56	(it)->MEL = (gt)->MEL;
57	(it)->MTAU1 = (gt)->MTAU1;
58	(it)->MZ2 = (gt)->MZ2;
59	(it)->MW2 = (gt)->MW2;
60	(it)->MHH = (gt)->MHH;
61	(it)->MDL = (gt)->MDL;
62	(it)->MT1 = (gt)->MT1;
63	(it)->MER = (gt)->MER;
64	(it)->MTAU2 = (gt)->MTAU2;
65	(it)->MZ3 = (gt)->MZ3;
66	(it)->MHA = (gt)->MHA;
67	(it)->MDR = (gt)->MDR;
68	(it)->MT2 = (gt)->MT2;
69	(it)->MZ4 = (gt)->MZ4;
70	(it)->MHp = (gt)->MHp;
71	match = true;
72	//std::cout <<"m0="<<(it)->Mzero<<", m1/2="<<(it)->Mhalf<< ", sq="<<(it)->MUL << ", gl="<<(it)->MGL <<std::endl;
73	}
74	}
75	//if (!match) std::cout << "No match for M0="<<(*it)->Mzero
76	// << ", M12="<<(*it)->Mhalf<<std::endl;
77	}
78	}
79
80	void TheLimits::FillGeneratorMasses(std::string file)
81	{
82	std::ifstream masses_file;
83	masses_file.open(file.c_str());
84	while (1) {
85	GeneratorMasses * p = new GeneratorMasses;
86	masses_file >> p->Mzero
87	>> p->Mhalf
88	>> p->TanBeta
89	>> p->Mu
90	>> p->Azero
91	>> p->Mtop
92	>> p->muQ
93	>> p->Q
94	>> p->M1
95	>> p->M2
96	>> p->M3
97	>> p->MGL
98	>> p->MUL
99	>> p->MB1
100	>> p->MSN
101	>> p->MNTAU
102	>> p->MZ1
103	>> p->MW1
104	>> p->MHL
105	>> p->MUR
106	>> p->MB2
107	>> p->MEL
108	>> p->MTAU1
109	>> p->MZ2
110	>> p->MW2
111	>> p->MHH
112	>> p->MDL
113	>> p->MT1
114	>> p->MER
115	>> p->MTAU2
116	>> p->MZ3
117	>> p->MHA
118	>> p->MDR
119	>> p->MT2
120	>> p->MZ4
121	>> p->MHp;
122
123	if (!masses_file.good()) break;
124	//std::cout << p->Mzero<<", "<<p->Mhalf<<", "<<p->TanBeta<<std::endl;
125	if (fabs(p->Mu)!=1.) {
126	std::cerr << "check lines near m0=" << p->Mzero << ", m1/2=" << p->Mhalf << std::endl;
127	break;
128	}
129	//std::cout <<"m0="<< p->Mzero<<", m1/2="<< p->Mhalf<< ", sq="<<p->MUL << ", gl="<<p->MGL <<std::endl;
130	_masses.push_back( p );
131	}
132	std::cout << "Read " << _masses.size() << " points from file " << file <<std::endl;
133
134	}
135
136
137	void TheLimits::OverwriteLimits(std::string flag)
138	{
139	std::cout<<"WARNING: OVERWRITING LIMITS!!!!"<<std::endl;
140	for (std::vector<SusyScan*>::iterator it=_scan.begin();it!=_scan.end();++it){
141	if (flag=="ABCD_MHT") {
142	(*it)->ExpNsigLimit = 17.85;
143	(*it)->PLExpNsigLimit = 12.558;
144	(*it)->FCExpNsigLimit = 15.75;
145	(*it)->MCMCExpNsigLimit = 15.8085;
146	(*it)->ObsNsigLimit = 10.0637;
147	(*it)->PLObsNsigLimit = 4.22124;
148	(*it)->FCObsNsigLimit = 0.75;
149	(*it)->MCMCObsNsigLimit = 8.75403;
150	}
151	else if (flag=="ABCD_HT") {
152	(*it)->ExpNsigLimit = 24.0438;
153	(*it)->PLExpNsigLimit = 17.8537;
154	(*it)->FCExpNsigLimit = 22.05;
155	(*it)->MCMCExpNsigLimit = 20.781;
156	(*it)->ObsNsigLimit = 18.24;
157	(*it)->PLObsNsigLimit = 9.70375;
158	(*it)->FCObsNsigLimit = 11.85;
159	(*it)->MCMCObsNsigLimit = 14.34;
160	}
161
162
163	(it)->ExpXsecLimit = (it)->ExpNsigLimit * (it)->Xsection / (it)->signal;
164	(it)->PLExpXsecLimit = (it)->PLExpNsigLimit * (it)->Xsection / (it)->signal;
165	(it)->FCExpXsecLimit = (it)->FCExpNsigLimit * (it)->Xsection / (it)->signal;
166	(it)->MCMCExpXsecLimit= (it)->MCMCExpNsigLimit * (it)->Xsection / (it)->signal;
167	(it)->ObsXsecLimit = (it)->ObsNsigLimit * (it)->Xsection / (it)->signal;
168	(it)->PLObsXsecLimit = (it)->PLObsNsigLimit * (it)->Xsection / (it)->signal;
169	(it)->FCObsXsecLimit = (it)->FCObsNsigLimit * (it)->Xsection / (it)->signal;
170	(it)->MCMCObsXsecLimit= (it)->MCMCObsNsigLimit * (it)->Xsection / (it)->signal;
171
172
173	}
174
175
176	}
177
178	void TheLimits::ExpandGrid(int s)
179	{
180	std::vector<SusyScan*> new_grid;
181	for (std::vector<SusyScan*>::iterator it=_scan.begin(); it!=_scan.end(); ++it){
182	double x=9999, y=9999;
183	std::vector<SusyScan*>::iterator next_x=_scan.end(), next_y=_scan.end();
184	for (std::vector<SusyScan*>::iterator nx=_scan.begin(); nx!=_scan.end(); ++nx)
185	if ((it)->Mzero<(nx)->Mzero && fabs((it)->Mzero-(nx)->Mzero)<x) {
186	x=fabs((it)->Mzero-(nx)->Mzero);
187	next_x=nx;
188	}
189	for (std::vector<SusyScan*>::iterator ny=_scan.begin(); ny!=_scan.end(); ++ny)
190	if ( (it)->Mzero==(ny)->Mzero && (it)->Mhalf<(ny)->Mhalf && fabs((it)->Mhalf-(ny)->Mhalf)<y) {
191	y=fabs((it)->Mhalf-(ny)->Mhalf);
192	next_y=ny;
193	}
194	if (next_y!=_scan.end()) {
195	SusyScan* n = new SusyScan((it + next_y) * 0.5);
196	new_grid.push_back( n );
197	//std::cout
198	//<< "l:"<<(it)->Mzero<<","<<(it)->Mhalf
199	//<< " <>"<< n->Mzero<<","<<n->Mhalf
200	//<< " <> r:"<<(next_y)->Mzero<<","<<(next_y)->Mhalf
201	//<<std::endl;
202	}
203
204	}
205	_scan.insert(_scan.end(), new_grid.begin(), new_grid.end());
206	}