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
#	User	Rev	Content
1	auterman	1.1	#include "TheLimits.h"
2			#include "SusyScan.h"
3			#include "GeneratorMasses.h"
4
5	auterman	1.4	#include <fstream>
6			#include <iostream>
7			#include <cmath>
8
9	auterman	1.6	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	auterman	1.10	std::cout << "Read " << _scan.size() << " points from argument list." <<std::endl;
16	auterman	1.6	}
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	auterman	1.10	std::cout << "Read " << _scan.size() << " points from file " << filelist <<std::endl;
30	auterman	1.6	masses_file.close();
31			}
32	auterman	1.4
33			void TheLimits::match()
34	auterman	1.1	{
35			for (std::vector<SusyScan*>::iterator it=_scan.begin();it!=_scan.end();++it){
36			bool match = false;
37	auterman	1.4	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	auterman	1.1	match = true;
72	auterman	1.10	//std::cout <<"m0="<<(it)->Mzero<<", m1/2="<<(it)->Mhalf<< ", sq="<<(it)->MUL << ", gl="<<(it)->MGL <<std::endl;
73	auterman	1.1	}
74			}
75			//if (!match) std::cout << "No match for M0="<<(*it)->Mzero
76			// << ", M12="<<(*it)->Mhalf<<std::endl;
77			}
78			}
79	auterman	1.4
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	auterman	1.8	//std::cout << p->Mzero<<", "<<p->Mhalf<<", "<<p->TanBeta<<std::endl;
125	auterman	1.4	if (fabs(p->Mu)!=1.) {
126			std::cerr << "check lines near m0=" << p->Mzero << ", m1/2=" << p->Mhalf << std::endl;
127			break;
128			}
129	auterman	1.10	//std::cout <<"m0="<< p->Mzero<<", m1/2="<< p->Mhalf<< ", sq="<<p->MUL << ", gl="<<p->MGL <<std::endl;
130	auterman	1.4	_masses.push_back( p );
131			}
132	auterman	1.10	std::cout << "Read " << _masses.size() << " points from file " << file <<std::endl;
133
134	auterman	1.4	}
135	auterman	1.5
136
137			void TheLimits::OverwriteLimits(std::string flag)
138			{
139	auterman	1.10	std::cout<<"WARNING: OVERWRITING LIMITS!!!!"<<std::endl;
140	auterman	1.5	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	auterman	1.9	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			}