Line data Source code
1 : #include <vector>
2 : #include <cmath>
3 : #include "PhotosBranch.h"
4 : #include "PhotosParticle.h"
5 : #include "PH_HEPEVT_Interface.h"
6 : #include "Log.h"
7 : using namespace std;
8 :
9 : namespace Photospp
10 : {
11 :
12 6 : vector<PhotosParticle*> PH_HEPEVT_Interface::m_particle_list;
13 : int PH_HEPEVT_Interface::ME_channel=0;
14 : int PH_HEPEVT_Interface::decay_idx=0;
15 :
16 : void PH_HEPEVT_Interface::clear(){
17 :
18 0 : m_particle_list.clear();
19 :
20 0 : hep.nevhep=0;
21 0 : hep.nhep=0;
22 :
23 :
24 : /** for(int i=0; i < NMXHEP; i++){
25 :
26 : hep.isthep[i]=0;
27 : hep.idhep[i]=0;
28 :
29 : for(int j=0; j<2; j++){
30 : hep.jmohep[i][j]=0;
31 : hep.jdahep[i][j]=0;
32 : }
33 :
34 : for(int j=0; j<5; j++)
35 : hep.phep[i][j]=0;
36 :
37 : for(int j=0; j<4; j++)
38 : hep.vhep[i][j]=0;
39 :
40 : ph_phoqed_.qedrad[i]=0;
41 :
42 : }**/
43 0 : }
44 :
45 : void PH_HEPEVT_Interface::add_particle(int i,PhotosParticle * particle,
46 : int first_mother, int last_mother,
47 : int first_daughter, int last_daughter){
48 :
49 0 : if(i>0)
50 0 : i--; //account for fortran indicies begining at 1
51 : else
52 0 : Log::Warning()<<"Index given to PH_HEPEVT_Interface::add_particle "
53 0 : <<"is too low (it must be > 0)."<<endl;
54 :
55 : //add to our internal list of pointer/index pairs
56 0 : m_particle_list.push_back(particle);
57 :
58 : //now set the element of PH_HEPEVT
59 0 : hep.nevhep=0; //dummy
60 0 : hep.nhep = hep.nhep + 1; // 1.II.2014: fix for gcc 4.8.1. Previously it was:
61 : // hep.nhep = hep.nhep++; which technically is an undefined operation
62 0 : hep.isthep[i]=particle->getStatus();
63 0 : hep.idhep[i]=particle->getPdgID();
64 :
65 0 : hep.jmohep[i][0]=first_mother;
66 0 : hep.jmohep[i][1]=last_mother;
67 :
68 0 : hep.jdahep[i][0]=first_daughter;
69 0 : hep.jdahep[i][1]=last_daughter;
70 :
71 0 : hep.phep[i][0]=particle->getPx();
72 0 : hep.phep[i][1]=particle->getPy();
73 0 : hep.phep[i][2]=particle->getPz();
74 0 : hep.phep[i][3]=particle->getE();
75 :
76 : // if massFrom4Vector=true (default) - get sqrt(e^2-p^2)
77 : // otherwise - get mass from event record
78 0 : if(!Photos::massFrom4Vector) hep.phep[i][4]=particle->getMass();
79 0 : else hep.phep[i][4]=particle->getVirtuality();
80 :
81 0 : int pdgid = abs(particle->getPdgID());
82 :
83 : // if 'forceMass' for this PDGID was used - overwrite mass
84 0 : if(Photos::forceMassList)
85 : {
86 0 : for(unsigned int j=0;j<Photos::forceMassList->size();j++)
87 : {
88 0 : if(pdgid == abs(Photos::forceMassList->at(j)->first))
89 : {
90 0 : double mass = Photos::forceMassList->at(j)->second;
91 :
92 : // when 'forceMass' is used the mass provided is larger than 0.0
93 : // when 'forceMassFromEventRecord' is used mass is -1.0
94 : // in this case - get mass from event record
95 0 : if(mass<0.0) mass = particle->getMass();
96 0 : hep.phep[i][4] = mass;
97 0 : }
98 : }
99 0 : }
100 :
101 0 : hep.vhep[i][0]=0;
102 0 : hep.vhep[i][1]=0;
103 0 : hep.vhep[i][2]=0;
104 0 : hep.vhep[i][3]=0;
105 :
106 0 : hep.qedrad[i]=1;
107 :
108 0 : }
109 :
110 : int PH_HEPEVT_Interface::set(PhotosBranch *branch)
111 : {
112 0 : PH_HEPEVT_Interface::clear();
113 : int idx=1;
114 :
115 : //get mothers
116 0 : vector<PhotosParticle *> mothers = branch->getMothers();
117 0 : int nmothers=mothers.size();
118 :
119 : //check if mid-particle exist
120 0 : decay_idx=0;
121 0 : PhotosParticle *decay_particle = branch->getDecayingParticle();
122 0 : if(decay_particle) decay_idx=nmothers+1;
123 :
124 : //get daughters
125 0 : vector<PhotosParticle *> daughters = branch->getDaughters();
126 0 : int ndaughters=daughters.size();
127 :
128 0 : for(int i=0;i<nmothers;i++)
129 : {
130 0 : if(decay_idx)
131 0 : add_particle(idx++,mothers.at(i),
132 : 0,0, //mothers
133 0 : decay_idx,decay_idx); //daughters
134 : else
135 0 : add_particle(idx++,mothers.at(i),
136 : 0,0, //mothers
137 0 : nmothers+1,nmothers+ndaughters); //daughters
138 : }
139 :
140 0 : if(decay_particle)
141 0 : add_particle(idx++,decay_particle,
142 : 1,nmothers, //mothers
143 0 : nmothers+2,nmothers+1+ndaughters); //daughters
144 :
145 0 : for(int i=0;i<ndaughters;i++)
146 : {
147 0 : if(decay_idx)
148 0 : add_particle(idx++,daughters.at(i),
149 0 : decay_idx,decay_idx, //mothers
150 : 0,0); //daughters
151 : else
152 0 : add_particle(idx++,daughters.at(i),
153 : 1,nmothers, //mothers
154 : 0,0); //daughters
155 : }
156 : //Log::RedirectOutput( phodmp_ , Log::Debug(1000) );
157 0 : Log::Debug(1000,false)<<"PH_HEPEVT returning: "<<( (decay_idx) ? decay_idx : 1 )<<" from "<<idx-1<<" particles."<<endl;
158 0 : return (decay_idx) ? decay_idx : 1;
159 0 : }
160 :
161 : void PH_HEPEVT_Interface::get(){
162 :
163 : int index = 0;
164 :
165 : //if no photons have been added to the event record, do nothing.
166 0 : if(hep.nhep == (int) m_particle_list.size())
167 0 : return;
168 :
169 : //phodmp_();
170 :
171 0 : int particle_count = m_particle_list.size();
172 0 : int daughters_start = hep.jmohep[hep.nhep-1][0];
173 0 : int photons = hep.nhep - m_particle_list.size();
174 0 : bool isPhotonCreated = (photons>0);
175 :
176 0 : std::vector<PhotosParticle*> photon_list; // list of added photons
177 : // which need kinematical treatment
178 : // in special case
179 :
180 : // we decipher daughters_start from last entry
181 : // that is last daughter in ph_hepevt_
182 : // another option of this functionality may be
183 : // hep.jdahep[ hep.jmohep[hep.nhep-1][0]-1][0];
184 : // Update daughters_start if there are two mothers
185 : // NOTE: daughters_start is index for C++ arrays, while hep.jmohep
186 : // contains indices for Fortran arrays.
187 0 : if(hep.jmohep[hep.nhep-1][1]>0)
188 0 : daughters_start = hep.jmohep[hep.nhep-1][1];
189 :
190 : index = particle_count;
191 :
192 : // Add extra photons
193 0 : for(;photons>0; photons--, index++){
194 :
195 0 : if(hep.idhep[index]!=PhotosParticle::GAMMA)
196 0 : Log::Fatal("PH_HEPEVT_Interface::get(): Extra particle added to the PH_HEPEVT common block in not a photon!",6);
197 :
198 : //create a new particle
199 0 : PhotosParticle * new_photon;
200 0 : new_photon = m_particle_list.at(0)->createNewParticle(hep.idhep[index],
201 0 : hep.isthep[index],
202 0 : hep.phep[index][4],
203 0 : hep.phep[index][0],
204 0 : hep.phep[index][1],
205 0 : hep.phep[index][2],
206 0 : hep.phep[index][3]);
207 :
208 : //add into the event record
209 : //get mother particle of photon
210 0 : PhotosParticle * mother = m_particle_list.at(hep.jmohep[index][0]-1);
211 0 : mother->addDaughter(new_photon);
212 :
213 : //add to list of photons
214 0 : photon_list.push_back(new_photon);
215 0 : }
216 :
217 : // Before we update particles, we check for special cases
218 : // At this step, particles are yet unmodified
219 : // but photons are already in the event record
220 : bool special=false;
221 : PhotosParticle *p1 = NULL;
222 : PhotosParticle *p2 = NULL;
223 :
224 0 : if( isPhotonCreated )
225 : {
226 0 : std::vector<PhotosParticle*> daughters;
227 :
228 : // in the following we create list of daughters,
229 : // later we calculate bool special which is true only if all
230 : // daughters self-decay
231 : // at peresent warning for mixed self-decay and not self decay
232 : // daughters is not printed.
233 :
234 0 : for(int i=daughters_start;i<particle_count;i++)
235 : {
236 0 : PhotosParticle *p = m_particle_list.at(i);
237 :
238 0 : daughters.push_back(p);
239 0 : }
240 :
241 : // Check if this is a special case
242 : special = true;
243 :
244 0 : if(daughters.size()==0) special = false;
245 :
246 : // special = false if there is a stable particle on the list
247 : // or there is a particle without self-decay
248 0 : for(unsigned int i=0;i<daughters.size();i++)
249 : {
250 0 : if(daughters[i]->getStatus()==1)
251 : {
252 : special = false;
253 0 : break;
254 : }
255 :
256 : // NOTE: We can use 'getDaughters' here, because vertices
257 : // of daughters are not being modified by Photos right now
258 : // (so there will be no caching)
259 0 : std::vector<PhotosParticle*> daughters2 = daughters[i]->getDaughters();
260 :
261 0 : if(daughters2.size()!=1 ||
262 0 : daughters2[0]->getPdgID() != daughters[i]->getPdgID() )
263 : {
264 : special = false;
265 0 : break;
266 : }
267 0 : }
268 :
269 0 : if( special )
270 : {
271 : double px1=0.0, py1=0.0, pz1=0.0, e1=0.0;
272 : double px2=0.0, py2=0.0, pz2=0.0, e2=0.0;
273 :
274 : // get sum of 4-momenta of unmodified particles
275 0 : for(unsigned int i=0;i<daughters.size();i++)
276 : {
277 : // ignore photons
278 0 : if(daughters[i]->getPdgID()==22) continue;
279 :
280 0 : px1+=daughters[i]->getPx();
281 0 : py1+=daughters[i]->getPy();
282 0 : pz1+=daughters[i]->getPz();
283 0 : e1 +=daughters[i]->getE();
284 0 : }
285 :
286 : // get sum of 4-momenta of particles in self-decay vertices
287 0 : for(unsigned int i=0;i<daughters.size();i++)
288 : {
289 : // ignore photons
290 0 : if(daughters[i]->getPdgID()==22) continue;
291 :
292 : // since 'allDaughtersSelfDecay()' is true
293 : // each of these particles has exactly one daughter
294 0 : px2 += daughters[i]->getDaughters().at(0)->getPx();
295 0 : py2 += daughters[i]->getDaughters().at(0)->getPy();
296 0 : pz2 += daughters[i]->getDaughters().at(0)->getPz();
297 0 : e2 += daughters[i]->getDaughters().at(0)->getE();
298 0 : }
299 :
300 : //cout<<"ORIG: "<<px1<<" "<<py1<<" "<<pz1<<" "<<e1<<endl;
301 : //cout<<"SELF: "<<px2<<" "<<py2<<" "<<pz2<<" "<<e2<<endl;
302 :
303 0 : p1 = m_particle_list.at(0)->createNewParticle(0,-1,0.0,px1,py1,pz1,e1);
304 0 : p2 = m_particle_list.at(0)->createNewParticle(0,-2,0.0,px2,py2,pz2,e2);
305 :
306 : // Finaly, boost photons to appropriate frame
307 0 : for(unsigned int i=0;i<photon_list.size();i++)
308 : {
309 0 : PhotosParticle *boosted = photon_list[i]->createNewParticle( 22, 1,
310 : 0.0,
311 0 : photon_list[i]->getPx(),
312 0 : photon_list[i]->getPy(),
313 0 : photon_list[i]->getPz(),
314 0 : photon_list[i]->getE() );
315 :
316 0 : boosted->boostToRestFrame(p1);
317 0 : boosted->boostFromRestFrame(p2);
318 :
319 0 : photon_list[i]->createSelfDecayVertex(boosted);
320 :
321 0 : delete boosted;
322 : }
323 :
324 0 : Log::Warning()<<"Hidden interaction, all daughters self decay."
325 0 : <<"Potentially over simplified solution applied."<<endl;
326 0 : }
327 0 : }
328 :
329 : //otherwise loop over particles which are already in the
330 : //event record and modify their 4 momentum
331 : //4.03.2012: Fix to prevent kinematical trap in vertex of simultaneous:
332 : // z-collinear and non-conservation pf E,p for dauthters of grandmothers
333 0 : for(index=daughters_start; index < particle_count && index < (int) m_particle_list.size(); index++){
334 :
335 0 : PhotosParticle * particle = m_particle_list.at(index);
336 :
337 0 : if(hep.idhep[index]!=particle->getPdgID())
338 0 : Log::Fatal("PH_HEPEVT_Interface::get(): Something is wrong with the PH_HEPEVT common block",5);
339 :
340 : // If photons were added - for each daughter create a history entry
341 0 : if(isPhotonCreated && Photos::isCreateHistoryEntries)
342 : {
343 0 : particle->createHistoryEntry();
344 : }
345 :
346 : //check to see if this particle's 4-momentum has been modified
347 : bool update=false;
348 :
349 : // don't update particle if difference lower than THRESHOLD * particle energy (default threshold = 10e-8)
350 0 : double threshold = NO_BOOST_THRESHOLD*hep.phep[index][3];
351 0 : if( fabs(hep.phep[index][0]-particle->getPx()) > threshold ||
352 0 : fabs(hep.phep[index][1]-particle->getPy()) > threshold ||
353 0 : fabs(hep.phep[index][2]-particle->getPz()) > threshold ||
354 0 : fabs(hep.phep[index][3]-particle->getE()) > threshold ) update=true;
355 :
356 0 : if(update)
357 : {
358 :
359 : //modify this particle's momentum and it's daughters momentum
360 : //Steps 1., 2. and 3. must be executed in order.
361 :
362 : //1. boost the particles daughters into it's (old) rest frame
363 0 : particle->boostDaughtersToRestFrame(particle);
364 :
365 : //2. change this particles 4 momentum
366 0 : particle->setPx(hep.phep[index][0]);
367 0 : particle->setPy(hep.phep[index][1]);
368 0 : particle->setPz(hep.phep[index][2]);
369 0 : particle->setE(hep.phep[index][3]);
370 :
371 : //3. boost the particles daughters back into the lab frame
372 0 : particle->boostDaughtersFromRestFrame(particle);
373 :
374 0 : if(special && particle->getDaughters().size()>0){
375 :
376 : // Algorithm for special case:
377 : // a. get self-daughter of 'particle'
378 0 : PhotosParticle *particled = particle->getDaughters().at(0);
379 :
380 : // b. boost 'particled' daughters to rest frame
381 0 : particled->boostDaughtersToRestFrame(particled);
382 :
383 : // c. copy four momentum of 'particle' into four momentum of
384 : // its self-daughter 'particled'
385 :
386 0 : particled->setPx( particle->getPx() );
387 0 : particled->setPy( particle->getPy() );
388 0 : particled->setPz( particle->getPz() );
389 0 : particled->setE ( particle->getE() );
390 :
391 : // d. boost self daughter to rest-frame of <e1>
392 : // boost self daughter from rest-frame of <e2>
393 :
394 0 : particled->boostToRestFrame(p1);
395 0 : particled->boostFromRestFrame(p2);
396 :
397 : // e. boost the 'particled' daughters back into the lab frame
398 0 : particled->boostDaughtersFromRestFrame(particled);
399 0 : }
400 :
401 : }
402 : }
403 :
404 : // cleanup
405 0 : if(p1) delete p1;
406 0 : if(p2) delete p2;
407 0 : }
408 :
409 : void PH_HEPEVT_Interface::prepare()
410 : {
411 0 : check_ME_channel();
412 0 : }
413 :
414 : void PH_HEPEVT_Interface::complete()
415 : {
416 :
417 0 : }
418 :
419 : void PH_HEPEVT_Interface::check_ME_channel()
420 : {
421 0 : ME_channel=0;
422 :
423 : // Check mothers:
424 :
425 0 : if(decay_idx==2) return; // Only one mother present
426 0 : if(hep.idhep[0]*hep.idhep[1]>0) return; // Mothers have same sign
427 :
428 0 : Log::Debug(900)<<"ME_channel: Mothers PDG: "<<hep.idhep[0]<<" "<<hep.idhep[1]<<endl;
429 0 : if(decay_idx)
430 0 : Log::Debug(900,false)<<" Intermediate: "<<hep.idhep[decay_idx-1]<<endl;
431 :
432 : int firstDaughter=3;
433 0 : if(decay_idx==0) firstDaughter=2; // if no intermediate particle - daughters start at idx 2
434 :
435 : // Are mothers in range +/- 1-6; +/- 11-16?
436 0 : int mother1 = abs(hep.idhep[0]);
437 0 : int mother2 = abs(hep.idhep[1]);
438 0 : if( mother1<1 || (mother1>6 && mother1<11) || mother1>16 ) return;
439 0 : if( mother2<1 || (mother2>6 && mother2<11) || mother2>16 ) return;
440 :
441 : //Check daughters
442 :
443 : // Z: check for pairs 11 -11 ; 13 -13 ; 15 -15
444 : // -------------------------------------------
445 : int firstPDG =0;
446 : int secondPDG=0;
447 0 : for(int i=firstDaughter; i<hep.nhep;i++)
448 : {
449 0 : int pdg = abs(hep.idhep[i]);
450 0 : if(pdg==11 || pdg==13 || pdg==15)
451 : {
452 0 : if(firstPDG==0) firstPDG=hep.idhep[i];
453 : else
454 : {
455 : secondPDG=hep.idhep[i];
456 : // Just in case two pairs are genereted - verify that we have a pair with oposite signs
457 0 : if(firstPDG*secondPDG>0) secondPDG=0;
458 0 : break;
459 : }
460 0 : }
461 0 : }
462 :
463 0 : if( ME_channel==0 && firstPDG!=0 && secondPDG!=0 &&
464 0 : firstPDG==-secondPDG ) ME_channel=1;
465 :
466 : // W: check for pairs 11 -12; -11 12; 13 -14; -13 14; 15 -16; -15 16
467 : // -----------------------------------------------------------------
468 : firstPDG =0;
469 : secondPDG=0;
470 0 : for(int i=firstDaughter; i<hep.nhep;i++)
471 : {
472 0 : int pdg = abs(hep.idhep[i]);
473 0 : if(pdg>=11 && pdg<=16)
474 : {
475 0 : if(firstPDG==0) firstPDG=hep.idhep[i];
476 : else
477 : {
478 : secondPDG=hep.idhep[i];
479 : // Just in case two pairs are genereted - verify that we have a pair with oposite signs
480 0 : if(firstPDG*secondPDG>0) secondPDG=0;
481 0 : break;
482 : }
483 0 : }
484 0 : }
485 :
486 0 : firstPDG =abs(firstPDG);
487 0 : secondPDG=abs(secondPDG);
488 :
489 0 : if( ME_channel==0 && firstPDG!=0 && secondPDG!=0 &&
490 0 : ( ( firstPDG==11 && secondPDG==12 ) || (firstPDG == 12 && secondPDG == 11) ||
491 0 : ( firstPDG==13 && secondPDG==14 ) || (firstPDG == 14 && secondPDG == 13) ||
492 0 : ( firstPDG==15 && secondPDG==16 ) || (firstPDG == 16 && secondPDG == 15)
493 : )
494 0 : ) ME_channel=2;
495 :
496 0 : Log::Debug(901)<<"ME_channel: Found ME_channel: "<<ME_channel<<endl;
497 :
498 : // Check intermediate particle (if exists):
499 :
500 : // Verify that intermediate particle PDG matches ME_channel found
501 0 : if(ME_channel>0 && decay_idx)
502 : {
503 0 : int pdg=hep.idhep[decay_idx-1];
504 :
505 0 : if(ME_channel==1 && !(pdg==22 || pdg==23) ) ME_channel=0; //gamma/Z
506 0 : if(ME_channel==2 && !(pdg==24 || pdg==-24)) ME_channel=0; //W+/W-
507 :
508 0 : if(ME_channel==0)
509 0 : Log::Debug(901,false)<<" but set to 0: wrong intermediate particle: "<<pdg<<endl;
510 0 : }
511 :
512 : // Check flags
513 :
514 0 : switch(ME_channel)
515 : {
516 : case 0: break; // Ok - no channel found
517 0 : case 1: if(!Photos::meCorrectionWtForZ) ME_channel=0; break;
518 0 : case 2: if(!Photos::meCorrectionWtForW) ME_channel=0; break;
519 0 : default: Log::Error()<<"Unimplemented ME channel: "<<ME_channel<<endl; break;
520 : }
521 0 : Log::Debug(902)<<"ME_channel: Finally, after flag check, ME_channel is: "<<ME_channel<<endl;
522 0 : }
523 :
524 :
525 : } // namespace Photospp
|
