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