PhotosBranch.cxx
1#include <vector>
2#include <list>
3#include "HEPEVT_struct.h"
4#include "PhotosParticle.h"
5#include "PhotosBranch.h"
6#include "Photos.h"
7#include "Log.h"
8using std::vector;
9using std::list;
10using std::endl;
11
12namespace Photospp
13{
14
16{
18
19 //Suppress if somehow got stable particle
20 if(daughters.size()==0)
21 {
22 Log::Debug(1)<<"Stable particle."<<endl;
23 suppression = 1;
24 forcing = 0;
25 particle = NULL;
26 return;
27 }
28 else if(daughters.at(0)->getMothers().size()==1)
29 {
30 // Regular case - one mother
31 Log::Debug(1)<<"Regular case."<<endl;
32 particle = p;
34 }
35 else
36 {
37 // Advanced case - branch with multiple mothers - no mid-particle
38 Log::Debug(1)<<"Advanced case."<<endl;
39 particle = NULL;
40 mothers = daughters.at(0)->getMothers();
41 }
42
43 //--------------------------------------------------
44 // Finalize suppression/forcing checks
45 // NOTE: if user forces decay of specific particle,
46 // this overrides any suppresion
47 //--------------------------------------------------
48
51 else suppression = 0;
52
53 // Even if forced or passed suppression check, we still have to check few things
54 if(!suppression)
55 {
56 // Check momentum conservation
58 if(suppression) Log::Warning()<<"Branching ignored due to 4-momentum non conservation"<<endl;
59
60 // Check if advanced case has only one daughter
61 if(!particle && daughters.size()==1) suppression=-1;
62
63 // If any of special cases is true, we're not forcing this branch
64 if(suppression) forcing=0;
65 }
66}
67
69{
70 Log::Debug(703)<<" Processing barcode: "<<( (particle) ? particle->getBarcode() : ( (mothers.size()) ? mothers.at(0)->getBarcode() : -1) )<<endl;
71 /*
72 cout<<"Particles send to photos (with barcodes in brackets):"<<endl;
73 vector<PhotosParticle *> get = getParticles();
74 for(int i=0;i<(int)get.size();i++) cout<<"ID: "<<get.at(i)->getPdgID()<<" ("<<get.at(i)->getBarcode()<<"), "; cout<<endl;
75 */
76 int index = HEPEVT_struct::set(this);
78 PHOTOS_MAKE_C(index);
82}
83
84vector<PhotosParticle *> PhotosBranch::getParticles()
85{
86 vector<PhotosParticle *> ret = mothers;
87 if(particle) ret.push_back(particle);
88 ret.insert(ret.end(),daughters.begin(),daughters.end());
89 return ret;
90}
91
93{
95 if(mothers.size()>0) return mothers.at(0)->checkMomentumConservation();
96 return true;
97}
98
99vector<PhotosBranch *> PhotosBranch::createBranches(vector<PhotosParticle *> particles)
100{
101 Log::Debug(700)<<"PhotosBranch::createBranches - filtering started"<<endl;
102 list<PhotosParticle *> list(particles.begin(),particles.end());
103 vector<PhotosBranch *> branches;
104
105 // First - add all forced decays
107 {
108 std::list<PhotosParticle *>::iterator it;
109 for(it=list.begin();it!=list.end();it++)
110 {
111 PhotosBranch *branch = new PhotosBranch(*it);
112 int forcing = branch->getForcingStatus();
113 if(forcing)
114 {
115 Log::Debug(701)<<" Forced: "<<(*it)->getPdgID()<<" (barcode: "<<(*it)->getBarcode()<<") with forcing status= "<<forcing<<endl;
116 branches.push_back(branch);
117 it = list.erase(it);
118 --it;
119 // If forcing consecutive decays
120 if(forcing==2)
121 {
122 PhotosParticle *p = branch->getDecayingParticle();
123 if(!p)
124 {
125 if(branch->getMothers().size()>0) p = branch->getMothers().at(0);
126 else continue;
127 }
128 vector<PhotosParticle *> tree = p->getDecayTree();
129 //Add branches for all particles from the list - max O(n*m)
130 std::list<PhotosParticle *>::iterator it2;
131 for(it2=list.begin();it2!=list.end();it2++)
132 {
133 for(int i=0;i<(int)tree.size();i++)
134 {
135 if(tree.at(i)->getBarcode()==(*it2)->getBarcode())
136 {
137 PhotosBranch *b = new PhotosBranch(*it2);
138 branches.push_back(b);
139 // If we were to delete our next particle in line
140 if(it==it2) --it;
141 it2 = list.erase(it2);
142 --it2;
143 break;
144 }
145 }
146 }
147 }
148 }
149 else delete branch;
150 }
151 }
152 // Quit if we're suppressing everything
153 if(Photos::isSuppressed) return branches;
154 // Now - check if remaining decays are suppressed
155 while(!list.empty())
156 {
157 PhotosParticle *particle = list.front();
158 list.pop_front();
159 if(!particle) continue;
160
161 PhotosBranch *branch = new PhotosBranch(particle);
162 int suppression = branch->getSuppressionStatus();
163 if(!suppression) branches.push_back(branch);
164 else
165 {
166 Log::Debug(702)<<" Suppressed: "<<particle->getPdgID()<<" (barcode: "<<particle->getBarcode()<<") with suppression status= "<<suppression<<endl;
167 //If suppressing consecutive decays
168 if(suppression==2)
169 {
170 PhotosParticle *p = branch->getDecayingParticle();
171 if(!p)
172 {
173 if(branch->getMothers().size()>0) p = branch->getMothers().at(0);
174 else continue;
175 }
176 vector<PhotosParticle *> tree = p->getDecayTree();
177 //Remove all particles from the list - max O(n*m)
178 std::list<PhotosParticle *>::iterator it;
179 for(it=list.begin();it!=list.end();it++)
180 {
181 for(int i=0;i<(int)tree.size();i++)
182 {
183 if(tree.at(i)->getBarcode()==(*it)->getBarcode())
184 {
185 it = list.erase(it);
186 --it;
187 break;
188 }
189 }
190 }
191 }
192 delete branch;
193 continue;
194 }
195
196 //In case we don't have mid-particle erase rest of the mothers from list
197 if(!branch->getDecayingParticle())
198 {
199 vector<PhotosParticle *> mothers = branch->getMothers();
200 for(int i=0;i<(int)mothers.size();i++)
201 {
202 PhotosParticle *m = mothers.at(i);
203 if(m->getBarcode()==particle->getBarcode()) continue;
204 std::list<PhotosParticle *>::iterator it;
205 for(it=list.begin();it!=list.end();it++)
206 if(m->getBarcode()==(*it)->getBarcode())
207 {
208 it = list.erase(it);
209 break;
210 }
211 }
212 }
213 }
214 return branches;
215}
216
217int PhotosBranch::checkList(bool forceOrSuppress)
218{
219 vector< vector<int>* > *list = (forceOrSuppress) ? Photos::forceBremList : Photos::supBremList;
220 if(!list) return 0;
221
222 // Can't check without pdgid
223 int motherID;
224 if(particle) motherID = particle->getPdgID();
225 else
226 {
227 if(mothers.size()==0) return 0;
228 motherID = mothers.at(0)->getPdgID();
229 }
230
231 // Create list of daughters
232 vector<int> dID;
233 for(int j=0;j<(int)daughters.size();j++) dID.push_back(daughters[j]->getPdgID());
234
235 vector< vector<int> *> &patternList = *list;
236
237 // Check if the mother and list of daughters matches any of the declared patterns
238 for(int j=0; j<(int)patternList.size();j++)
239 {
240 // Skip patterns that don't have our mother
241 if(motherID!=(*patternList[j])[0]) continue;
242
243 // Compare decay daughters with pattern - max O(n*m)
244 vector<int> &pattern = *patternList[j];
245 bool fullMatch=true;
246 for(int k = 1; k<(int)pattern.size()-1; k++)
247 {
248 bool oneMatch=false;
249 for(int l=0;l<(int)dID.size(); l++)
250 if(pattern[k]==dID[l]) { oneMatch=true; break; }
251 if(!oneMatch) { fullMatch=false; break; }
252 }
253 // Check if the matching pattern is set for consecutive suppression
254 /*
255 Currently minimal matching is used.
256 If e.g. 25 -> -15 is suppressed, then 25 -> 15,-15 and 25 -> 15,-15,22 etc. is suppressed too
257 For exact matching (then suppress 25 -> 15,-15 ; 25 -> 15,-15,22 etc. must be done separately) uncoment line ...:
258 */
259 // if(pattern.size()<=2 || (fullMatch && dID.size()==pattern.size()-2) )
260 // ...and comment out line:
261 if(pattern.size()<=2 || fullMatch)
262 return (pattern.back()==1) ? 2 : 1;
263 }
264 return 0;
265}
266
267} // namespace Photospp
static int set(PhotosBranch *branch)
static ostream & Debug(unsigned short int code=0, bool count=true)
Definition Log.cxx:33
static vector< PhotosBranch * > createBranches(vector< PhotosParticle * > particles)
int checkList(bool forceOrSuppress)
vector< PhotosParticle * > getMothers()
vector< PhotosParticle * > mothers
vector< PhotosParticle * > getParticles()
vector< PhotosParticle * > daughters
PhotosBranch(PhotosParticle *p)
virtual std::vector< PhotosParticle * > getDaughters()=0
virtual bool checkMomentumConservation()=0
virtual int getBarcode()=0
virtual int getPdgID()=0
std::vector< PhotosParticle * > getDecayTree()
std::vector< PhotosParticle * > findProductionMothers()
static vector< vector< int > * > * supBremList
static vector< vector< int > * > * forceBremList