| 1 |
csander |
1.1 |
////////////////////////////////////////////////////////////////
|
| 2 |
|
|
// IsaFialStateAnalyzer v1.0
|
| 3 |
|
|
// date 15.06.007
|
| 4 |
|
|
// original author: Christian Sander
|
| 5 |
|
|
//
|
| 6 |
|
|
// This program reads in the ascii decay table which is output
|
| 7 |
|
|
// from the RGE code ISASUGRA. Then the branching ratio of any
|
| 8 |
|
|
// particle into a wanted final state defined by number of
|
| 9 |
|
|
// leptons, quarks, taus, gluons can be calculated.
|
| 10 |
|
|
//
|
| 11 |
|
|
// To do:
|
| 12 |
|
|
// - direct interface to ISASUGRA (done by shell script)
|
| 13 |
|
|
// - allow decay of taus and gluons
|
| 14 |
|
|
// - allow definition final state with number of Ws, Zs and or
|
| 15 |
|
|
// Higgses (almost done)
|
| 16 |
|
|
// - develop routines to scan over parameter space and fill
|
| 17 |
|
|
// histograms or ntuple (will be done externally)
|
| 18 |
|
|
// - include cross section to create primary particle to get
|
| 19 |
|
|
// estimates of total number of events (almost done)
|
| 20 |
|
|
// - split into several files
|
| 21 |
|
|
////////////////////////////////////////////////////////////////
|
| 22 |
|
|
|
| 23 |
|
|
#include<string>
|
| 24 |
|
|
#include<iostream>
|
| 25 |
|
|
#include<fstream>
|
| 26 |
|
|
#include<vector>
|
| 27 |
|
|
#include<map>
|
| 28 |
|
|
#include<sstream>
|
| 29 |
|
|
#include<stdlib.h>
|
| 30 |
|
|
|
| 31 |
|
|
using namespace std;
|
| 32 |
|
|
|
| 33 |
|
|
////////////////////////////////////////////////////////////////
|
| 34 |
|
|
// This class describes a particle
|
| 35 |
|
|
////////////////////////////////////////////////////////////////
|
| 36 |
|
|
class particleClass
|
| 37 |
|
|
{
|
| 38 |
|
|
private:
|
| 39 |
|
|
string Pid;
|
| 40 |
|
|
map<string, string> APid;
|
| 41 |
|
|
|
| 42 |
|
|
public:
|
| 43 |
|
|
string getPid(){
|
| 44 |
|
|
return Pid;
|
| 45 |
|
|
};
|
| 46 |
|
|
string getAPid(){
|
| 47 |
|
|
return APid[Pid];
|
| 48 |
|
|
};
|
| 49 |
|
|
void setPid(string input){
|
| 50 |
|
|
Pid.clear();
|
| 51 |
|
|
Pid.append(input);
|
| 52 |
|
|
};
|
| 53 |
|
|
bool isQuark(){
|
| 54 |
|
|
string quarks = "UP UB DN DB ST SB CH CB BT BB";
|
| 55 |
|
|
string::size_type loc = quarks.find( Pid, 0 );
|
| 56 |
|
|
if( loc != string::npos ) {
|
| 57 |
|
|
return true;
|
| 58 |
|
|
} else {
|
| 59 |
|
|
return false;
|
| 60 |
|
|
}
|
| 61 |
|
|
};
|
| 62 |
|
|
bool isTop(){
|
| 63 |
|
|
string tops = "TP TB";
|
| 64 |
|
|
string::size_type loc = tops.find( Pid, 0 );
|
| 65 |
|
|
if( loc != string::npos ) {
|
| 66 |
|
|
return true;
|
| 67 |
|
|
} else {
|
| 68 |
|
|
return false;
|
| 69 |
|
|
}
|
| 70 |
|
|
};
|
| 71 |
|
|
bool isLepton(){
|
| 72 |
|
|
string leptons = "E+ E- MU+ MU-";
|
| 73 |
|
|
string::size_type loc = leptons.find( Pid, 0 );
|
| 74 |
|
|
if( loc != string::npos ) {
|
| 75 |
|
|
return true;
|
| 76 |
|
|
} else {
|
| 77 |
|
|
return false;
|
| 78 |
|
|
}
|
| 79 |
|
|
};
|
| 80 |
|
|
// Taus are treated separatly since they could decay to hadronically or leptonically
|
| 81 |
|
|
bool isTau(){
|
| 82 |
|
|
string taus = "TAU+ TAU-";
|
| 83 |
|
|
string::size_type loc = taus.find( Pid, 0 );
|
| 84 |
|
|
if( loc != string::npos ) {
|
| 85 |
|
|
return true;
|
| 86 |
|
|
} else {
|
| 87 |
|
|
return false;
|
| 88 |
|
|
}
|
| 89 |
|
|
};
|
| 90 |
|
|
bool isSparticle(){
|
| 91 |
|
|
string sparticles = "UPL UBL DNL DBL STL SBL CHL CBL BT1 BB1 TP1 TB1 UPR UBR DNR DBR STR SBR CHR CBR BT2 BB2 TP2 TB2 NUEL ANUEL EL- EL+ NUML ANUML MUL- MUL+ NUTL ANUTL TAU1- TAU1+ NUER ANUER ER- ER+ NUMR ANUMR MUR- MUR+ NUTR ANUTR TAU2- TAU2+ GLSS Z1SS W1SS+ W1SS- Z2SS W2SS+ W2SS- Z3SS Z4SS GVSS GRAV";
|
| 92 |
|
|
string::size_type loc = sparticles.find( Pid, 0 );
|
| 93 |
|
|
if( loc != string::npos ) {
|
| 94 |
|
|
return true;
|
| 95 |
|
|
} else {
|
| 96 |
|
|
return false;
|
| 97 |
|
|
}
|
| 98 |
|
|
};
|
| 99 |
|
|
bool isHiggs(){
|
| 100 |
|
|
string higgs = "HL0 HH0 HA0 H+ H-";
|
| 101 |
|
|
string::size_type loc = higgs.find( Pid, 0 );
|
| 102 |
|
|
if( loc != string::npos ) {
|
| 103 |
|
|
return true;
|
| 104 |
|
|
} else {
|
| 105 |
|
|
return false;
|
| 106 |
|
|
}
|
| 107 |
|
|
};
|
| 108 |
|
|
bool isW(){
|
| 109 |
|
|
string Ws = "W+ W-";
|
| 110 |
|
|
string::size_type loc = Ws.find( Pid, 0 );
|
| 111 |
|
|
if( loc != string::npos ) {
|
| 112 |
|
|
return true;
|
| 113 |
|
|
} else {
|
| 114 |
|
|
return false;
|
| 115 |
|
|
}
|
| 116 |
|
|
};
|
| 117 |
|
|
bool isZ(){
|
| 118 |
|
|
string Z = "Z0";
|
| 119 |
|
|
string::size_type loc = Z.find( Pid, 0 );
|
| 120 |
|
|
if( loc != string::npos ) {
|
| 121 |
|
|
return true;
|
| 122 |
|
|
} else {
|
| 123 |
|
|
return false;
|
| 124 |
|
|
}
|
| 125 |
|
|
};
|
| 126 |
|
|
bool isBoson(){
|
| 127 |
|
|
string Bo = "W+ W- Z0 HL0";
|
| 128 |
|
|
string::size_type loc = Bo.find( Pid, 0 );
|
| 129 |
|
|
if( loc != string::npos ) {
|
| 130 |
|
|
return true;
|
| 131 |
|
|
} else {
|
| 132 |
|
|
return false;
|
| 133 |
|
|
}
|
| 134 |
|
|
};
|
| 135 |
|
|
bool isGl(){
|
| 136 |
|
|
string Gl = "GL";
|
| 137 |
|
|
string::size_type loc = Gl.find( Pid, 0 );
|
| 138 |
|
|
if( loc != string::npos ) {
|
| 139 |
|
|
return true;
|
| 140 |
|
|
} else {
|
| 141 |
|
|
return false;
|
| 142 |
|
|
}
|
| 143 |
|
|
};
|
| 144 |
|
|
bool isJet(){
|
| 145 |
|
|
string jet = "UP UB DN DB ST SB CH CB BT BB GL TAU+ TAU-";
|
| 146 |
|
|
string::size_type loc = jet.find( Pid, 0 );
|
| 147 |
|
|
if( loc != string::npos ) {
|
| 148 |
|
|
return true;
|
| 149 |
|
|
} else {
|
| 150 |
|
|
return false;
|
| 151 |
|
|
}
|
| 152 |
|
|
};
|
| 153 |
|
|
//constructor
|
| 154 |
|
|
particleClass(string input = "PP"){
|
| 155 |
|
|
setPid(input);
|
| 156 |
|
|
// define pairs of particle and antiparticle because some decays are not
|
| 157 |
|
|
// in the decay table since the corresponding antiparticle process has
|
| 158 |
|
|
// the same branching ratio
|
| 159 |
|
|
APid["PP"] = "PP";
|
| 160 |
|
|
APid["UP"] = "UB";
|
| 161 |
|
|
APid["UB"] = "UP";
|
| 162 |
|
|
APid["UPL"] = "UBL";
|
| 163 |
|
|
APid["UBL"] = "UPL";
|
| 164 |
|
|
APid["UPR"] = "UBR";
|
| 165 |
|
|
APid["UBR"] = "UPR";
|
| 166 |
|
|
APid["DN"] = "DB";
|
| 167 |
|
|
APid["DB"] = "DN";
|
| 168 |
|
|
APid["DNL"] = "DBL";
|
| 169 |
|
|
APid["DBL"] = "DNL";
|
| 170 |
|
|
APid["DNR"] = "DBR";
|
| 171 |
|
|
APid["DBR"] = "DNR";
|
| 172 |
|
|
APid["CH"] = "CB";
|
| 173 |
|
|
APid["CB"] = "CH";
|
| 174 |
|
|
APid["CHL"] = "CBL";
|
| 175 |
|
|
APid["CBL"] = "CHL";
|
| 176 |
|
|
APid["CHR"] = "CBR";
|
| 177 |
|
|
APid["CBR"] = "CHR";
|
| 178 |
|
|
APid["ST"] = "SB";
|
| 179 |
|
|
APid["SB"] = "ST";
|
| 180 |
|
|
APid["STL"] = "SBL";
|
| 181 |
|
|
APid["SBL"] = "STL";
|
| 182 |
|
|
APid["STR"] = "SBR";
|
| 183 |
|
|
APid["SBR"] = "STR";
|
| 184 |
|
|
APid["TP"] = "TB";
|
| 185 |
|
|
APid["TB"] = "TP";
|
| 186 |
|
|
APid["TP1"] = "TB1";
|
| 187 |
|
|
APid["TB1"] = "TP1";
|
| 188 |
|
|
APid["TP2"] = "TB2";
|
| 189 |
|
|
APid["TB2"] = "TP2";
|
| 190 |
|
|
APid["BT"] = "BB";
|
| 191 |
|
|
APid["BB"] = "BT";
|
| 192 |
|
|
APid["BT1"] = "BB1";
|
| 193 |
|
|
APid["BB1"] = "BT1";
|
| 194 |
|
|
APid["BT2"] = "BB2";
|
| 195 |
|
|
APid["BB2"] = "BT2";
|
| 196 |
|
|
APid["E+"] = "E-";
|
| 197 |
|
|
APid["E-"] = "E+";
|
| 198 |
|
|
APid["ER+"] = "ER-";
|
| 199 |
|
|
APid["ER-"] = "ER+";
|
| 200 |
|
|
APid["EL+"] = "EL-";
|
| 201 |
|
|
APid["EL-"] = "EL+";
|
| 202 |
|
|
APid["NUE"] = "ANUE";
|
| 203 |
|
|
APid["ANUE"] = "NUE";
|
| 204 |
|
|
APid["NUEL"] = "ANUEL";
|
| 205 |
|
|
APid["ANUEL"] = "NUEL";
|
| 206 |
|
|
APid["MU+"] = "MU-";
|
| 207 |
|
|
APid["MU-"] = "MU+";
|
| 208 |
|
|
APid["MUR+"] = "MUR-";
|
| 209 |
|
|
APid["MUR-"] = "MUR+";
|
| 210 |
|
|
APid["MUL+"] = "MUL-";
|
| 211 |
|
|
APid["MUL-"] = "MUL+";
|
| 212 |
|
|
APid["NUM"] = "ANUM";
|
| 213 |
|
|
APid["ANUM"] = "NUM";
|
| 214 |
|
|
APid["NUML"] = "ANUML";
|
| 215 |
|
|
APid["ANUML"] = "NUML";
|
| 216 |
|
|
APid["TAU+"] = "TAU-";
|
| 217 |
|
|
APid["TAU-"] = "TAU+";
|
| 218 |
|
|
APid["TAU1+"] = "TAU1-";
|
| 219 |
|
|
APid["TAU1-"] = "TAU1+";
|
| 220 |
|
|
APid["TAU2+"] = "TAU2-";
|
| 221 |
|
|
APid["TAU2-"] = "TAU2+";
|
| 222 |
|
|
APid["NUT"] = "ANUT";
|
| 223 |
|
|
APid["ANUT"] = "NUT";
|
| 224 |
|
|
APid["NUTL"] = "ANUTL";
|
| 225 |
|
|
APid["ANUTL"] = "NUTL";
|
| 226 |
|
|
APid["W+"] = "W-";
|
| 227 |
|
|
APid["W-"] = "W+";
|
| 228 |
|
|
APid["W1SS+"] = "W1SS-";
|
| 229 |
|
|
APid["W1SS-"] = "W1SS+";
|
| 230 |
|
|
APid["W2SS+"] = "W2SS-";
|
| 231 |
|
|
APid["W2SS-"] = "W2SS+";
|
| 232 |
|
|
APid["Z0"] = "Z0";
|
| 233 |
|
|
APid["Z1SS"] = "Z1SS";
|
| 234 |
|
|
APid["Z2SS"] = "Z2SS";
|
| 235 |
|
|
APid["Z3SS"] = "Z3SS";
|
| 236 |
|
|
APid["Z4SS"] = "Z4SS";
|
| 237 |
|
|
APid["GLSS"] = "GLSS";
|
| 238 |
|
|
APid["HL0"] = "HL0";
|
| 239 |
|
|
APid["HH0"] = "HH0";
|
| 240 |
|
|
APid["HA0"] = "HA0";
|
| 241 |
|
|
APid["H+"] = "H-";
|
| 242 |
|
|
APid["H-"] = "H+";
|
| 243 |
|
|
APid["GVSS"] = "GVSS";
|
| 244 |
|
|
APid["GRAV"] = "GRAV";
|
| 245 |
|
|
APid["GL"] = "GL";
|
| 246 |
|
|
APid["GM"] = "GM";
|
| 247 |
|
|
|
| 248 |
|
|
};
|
| 249 |
|
|
//destructor
|
| 250 |
|
|
~particleClass(){};
|
| 251 |
|
|
};
|
| 252 |
|
|
////////////////////////////////////////////////////////////////
|
| 253 |
|
|
|
| 254 |
|
|
////////////////////////////////////////////////////////////////
|
| 255 |
|
|
// This class describes a final state
|
| 256 |
|
|
////////////////////////////////////////////////////////////////
|
| 257 |
|
|
class finalStateClass
|
| 258 |
|
|
{
|
| 259 |
|
|
private:
|
| 260 |
|
|
particleClass mother;
|
| 261 |
|
|
int Nleptons;
|
| 262 |
|
|
int Nquarks;
|
| 263 |
|
|
int Ntops;
|
| 264 |
|
|
int Ntaus;
|
| 265 |
|
|
int Ngluons;
|
| 266 |
|
|
int NWs;
|
| 267 |
|
|
int NZs;
|
| 268 |
|
|
int Nbosons;
|
| 269 |
|
|
int Njets;
|
| 270 |
|
|
double BrRatio;
|
| 271 |
|
|
|
| 272 |
|
|
public:
|
| 273 |
|
|
// returns true is final state has same mother and same decay products as the wanted one
|
| 274 |
|
|
bool isFinalState(finalStateClass *wantedState){
|
| 275 |
|
|
int Nl=wantedState->getNleptons();
|
| 276 |
|
|
int Nq=wantedState->getNquarks();
|
| 277 |
|
|
int Ntps=wantedState->getNtops();
|
| 278 |
|
|
int Nt=wantedState->getNtaus();
|
| 279 |
|
|
int Ng=wantedState->getNgluons();
|
| 280 |
|
|
int NW=wantedState->getNWs();
|
| 281 |
|
|
int NZ=wantedState->getNZs();
|
| 282 |
|
|
int Nboson=wantedState->getNbosons();
|
| 283 |
|
|
int Njet=wantedState->getNjets();
|
| 284 |
|
|
if ((Nl==Nleptons || Nl<0) && (Nq==Nquarks || Nq<0) && (Ntps==Ntops || Ntps<0) && (Nt==Ntaus || Nt<0) && (Ng==Ngluons || Ng<0) && (NW==NWs || NW<0) && (NZ==NZs || NZ<0) && (Nboson==Nbosons || Nboson<0)&& (Njet==Njets || Njet<0) && wantedState->getMother().getPid()==mother.getPid()){
|
| 285 |
|
|
return true;
|
| 286 |
|
|
} else {
|
| 287 |
|
|
return false;
|
| 288 |
|
|
}
|
| 289 |
|
|
};
|
| 290 |
|
|
// returns true is final state has same mother and same decay products as the wanted one
|
| 291 |
|
|
bool isFinalStateExact(finalStateClass *wantedState){
|
| 292 |
|
|
int Nl=wantedState->getNleptons();
|
| 293 |
|
|
int Nq=wantedState->getNquarks();
|
| 294 |
|
|
int Ntps=wantedState->getNtops();
|
| 295 |
|
|
int Nt=wantedState->getNtaus();
|
| 296 |
|
|
int Ng=wantedState->getNgluons();
|
| 297 |
|
|
int NW=wantedState->getNWs();
|
| 298 |
|
|
int NZ=wantedState->getNZs();
|
| 299 |
|
|
int Nboson=wantedState->getNbosons();
|
| 300 |
|
|
int Njet=wantedState->getNjets();
|
| 301 |
|
|
if ((Nl==Nleptons) && (Nq==Nquarks) && (Ntps==Ntops) && (Nt==Ntaus) && (Ng==Ngluons) && (NW==NWs) && (NZ==NZs) && (Nboson==Nbosons) && (Njet==Njets) && wantedState->getMother().getPid()==mother.getPid()){
|
| 302 |
|
|
return true;
|
| 303 |
|
|
} else {
|
| 304 |
|
|
return false;
|
| 305 |
|
|
}
|
| 306 |
|
|
};
|
| 307 |
|
|
// returns true is final state has still the chance to become the wanted one
|
| 308 |
|
|
bool isPossible(finalStateClass *wantedState){
|
| 309 |
|
|
int Nl=wantedState->getNleptons();
|
| 310 |
|
|
int Nq=wantedState->getNquarks();
|
| 311 |
|
|
int Ntps=wantedState->getNtops();
|
| 312 |
|
|
int Nt=wantedState->getNtaus();
|
| 313 |
|
|
int Ng=wantedState->getNgluons();
|
| 314 |
|
|
int NW=wantedState->getNWs();
|
| 315 |
|
|
int NZ=wantedState->getNZs();
|
| 316 |
|
|
int Nboson=wantedState->getNbosons();
|
| 317 |
|
|
int Njet=wantedState->getNjets();
|
| 318 |
|
|
if ((Nl>=Nleptons || Nl<0) && (Nq>=Nquarks || Nq<0) && (Ntps>=Ntops || Ntps<0)&& (Nt>=Ntaus || Nt<0) && (Ng>=Ngluons || Ng<0) && (NW>=NWs || NW<0) && (NZ>=NZs || NZ<0) && (Nboson>=Nbosons || Nboson<0) && (Njet>=Njets || Njet<0) && wantedState->getMother().getPid()==mother.getPid()){
|
| 319 |
|
|
return true;
|
| 320 |
|
|
} else {
|
| 321 |
|
|
return false;
|
| 322 |
|
|
}
|
| 323 |
|
|
};
|
| 324 |
|
|
void incNleptons(){
|
| 325 |
|
|
++Nleptons;
|
| 326 |
|
|
};
|
| 327 |
|
|
void incNquarks(){
|
| 328 |
|
|
++Nquarks;
|
| 329 |
|
|
};
|
| 330 |
|
|
void incNtops(){
|
| 331 |
|
|
++Ntops;
|
| 332 |
|
|
};
|
| 333 |
|
|
void incNtaus(){
|
| 334 |
|
|
++Ntaus;
|
| 335 |
|
|
};
|
| 336 |
|
|
void incNgluons(){
|
| 337 |
|
|
++Ngluons;
|
| 338 |
|
|
};
|
| 339 |
|
|
void incNWs(){
|
| 340 |
|
|
++NWs;
|
| 341 |
|
|
};
|
| 342 |
|
|
void incNZs(){
|
| 343 |
|
|
++NZs;
|
| 344 |
|
|
};
|
| 345 |
|
|
void incNbosons(){
|
| 346 |
|
|
++Nbosons;
|
| 347 |
|
|
};
|
| 348 |
|
|
void incNjets(){
|
| 349 |
|
|
++Njets;
|
| 350 |
|
|
};
|
| 351 |
|
|
void setNleptons(int input){
|
| 352 |
|
|
Nleptons=input;
|
| 353 |
|
|
};
|
| 354 |
|
|
void setNquarks(int input){
|
| 355 |
|
|
Nquarks=input;
|
| 356 |
|
|
};
|
| 357 |
|
|
void setNtops(int input){
|
| 358 |
|
|
Ntops=input;
|
| 359 |
|
|
};
|
| 360 |
|
|
void setNtaus(int input){
|
| 361 |
|
|
Ntaus=input;
|
| 362 |
|
|
};
|
| 363 |
|
|
void setNgluons(int input){
|
| 364 |
|
|
Ngluons=input;
|
| 365 |
|
|
};
|
| 366 |
|
|
void setNWs(int input){
|
| 367 |
|
|
NWs=input;
|
| 368 |
|
|
};
|
| 369 |
|
|
void setNZs(int input){
|
| 370 |
|
|
NZs=input;
|
| 371 |
|
|
};
|
| 372 |
|
|
void setNbosons(int input){
|
| 373 |
|
|
Nbosons=input;
|
| 374 |
|
|
};
|
| 375 |
|
|
void setNjets(int input){
|
| 376 |
|
|
Njets=input;
|
| 377 |
|
|
};
|
| 378 |
|
|
int getNleptons(){
|
| 379 |
|
|
return Nleptons;
|
| 380 |
|
|
};
|
| 381 |
|
|
int getNquarks(){
|
| 382 |
|
|
return Nquarks;
|
| 383 |
|
|
};
|
| 384 |
|
|
int getNtops(){
|
| 385 |
|
|
return Ntops;
|
| 386 |
|
|
};
|
| 387 |
|
|
int getNtaus(){
|
| 388 |
|
|
return Ntaus;
|
| 389 |
|
|
};
|
| 390 |
|
|
int getNgluons(){
|
| 391 |
|
|
return Ngluons;
|
| 392 |
|
|
};
|
| 393 |
|
|
int getNWs(){
|
| 394 |
|
|
return NWs;
|
| 395 |
|
|
};
|
| 396 |
|
|
int getNZs(){
|
| 397 |
|
|
return NZs;
|
| 398 |
|
|
};
|
| 399 |
|
|
int getNbosons(){
|
| 400 |
|
|
return Nbosons;
|
| 401 |
|
|
};
|
| 402 |
|
|
int getNjets(){
|
| 403 |
|
|
return Njets;
|
| 404 |
|
|
};
|
| 405 |
|
|
double getBrRatio(){
|
| 406 |
|
|
return BrRatio;
|
| 407 |
|
|
};
|
| 408 |
|
|
void setBrRatio(double input){
|
| 409 |
|
|
BrRatio=input;
|
| 410 |
|
|
};
|
| 411 |
|
|
void setMotherId(string input){
|
| 412 |
|
|
mother.setPid(input);
|
| 413 |
|
|
};
|
| 414 |
|
|
particleClass getMother(){
|
| 415 |
|
|
return mother;
|
| 416 |
|
|
};
|
| 417 |
|
|
// operator to subtract final states from each other
|
| 418 |
|
|
// the mother of the resulting FS is set to the one of the second summand
|
| 419 |
|
|
finalStateClass operator-(finalStateClass FS2){
|
| 420 |
|
|
int Nl2=FS2.getNleptons();
|
| 421 |
|
|
int Nq2=FS2.getNquarks();
|
| 422 |
|
|
int Ntps2=FS2.getNtops();
|
| 423 |
|
|
int Nt2=FS2.getNtaus();
|
| 424 |
|
|
int Ng2=FS2.getNgluons();
|
| 425 |
|
|
int NW2=FS2.getNWs();
|
| 426 |
|
|
int NZ2=FS2.getNZs();
|
| 427 |
|
|
int Nboson2=FS2.getNbosons();
|
| 428 |
|
|
int Njet2=FS2.getNjets();
|
| 429 |
|
|
finalStateClass FSresult;
|
| 430 |
|
|
|
| 431 |
|
|
if (Nleptons>0){
|
| 432 |
|
|
FSresult.setNleptons(Nleptons-Nl2);
|
| 433 |
|
|
} else {
|
| 434 |
|
|
FSresult.setNleptons(Nleptons);
|
| 435 |
|
|
};
|
| 436 |
|
|
if (Nquarks>0){
|
| 437 |
|
|
FSresult.setNquarks(Nquarks-Nq2);
|
| 438 |
|
|
} else {
|
| 439 |
|
|
FSresult.setNquarks(Nquarks);
|
| 440 |
|
|
};
|
| 441 |
|
|
if (Ntops>0){
|
| 442 |
|
|
FSresult.setNtops(Ntops-Ntps2);
|
| 443 |
|
|
} else {
|
| 444 |
|
|
FSresult.setNtops(Ntops);
|
| 445 |
|
|
};
|
| 446 |
|
|
if (Ntaus>0){
|
| 447 |
|
|
FSresult.setNtaus(Ntaus-Nt2);
|
| 448 |
|
|
} else {
|
| 449 |
|
|
FSresult.setNtaus(Ntaus);
|
| 450 |
|
|
};
|
| 451 |
|
|
if (Ngluons>0){
|
| 452 |
|
|
FSresult.setNgluons(Ngluons-Nl2);
|
| 453 |
|
|
} else {
|
| 454 |
|
|
FSresult.setNgluons(Ngluons);
|
| 455 |
|
|
};
|
| 456 |
|
|
if (NWs>0){
|
| 457 |
|
|
FSresult.setNWs(NWs-NW2);
|
| 458 |
|
|
} else {
|
| 459 |
|
|
FSresult.setNWs(NWs);
|
| 460 |
|
|
};
|
| 461 |
|
|
if (NZs>0){
|
| 462 |
|
|
FSresult.setNZs(NZs-NZ2);
|
| 463 |
|
|
} else {
|
| 464 |
|
|
FSresult.setNZs(NZs);
|
| 465 |
|
|
};
|
| 466 |
|
|
if (Nbosons>0){
|
| 467 |
|
|
FSresult.setNbosons(Nbosons-Nboson2);
|
| 468 |
|
|
} else {
|
| 469 |
|
|
FSresult.setNbosons(Nbosons);
|
| 470 |
|
|
};
|
| 471 |
|
|
if (Njets>0){
|
| 472 |
|
|
FSresult.setNjets(Njets-Njet2);
|
| 473 |
|
|
} else {
|
| 474 |
|
|
FSresult.setNjets(Njets);
|
| 475 |
|
|
};
|
| 476 |
|
|
|
| 477 |
|
|
FSresult.setMotherId(FS2.getMother().getPid());
|
| 478 |
|
|
FSresult.setBrRatio(FS2.getBrRatio());
|
| 479 |
|
|
return FSresult;
|
| 480 |
|
|
};
|
| 481 |
|
|
void printFinalState(){
|
| 482 |
|
|
cout<<mother.getPid()<<" --> Nls= "<<Nleptons<<" Nqs= "<<Nquarks<<" Ntops= "<<Ntops<<" Ntau= "<<Ntaus<<" Ngl= "<<Ngluons<<" NWs= "<<NWs<<" NZs= "<<NZs<<" Nbosons= "<<Nbosons<<" Njets= "<<Njets<<" Br= "<<BrRatio<<endl;
|
| 483 |
|
|
};
|
| 484 |
|
|
//constructor
|
| 485 |
|
|
finalStateClass(string firstMotherId="GLSS", int Nl=0, int Nq=0, int Ntps=0, int Nt=0, int Ng=0, int NW=0, int NZ=0, int Nboson=0, int Njet=0){
|
| 486 |
|
|
mother.setPid(firstMotherId);
|
| 487 |
|
|
BrRatio=1.; Nleptons=Nl; Nquarks=Nq; Ntops=Ntps; Ntaus=Nt; Ngluons=Ng; NWs=NW; NZs=NZ; Nbosons=Nboson; Njets=Njet;
|
| 488 |
|
|
};
|
| 489 |
|
|
//destructor
|
| 490 |
|
|
~finalStateClass(){};
|
| 491 |
|
|
};
|
| 492 |
|
|
////////////////////////////////////////////////////////////////
|
| 493 |
|
|
|
| 494 |
|
|
////////////////////////////////////////////////////////////////
|
| 495 |
|
|
// This class describes a decay
|
| 496 |
|
|
////////////////////////////////////////////////////////////////
|
| 497 |
|
|
class decayClass
|
| 498 |
|
|
{
|
| 499 |
|
|
private:
|
| 500 |
|
|
particleClass mother;
|
| 501 |
|
|
vector<particleClass> daughter;
|
| 502 |
|
|
double BrRatio;
|
| 503 |
|
|
|
| 504 |
|
|
public:
|
| 505 |
|
|
particleClass getMother(){
|
| 506 |
|
|
return mother;
|
| 507 |
|
|
};
|
| 508 |
|
|
vector<particleClass> getDaughter(){
|
| 509 |
|
|
return daughter;
|
| 510 |
|
|
};
|
| 511 |
|
|
void addDaughter(string Did){
|
| 512 |
|
|
particleClass tmp(Did);
|
| 513 |
|
|
daughter.push_back(tmp);
|
| 514 |
|
|
};
|
| 515 |
|
|
double getBrRatio(){
|
| 516 |
|
|
return BrRatio;
|
| 517 |
|
|
};
|
| 518 |
|
|
void printDecay(){
|
| 519 |
|
|
cout<<getMother().getPid()<<" --> ";
|
| 520 |
|
|
vector<particleClass>::iterator iter;
|
| 521 |
|
|
for (iter=daughter.begin(); iter!=daughter.end();++iter){
|
| 522 |
|
|
cout<<(*iter).getPid()<<" ";
|
| 523 |
|
|
}
|
| 524 |
|
|
cout<<" branching ratio: "<<getBrRatio()<<endl;
|
| 525 |
|
|
};
|
| 526 |
|
|
//constructor
|
| 527 |
|
|
decayClass(string inputMother = "GLSS", double inputBR=0.){
|
| 528 |
|
|
mother.setPid(inputMother);
|
| 529 |
|
|
BrRatio=inputBR;
|
| 530 |
|
|
};
|
| 531 |
|
|
//destructor
|
| 532 |
|
|
~decayClass(){};
|
| 533 |
|
|
};
|
| 534 |
|
|
////////////////////////////////////////////////////////////////
|
| 535 |
|
|
|
| 536 |
|
|
////////////////////////////////////////////////////////////////
|
| 537 |
|
|
// This class describes a decay table
|
| 538 |
|
|
////////////////////////////////////////////////////////////////
|
| 539 |
|
|
class decayTableClass
|
| 540 |
|
|
{
|
| 541 |
|
|
private:
|
| 542 |
|
|
static const double eps=1.e-3; //smaller Br are neglected
|
| 543 |
|
|
vector<decayClass> decayTable;
|
| 544 |
|
|
vector<finalStateClass> availableFinalStates;
|
| 545 |
|
|
// returns a list of decays with same mother particle
|
| 546 |
|
|
vector<decayClass> getDecaysByMotherId(particleClass mother){
|
| 547 |
|
|
vector<decayClass> decayTableByMother;
|
| 548 |
|
|
// exclude particle which should not decay further
|
| 549 |
|
|
if (!mother.isBoson() && !mother.isTop()){
|
| 550 |
|
|
vector<decayClass>::iterator iter;
|
| 551 |
|
|
for (iter=decayTable.begin(); iter!=decayTable.end();++iter){
|
| 552 |
|
|
if ((*iter).getMother().getPid()==mother.getPid()){
|
| 553 |
|
|
decayTableByMother.push_back(*iter);
|
| 554 |
|
|
}
|
| 555 |
|
|
}
|
| 556 |
|
|
if (decayTableByMother.empty()){
|
| 557 |
|
|
for (iter=decayTable.begin(); iter!=decayTable.end();++iter){
|
| 558 |
|
|
if ((*iter).getMother().getAPid()==mother.getPid()){
|
| 559 |
|
|
decayTableByMother.push_back(*iter);
|
| 560 |
|
|
}
|
| 561 |
|
|
}
|
| 562 |
|
|
}
|
| 563 |
|
|
}
|
| 564 |
|
|
return decayTableByMother;
|
| 565 |
|
|
};
|
| 566 |
|
|
// returns a list of final states which has each final object
|
| 567 |
|
|
// number less or equal
|
| 568 |
|
|
vector<finalStateClass> smallerFinalStates(finalStateClass wantedFS){
|
| 569 |
|
|
vector<finalStateClass> FinalStates;
|
| 570 |
|
|
string motherId = wantedFS.getMother().getPid();
|
| 571 |
|
|
for (int i=0;i<=max(0,wantedFS.getNleptons());++i){
|
| 572 |
|
|
for (int j=0;j<=max(0,wantedFS.getNquarks());++j){
|
| 573 |
|
|
for (int k=0;k<=max(0,wantedFS.getNtaus());++k){
|
| 574 |
|
|
for (int l=0;l<=max(0,wantedFS.getNgluons());++l){
|
| 575 |
|
|
for (int m=0;m<=max(0,wantedFS.getNWs());++m){
|
| 576 |
|
|
for (int n=0;n<=max(0,wantedFS.getNZs());++n){
|
| 577 |
|
|
for (int o=0;o<=max(0,wantedFS.getNbosons());++o){
|
| 578 |
|
|
for (int p=0;p<=max(0,wantedFS.getNjets());++p){
|
| 579 |
|
|
for (int q=0;q<=max(0,wantedFS.getNtops());++q){
|
| 580 |
|
|
int wleptons=i;
|
| 581 |
|
|
if (wantedFS.getNleptons()<0) wleptons=-1;
|
| 582 |
|
|
int wquarks=j;
|
| 583 |
|
|
if (wantedFS.getNquarks()<0) wquarks=-1;
|
| 584 |
|
|
int wtaus=k;
|
| 585 |
|
|
if (wantedFS.getNtaus()<0) wtaus=-1;
|
| 586 |
|
|
int wgluons=l;
|
| 587 |
|
|
if (wantedFS.getNgluons()<0) wgluons=-1;
|
| 588 |
|
|
int wWs=m;
|
| 589 |
|
|
if (wantedFS.getNWs()<0) wWs=-1;
|
| 590 |
|
|
int wZs=n;
|
| 591 |
|
|
if (wantedFS.getNZs()<0) wZs=-1;
|
| 592 |
|
|
int wbosons=o;
|
| 593 |
|
|
if (wantedFS.getNbosons()<0) wbosons=-1;
|
| 594 |
|
|
int wjets=p;
|
| 595 |
|
|
if (wantedFS.getNjets()<0) wjets=-1;
|
| 596 |
|
|
int wtops=q;
|
| 597 |
|
|
if (wantedFS.getNtops()<0) wtops=-1;
|
| 598 |
|
|
|
| 599 |
|
|
finalStateClass FStmp(motherId, wleptons, wquarks, wtops, wtaus, wgluons, wWs, wZs, wbosons, wjets);
|
| 600 |
|
|
FinalStates.push_back(FStmp);
|
| 601 |
|
|
}
|
| 602 |
|
|
}
|
| 603 |
|
|
}
|
| 604 |
|
|
}
|
| 605 |
|
|
}
|
| 606 |
|
|
}
|
| 607 |
|
|
}
|
| 608 |
|
|
}
|
| 609 |
|
|
}
|
| 610 |
|
|
return FinalStates;
|
| 611 |
|
|
};
|
| 612 |
|
|
// returns branching ratio of wanted final state
|
| 613 |
|
|
double BrRatioByFinalState(finalStateClass wantedFS){
|
| 614 |
|
|
//cout<<"wanted FinalState: ";
|
| 615 |
|
|
//wantedFS.printFinalState();
|
| 616 |
|
|
double returnValue;
|
| 617 |
|
|
if (getBRbyFinalState(wantedFS)){
|
| 618 |
|
|
returnValue=wantedFS.getBrRatio();
|
| 619 |
|
|
} else {
|
| 620 |
|
|
string MotherId=wantedFS.getMother().getPid();
|
| 621 |
|
|
// BrRatio of input final state
|
| 622 |
|
|
double BR=0.;
|
| 623 |
|
|
vector<decayClass> decayTableByMother = getDecaysByMotherId(wantedFS.getMother());
|
| 624 |
|
|
// if there exists at least one decay mode
|
| 625 |
|
|
if (!decayTableByMother.empty()){
|
| 626 |
|
|
vector<decayClass>::iterator i;
|
| 627 |
|
|
// loop over decay modes
|
| 628 |
|
|
for (i=decayTableByMother.begin(); i!=decayTableByMother.end();++i){
|
| 629 |
|
|
vector<particleClass> daughters=(*i).getDaughter();
|
| 630 |
|
|
vector<particleClass>::iterator j;
|
| 631 |
|
|
finalStateClass FSnew;
|
| 632 |
|
|
FSnew.setMotherId(MotherId);
|
| 633 |
|
|
// loop over daughters of each decay mode
|
| 634 |
|
|
bool onlyFinals = true;
|
| 635 |
|
|
for (j=daughters.begin(); j!=daughters.end(); ++j){
|
| 636 |
|
|
vector<decayClass> decayTableOfDaughter = getDecaysByMotherId(*j);
|
| 637 |
|
|
if (decayTableOfDaughter.empty()){
|
| 638 |
|
|
// increment final state counters
|
| 639 |
|
|
if ((*j).isLepton()) {FSnew.incNleptons();}
|
| 640 |
|
|
if ((*j).isQuark()) {FSnew.incNquarks();}
|
| 641 |
|
|
if ((*j).isTau()) {FSnew.incNtaus();}
|
| 642 |
|
|
if ((*j).isGl()) {FSnew.incNgluons();}
|
| 643 |
|
|
if ((*j).isJet()) {FSnew.incNjets();}
|
| 644 |
|
|
// treated as stable particles (see getDecaysByMotherId)
|
| 645 |
|
|
if ((*j).isTop()) {FSnew.incNtops();}
|
| 646 |
|
|
if ((*j).isW()) {FSnew.incNWs();}
|
| 647 |
|
|
if ((*j).isZ()) {FSnew.incNZs();}
|
| 648 |
|
|
if ((*j).isBoson()) {FSnew.incNbosons();}
|
| 649 |
|
|
} else {
|
| 650 |
|
|
// treated as unstable particles (see getDecaysByMotherId)
|
| 651 |
|
|
//if ((*j).isTop()) {FSnew.incNtops();}
|
| 652 |
|
|
//if ((*j).isW()) {FSnew.incNWs();}
|
| 653 |
|
|
//if ((*j).isZ()) {FSnew.incNZs();}
|
| 654 |
|
|
//if ((*j).isBoson()) {FSnew.incNbosons();}
|
| 655 |
|
|
onlyFinals=false;
|
| 656 |
|
|
}
|
| 657 |
|
|
}
|
| 658 |
|
|
if (onlyFinals){
|
| 659 |
|
|
// check if final state is wanted final state
|
| 660 |
|
|
if (FSnew.isFinalState(&wantedFS)){
|
| 661 |
|
|
BR+=(*i).getBrRatio();
|
| 662 |
|
|
}
|
| 663 |
|
|
} else {
|
| 664 |
|
|
// check if wanted final state is still possible
|
| 665 |
|
|
if (FSnew.isPossible(&wantedFS)){
|
| 666 |
|
|
// MotherId of wantedFStmp has not to be set since this is done in BrRatioByFinalStateFromParticleVector
|
| 667 |
|
|
if (!daughters.empty()){
|
| 668 |
|
|
// calculate BrRatio to get wanted final state from all daughters and multiply with BrRatio of current decay mode
|
| 669 |
|
|
finalStateClass wantedFSnew = wantedFS;
|
| 670 |
|
|
wantedFSnew=wantedFS-FSnew;
|
| 671 |
|
|
double BRtmp = (*i).getBrRatio();
|
| 672 |
|
|
if (BRtmp>eps){
|
| 673 |
|
|
BR+=BRtmp*BrRatioByFinalStateFromParticleVector(wantedFSnew, &daughters);
|
| 674 |
|
|
}
|
| 675 |
|
|
}
|
| 676 |
|
|
}
|
| 677 |
|
|
}
|
| 678 |
|
|
}
|
| 679 |
|
|
wantedFS.setBrRatio(BR);
|
| 680 |
|
|
addFinalState(wantedFS);
|
| 681 |
|
|
returnValue=BR;
|
| 682 |
|
|
} else {
|
| 683 |
|
|
// if there exist no decay mode, check if particle is already wanted final state
|
| 684 |
|
|
finalStateClass FSsingle;
|
| 685 |
|
|
FSsingle.setMotherId(wantedFS.getMother().getPid());
|
| 686 |
|
|
// only to be used, if first mother id stable
|
| 687 |
|
|
// this should never be the case
|
| 688 |
|
|
//if (wantedFS.getMother().isLepton()) {FSsingle.incNleptons();}
|
| 689 |
|
|
//if (wantedFS.getMother().isQuark()) {FSsingle.incNquarks();}
|
| 690 |
|
|
//if (wantedFS.getMother().isTau()) {FSsingle.incNtaus();}
|
| 691 |
|
|
//if (wantedFS.getMother().isGl()) {FSsingle.incNgluons();}
|
| 692 |
|
|
//if (wantedFS.getMother().isJet()) {FSsingle.incNjets();}
|
| 693 |
|
|
//if (wantedFS.getMother().isTop()) {FSsingle.incNtops();}
|
| 694 |
|
|
//if (wantedFS.getMother().isW()) {FSsingle.incNWs();}
|
| 695 |
|
|
//if (wantedFS.getMother().isZ()) {FSsingle.incNZs();}
|
| 696 |
|
|
//if (wantedFS.getMother().isBoson()) {FSsingle.incNbosons();}
|
| 697 |
|
|
if (FSsingle.isFinalState(&wantedFS)){
|
| 698 |
|
|
BR=1.;
|
| 699 |
|
|
}
|
| 700 |
|
|
returnValue=BR;
|
| 701 |
|
|
}
|
| 702 |
|
|
}
|
| 703 |
|
|
return returnValue;
|
| 704 |
|
|
};
|
| 705 |
|
|
double BrRatioByFinalStateFromParticleVector(finalStateClass wantedFS, vector<particleClass> *daughters){
|
| 706 |
|
|
double BR=0.;
|
| 707 |
|
|
if (!daughters->empty()){
|
| 708 |
|
|
// take the last daughter and remove it from the list of the remaining daughter
|
| 709 |
|
|
particleClass Part = daughters->back();
|
| 710 |
|
|
wantedFS.setMotherId(Part.getPid());
|
| 711 |
|
|
vector<particleClass> daughtersNew = *daughters;
|
| 712 |
|
|
daughtersNew.pop_back();
|
| 713 |
|
|
if (!daughtersNew.empty()){
|
| 714 |
|
|
vector<finalStateClass> FinalStates = smallerFinalStates(wantedFS);
|
| 715 |
|
|
if (!FinalStates.empty()){
|
| 716 |
|
|
vector<finalStateClass>::iterator iter;
|
| 717 |
|
|
// loop over all combinations of final states of single daughter and list of remaining daughters
|
| 718 |
|
|
// which gives as a sum the wanted FS
|
| 719 |
|
|
for (iter=FinalStates.begin(); iter!=FinalStates.end(); ++iter){
|
| 720 |
|
|
finalStateClass wantedFSnew = wantedFS;
|
| 721 |
|
|
wantedFSnew=wantedFS-(*iter);
|
| 722 |
|
|
// MotherId of wantedFStmp has not to be set since this is done in BrRatioByFinalStateFromParticleVector
|
| 723 |
|
|
double BRtmp=BrRatioByFinalState((*iter));
|
| 724 |
|
|
if (BRtmp>eps){
|
| 725 |
|
|
BR+=BRtmp*BrRatioByFinalStateFromParticleVector(wantedFSnew, &daughtersNew);
|
| 726 |
|
|
}
|
| 727 |
|
|
}
|
| 728 |
|
|
}
|
| 729 |
|
|
} else {
|
| 730 |
|
|
BR=BrRatioByFinalState(wantedFS);
|
| 731 |
|
|
}
|
| 732 |
|
|
} else {
|
| 733 |
|
|
// no daughters available should never be the case
|
| 734 |
|
|
BR=1.;
|
| 735 |
|
|
cout<<"You should never read this!"<<endl;
|
| 736 |
|
|
}
|
| 737 |
|
|
return BR;
|
| 738 |
|
|
};
|
| 739 |
|
|
void addFinalState(finalStateClass input){
|
| 740 |
|
|
//cout<<"Add FinalState: ";
|
| 741 |
|
|
//input.printFinalState();
|
| 742 |
|
|
availableFinalStates.push_back(input);
|
| 743 |
|
|
};
|
| 744 |
|
|
bool getBRbyFinalState(finalStateClass &wantedFS){
|
| 745 |
|
|
vector<finalStateClass>::iterator iter;
|
| 746 |
|
|
bool returnValue=false;
|
| 747 |
|
|
for (iter=availableFinalStates.begin(); iter!=availableFinalStates.end();++iter){
|
| 748 |
|
|
if ((*iter).isFinalStateExact(&wantedFS)){
|
| 749 |
|
|
wantedFS.setBrRatio((*iter).getBrRatio());
|
| 750 |
|
|
//cout<<"Is avaliable: ";
|
| 751 |
|
|
//wantedFS.printFinalState();
|
| 752 |
|
|
returnValue=true;
|
| 753 |
|
|
}
|
| 754 |
|
|
}
|
| 755 |
|
|
return returnValue;
|
| 756 |
|
|
};
|
| 757 |
|
|
|
| 758 |
|
|
public:
|
| 759 |
|
|
void addDecay(decayClass decay){
|
| 760 |
|
|
decayTable.push_back(decay);
|
| 761 |
|
|
}
|
| 762 |
|
|
double BrRatioByDecayProducts(string firstMotherId, int Nleptons, int Nquarks, int Ntops, int Ntaus, int Ngluons, int NWs, int NZs, int Nbosons, int Njets){
|
| 763 |
|
|
// set wanted final state
|
| 764 |
|
|
finalStateClass wantedFS(firstMotherId, Nleptons, Nquarks, Ntops, Ntaus, Ngluons, NWs, NZs, Nbosons, Njets);
|
| 765 |
|
|
return BrRatioByFinalState(wantedFS);
|
| 766 |
|
|
};
|
| 767 |
|
|
void printDecayTable(){
|
| 768 |
|
|
vector<decayClass>::iterator iter;
|
| 769 |
|
|
for (iter=decayTable.begin(); iter!=decayTable.end();++iter){
|
| 770 |
|
|
(*iter).printDecay();
|
| 771 |
|
|
}
|
| 772 |
|
|
};
|
| 773 |
|
|
void printFinalStateTable(){
|
| 774 |
|
|
vector<finalStateClass>::iterator iter;
|
| 775 |
|
|
for (iter=availableFinalStates.begin(); iter!=availableFinalStates.end();++iter){
|
| 776 |
|
|
(*iter).printFinalState();
|
| 777 |
|
|
};
|
| 778 |
|
|
};
|
| 779 |
|
|
decayTableClass(){}; // constructor
|
| 780 |
|
|
~decayTableClass(){}; // destructor
|
| 781 |
|
|
};
|
| 782 |
|
|
////////////////////////////////////////////////////////////////
|
| 783 |
|
|
|
| 784 |
|
|
////////////////////////////////////////////////////////////////
|
| 785 |
|
|
// This class describes a subprocess table
|
| 786 |
|
|
////////////////////////////////////////////////////////////////
|
| 787 |
|
|
class processTableClass
|
| 788 |
|
|
{
|
| 789 |
|
|
private:
|
| 790 |
|
|
map<const int, const double> processMap;
|
| 791 |
|
|
|
| 792 |
|
|
public:
|
| 793 |
|
|
void addProcess(int id=0, double xs=1.0){
|
| 794 |
|
|
processMap.insert(make_pair(id,xs));
|
| 795 |
|
|
}
|
| 796 |
|
|
double getProcessWeight(int id=0){
|
| 797 |
|
|
double result=0.;
|
| 798 |
|
|
map<const int, const double>::iterator iter=processMap.find(id);
|
| 799 |
|
|
if (iter != processMap.end()){
|
| 800 |
|
|
result=processMap.find(id)->second;
|
| 801 |
|
|
result/=processMap.find(0)->second;
|
| 802 |
|
|
}
|
| 803 |
|
|
return result;
|
| 804 |
|
|
}
|
| 805 |
|
|
void printProcessTable(){
|
| 806 |
|
|
map<const int, const double>::iterator iter;
|
| 807 |
|
|
for (iter=processMap.begin(); iter!=processMap.end();++iter){
|
| 808 |
|
|
int pid=(*iter).first;
|
| 809 |
|
|
double wq=(*iter).second;
|
| 810 |
|
|
cout<<" Subprocess Id: "<<pid<<" Xsec (mb) = "<<wq<<endl;
|
| 811 |
|
|
}
|
| 812 |
|
|
}
|
| 813 |
|
|
double getParticlesPerEvent(string particle = "GLSS"){
|
| 814 |
|
|
vector<int> subprocesses;
|
| 815 |
|
|
|
| 816 |
|
|
//define processes (pythia id) which produces wanted particle
|
| 817 |
|
|
//listed twice if pair production
|
| 818 |
|
|
|
| 819 |
|
|
if (particle=="EL+" || particle == "EL-" || particle == "EL"){
|
| 820 |
|
|
int processes[] = {201,201,203,210};
|
| 821 |
|
|
subprocesses.assign(processes,processes+4);
|
| 822 |
|
|
}
|
| 823 |
|
|
if (particle=="ER+" || particle == "ER-" || particle == "ER"){
|
| 824 |
|
|
int processes[] = {202,202,203};
|
| 825 |
|
|
subprocesses.assign(processes,processes+3);
|
| 826 |
|
|
}
|
| 827 |
|
|
if (particle=="MUL+" || particle == "MUL-" || particle == "MUL"){
|
| 828 |
|
|
int processes[] = {204,204,206,210};
|
| 829 |
|
|
subprocesses.assign(processes,processes+4);
|
| 830 |
|
|
}
|
| 831 |
|
|
if (particle=="MUR+" || particle == "MUR-" || particle == "MUR"){
|
| 832 |
|
|
int processes[] = {205,205,206};
|
| 833 |
|
|
subprocesses.assign(processes,processes+3);
|
| 834 |
|
|
}
|
| 835 |
|
|
if (particle=="TAU1+" || particle == "TAU1-" || particle == "TAU1"){
|
| 836 |
|
|
int processes[] = {207,207,209};
|
| 837 |
|
|
subprocesses.assign(processes,processes+3);
|
| 838 |
|
|
}
|
| 839 |
|
|
if (particle=="TAU2+" || particle == "TAU2-" || particle == "TAU2"){
|
| 840 |
|
|
int processes[] = {208,208,209};
|
| 841 |
|
|
subprocesses.assign(processes,processes+3);
|
| 842 |
|
|
}
|
| 843 |
|
|
if (particle=="NUEL" || particle == "ANUEL" ||
|
| 844 |
|
|
particle=="NUER" || particle == "ANUER" ||
|
| 845 |
|
|
particle=="NUML" || particle == "ANUML" ||
|
| 846 |
|
|
particle=="NUMR" || particle == "ANUMR" || particle == "NULEP"){
|
| 847 |
|
|
int processes[] = {210,213,213};
|
| 848 |
|
|
subprocesses.assign(processes,processes+3);
|
| 849 |
|
|
}
|
| 850 |
|
|
if (particle=="NUTL" || particle == "ANUTL" ||
|
| 851 |
|
|
particle=="NUTR" || particle == "ANUTR" || particle == "NUTAU"){
|
| 852 |
|
|
int processes[] = {211,212,214,214};
|
| 853 |
|
|
subprocesses.assign(processes,processes+4);
|
| 854 |
|
|
}
|
| 855 |
|
|
if (particle=="GLSS"){
|
| 856 |
|
|
int processes[] = {237,238,239,240,241,242,243,243,244,244,
|
| 857 |
|
|
258,259,294,295};
|
| 858 |
|
|
subprocesses.assign(processes,processes+14);
|
| 859 |
|
|
}
|
| 860 |
|
|
if (particle=="Z1SS"){
|
| 861 |
|
|
int processes[] = {216,216,220,221,222,229,233,237,246,247};
|
| 862 |
|
|
subprocesses.assign(processes,processes+10);
|
| 863 |
|
|
}
|
| 864 |
|
|
if (particle=="Z2SS"){
|
| 865 |
|
|
int processes[] = {217,217,220,223,224,230,234,238,248,249};
|
| 866 |
|
|
subprocesses.assign(processes,processes+10);
|
| 867 |
|
|
}
|
| 868 |
|
|
if (particle=="Z3SS"){
|
| 869 |
|
|
int processes[] = {218,218,221,223,225,231,235,239,250,251};
|
| 870 |
|
|
subprocesses.assign(processes,processes+10);
|
| 871 |
|
|
}
|
| 872 |
|
|
if (particle=="Z4SS"){
|
| 873 |
|
|
int processes[] = {219,219,222,224,225,232,236,240,252,253};
|
| 874 |
|
|
subprocesses.assign(processes,processes+10);
|
| 875 |
|
|
}
|
| 876 |
|
|
if (particle=="W1SS+" || particle=="W1SS-" || particle=="W1SS"){
|
| 877 |
|
|
int processes[] = {226,226,228,229,230,231,232,241,254};
|
| 878 |
|
|
subprocesses.assign(processes,processes+9);
|
| 879 |
|
|
}
|
| 880 |
|
|
if (particle=="W2SS+" || particle=="W2SS-" || particle=="W2SS"){
|
| 881 |
|
|
int processes[] = {227,227,228,233,234,235,236,242,256};
|
| 882 |
|
|
subprocesses.assign(processes,processes+9);
|
| 883 |
|
|
}
|
| 884 |
|
|
if (particle=="UPL" || particle=="UBL" ||
|
| 885 |
|
|
particle=="DNL" || particle=="DBL" ||
|
| 886 |
|
|
particle=="STL" || particle=="SBL" ||
|
| 887 |
|
|
particle=="CHL" || particle=="CBL" || particle=="SQL"){
|
| 888 |
|
|
int processes[] = {246,248,250,252,254,256,258,271,271,273,
|
| 889 |
|
|
274,274,276,277,277,279,279,281,284};
|
| 890 |
|
|
subprocesses.assign(processes,processes+19);
|
| 891 |
|
|
}
|
| 892 |
|
|
if (particle=="BT1" || particle=="BB1" || particle=="B1"){
|
| 893 |
|
|
int processes[] = {281,283,284,286,287,287,289,289,291,291,
|
| 894 |
|
|
293,294,296};
|
| 895 |
|
|
subprocesses.assign(processes,processes+13);
|
| 896 |
|
|
}
|
| 897 |
|
|
if (particle=="BT2" || particle=="BB2" || particle=="B2"){
|
| 898 |
|
|
int processes[] = {282,285,288,288,290,290,292,292,293,295,296};
|
| 899 |
|
|
subprocesses.assign(processes,processes+11);
|
| 900 |
|
|
}
|
| 901 |
|
|
if (particle=="TP1" || particle=="TB1" || particle=="T1"){
|
| 902 |
|
|
int processes[] = {261,261,263,264,264};
|
| 903 |
|
|
subprocesses.assign(processes,processes+5);
|
| 904 |
|
|
}
|
| 905 |
|
|
if (particle=="TP2" || particle=="TB2" || particle=="T2"){
|
| 906 |
|
|
int processes[] = {262,262,263,265,265};
|
| 907 |
|
|
subprocesses.assign(processes,processes+5);
|
| 908 |
|
|
}
|
| 909 |
|
|
if (particle=="UPR" || particle=="UBR" ||
|
| 910 |
|
|
particle=="DNR" || particle=="DBR" ||
|
| 911 |
|
|
particle=="STR" || particle=="SBR" ||
|
| 912 |
|
|
particle=="CHR" || particle=="CBR" || particle=="SQR"){
|
| 913 |
|
|
int processes[] = {247,249,251,253,259,272,272,273,275,275,
|
| 914 |
|
|
276,278,278,280,280,282,283,285,286};
|
| 915 |
|
|
subprocesses.assign(processes,processes+19);
|
| 916 |
|
|
}
|
| 917 |
|
|
if (particle=="HL0"){
|
| 918 |
|
|
int processes[] = {297,299};
|
| 919 |
|
|
subprocesses.assign(processes,processes+2);
|
| 920 |
|
|
}
|
| 921 |
|
|
if (particle=="HH0"){
|
| 922 |
|
|
int processes[] = {298,300};
|
| 923 |
|
|
subprocesses.assign(processes,processes+2);
|
| 924 |
|
|
}
|
| 925 |
|
|
if (particle=="HA0"){
|
| 926 |
|
|
int processes[] = {299,300};
|
| 927 |
|
|
subprocesses.assign(processes,processes+2);
|
| 928 |
|
|
}
|
| 929 |
|
|
if (particle=="H+" || particle=="H-" || particle=="H+-"){
|
| 930 |
|
|
int processes[] = {297,298,301,301};
|
| 931 |
|
|
subprocesses.assign(processes,processes+4);
|
| 932 |
|
|
}
|
| 933 |
|
|
|
| 934 |
|
|
//calculate number of wanted particles per event
|
| 935 |
|
|
double result=0;
|
| 936 |
|
|
double result2=0;
|
| 937 |
|
|
for (vector<int>::iterator iter=subprocesses.begin();
|
| 938 |
|
|
iter != subprocesses.end(); ++ iter){
|
| 939 |
|
|
result+=processMap.find(*iter)->second;
|
| 940 |
|
|
}
|
| 941 |
|
|
result /= processMap.find(0)->second;
|
| 942 |
|
|
|
| 943 |
|
|
if (particle=="EL+" || particle == "EL-" ||
|
| 944 |
|
|
particle=="ER+" || particle == "ER-" ||
|
| 945 |
|
|
particle=="MUL+" || particle == "MUL-" ||
|
| 946 |
|
|
particle=="MUR+" || particle == "MUR-" ||
|
| 947 |
|
|
particle=="TAU1+" || particle == "TAU1-" ||
|
| 948 |
|
|
particle=="TAU2+" || particle == "TAU2-" ||
|
| 949 |
|
|
particle=="W1SS+" || particle=="W1SS-" ||
|
| 950 |
|
|
particle=="W2SS+" || particle=="W2SS-" ||
|
| 951 |
|
|
particle=="H+" || particle=="H-" ||
|
| 952 |
|
|
particle=="BT1" || particle=="BB1" ||
|
| 953 |
|
|
particle=="BT2" || particle=="BB2" ||
|
| 954 |
|
|
particle=="TP1" || particle=="TB1" ||
|
| 955 |
|
|
particle=="TP2" || particle=="TB2"){
|
| 956 |
|
|
result /= 2.;
|
| 957 |
|
|
}
|
| 958 |
|
|
if (particle=="NUEL" || particle == "ANUEL" ||
|
| 959 |
|
|
particle=="NUER" || particle == "ANUER" ||
|
| 960 |
|
|
particle=="NUML" || particle == "ANUML" ||
|
| 961 |
|
|
particle=="NUMR" || particle == "ANUMR"){
|
| 962 |
|
|
result /= 8.;
|
| 963 |
|
|
}
|
| 964 |
|
|
if (particle=="NUTL" || particle == "ANUTL" ||
|
| 965 |
|
|
particle=="NUTR" || particle == "ANUTR"){
|
| 966 |
|
|
result /= 4.;
|
| 967 |
|
|
}
|
| 968 |
|
|
if (particle=="UPL" || particle=="UBL" ||
|
| 969 |
|
|
particle=="DNL" || particle=="DBL" ||
|
| 970 |
|
|
particle=="STL" || particle=="SBL" ||
|
| 971 |
|
|
particle=="CHL" || particle=="CBL" ||
|
| 972 |
|
|
particle=="UPR" || particle=="UBR" ||
|
| 973 |
|
|
particle=="DNR" || particle=="DBR" ||
|
| 974 |
|
|
particle=="STR" || particle=="SBR" ||
|
| 975 |
|
|
particle=="CHR" || particle=="CBR"){
|
| 976 |
|
|
result /= 8.;
|
| 977 |
|
|
}
|
| 978 |
|
|
return result;
|
| 979 |
|
|
}
|
| 980 |
|
|
processTableClass(){}; // constructor
|
| 981 |
|
|
~processTableClass(){}; // destructor
|
| 982 |
|
|
};
|
| 983 |
|
|
|
| 984 |
|
|
////////////////////////////////////////////////////////////////
|
| 985 |
|
|
// small function to remove leading and trailing spaces from
|
| 986 |
|
|
// a string
|
| 987 |
|
|
////////////////////////////////////////////////////////////////
|
| 988 |
|
|
string trim(string str)
|
| 989 |
|
|
{
|
| 990 |
|
|
char const* delims = " ";
|
| 991 |
|
|
int pos;
|
| 992 |
|
|
|
| 993 |
|
|
// trim leading whitespace
|
| 994 |
|
|
string::size_type notwhite = str.find_first_not_of(delims);
|
| 995 |
|
|
str.erase(0,notwhite);
|
| 996 |
|
|
|
| 997 |
|
|
// trim trailing whitespace
|
| 998 |
|
|
notwhite = str.find_last_not_of(delims);
|
| 999 |
|
|
str.erase(notwhite+1);
|
| 1000 |
|
|
|
| 1001 |
|
|
return str;
|
| 1002 |
|
|
}
|
| 1003 |
|
|
////////////////////////////////////////////////////////////////
|
| 1004 |
|
|
|
| 1005 |
|
|
////////////////////////////////////////////////////////////////
|
| 1006 |
|
|
int main(){
|
| 1007 |
|
|
|
| 1008 |
|
|
//read in partial cross sections
|
| 1009 |
|
|
//define input file (output of Pythia)
|
| 1010 |
|
|
ifstream fin2("./pythia.out");
|
| 1011 |
|
|
string line2;
|
| 1012 |
|
|
processTableClass processTable;
|
| 1013 |
|
|
while( getline(fin2,line2) ) {
|
| 1014 |
|
|
string::size_type loc = line2.find( "All included subprocesses", 0 );
|
| 1015 |
|
|
if( loc != string::npos ) {
|
| 1016 |
|
|
string tmpPid = trim(line2.substr(3,3));
|
| 1017 |
|
|
string tmpWq = trim(line2.substr(68,9));
|
| 1018 |
|
|
int convtmpPid;
|
| 1019 |
|
|
double convtmpWq;
|
| 1020 |
|
|
istringstream ins; // Declare an input string stream.
|
| 1021 |
|
|
ins.str(tmpPid); // Specify string to read.
|
| 1022 |
|
|
ins >> convtmpPid; // Reads the integer from the string.
|
| 1023 |
|
|
istringstream ins2; // Declare an input string stream.
|
| 1024 |
|
|
ins2.str(tmpWq); // Specify string to read.
|
| 1025 |
|
|
ins2 >> convtmpWq; // Reads the double from the string.
|
| 1026 |
|
|
processTable.addProcess(convtmpPid,convtmpWq);
|
| 1027 |
|
|
}
|
| 1028 |
|
|
loc = line2.find( "->", 0 );
|
| 1029 |
|
|
if( loc != string::npos ) {
|
| 1030 |
|
|
if (line2.at(1)=='I'){
|
| 1031 |
|
|
string tmpPid = trim(line2.substr(3,3));
|
| 1032 |
|
|
string tmpWq = trim(line2.substr(68,9));
|
| 1033 |
|
|
int convtmpPid;
|
| 1034 |
|
|
double convtmpWq;
|
| 1035 |
|
|
istringstream ins; // Declare an input string stream.
|
| 1036 |
|
|
ins.str(tmpPid); // Specify string to read.
|
| 1037 |
|
|
ins >> convtmpPid; // Reads the integer from the string.
|
| 1038 |
|
|
istringstream ins2; // Declare an input string stream.
|
| 1039 |
|
|
ins2.str(tmpWq); // Specify string to read.
|
| 1040 |
|
|
ins2 >> convtmpWq; // Reads the double from the string.
|
| 1041 |
|
|
processTable.addProcess(convtmpPid,convtmpWq);
|
| 1042 |
|
|
}
|
| 1043 |
|
|
}
|
| 1044 |
|
|
}
|
| 1045 |
|
|
//processTable.printProcessTable();
|
| 1046 |
|
|
|
| 1047 |
|
|
//read in decay table
|
| 1048 |
|
|
//define input file (output of IsaJet)
|
| 1049 |
|
|
//ifstream fin("./LM4.out");
|
| 1050 |
|
|
ifstream fin("./isa.out");
|
| 1051 |
|
|
string line;
|
| 1052 |
|
|
decayTableClass decayTable;
|
| 1053 |
|
|
while( getline(fin,line) ) {
|
| 1054 |
|
|
string::size_type loc = line.find( " --> ", 0 );
|
| 1055 |
|
|
if( loc != string::npos ) {
|
| 1056 |
|
|
string tmpMother = trim(line.substr(0,8));
|
| 1057 |
|
|
string tmpBR = trim(line.substr(60,13));
|
| 1058 |
|
|
double convtmpBR;
|
| 1059 |
|
|
istringstream ins; // Declare an input string stream.
|
| 1060 |
|
|
ins.str(tmpBR); // Specify string to read.
|
| 1061 |
|
|
ins >> convtmpBR; // Reads the double from the string.
|
| 1062 |
|
|
if (tmpMother != "TP" && tmpMother != "TB"){
|
| 1063 |
|
|
decayClass decay(tmpMother,convtmpBR);
|
| 1064 |
|
|
string tmpDaughter;
|
| 1065 |
|
|
for (int i=0; i<4; ++i){
|
| 1066 |
|
|
tmpDaughter= trim(line.substr(12+i*7,7));
|
| 1067 |
|
|
if (tmpDaughter.size()>0){
|
| 1068 |
|
|
decay.addDaughter(tmpDaughter);
|
| 1069 |
|
|
}
|
| 1070 |
|
|
}
|
| 1071 |
|
|
decayTable.addDecay(decay);
|
| 1072 |
|
|
}
|
| 1073 |
|
|
}
|
| 1074 |
|
|
}
|
| 1075 |
|
|
|
| 1076 |
|
|
// Add top, Z and W BrRatios by hand
|
| 1077 |
|
|
string tmpMother = "TP";
|
| 1078 |
|
|
double tmpBR = 1.0;
|
| 1079 |
|
|
decayClass decay0(tmpMother,tmpBR);
|
| 1080 |
|
|
decay0.addDaughter("W+");
|
| 1081 |
|
|
decay0.addDaughter("BT");
|
| 1082 |
|
|
decayTable.addDecay(decay0);
|
| 1083 |
|
|
|
| 1084 |
|
|
tmpMother = "Z0";
|
| 1085 |
|
|
tmpBR = 0.03658;
|
| 1086 |
|
|
decayClass decay1(tmpMother,tmpBR);
|
| 1087 |
|
|
decay1.addDaughter("E+");
|
| 1088 |
|
|
decay1.addDaughter("E-");
|
| 1089 |
|
|
decayTable.addDecay(decay1);
|
| 1090 |
|
|
|
| 1091 |
|
|
tmpMother = "Z0";
|
| 1092 |
|
|
tmpBR = 0.06342;
|
| 1093 |
|
|
decayClass decay1a(tmpMother,tmpBR);
|
| 1094 |
|
|
decay1a.addDaughter("NUE");
|
| 1095 |
|
|
decay1a.addDaughter("NUE");
|
| 1096 |
|
|
decayTable.addDecay(decay1a);
|
| 1097 |
|
|
|
| 1098 |
|
|
tmpMother = "Z0";
|
| 1099 |
|
|
tmpBR = 0.03658;
|
| 1100 |
|
|
decayClass decay2(tmpMother,tmpBR);
|
| 1101 |
|
|
decay2.addDaughter("MU+");
|
| 1102 |
|
|
decay2.addDaughter("MU-");
|
| 1103 |
|
|
decayTable.addDecay(decay2);
|
| 1104 |
|
|
|
| 1105 |
|
|
tmpMother = "Z0";
|
| 1106 |
|
|
tmpBR = 0.06342;
|
| 1107 |
|
|
decayClass decay1b(tmpMother,tmpBR);
|
| 1108 |
|
|
decay1b.addDaughter("NUM");
|
| 1109 |
|
|
decay1b.addDaughter("NUM");
|
| 1110 |
|
|
decayTable.addDecay(decay1b);
|
| 1111 |
|
|
|
| 1112 |
|
|
tmpMother = "Z0";
|
| 1113 |
|
|
tmpBR = 0.03658;
|
| 1114 |
|
|
decayClass decay3(tmpMother,tmpBR);
|
| 1115 |
|
|
decay3.addDaughter("TAU+");
|
| 1116 |
|
|
decay3.addDaughter("TAU-");
|
| 1117 |
|
|
decayTable.addDecay(decay3);
|
| 1118 |
|
|
|
| 1119 |
|
|
tmpMother = "Z0";
|
| 1120 |
|
|
tmpBR = 0.06342;
|
| 1121 |
|
|
decayClass decay1c(tmpMother,tmpBR);
|
| 1122 |
|
|
decay1c.addDaughter("NUT");
|
| 1123 |
|
|
decay1c.addDaughter("NUT");
|
| 1124 |
|
|
decayTable.addDecay(decay1c);
|
| 1125 |
|
|
|
| 1126 |
|
|
tmpMother = "Z0";
|
| 1127 |
|
|
tmpBR = 0.101;
|
| 1128 |
|
|
decayClass decay4(tmpMother,tmpBR);
|
| 1129 |
|
|
decay4.addDaughter("UP");
|
| 1130 |
|
|
decay4.addDaughter("UB");
|
| 1131 |
|
|
decayTable.addDecay(decay4);
|
| 1132 |
|
|
|
| 1133 |
|
|
tmpMother = "Z0";
|
| 1134 |
|
|
tmpBR = 0.101;
|
| 1135 |
|
|
decayClass decay5(tmpMother,tmpBR);
|
| 1136 |
|
|
decay5.addDaughter("CH");
|
| 1137 |
|
|
decay5.addDaughter("CB");
|
| 1138 |
|
|
decayTable.addDecay(decay5);
|
| 1139 |
|
|
|
| 1140 |
|
|
tmpMother = "Z0";
|
| 1141 |
|
|
tmpBR = 0.166;
|
| 1142 |
|
|
decayClass decay6(tmpMother,tmpBR);
|
| 1143 |
|
|
decay6.addDaughter("DN");
|
| 1144 |
|
|
decay6.addDaughter("DB");
|
| 1145 |
|
|
decayTable.addDecay(decay6);
|
| 1146 |
|
|
|
| 1147 |
|
|
tmpMother = "Z0";
|
| 1148 |
|
|
tmpBR = 0.166;
|
| 1149 |
|
|
decayClass decay7(tmpMother,tmpBR);
|
| 1150 |
|
|
decay7.addDaughter("ST");
|
| 1151 |
|
|
decay7.addDaughter("SB");
|
| 1152 |
|
|
decayTable.addDecay(decay7);
|
| 1153 |
|
|
|
| 1154 |
|
|
tmpMother = "Z0";
|
| 1155 |
|
|
tmpBR = 0.166;
|
| 1156 |
|
|
decayClass decay8(tmpMother,tmpBR);
|
| 1157 |
|
|
decay8.addDaughter("BT");
|
| 1158 |
|
|
decay8.addDaughter("BB");
|
| 1159 |
|
|
decayTable.addDecay(decay8);
|
| 1160 |
|
|
|
| 1161 |
|
|
tmpMother = "W+";
|
| 1162 |
|
|
tmpBR = 0.1068;
|
| 1163 |
|
|
decayClass decay9(tmpMother,tmpBR);
|
| 1164 |
|
|
decay9.addDaughter("E+");
|
| 1165 |
|
|
decay9.addDaughter("NUE");
|
| 1166 |
|
|
decayTable.addDecay(decay9);
|
| 1167 |
|
|
|
| 1168 |
|
|
tmpMother = "W+";
|
| 1169 |
|
|
tmpBR = 0.1068;
|
| 1170 |
|
|
decayClass decay10(tmpMother,tmpBR);
|
| 1171 |
|
|
decay10.addDaughter("MU+");
|
| 1172 |
|
|
decay10.addDaughter("NUM");
|
| 1173 |
|
|
decayTable.addDecay(decay10);
|
| 1174 |
|
|
|
| 1175 |
|
|
tmpMother = "W+";
|
| 1176 |
|
|
tmpBR = 0.1068;
|
| 1177 |
|
|
decayClass decay11(tmpMother,tmpBR);
|
| 1178 |
|
|
decay11.addDaughter("TAU+");
|
| 1179 |
|
|
decay11.addDaughter("NUT");
|
| 1180 |
|
|
decayTable.addDecay(decay11);
|
| 1181 |
|
|
|
| 1182 |
|
|
tmpMother = "W+";
|
| 1183 |
|
|
tmpBR = 0.3398;
|
| 1184 |
|
|
decayClass decay12(tmpMother,tmpBR);
|
| 1185 |
|
|
decay12.addDaughter("UP");
|
| 1186 |
|
|
decay12.addDaughter("DB");
|
| 1187 |
|
|
decayTable.addDecay(decay12);
|
| 1188 |
|
|
|
| 1189 |
|
|
tmpMother = "W+";
|
| 1190 |
|
|
tmpBR = 0.3398;
|
| 1191 |
|
|
decayClass decay13(tmpMother,tmpBR);
|
| 1192 |
|
|
decay13.addDaughter("CH");
|
| 1193 |
|
|
decay13.addDaughter("CB");
|
| 1194 |
|
|
decayTable.addDecay(decay13);
|
| 1195 |
|
|
|
| 1196 |
|
|
//add all the production processes to the decay table
|
| 1197 |
|
|
//read in subprocess table
|
| 1198 |
|
|
ifstream fin3("./processTable.dat");
|
| 1199 |
|
|
string line3;
|
| 1200 |
|
|
double totalBR=0.;
|
| 1201 |
|
|
while( getline(fin3,line3) ) {
|
| 1202 |
|
|
string::size_type loc = line3.find( "PP", 0 );
|
| 1203 |
|
|
if( loc != string::npos ) {
|
| 1204 |
|
|
string tmpMotherPP = trim(line3.substr(0,9));
|
| 1205 |
|
|
string tmpDaughterPP1 = trim(line3.substr(10,9));
|
| 1206 |
|
|
string tmpDaughterPP2 = trim(line3.substr(20,9));
|
| 1207 |
|
|
string tmpPythiaID = trim(line3.substr(30,9));
|
| 1208 |
|
|
string tmpWeight = trim(line3.substr(40,9));
|
| 1209 |
|
|
int convID;
|
| 1210 |
|
|
istringstream ins; // Declare an input string stream.
|
| 1211 |
|
|
ins.str(tmpPythiaID); // Specify string to read.
|
| 1212 |
|
|
ins >> convID; // Reads the double from the string.
|
| 1213 |
|
|
double convWeight;
|
| 1214 |
|
|
istringstream ins2; // Declare an input string stream.
|
| 1215 |
|
|
ins2.str(tmpWeight); // Specify string to read.
|
| 1216 |
|
|
ins2 >> convWeight; // Reads the double from the string.
|
| 1217 |
|
|
double tmpBRPP = processTable.getProcessWeight(convID)*convWeight;
|
| 1218 |
|
|
totalBR+=tmpBRPP;
|
| 1219 |
|
|
decayClass decayPP(tmpMotherPP,tmpBRPP);
|
| 1220 |
|
|
decayPP.addDaughter(tmpDaughterPP1);
|
| 1221 |
|
|
decayPP.addDaughter(tmpDaughterPP2);
|
| 1222 |
|
|
decayTable.addDecay(decayPP);
|
| 1223 |
|
|
}
|
| 1224 |
|
|
}
|
| 1225 |
|
|
|
| 1226 |
|
|
cout<<"Total pp branching ratio (should be exactly 1 !!!) = "<<totalBR<<endl;
|
| 1227 |
|
|
|
| 1228 |
|
|
//decayTable.printDecayTable();
|
| 1229 |
|
|
|
| 1230 |
|
|
//----------------------------------------------------------------------------
|
| 1231 |
|
|
//calculate how often a particular decay chain can be found in a pp event
|
| 1232 |
|
|
double tot_prod = 0.;
|
| 1233 |
|
|
double EL_prod = processTable.getParticlesPerEvent("EL");
|
| 1234 |
|
|
tot_prod += EL_prod;
|
| 1235 |
|
|
double ER_prod = processTable.getParticlesPerEvent("ER");
|
| 1236 |
|
|
tot_prod += ER_prod;
|
| 1237 |
|
|
double MUL_prod = processTable.getParticlesPerEvent("MUL");
|
| 1238 |
|
|
tot_prod += MUL_prod;
|
| 1239 |
|
|
double MUR_prod = processTable.getParticlesPerEvent("MUR");
|
| 1240 |
|
|
tot_prod += MUR_prod;
|
| 1241 |
|
|
double TAU1_prod = processTable.getParticlesPerEvent("TAU1");
|
| 1242 |
|
|
tot_prod += TAU1_prod;
|
| 1243 |
|
|
double TAU2_prod = processTable.getParticlesPerEvent("TAU2");
|
| 1244 |
|
|
tot_prod += TAU2_prod;
|
| 1245 |
|
|
double NULEP_prod = processTable.getParticlesPerEvent("NULEP");
|
| 1246 |
|
|
tot_prod += NULEP_prod;
|
| 1247 |
|
|
double NUTAU_prod = processTable.getParticlesPerEvent("NUTAU");
|
| 1248 |
|
|
tot_prod += NUTAU_prod;
|
| 1249 |
|
|
double SQL_prod = processTable.getParticlesPerEvent("SQL");
|
| 1250 |
|
|
tot_prod += SQL_prod;
|
| 1251 |
|
|
double SQR_prod = processTable.getParticlesPerEvent("SQR");
|
| 1252 |
|
|
tot_prod += SQR_prod;
|
| 1253 |
|
|
double BT1_prod = processTable.getParticlesPerEvent("B1");
|
| 1254 |
|
|
tot_prod += BT1_prod;
|
| 1255 |
|
|
double BT2_prod = processTable.getParticlesPerEvent("B2");
|
| 1256 |
|
|
tot_prod += BT2_prod;
|
| 1257 |
|
|
double TP1_prod = processTable.getParticlesPerEvent("T1");
|
| 1258 |
|
|
tot_prod += TP1_prod;
|
| 1259 |
|
|
double TP2_prod = processTable.getParticlesPerEvent("T2");
|
| 1260 |
|
|
tot_prod += TP2_prod;
|
| 1261 |
|
|
double GLSS_prod = processTable.getParticlesPerEvent("GLSS");
|
| 1262 |
|
|
tot_prod += GLSS_prod;
|
| 1263 |
|
|
double Z1SS_prod = processTable.getParticlesPerEvent("Z1SS");
|
| 1264 |
|
|
tot_prod += Z1SS_prod;
|
| 1265 |
|
|
double Z2SS_prod = processTable.getParticlesPerEvent("Z2SS");
|
| 1266 |
|
|
tot_prod += Z2SS_prod;
|
| 1267 |
|
|
double Z3SS_prod = processTable.getParticlesPerEvent("Z3SS");
|
| 1268 |
|
|
tot_prod += Z3SS_prod;
|
| 1269 |
|
|
double Z4SS_prod = processTable.getParticlesPerEvent("Z4SS");
|
| 1270 |
|
|
tot_prod += Z4SS_prod;
|
| 1271 |
|
|
double W1SS_prod = processTable.getParticlesPerEvent("W1SS");
|
| 1272 |
|
|
tot_prod += W1SS_prod;
|
| 1273 |
|
|
double W2SS_prod = processTable.getParticlesPerEvent("W2SS");
|
| 1274 |
|
|
tot_prod += W2SS_prod;
|
| 1275 |
|
|
double HL0_prod = processTable.getParticlesPerEvent("HL0");
|
| 1276 |
|
|
tot_prod += HL0_prod;
|
| 1277 |
|
|
double HH0_prod = processTable.getParticlesPerEvent("HH0");
|
| 1278 |
|
|
tot_prod += HH0_prod;
|
| 1279 |
|
|
double HA0_prod = processTable.getParticlesPerEvent("HA0");
|
| 1280 |
|
|
tot_prod += HA0_prod;
|
| 1281 |
|
|
double Hpn_prod = processTable.getParticlesPerEvent("H+-");
|
| 1282 |
|
|
tot_prod += Hpn_prod;
|
| 1283 |
|
|
|
| 1284 |
|
|
/*
|
| 1285 |
|
|
|
| 1286 |
|
|
cout<<"EL per event = "<<EL_prod<<endl;
|
| 1287 |
|
|
cout<<"ER per event = "<<ER_prod<<endl;
|
| 1288 |
|
|
cout<<"MUL per event = "<<MUL_prod<<endl;
|
| 1289 |
|
|
cout<<"MUR per event = "<<MUR_prod<<endl;
|
| 1290 |
|
|
cout<<"TAU1 per event = "<<TAU1_prod<<endl;
|
| 1291 |
|
|
cout<<"TAU2 per event = "<<TAU2_prod<<endl;
|
| 1292 |
|
|
cout<<"NULEP per event = "<<NULEP_prod<<endl;
|
| 1293 |
|
|
cout<<"NUTAU per event = "<<NUTAU_prod<<endl;
|
| 1294 |
|
|
|
| 1295 |
|
|
cout<<"SQL per event = "<<SQL_prod<<endl;
|
| 1296 |
|
|
cout<<"SQR per event = "<<SQR_prod<<endl;
|
| 1297 |
|
|
cout<<"BT1 per event = "<<BT1_prod<<endl;
|
| 1298 |
|
|
cout<<"BT2 per event = "<<BT2_prod<<endl;
|
| 1299 |
|
|
cout<<"TP1 per event = "<<TP1_prod<<endl;
|
| 1300 |
|
|
cout<<"TP2 per event = "<<TP2_prod<<endl;
|
| 1301 |
|
|
|
| 1302 |
|
|
cout<<"GLSS per event = "<<GLSS_prod<<endl;
|
| 1303 |
|
|
cout<<"Z1SS per event = "<<Z1SS_prod<<endl;
|
| 1304 |
|
|
cout<<"Z2SS per event = "<<Z2SS_prod<<endl;
|
| 1305 |
|
|
cout<<"Z3SS per event = "<<Z3SS_prod<<endl;
|
| 1306 |
|
|
cout<<"Z4SS per event = "<<Z4SS_prod<<endl;
|
| 1307 |
|
|
cout<<"W1SS per event = "<<W1SS_prod<<endl;
|
| 1308 |
|
|
cout<<"W2SS per event = "<<W2SS_prod<<endl;
|
| 1309 |
|
|
|
| 1310 |
|
|
cout<<"HL0 per event = "<<HL0_prod<<endl;
|
| 1311 |
|
|
cout<<"HH0 per event = "<<HH0_prod<<endl;
|
| 1312 |
|
|
cout<<"HA0 per event = "<<HA0_prod<<endl;
|
| 1313 |
|
|
cout<<"H+- per event = "<<Hpn_prod<<endl;
|
| 1314 |
|
|
|
| 1315 |
|
|
*/
|
| 1316 |
|
|
|
| 1317 |
|
|
cout<<"total per event (should be exactly 2 !!!) = "<<tot_prod<<endl;
|
| 1318 |
|
|
|
| 1319 |
|
|
int i, j, k, l, m, n, o, p, q;
|
| 1320 |
|
|
double final=0., finalWAv=0., finalW=0., finalTopAv=0., finalTop=0.;
|
| 1321 |
|
|
|
| 1322 |
|
|
/*
|
| 1323 |
|
|
|
| 1324 |
|
|
i=-1; //Leptons e,mu
|
| 1325 |
|
|
j=-1; //Quarks u,d,c,s,b
|
| 1326 |
|
|
k=-1; //Gluons
|
| 1327 |
|
|
l=-1; //taus
|
| 1328 |
|
|
m=-1; //Ws
|
| 1329 |
|
|
n=-1; //Zs
|
| 1330 |
|
|
o=1; //W,Z,h
|
| 1331 |
|
|
p=-1; //jets (incl taus)
|
| 1332 |
|
|
q=0; //tops
|
| 1333 |
|
|
|
| 1334 |
|
|
double final=0;
|
| 1335 |
|
|
final += decayTable.BrRatioByDecayProducts("EL+", i, j, q, k, l, m, n, o, p)*EL_prod;
|
| 1336 |
|
|
final += decayTable.BrRatioByDecayProducts("ER+", i, j, q, k, l, m, n, o, p)*ER_prod;
|
| 1337 |
|
|
final += decayTable.BrRatioByDecayProducts("MUL+", i, j, q, k, l, m, n, o, p)*MUL_prod;
|
| 1338 |
|
|
final += decayTable.BrRatioByDecayProducts("MUR+", i, j, q, k, l, m, n, o, p)*MUR_prod;
|
| 1339 |
|
|
final += decayTable.BrRatioByDecayProducts("TAU1+", i, j, q, k, l, m, n, o, p)*TAU1_prod;
|
| 1340 |
|
|
final += decayTable.BrRatioByDecayProducts("TAU2+", i, j, q, k, l, m, n, o, p)*TAU2_prod;
|
| 1341 |
|
|
final += decayTable.BrRatioByDecayProducts("NUEL", i, j, q, k, l, m, n, o, p)*NULEP_prod;
|
| 1342 |
|
|
final += decayTable.BrRatioByDecayProducts("NUTL", i, j, q, k, l, m, n, o, p)*NUTAU_prod;
|
| 1343 |
|
|
final += decayTable.BrRatioByDecayProducts("GLSS", i, j, q, k, l, m, n, o, p)*GLSS_prod;
|
| 1344 |
|
|
final += decayTable.BrRatioByDecayProducts("Z1SS", i, j, q, k, l, m, n, o, p)*Z1SS_prod;
|
| 1345 |
|
|
final += decayTable.BrRatioByDecayProducts("Z2SS", i, j, q, k, l, m, n, o, p)*Z2SS_prod;
|
| 1346 |
|
|
final += decayTable.BrRatioByDecayProducts("Z3SS", i, j, q, k, l, m, n, o, p)*Z3SS_prod;
|
| 1347 |
|
|
final += decayTable.BrRatioByDecayProducts("Z4SS", i, j, q, k, l, m, n, o, p)*Z4SS_prod;
|
| 1348 |
|
|
final += decayTable.BrRatioByDecayProducts("W1SS+", i, j, q, k, l, m, n, o, p)*W1SS_prod;
|
| 1349 |
|
|
final += decayTable.BrRatioByDecayProducts("W2SS+", i, j, q, k, l, m, n, o, p)*W2SS_prod;
|
| 1350 |
|
|
final += decayTable.BrRatioByDecayProducts("UPL", i, j, q, k, l, m, n, o, p)*SQL_prod/4.;
|
| 1351 |
|
|
final += decayTable.BrRatioByDecayProducts("UPR", i, j, q, k, l, m, n, o, p)*SQR_prod/4.;
|
| 1352 |
|
|
final += decayTable.BrRatioByDecayProducts("DNL", i, j, q, k, l, m, n, o, p)*SQL_prod/4.;
|
| 1353 |
|
|
final += decayTable.BrRatioByDecayProducts("DNR", i, j, q, k, l, m, n, o, p)*SQR_prod/4.;
|
| 1354 |
|
|
final += decayTable.BrRatioByDecayProducts("CHL", i, j, q, k, l, m, n, o, p)*SQL_prod/4.;
|
| 1355 |
|
|
final += decayTable.BrRatioByDecayProducts("CHR", i, j, q, k, l, m, n, o, p)*SQR_prod/4.;
|
| 1356 |
|
|
final += decayTable.BrRatioByDecayProducts("STL", i, j, q, k, l, m, n, o, p)*SQL_prod/4.;
|
| 1357 |
|
|
final += decayTable.BrRatioByDecayProducts("STR", i, j, q, k, l, m, n, o, p)*SQR_prod/4.;
|
| 1358 |
|
|
final += decayTable.BrRatioByDecayProducts("BT1", i, j, q, k, l, m, n, o, p)*BT1_prod;
|
| 1359 |
|
|
final += decayTable.BrRatioByDecayProducts("BT2", i, j, q, k, l, m, n, o, p)*BT2_prod;
|
| 1360 |
|
|
final += decayTable.BrRatioByDecayProducts("TP1", i, j, q, k, l, m, n, o, p)*TP1_prod;
|
| 1361 |
|
|
final += decayTable.BrRatioByDecayProducts("TP2", i, j, q, k, l, m, n, o, p)*TP2_prod;
|
| 1362 |
|
|
final += decayTable.BrRatioByDecayProducts("HL0", i, j, q, k, l, m, n, o, p)*HL0_prod;
|
| 1363 |
|
|
final += decayTable.BrRatioByDecayProducts("HH0", i, j, q, k, l, m, n, o, p)*HH0_prod;
|
| 1364 |
|
|
final += decayTable.BrRatioByDecayProducts("HA0", i, j, q, k, l, m, n, o, p)*HA0_prod;
|
| 1365 |
|
|
final += decayTable.BrRatioByDecayProducts("H+", i, j, q, k, l, m, n, o, p)*Hpn_prod;
|
| 1366 |
|
|
|
| 1367 |
|
|
cout<<"how often appears final state in one cascade of pp: leptons:"<<i
|
| 1368 |
|
|
<<" quarks:"<<j
|
| 1369 |
|
|
<<" tops: "<<q
|
| 1370 |
|
|
<<" gluons:"<<k
|
| 1371 |
|
|
<<" taus:"<<l
|
| 1372 |
|
|
<<" Ws: "<<m
|
| 1373 |
|
|
<<" Zs: "<<n
|
| 1374 |
|
|
<<" Bosons: "<<o
|
| 1375 |
|
|
<<" Jets: "<<p
|
| 1376 |
|
|
<<" Nevents: "<<final<<endl;
|
| 1377 |
|
|
//----------------------------------------------------------------------------
|
| 1378 |
|
|
*/
|
| 1379 |
|
|
|
| 1380 |
|
|
//----------------------------------------------------------------------------
|
| 1381 |
|
|
//calculate how often a particular final state can be found in a pp event
|
| 1382 |
|
|
// for (int oo=0; oo<5; ++oo){
|
| 1383 |
|
|
// i=-1;
|
| 1384 |
|
|
// j=-1;
|
| 1385 |
|
|
// k=-1;
|
| 1386 |
|
|
// l=-1;
|
| 1387 |
|
|
// m=-1;
|
| 1388 |
|
|
// n=-1;
|
| 1389 |
|
|
// o=oo;
|
| 1390 |
|
|
// p=-1;
|
| 1391 |
|
|
// q=-1;
|
| 1392 |
|
|
// final = decayTable.BrRatioByDecayProducts("PP", i, j, q, k, l, m, n, o, p);
|
| 1393 |
|
|
// cout<<"pp --> leptons:"<<i
|
| 1394 |
|
|
// <<" quarks:"<<j
|
| 1395 |
|
|
// <<" tops: "<<q
|
| 1396 |
|
|
// <<" gluons:"<<k
|
| 1397 |
|
|
// <<" taus:"<<l
|
| 1398 |
|
|
// <<" Ws: "<<m
|
| 1399 |
|
|
// <<" Zs: "<<n
|
| 1400 |
|
|
// <<" Bosons: "<<o
|
| 1401 |
|
|
// <<" Jets: "<<p
|
| 1402 |
|
|
// <<" Nevents: "<<final<<endl;
|
| 1403 |
|
|
// finalWAv += (final*o);
|
| 1404 |
|
|
// finalW += final;
|
| 1405 |
|
|
// }
|
| 1406 |
|
|
// cout<<"Bosons per event: "<<finalWAv<<endl;
|
| 1407 |
|
|
// cout<<"Events with bosons: "<<finalW<<endl;
|
| 1408 |
|
|
|
| 1409 |
|
|
for (int qq=0; qq<5; ++qq){
|
| 1410 |
|
|
i=-1;
|
| 1411 |
|
|
j=-1;
|
| 1412 |
|
|
k=-1;
|
| 1413 |
|
|
l=-1;
|
| 1414 |
|
|
m=-1;
|
| 1415 |
|
|
n=-1;
|
| 1416 |
|
|
o=-1;
|
| 1417 |
|
|
p=-1;
|
| 1418 |
|
|
q=qq;
|
| 1419 |
|
|
final = decayTable.BrRatioByDecayProducts("PP", i, j, q, k, l, m, n, o, p);
|
| 1420 |
|
|
cout<<"pp --> leptons:"<<i
|
| 1421 |
|
|
<<" quarks:"<<j
|
| 1422 |
|
|
<<" tops: "<<q
|
| 1423 |
|
|
<<" gluons:"<<k
|
| 1424 |
|
|
<<" taus:"<<l
|
| 1425 |
|
|
<<" Ws: "<<m
|
| 1426 |
|
|
<<" Zs: "<<n
|
| 1427 |
|
|
<<" Bosons: "<<o
|
| 1428 |
|
|
<<" Jets: "<<p
|
| 1429 |
|
|
<<" Nevents: "<<final<<endl;
|
| 1430 |
|
|
finalTopAv += (final*qq);
|
| 1431 |
|
|
finalTop += final;
|
| 1432 |
|
|
}
|
| 1433 |
|
|
cout<<"Tops per event: "<<finalTopAv<<endl;
|
| 1434 |
|
|
cout<<"Events with topst: "<<finalTop<<endl;
|
| 1435 |
|
|
|
| 1436 |
|
|
//decayTable.printFinalStateTable();
|
| 1437 |
|
|
|
| 1438 |
|
|
return 0;
|
| 1439 |
|
|
}
|
