Line data Source code
1 : //--------------------------------------------------------------------------
2 : //
3 : // Environment:
4 : // This software is part of the EvtGen package developed jointly
5 : // for the BaBar and CLEO collaborations. If you use all or part
6 : // of it, please give an appropriate acknowledgement.
7 : //
8 : // Copyright Information: See EvtGen/COPYRIGHT
9 : // Copyright (C) 1998 Caltech, UCSB
10 : //
11 : // Module: EvtGen/EvtGenericDalitz.hh
12 : //
13 : // Description: Model to describe a generic dalitz decay
14 : //
15 : // Modification history:
16 : //
17 : // DCC 16 December, 2011 Module created
18 : //
19 : //------------------------------------------------------------------------
20 :
21 : #include "EvtGenModels/EvtDalitzTable.hh"
22 : #include "EvtGenBase/EvtReport.hh"
23 : #include "EvtGenBase/EvtParserXml.hh"
24 : #include "EvtGenBase/EvtPDL.hh"
25 : #include "EvtGenBase/EvtSpinType.hh"
26 : #include "EvtGenBase/EvtDalitzPlot.hh"
27 : #include "EvtGenBase/EvtCyclic3.hh"
28 :
29 : #include <stdlib.h>
30 : #include <sstream>
31 :
32 : using std::endl;
33 : using std::fstream;
34 : using std::ifstream;
35 :
36 0 : EvtDalitzTable::EvtDalitzTable() {
37 0 : _dalitztable.clear();
38 0 : _readFiles.clear();
39 0 : }
40 :
41 0 : EvtDalitzTable::~EvtDalitzTable() {
42 0 : _dalitztable.clear();
43 0 : _readFiles.clear();
44 0 : }
45 :
46 : EvtDalitzTable* EvtDalitzTable::getInstance(const std::string dec_name, bool verbose) {
47 :
48 : static EvtDalitzTable* theDalitzTable = 0;
49 :
50 0 : if(theDalitzTable == 0) {
51 0 : theDalitzTable = new EvtDalitzTable();
52 0 : }
53 :
54 0 : if(!theDalitzTable->fileHasBeenRead(dec_name)) {
55 0 : theDalitzTable->readXMLDecayFile(dec_name,verbose);
56 0 : }
57 :
58 0 : return theDalitzTable;
59 :
60 0 : }
61 :
62 : bool EvtDalitzTable::fileHasBeenRead(const std::string dec_name) {
63 0 : std::vector<std::string>::iterator i = _readFiles.begin();
64 0 : for( ; i!=_readFiles.end(); i++) {
65 0 : if((*i).compare(dec_name) == 0) {
66 0 : return true;
67 : }
68 : }
69 0 : return false;
70 0 : }
71 :
72 : void EvtDalitzTable::readXMLDecayFile(const std::string dec_name, bool verbose){
73 :
74 0 : if (verbose) {
75 0 : report(Severity::Info,"EvtGen")<<"EvtDalitzTable: Reading in xml parameter file "<<dec_name<<endl;
76 0 : }
77 :
78 0 : _readFiles.push_back(dec_name);
79 :
80 : EvtDalitzDecayInfo* dalitzDecay = 0;
81 : double probMax = 0;
82 0 : EvtId ipar;
83 0 : std::string decayParent = "";
84 0 : std::string daugStr = "";
85 0 : EvtId daughter[3];
86 :
87 0 : EvtDalitzPlot dp;
88 0 : EvtComplex cAmp;
89 0 : std::vector< std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> > angAndResPairs;
90 0 : std::string shape("");
91 : EvtSpinType::spintype spinType(EvtSpinType::SCALAR);
92 : double mass(0.), width(0.), FFp(0.), FFr(0.);
93 0 : std::vector<EvtFlatteParam> flatteParams;
94 : //Nonres parameters
95 : double alpha(0.);
96 : //LASS parameters
97 : double aLass(0.), rLass(0.), BLass(0.), phiBLass(0.), RLass(0.), phiRLass(0.), cutoffLass(-1.);
98 :
99 0 : EvtParserXml parser;
100 0 : parser.open(dec_name);
101 :
102 : bool endReached = false;
103 :
104 0 : while(parser.readNextTag()) {
105 : //TAGS FOUND UNDER DATA
106 0 : if(parser.getParentTagTitle() == "data") {
107 0 : if(parser.getTagTitle() == "dalitzDecay") {
108 : int nDaughters = 0;
109 :
110 0 : decayParent = parser.readAttribute("particle");
111 0 : daugStr = parser.readAttribute("daughters");
112 0 : probMax = parser.readAttributeDouble("probMax",-1);
113 :
114 0 : checkParticle(decayParent);
115 0 : ipar=EvtPDL::getId(decayParent);
116 :
117 0 : std::istringstream daugStream(daugStr);
118 :
119 0 : std::string daugh;
120 0 : while(std::getline(daugStream, daugh, ' ')) {
121 0 : checkParticle(daugh);
122 0 : daughter[nDaughters++] = EvtPDL::getId(daugh);
123 : }
124 :
125 0 : if(nDaughters!=3) {
126 0 : report(Severity::Error,"EvtGen") <<
127 0 : "Expected to find three daughters for dalitzDecay of "<<decayParent<<" near line "
128 0 : <<parser.getLineNumber()<<", "<<"found "<<nDaughters<<endl;
129 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
130 0 : ::abort();
131 : }
132 :
133 0 : double m_d1 = EvtPDL::getMass(daughter[0]), m_d2 = EvtPDL::getMass(daughter[1]), m_d3 = EvtPDL::getMass(daughter[2]), M = EvtPDL::getMass(ipar);
134 :
135 0 : dp = EvtDalitzPlot( m_d1, m_d2, m_d3, M );
136 :
137 0 : dalitzDecay = new EvtDalitzDecayInfo(daughter[0],daughter[1],daughter[2]);
138 :
139 0 : } else if(parser.getTagTitle() == "copyDalitz") {
140 : int nDaughters = 0;
141 0 : EvtId daughter[3];
142 : int nCopyDaughters = 0;
143 0 : EvtId copyDaughter[3];
144 :
145 0 : decayParent = parser.readAttribute("particle");
146 0 : daugStr = parser.readAttribute("daughters");
147 :
148 0 : std::string copyParent = parser.readAttribute("copy");
149 0 : std::string copyDaugStr = parser.readAttribute("copyDaughters");
150 :
151 0 : checkParticle(decayParent);
152 0 : ipar=EvtPDL::getId(decayParent);
153 :
154 0 : checkParticle(copyParent);
155 0 : EvtId copypar=EvtPDL::getId(copyParent);
156 :
157 0 : std::istringstream daugStream(daugStr);
158 0 : std::istringstream copyDaugStream(copyDaugStr);
159 :
160 0 : std::string daugh;
161 0 : while(std::getline(daugStream, daugh, ' ')) {
162 0 : checkParticle(daugh);
163 0 : daughter[nDaughters++] = EvtPDL::getId(daugh);
164 : }
165 0 : while(std::getline(copyDaugStream, daugh, ' ')) {
166 0 : checkParticle(daugh);
167 0 : copyDaughter[nCopyDaughters++] = EvtPDL::getId(daugh);
168 : }
169 :
170 0 : if(nDaughters!=3 || nCopyDaughters!=3) {
171 0 : report(Severity::Error,"EvtGen") <<
172 0 : "Expected to find three daughters for copyDecay of "<<decayParent<<
173 0 : " from "<<copyParent<<" near line "<<parser.getLineNumber()<<
174 0 : ", "<<"found "<<nDaughters<<" and "<<nCopyDaughters<<endl;
175 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
176 0 : ::abort();
177 : }
178 :
179 0 : copyDecay(ipar, daughter, copypar, copyDaughter);
180 :
181 0 : } else if(parser.getTagTitle() == "/data") { //end of data
182 : endReached = true;
183 0 : parser.close();
184 : break;
185 : }
186 : //TAGS FOUND UNDER DALITZDECAY
187 0 : } else if(parser.getParentTagTitle() == "dalitzDecay") {
188 0 : if(parser.getTagTitle() == "resonance") {
189 :
190 0 : flatteParams.clear();
191 :
192 : //Amplitude
193 0 : EvtComplex ampFactor(parser.readAttributeDouble("ampFactorReal",1.),
194 0 : parser.readAttributeDouble("ampFactorImag",0.));
195 0 : double mag = parser.readAttributeDouble("mag",-999.);
196 0 : double phase = parser.readAttributeDouble("phase",-999.);
197 0 : double real = parser.readAttributeDouble("real",-999.);
198 0 : double imag = parser.readAttributeDouble("imag",-999.);
199 :
200 0 : if((real!=-999. || imag!=-999.) && mag==-999. && phase==-999.) {
201 0 : if(real==-999.) { real = 0; }
202 0 : if(imag==-999.) { imag = 0; }
203 0 : mag = sqrt(real*real + imag*imag);
204 0 : phase = atan2(imag,real) * EvtConst::radToDegrees;
205 0 : }
206 0 : if( mag==-999. ) {
207 : mag = 1.;
208 0 : }
209 0 : if( phase==-999. ) {
210 : phase = 0.;
211 0 : }
212 :
213 0 : cAmp = ampFactor*EvtComplex(cos(phase*1.0/EvtConst::radToDegrees),sin(phase*1.0/EvtConst::radToDegrees))*mag;
214 :
215 : //Resonance particle properties
216 : mass = 0.;
217 : width = 0.;
218 : spinType = EvtSpinType::SCALAR;
219 :
220 0 : std::string particle = parser.readAttribute("particle");
221 0 : if(particle != "") {
222 0 : EvtId resId = EvtPDL::getId(particle);
223 0 : if(resId == EvtId(-1,-1)) {
224 0 : report(Severity::Error,"EvtGen") <<"Unknown particle name:"<<particle.c_str()<<endl;
225 0 : report(Severity::Error,"EvtGen") <<"Will terminate execution!"<<endl;
226 0 : ::abort();
227 : } else {
228 0 : mass = EvtPDL::getMeanMass(resId);
229 0 : width = EvtPDL::getWidth(resId);
230 0 : spinType = EvtPDL::getSpinType(resId);
231 : }
232 0 : }
233 :
234 0 : width = parser.readAttributeDouble("width",width);
235 0 : mass = parser.readAttributeDouble("mass",mass);
236 0 : switch(parser.readAttributeInt("spin",-1)) {
237 : case -1://not set here
238 : break;
239 : case 0:
240 : spinType = EvtSpinType::SCALAR;
241 0 : break;
242 : case 1:
243 : spinType = EvtSpinType::VECTOR;
244 0 : break;
245 : case 2:
246 : spinType = EvtSpinType::TENSOR;
247 0 : break;
248 : default:
249 0 : report(Severity::Error,"EvtGen") << "Unsupported spin near line "<<parser.getLineNumber()<<" of XML file."<<endl;
250 0 : ::abort();
251 : }
252 :
253 : //Shape and form factors
254 0 : shape = parser.readAttribute("shape");
255 0 : FFp = parser.readAttributeDouble("BlattWeisskopfFactorParent",0.0);
256 0 : FFr = parser.readAttributeDouble("BlattWeisskopfFactorResonance",1.5);
257 :
258 : //Shape specific attributes
259 0 : if(shape=="NonRes_Exp") {
260 0 : alpha = parser.readAttributeDouble("alpha",0.0);
261 0 : }
262 0 : if(shape=="LASS") {
263 0 : aLass = parser.readAttributeDouble("a",2.07);
264 0 : rLass = parser.readAttributeDouble("r",3.32);
265 0 : BLass = parser.readAttributeDouble("B",1.0);
266 0 : phiBLass = parser.readAttributeDouble("phiB",0.0);
267 0 : RLass = parser.readAttributeDouble("R",1.0);
268 0 : phiRLass = parser.readAttributeDouble("phiR",0.0);
269 0 : cutoffLass = parser.readAttributeDouble("cutoff",-1.0);
270 0 : }
271 :
272 : //Daughter pairs for resonance
273 0 : angAndResPairs.clear();
274 :
275 0 : std::string resDaugStr = parser.readAttribute("resDaughters");
276 0 : if(resDaugStr != "") {
277 0 : std::istringstream resDaugStream(resDaugStr);
278 0 : std::string resDaug;
279 : int nResDaug(0);
280 0 : EvtId resDaughter[2];
281 0 : while(std::getline(resDaugStream, resDaug, ' ')) {
282 0 : checkParticle(resDaug);
283 0 : resDaughter[nResDaug++] = EvtPDL::getId(resDaug);
284 : }
285 0 : if(nResDaug != 2) {
286 0 : report(Severity::Error,"EvtGen") << "Resonance must have exactly 2 daughters near line "<<parser.getLineNumber()<<" of XML file."<<endl;
287 0 : ::abort();
288 : }
289 0 : int nRes = getDaughterPairs(resDaughter, daughter, angAndResPairs);
290 0 : if(nRes==0) {
291 0 : report(Severity::Error,"EvtGen") << "Resonance daughters must match decay daughters near line "<<parser.getLineNumber()<<" of XML file."<<endl;
292 0 : ::abort();
293 : }
294 0 : if(parser.readAttributeBool("normalise",true)) cAmp /= sqrt(nRes);
295 0 : }
296 0 : if(angAndResPairs.empty()) {
297 0 : switch(parser.readAttributeInt("daughterPair")) {
298 : case 1:
299 0 : angAndResPairs.push_back(std::make_pair(EvtCyclic3::BC,EvtCyclic3::AB));
300 : break;
301 : case 2:
302 0 : angAndResPairs.push_back(std::make_pair(EvtCyclic3::CA,EvtCyclic3::BC));
303 : break;
304 : case 3:
305 0 : angAndResPairs.push_back(std::make_pair(EvtCyclic3::AB,EvtCyclic3::CA));
306 : break;
307 : default:
308 0 : if(shape == "NonRes") { //We don't expect a pair for non-resonant terms but add dummy values for convenience
309 0 : angAndResPairs.push_back(std::make_pair(EvtCyclic3::AB,EvtCyclic3::AB));
310 : break;
311 : }
312 0 : report(Severity::Error,"EvtGen") << "Daughter pair must be 1, 2 or 3 near line "<<parser.getLineNumber()<<" of XML file."<<endl;
313 0 : ::abort();
314 : }
315 : }
316 :
317 0 : if(parser.isTagInline()) {
318 0 : std::vector< std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> >::iterator it = angAndResPairs.begin();
319 0 : for( ; it != angAndResPairs.end(); it++) {
320 0 : std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> pairs = *it;
321 0 : EvtDalitzReso resonance = getResonance(shape, dp, pairs.first, pairs.second, spinType, mass, width, FFp, FFr, alpha, aLass, rLass, BLass, phiBLass, RLass, phiRLass, cutoffLass);
322 0 : dalitzDecay->addResonance(cAmp,resonance);
323 0 : }
324 0 : }
325 0 : } else if(parser.getTagTitle() == "/dalitzDecay") {
326 0 : if(probMax < 0) {
327 0 : report(Severity::Info,"EvtGen") << "probMax is not defined for " << decayParent << " -> " << daugStr << endl;
328 0 : report(Severity::Info,"EvtGen") << "Will now estimate probMax. This may take a while. Once probMax is calculated, update the XML file to skip this step in future." << endl;
329 0 : probMax = calcProbMax(dp,dalitzDecay);
330 0 : }
331 0 : dalitzDecay->setProbMax(probMax);
332 0 : addDecay(ipar, *dalitzDecay);
333 0 : delete dalitzDecay;
334 : dalitzDecay = 0;
335 0 : } else if(verbose) {
336 0 : report(Severity::Info,"EvtGen") << "Unexpected tag "<<parser.getTagTitle()
337 0 : <<" found in XML decay file near line "
338 0 : <<parser.getLineNumber()<<". Tag will be ignored."<<endl;
339 0 : }
340 : //TAGS FOUND UNDER RESONANCE
341 0 : } else if(parser.getParentTagTitle() == "resonance"){
342 0 : if(parser.getTagTitle() == "flatteParam") {
343 0 : EvtFlatteParam param(parser.readAttributeDouble("mass1"),
344 0 : parser.readAttributeDouble("mass2"),
345 0 : parser.readAttributeDouble("g"));
346 0 : flatteParams.push_back(param);
347 0 : } else if(parser.getTagTitle() == "/resonance") {
348 0 : std::vector< std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> >::iterator it = angAndResPairs.begin();
349 0 : for( ; it != angAndResPairs.end(); it++) {
350 0 : std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> pairs = *it;
351 0 : EvtDalitzReso resonance = getResonance(shape, dp, pairs.first, pairs.second, spinType, mass, width, FFp, FFr, alpha, aLass, rLass, BLass, phiBLass, RLass, phiRLass, cutoffLass);
352 :
353 0 : std::vector<EvtFlatteParam>::iterator flatteIt = flatteParams.begin();
354 0 : for( ; flatteIt != flatteParams.end(); flatteIt++) {
355 0 : resonance.addFlatteParam((*flatteIt));
356 : }
357 :
358 0 : dalitzDecay->addResonance(cAmp,resonance);
359 0 : }
360 0 : }
361 : }
362 : }
363 :
364 0 : if(!endReached) {
365 0 : report(Severity::Error,"EvtGen") << "Either the decay file ended prematurely or the file is badly formed.\n"
366 0 : <<"Error occured near line"<<parser.getLineNumber()<<endl;
367 0 : ::abort();
368 : }
369 0 : }
370 :
371 : void EvtDalitzTable::checkParticle(std::string particle) {
372 0 : if (EvtPDL::getId(particle)==EvtId(-1,-1)) {
373 0 : report(Severity::Error,"EvtGen") <<"Unknown particle name:"<<particle.c_str()<<endl;
374 0 : report(Severity::Error,"EvtGen") <<"Will terminate execution!"<<endl;
375 0 : ::abort();
376 : }
377 0 : }
378 :
379 : void EvtDalitzTable::addDecay(EvtId parent, const EvtDalitzDecayInfo& dec) {
380 0 : if(_dalitztable.find(parent)!=_dalitztable.end()) {
381 0 : _dalitztable[parent].push_back(dec);
382 : } else {
383 : _dalitztable[parent].push_back(dec);
384 : }
385 0 : }
386 :
387 : void EvtDalitzTable::copyDecay(EvtId parent, EvtId* daughters,
388 : EvtId copy, EvtId* copyd) {
389 0 : EvtDalitzDecayInfo decay(daughters[0],daughters[1],daughters[2]);
390 0 : std::vector<EvtDalitzDecayInfo> copyTable = getDalitzTable(copy);
391 0 : std::vector<EvtDalitzDecayInfo>::iterator i = copyTable.begin();
392 0 : for( ; i != copyTable.end(); i++) {
393 0 : EvtId daughter1 = (*i).daughter1();
394 0 : EvtId daughter2 = (*i).daughter2();
395 0 : EvtId daughter3 = (*i).daughter3();
396 :
397 0 : if((copyd[0] == daughter1 && copyd[1] == daughter2 && copyd[2] == daughter3) ||
398 0 : (copyd[0] == daughter1 && copyd[1] == daughter3 && copyd[2] == daughter2) ||
399 0 : (copyd[0] == daughter2 && copyd[1] == daughter1 && copyd[2] == daughter3) ||
400 0 : (copyd[0] == daughter2 && copyd[1] == daughter3 && copyd[2] == daughter1) ||
401 0 : (copyd[0] == daughter3 && copyd[1] == daughter1 && copyd[2] == daughter2) ||
402 0 : (copyd[0] == daughter3 && copyd[1] == daughter2 && copyd[2] == daughter1)) {
403 0 : decay.setProbMax((*i).getProbMax());
404 0 : std::vector<std::pair<EvtComplex, EvtDalitzReso> >::const_iterator j = (*i).getResonances().begin();
405 0 : for( ; j != (*i).getResonances().end(); j++) {
406 0 : decay.addResonance((*j));
407 : }
408 0 : addDecay(parent,decay);
409 : return;
410 0 : }
411 0 : }
412 : //if we get here then there was no match
413 0 : report(Severity::Error,"EvtGen") << "Did not find dalitz decays for particle:"
414 0 : <<copy<<"\n";
415 0 : }
416 :
417 : std::vector<EvtDalitzDecayInfo> EvtDalitzTable::getDalitzTable(const EvtId& parent) {
418 0 : std::vector<EvtDalitzDecayInfo> table;
419 0 : if ( _dalitztable.find(parent)!=_dalitztable.end() ) {
420 0 : table=_dalitztable[parent];
421 : }
422 :
423 0 : if (table.empty()){
424 0 : report(Severity::Error,"EvtGen") << "Did not find dalitz decays for particle:"
425 0 : <<parent<<"\n";
426 : }
427 :
428 : return table;
429 0 : }
430 :
431 :
432 : EvtDalitzReso EvtDalitzTable::getResonance(std::string shape, EvtDalitzPlot dp, EvtCyclic3::Pair angPair, EvtCyclic3::Pair resPair,
433 : EvtSpinType::spintype spinType, double mass, double width, double FFp, double FFr, double alpha,
434 : double aLass, double rLass, double BLass, double phiBLass, double RLass, double phiRLass, double cutoffLass) {
435 0 : if( shape=="RBW" || shape=="RBW_CLEO") {
436 0 : return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::RBW_CLEO, FFp, FFr );
437 0 : } else if( shape=="RBW_CLEO_ZEMACH" ) {
438 0 : return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::RBW_CLEO_ZEMACH, FFp, FFr );
439 0 : } else if( shape=="GS" || shape=="GS_CLEO" ) {
440 0 : return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::GS_CLEO, FFp, FFr );
441 0 : } else if( shape=="GS_CLEO_ZEMACH" ) {
442 0 : return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::GS_CLEO_ZEMACH, FFp, FFr );
443 0 : } else if( shape=="GAUSS" || shape=="GAUSS_CLEO" ) {
444 0 : return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::GAUSS_CLEO, FFp, FFr );
445 0 : } else if( shape=="GAUSS_CLEO_ZEMACH" ) {
446 0 : return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::GAUSS_CLEO_ZEMACH, FFp, FFr );
447 0 : } else if( shape=="Flatte" ) {
448 0 : return EvtDalitzReso( dp, resPair, mass );
449 0 : } else if( shape=="LASS" ) {
450 0 : return EvtDalitzReso( dp, resPair, mass, width, aLass, rLass, BLass, phiBLass, RLass, phiRLass, cutoffLass, true );
451 0 : } else if( shape=="NonRes" ) {
452 0 : return EvtDalitzReso( );
453 0 : } else if( shape=="NonRes_Linear" ) {
454 0 : return EvtDalitzReso( dp, resPair, EvtDalitzReso::NON_RES_LIN );
455 0 : } else if( shape=="NonRes_Exp" ) {
456 0 : return EvtDalitzReso( dp, resPair, EvtDalitzReso::NON_RES_EXP, alpha );
457 : } else { //NBW
458 0 : if( shape!="NBW") report(Severity::Warning,"EvtGen")<<"EvtDalitzTable: shape "<<shape<<" is unknown. Defaulting to NBW."<<endl;
459 0 : return EvtDalitzReso( dp, angPair, resPair, spinType, mass, width, EvtDalitzReso::NBW, FFp, FFr );
460 : }
461 0 : }
462 :
463 : int EvtDalitzTable::getDaughterPairs(EvtId* resDaughter, EvtId* daughter, std::vector< std::pair<EvtCyclic3::Pair,EvtCyclic3::Pair> >& angAndResPairs) {
464 : int n(0);
465 0 : if(resDaughter[0]==daughter[0] && resDaughter[1]==daughter[1]) {
466 0 : angAndResPairs.push_back(std::make_pair(EvtCyclic3::BC,EvtCyclic3::AB)); n++;
467 0 : } else if(resDaughter[0]==daughter[1] && resDaughter[1]==daughter[0]) {
468 0 : angAndResPairs.push_back(std::make_pair(EvtCyclic3::CA,EvtCyclic3::AB)); n++;
469 0 : }
470 :
471 0 : if(resDaughter[0]==daughter[1] && resDaughter[1]==daughter[2]) {
472 0 : angAndResPairs.push_back(std::make_pair(EvtCyclic3::CA,EvtCyclic3::BC)); n++;
473 0 : } else if(resDaughter[0]==daughter[2] && resDaughter[1]==daughter[1]) {
474 0 : angAndResPairs.push_back(std::make_pair(EvtCyclic3::AB,EvtCyclic3::BC)); n++;
475 0 : }
476 :
477 0 : if(resDaughter[0]==daughter[2] && resDaughter[1]==daughter[0]) {
478 0 : angAndResPairs.push_back(std::make_pair(EvtCyclic3::AB,EvtCyclic3::CA)); n++;
479 0 : } else if(resDaughter[0]==daughter[0] && resDaughter[1]==daughter[2]) {
480 0 : angAndResPairs.push_back(std::make_pair(EvtCyclic3::BC,EvtCyclic3::CA)); n++;
481 0 : }
482 :
483 0 : return n;
484 : }
485 :
486 : double EvtDalitzTable::calcProbMax(EvtDalitzPlot dp, EvtDalitzDecayInfo* model) {
487 :
488 : double factor = 1.2; //factor to increase our final answer by
489 : int nStep(1000); //number of steps - total points will be 3*nStep*nStep
490 :
491 : double maxProb(0);
492 : double min(0), max(0), step(0), min2(0), max2(0), step2(0);
493 :
494 : //first do AB, BC
495 0 : min = dp.qAbsMin(EvtCyclic3::AB);
496 0 : max = dp.qAbsMax(EvtCyclic3::AB);
497 0 : step = (max-min)/nStep;
498 0 : for(int i=0; i<nStep; ++i) {
499 0 : double qAB = min + i*step;
500 0 : min2 = dp.qMin(EvtCyclic3::BC,EvtCyclic3::AB,qAB);
501 0 : max2 = dp.qMax(EvtCyclic3::BC,EvtCyclic3::AB,qAB);
502 0 : step2 = (max2-min2)/nStep;
503 0 : for(int j=0; j<nStep; ++j) {
504 0 : double qBC = min2+ j*step2;
505 0 : EvtDalitzCoord coord(EvtCyclic3::AB,qAB,EvtCyclic3::BC,qBC);
506 0 : EvtDalitzPoint point(dp,coord);
507 0 : double prob = calcProb(point,model);
508 0 : if(prob > maxProb) maxProb = prob;
509 0 : }
510 : }
511 :
512 : //next do BC, CA
513 0 : min = dp.qAbsMin(EvtCyclic3::BC);
514 0 : max = dp.qAbsMax(EvtCyclic3::BC);
515 0 : step = (max-min)/nStep;
516 0 : for(int i=0; i<nStep; ++i) {
517 0 : double qBC = min + i*step;
518 0 : min2 = dp.qMin(EvtCyclic3::CA,EvtCyclic3::BC,qBC);
519 0 : max2 = dp.qMax(EvtCyclic3::CA,EvtCyclic3::BC,qBC);
520 0 : step2 = (max2-min2)/nStep;
521 0 : for(int j=0; j<nStep; ++j) {
522 0 : double qCA = min2+ j*step2;
523 0 : EvtDalitzCoord coord(EvtCyclic3::BC,qBC,EvtCyclic3::CA,qCA);
524 0 : EvtDalitzPoint point(dp,coord);
525 0 : double prob = calcProb(point,model);
526 0 : if(prob > maxProb) maxProb = prob;
527 0 : }
528 : }
529 :
530 : //finally do CA, AB
531 0 : min = dp.qAbsMin(EvtCyclic3::CA);
532 0 : max = dp.qAbsMax(EvtCyclic3::CA);
533 0 : step = (max-min)/nStep;
534 0 : for(int i=0; i<nStep; ++i) {
535 0 : double qCA = min + i*step;
536 0 : min2 = dp.qMin(EvtCyclic3::AB,EvtCyclic3::CA,qCA);
537 0 : max2 = dp.qMax(EvtCyclic3::AB,EvtCyclic3::CA,qCA);
538 0 : step2 = (max2-min2)/nStep;
539 0 : for(int j=0; j<nStep; ++j) {
540 0 : double qAB = min2+ j*step2;
541 0 : EvtDalitzCoord coord(EvtCyclic3::CA,qCA,EvtCyclic3::AB,qAB);
542 0 : EvtDalitzPoint point(dp,coord);
543 0 : double prob = calcProb(point,model);
544 0 : if(prob > maxProb) maxProb = prob;
545 0 : }
546 : }
547 0 : report(Severity::Info,"EvtGen") << "Largest probability found was " << maxProb << endl;
548 0 : report(Severity::Info,"EvtGen") << "Setting probMax to " << factor*maxProb << endl;
549 0 : return factor*maxProb;
550 0 : }
551 :
552 : double EvtDalitzTable::calcProb(EvtDalitzPoint point, EvtDalitzDecayInfo* model) {
553 :
554 0 : std::vector<std::pair<EvtComplex,EvtDalitzReso> > resonances = model->getResonances();
555 :
556 0 : EvtComplex amp(0,0);
557 0 : std::vector<std::pair<EvtComplex,EvtDalitzReso> >::iterator i = resonances.begin();
558 0 : for( ; i!= resonances.end(); i++) {
559 0 : std::pair<EvtComplex,EvtDalitzReso> res = (*i);
560 0 : amp += res.first * res.second.evaluate( point );
561 0 : }
562 0 : return abs2(amp);
563 0 : }
|
