pairs.cxx
1#include "Photos.h"
2#include "pairs.h"
3
4#include <cmath>
5#include <stdio.h>
6#include <stdlib.h>
7using namespace Photospp;
8
9namespace Photospp {
10
11//
12 inline double xlam(double A,double B,double C){return sqrt((A-B-C)*(A-B-C)-4.0*B*C);}
13
14 inline double max(double a, double b)
15 {
16 return (a > b) ? a : b;
17 }
18 //
19 //extern "C" void varran_( double RRR[], int *N);
20 double angfi(double X,double Y){
21 double THE;
22 const double PI=3.141592653589793238462643;
23
24 if(X==0.0 && Y==0.0) return 0.0;
25 if(fabs(Y) < fabs(X)){
26 THE=atan(fabs(Y/X));
27 if(X<=0.0) THE=PI-THE;}
28 else{
29 THE=acos(X/sqrt(X*X +Y*Y));
30 }
31 if(Y<0.0) THE=2*PI-THE;
32 return THE;
33}
34
35 double angxy(double X,double Y){
36 double THE;
37 const double PI=3.141592653589793238462643;
38
39 if(X==0.0 && Y==0.0) return 0.0;
40
41 if(fabs(Y) < fabs(X)){
42 THE=atan(fabs(Y/X));
43 if(X<=0.0) THE=PI-THE;}
44 else{
45 THE=acos(X/sqrt(X*X +Y*Y));
46 }
47 return THE;
48}
49
50void bostd3(double EXE,double PVEC[4],double QVEC[4]){
51 // ----------------------------------------------------------------------
52 // BOOST ALONG Z AXIS, EXE=EXP(ETA), ETA= HIPERBOLIC VELOCITY.
53 //
54 // USED BY : KORALZ RADKOR
55 /// ----------------------------------------------------------------------
56 int j=1; // convention of indices of Riemann space must be preserved.
57 double RPL,RMI,QPL,QMI;
58 double RVEC[4];
59
60
61 RVEC[1-j]=PVEC[1-j];
62 RVEC[2-j]=PVEC[2-j];
63 RVEC[3-j]=PVEC[3-j];
64 RVEC[4-j]=PVEC[4-j];
65 RPL=RVEC[4-j]+RVEC[3-j];
66 RMI=RVEC[4-j]-RVEC[3-j];
67 QPL=RPL*EXE;
68 QMI=RMI/EXE;
69 QVEC[1-j]=RVEC[1-j];
70 QVEC[2-j]=RVEC[2-j];
71 QVEC[3-j]=(QPL-QMI)/2;
72 QVEC[4-j]=(QPL+QMI)/2;
73}
74
75// after investigations PHORO3 of PhotosUtilities.cxx will be used instead
76// but it must be checked first if it works
77
78void rotod3(double ANGLE,double PVEC[4],double QVEC[4]){
79
80
81 int j=1; // convention of indices of Riemann space must be preserved.
82 double CS,SN;
83 // printf ("%5.2f\n",cos(ANGLE));
84 CS=cos(ANGLE)*PVEC[1-j]-sin(ANGLE)*PVEC[2-j];
85 SN=sin(ANGLE)*PVEC[1-j]+cos(ANGLE)*PVEC[2-j];
86
87 QVEC[1-j]=CS;
88 QVEC[2-j]=SN;
89 QVEC[3-j]=PVEC[3-j];
90 QVEC[4-j]=PVEC[4-j];
91}
92
93
94
95void rotod2(double PHI,double PVEC[4],double QVEC[4]){
96
97 double RVEC[4];
98 int j=1; // convention of indices of Riemann space must be preserved.
99 double CS,SN;
100
101 CS=cos(PHI);
102 SN=sin(PHI);
103
104 RVEC[1-j]=PVEC[1-j];
105 RVEC[2-j]=PVEC[2-j];
106 RVEC[3-j]=PVEC[3-j];
107 RVEC[4-j]=PVEC[4-j];
108
109 QVEC[1-j]= CS*RVEC[1-j]+SN*RVEC[3-j];
110 QVEC[2-j]=RVEC[2-j];
111 QVEC[3-j]=-SN*RVEC[1-j]+CS*RVEC[3-j];
112 QVEC[4-j]=RVEC[4-j];
113 // printf ("%15.12f %15.12f %15.12f %15.12f \n",QVEC[0],QVEC[1],QVEC[2],QVEC[3]);
114 // exit(-1);
115}
116
117void lortra(int KEY,double PRM,double PNEUTR[4],double PNU[4],double PAA[4],double PP[4],double PE[4]){
118 // ---------------------------------------------------------------------
119 // THIS ROUTINE PERFORMS LORENTZ TRANSFORMATION ON MANY 4-VECTORS
120 // KEY =1 BOOST ALONG 3RD AXIS
121 // =2 ROTATION AROUND 2ND AXIS
122 // =3 ROTATION AROUND 3RD AXIS
123 // PRM TRANSFORMATION PARAMETER - ANGLE OR EXP(HIPERANGLE).
124 //
125 // called by : RADCOR
126 // ---------------------------------------------------------------------
127 if(KEY==1){
128 bostd3(PRM,PNEUTR,PNEUTR);
129 bostd3(PRM,PNU ,PNU );
130 bostd3(PRM,PAA ,PAA );
131 bostd3(PRM,PE ,PE );
132 bostd3(PRM,PP ,PP );
133 }
134 else if (KEY==2){
135 rotod2(PRM,PNEUTR,PNEUTR);
136 rotod2(PRM,PNU ,PNU );
137 rotod2(PRM,PAA ,PAA );
138 rotod2(PRM,PE ,PE );
139 rotod2(PRM,PP ,PP );
140 }
141 else if (KEY==3){
142 rotod3(PRM,PNEUTR,PNEUTR);
143 rotod3(PRM,PNU ,PNU );
144 rotod3(PRM,PAA ,PAA );
145 rotod3(PRM,PE ,PE );
146 rotod3(PRM,PP ,PP );
147 }
148 else{
149 printf (" STOP IN LOTRA. WRONG KEYTRA");
150 exit(-1);
151 }
152}
153
154double amast(double VEC[4]){
155 int i=1; // convention of indices of Riemann space must be preserved
156 double ama= VEC[4-i]*VEC[4-i]-VEC[1-i]*VEC[1-i]-VEC[2-i]*VEC[2-i]-VEC[3-i]*VEC[3-i];
157 ama=sqrt(fabs(ama));
158 return ama;
159}
160
161void spaj(int KUDA,double P2[4],double Q2[4],double PP[4],double PE[4]){
162 // **********************
163 // THIS PRINTS OUT FOUR MOMENTA OF PHOTONS
164 // ON OUTPUT UNIT NOUT
165
166 double SUM[4];
167 const int KLUCZ=0;
168 if (KLUCZ==0) return;
169
170 printf (" %10i =====================SPAJ==================== \n", KUDA);
171 printf (" P2 %18.13f %18.13f %18.13f %18.13f \n",P2[0],P2[1],P2[2],P2[3]);
172 printf (" Q2 %18.13f %18.13f %18.13f %18.13f \n",Q2[0],Q2[1],Q2[2],Q2[3]);
173 printf (" PE %18.13f %18.13f %18.13f %18.13f \n",PE[0],PE[1],PE[2],PE[3]);
174 printf (" PP %18.13f %18.13f %18.13f %18.13f \n",PP[0],PP[1],PP[2],PP[3]);
175
176 for( int k=0;k<=3;k++) SUM[k]=P2[k]+Q2[k]+PE[k]+PP[k];
177
178 printf ("SUM %18.13f %18.13f %18.13f %18.13f \n",SUM[0],SUM[1],SUM[2],SUM[3]);
179}
180
181//extern "C" {
182 struct PARKIN {
183 double fi0; // FI0
184 double fi1; // FI1
185 double fi2; // FI2
186 double fi3; // FI3
187 double fi4; // FI4
188 double fi5; // FI5
189 double th0; // TH0
190 double th1; // TH1
191 double th3; // TH3
192 double th4; // TH4
193 double parneu; // PARNEU
194 double parch; // PARCH
195 double bpar; // BPAR
196 double bsta; // BSTA
197 double bstb; // BSTB
198 } parkin;
199//}
200
201//struct PARKIN parkin;
202
203void partra(int IBRAN,double PHOT[4]){
204
205
206
207 rotod3(-parkin.fi0,PHOT,PHOT);
208 rotod2(-parkin.th0,PHOT,PHOT);
209 bostd3(parkin.bsta,PHOT,PHOT);
210 rotod3(-parkin.fi1,PHOT,PHOT);
211 rotod2(-parkin.th1,PHOT,PHOT);
212 rotod3(-parkin.fi2,PHOT,PHOT);
213
214 if(IBRAN==-1){
215 bostd3(parkin.parneu,PHOT,PHOT);
216 }
217 else{
218 bostd3(parkin.parch,PHOT,PHOT);
219 }
220
221 rotod3(-parkin.fi3,PHOT,PHOT);
222 rotod2(-parkin.th3,PHOT,PHOT);
223 bostd3(parkin.bpar,PHOT,PHOT);
224 rotod3( parkin.fi4,PHOT,PHOT);
225 rotod2(-parkin.th4,PHOT,PHOT);
226 rotod3(-parkin.fi5,PHOT,PHOT);
227 rotod3( parkin.fi2,PHOT,PHOT);
228 rotod2( parkin.th1,PHOT,PHOT);
229 rotod3( parkin.fi1,PHOT,PHOT);
230 bostd3(parkin.bstb,PHOT,PHOT);
231 rotod2( parkin.th0,PHOT,PHOT);
232 rotod3( parkin.fi0,PHOT,PHOT);
233
234}
235
236
237 void trypar(bool *JESLI,double *pSTRENG,double AMCH, double AMEL, double PA[4],double PB[4],double PE[4],double PP[4],bool *sameflav){
238 double &STRENG = *pSTRENG;
239 // COMMON /PARKIN/
240 double &FI0=parkin.fi0;
241 double &FI1=parkin.fi1;
242 double &FI2=parkin.fi2;
243 double &FI3=parkin.fi3;
244 double &FI4=parkin.fi4;
245 double &FI5=parkin.fi5;
246 double &TH0=parkin.th0;
247 double &TH1=parkin.th1;
248 double &TH3=parkin.th3;
249 double &TH4=parkin.th4;
250 double &PARNEU=parkin.parneu;
251 double &PARCH=parkin.parch;
252 double &BPAR=parkin.bpar;
253 double &BSTA=parkin.bsta;
254 double &BSTB=parkin.bstb;
255
256 double PNEUTR[4],PAA[4],PHOT[4],PSUM[4];
257 double VEC[4];
258 double RRR[8];
259 bool JESLIK;
260 const double PI=3.141592653589793238462643;
261 const double ALFINV= 137.01;
262 const int j=1; // convention of indices of Riemann space must be preserved.
263
264 PA[4-j]=max(PA[4-j],sqrt(PA[1-j]*PA[1-j]+PA[2-j]*PA[2-j]+PA[3-j]*PA[3-j]));
265 PB[4-j]=max(PB[4-j],sqrt(PB[1-j]*PB[1-j]+PB[2-j]*PB[2-j]+PB[3-j]*PB[3-j]));
266
267 // 4-MOMENTUM OF THE NEUTRAL SYSTEM
268 for( int k=0;k<=3;k++){
269 PE[k] =0.0;
270 PP[k] =0.0;
271 PSUM[k] =PA[k]+PB[k];
272 PAA[k] =PA[k];
273 PNEUTR[k]=PB[k];
274 }
275 if((PAA[4-j]+PNEUTR[4-j])< 0.01 ){
276printf (" too small energy to emit %10.7f\n",PAA[4-j]+PNEUTR[4-j]);
277 *JESLI=false;
278 return;
279 }
280
281 // MASSES OF THE NEUTRAL AND CHARGED SYSTEMS AND OVERALL MASS
282 // FIRST WE HAVE TO GO TO THE RESTFRAME TO GET RID OF INSTABILITIES
283 // FROM BHLUMI OR ANYTHING ELSE
284 // THIRD AXIS ALONG PNEUTR
285 double X1 = PSUM[1-j];
286 double X2 = PSUM[2-j];
287 FI0 =angfi(X1,X2);
288 X1 = PSUM[3-j];
289 X2 = sqrt(PSUM[1-j]*PSUM[1-j]+PSUM[2-j]*PSUM[2-j]);
290 TH0 =angxy(X1,X2) ;
291 spaj(-2,PNEUTR,PAA,PP,PE);
292 lortra(3,-FI0,PNEUTR,VEC,PAA,PP,PE);
293 lortra(2,-TH0,PNEUTR,VEC,PAA,PP,PE);
294 rotod3(-FI0,PSUM,PSUM);
295 rotod2(-TH0,PSUM,PSUM);
296 BSTA=(PSUM[4-j]-PSUM[3-j])/sqrt(PSUM[4-j]*PSUM[4-j]-PSUM[3-j]*PSUM[3-j]);
297 BSTB=(PSUM[4-j]+PSUM[3-j])/sqrt(PSUM[4-j]*PSUM[4-j]-PSUM[3-j]*PSUM[3-j]);
298 lortra(1,BSTA,PNEUTR,VEC,PAA,PP,PE);
299 spaj(-1,PNEUTR,PAA,PP,PE);
300 double AMNE=amast(PNEUTR);
301 AMCH=amast(PAA); // to be improved. May be dangerous because of rounding error
302 if(AMCH<0.0) AMCH=AMEL;
303 if (AMNE<0.0) AMNE=0.0;
304 double AMTO =PAA[4-j]+PNEUTR[4-j];
305
306
307 for( int k=0;k<=7;k++) RRR[k]=Photos::randomDouble();
308
309 if(STRENG==0.0) {*JESLI=false; return;}
310
311 double PRHARD;
312 PRHARD= STRENG // NOTE: logs from phase space presamplers not MEs
313 *0.5*(1.0/PI/ALFINV)*(1.0/PI/ALFINV) // normalization of triple log 1/36 from
314 // journals.aps.org/prd/pdf/10.1103/PhysRevD.49.1178
315 // must come from rejection
316 // 0.5 is because it is for 1-leg only
317 // STRENG=0,5 because of calls before and after photons
318 // other logs should come from rejection
319 *2*log(AMTO/AMEL/2.0) // virtuality
320 *log(AMTO/AMEL/2.0) // soft
321 *log((AMTO*AMTO+2*AMCH*AMCH)/2.0/AMCH/AMCH);// collinear
322 double FREJECT=2.; // to make room for interference second pair posiblty.
323 PRHARD=PRHARD*FREJECT;
324 // PRHARD=PRHARD*50; // to increase number of pairs in test of mu mu from mu
325 // fror mumuee set *15
326 // enforces hard pairs to be generated 'always'
327 // for the sake of tests with high statistics, also for flat phase space.
328 // PRHARD=0.99* STRENG*2;
329 // STRENG=0.0;
330 if (PRHARD>1.0) {
331 printf(" stop from Photos pairs.cxx PRHARD= %18.13f \n",PRHARD);
332 exit(0);
333 }
334 // delta is for tests with PhysRevD.49.1178, default is AMTO*2 no restriction on pair phase space
335 double delta=AMTO*2; //5;//.125; //AMTO*2; //.125; //AMTO*2; ;0.25;
336
337
338 if (RRR[7-j]>PRHARD){ // compensate crude probablilities; for pairs from consecutive sources
339 STRENG=STRENG/(1.0-PRHARD);
340 *JESLI=false;
341 return;
342 }
343 else STRENG=0.0;
344
345
346
347
348
349 //
350
351 //virtuality of lepton pair
352 double XMP=2.0*AMEL*exp(RRR[1-j]*log(AMTO/2.0/AMEL));
353 // XMP=2.0*AMEL*2.0*AMEL+RRR[1-j]*(AMTO-2.0*AMEL)*(AMTO-2.0*AMEL); XMP=sqrt(XMP); // option of no presampler
354
355
356 // energy of lepton pair replace virtuality of CHAR+NEU system in phase space parametrization
357 double XPmin=2.0*AMEL;
358 double XPdelta=AMTO-XPmin;
359 double XP= XPmin*exp(RRR[2-j]*log((XPdelta+XPmin)/XPmin));
360 // XP= XPmin +RRR[2-j]*XPdelta; // option of no presampler
361 double XMK2=(AMTO*AMTO+XMP*XMP)-2.0*AMTO*XP;
362
363 // angles of lepton pair
364 double eps=4.0*AMCH*AMCH/AMTO/AMTO;
365 double C1 =1.0+eps -eps*exp(RRR[3-j]*log((2+eps)/eps));
366 // C1=1.0-2.0*RRR[3-j]; // option of no presampler
367 double FIX1=2.0*PI*RRR[4-j];
368
369 // angles of lepton in restframe of lepton pair
370 double C2 =1.0-2.0*RRR[5-j];
371 double FIX2=2.0*PI*RRR[6-j];
372
373
374
375 // histograming .......................
376 JESLIK= (XP <((AMTO*AMTO+XMP*XMP-(AMCH+AMNE)*(AMCH+AMNE))/2.0/AMTO));
377 double WTA=0.0;
378 WTA=WTA*5;
379 if(JESLIK) WTA=1.0;
380 // GMONIT( 0,101 ,WTA,1D0,0D0)
381 JESLIK= (XMP<(AMTO-AMNE-AMCH));
382 WTA=0.0;
383 if(JESLIK) WTA=1.0;
384 // GMONIT( 0,102 ,WTA,1D0,0D0)
385 JESLIK= (XMP<(AMTO-AMNE-AMCH))&& (XP >XMP);
386
387 WTA=0.0;
388 if(JESLIK) WTA=1.0;
389 // GMONIT( 0,103 ,WTA,1D0,0D0)
390 JESLIK=
391 (XMP<(AMTO-AMNE-AMCH))&&
392 (XP >XMP) &&
393 (XP <((AMTO*AMTO+XMP*XMP-(AMCH+AMNE)*(AMCH+AMNE))/2.0/AMTO));
394 WTA=0.0;
395 if (JESLIK) WTA=1.0;
396 // GMONIT( 0,104 ,WTA,1D0,0D0)
397 // end of histograming ................
398
399 // phase space limits rejection variable
400 *JESLI=(RRR[7-j]<PRHARD) &&
401 (XMP<(AMTO-AMNE-AMCH)) &&
402 (XP >XMP) &&
403 (XP <((AMTO*AMTO+XMP*XMP-(AMCH+AMNE)*(AMCH+AMNE))/2.0/AMTO));
404
405
406 // rejection for phase space restriction: for tests with PhysRevD.49.1178
407 *JESLI= *JESLI && XP< delta;
408 if (!*JESLI) return;
409
410 // for events in phase: jacobians weights etc.
411
412 // virtuality of added lepton pair
413 double F= (AMTO*AMTO-4.0*AMEL*AMEL) // flat phase space
414 /(AMTO*AMTO-4.0*AMEL*AMEL) *XMP*XMP; // use this when presampler is on (log moved to PRHARD)
415
416
417 // Energy of added lepton pair represented by virtuality of CH+N pair
418 double G= 2*AMTO*XPdelta // flat phase space
419 /(2*AMTO*XPdelta) *2*AMTO*XP; // use this when presampler is on (log moved to PRHARD)
420
421
422 // scattering angle of emitted lepton pair (also flat factors for other angles)
423 double H= 2.0 // flat phase space
424 /2.0 *(1.0+eps-C1)/2.0; // use this when presampler is on (log moved to PRHARD)
425
426 double H1=2.0 // for other generation arm: char neutr replaced
427 /2.0 *(1.0+eps-C1);
428 double H2=2.0
429 /2.0 *(1.0+eps+C1);
430 H=1./(0.5/H1+0.5/H2);
431
432 //*2*PI*4*PI /2/PI/4/PI; // other angles normalization of transformation to random numbers.
433
434 double XPMAX=(AMTO*AMTO+XMP*XMP-(AMCH+AMNE)*(AMCH+AMNE))/2.0/AMTO;
435
436 double YOT3=F*G*H; // jacobians for phase space variables
437 double YOT2= // lambda factors:
438 xlam(1.0,AMEL*AMEL/XMP/XMP, AMEL*AMEL/XMP/XMP)*
439 xlam(1.0,XMK2/AMTO/AMTO,XMP*XMP/AMTO/AMTO)*
440 xlam(1.0,AMCH*AMCH/XMK2,AMNE*AMNE/XMK2)
441 / xlam(1.0,AMCH*AMCH/AMTO/AMTO,AMNE*AMNE/AMTO/AMTO);
442
443
444 //C histograming .......................
445 // GMONIT( 0,105 ,WT ,1D0,0D0)
446 // GMONIT( 0,106 ,YOT1,1D0,0D0)
447 // GMONIT( 0,107 ,YOT2,1D0,0D0)
448 // GMONIT( 0,108 ,YOT3,1D0,0D0)
449 // GMONIT( 0,109 ,YOT4,1D0,0D0)
450 // end of histograming ................
451
452
453 //
454 //
455 // FRAGMENTATION COMES !!
456 //
457 // THIRD AXIS ALONG PNEUTR
458 X1 = PNEUTR[1-j];
459 X2 = PNEUTR[2-j];
460 FI1 =angfi(X1,X2);
461 X1 = PNEUTR[3-j];
462 X2 = sqrt(PNEUTR[1-j]*PNEUTR[1-j]+PNEUTR[2-j]*PNEUTR[2-j]) ;
463 TH1 =angxy(X1,X2);
464 spaj(0,PNEUTR,PAA,PP,PE);
465 lortra(3,-FI1,PNEUTR,VEC,PAA,PP,PE);
466 lortra(2,-TH1,PNEUTR,VEC,PAA,PP,PE);
467 VEC[4-j]=0.0;
468 VEC[3-j]=0.0;
469 VEC[2-j]=0.0;
470 VEC[1-j]=1.0;
471 FI2=angfi(VEC[1-j],VEC[2-j]);
472 lortra(3,-FI2,PNEUTR,VEC,PAA,PP,PE);
473 spaj(1,PNEUTR,PAA,PP,PE);
474
475 // STEALING FROM PAA AND PNEUTR ENERGY FOR THE pair
476 // ====================================================
477 // NEW MOMENTUM OF PAA AND PNEUTR (IN THEIR VERY REST FRAME)
478 // 1) PARAMETERS.....
479 double AMCH2=AMCH*AMCH;
480 double AMNE2=AMNE*AMNE;
481 double AMTOST=XMK2;
482 double QNEW=xlam(AMTOST,AMNE2,AMCH2)/sqrt(AMTOST)/2.0;
483 double QOLD=PNEUTR[3-j];
484 double GCHAR=(QNEW*QNEW+QOLD*QOLD+AMCH*AMCH)/
485 (QNEW*QOLD+sqrt((QNEW*QNEW+AMCH*AMCH)*(QOLD*QOLD+AMCH*AMCH)));
486 double GNEU= (QNEW*QNEW+QOLD*QOLD+AMNE*AMNE)/
487 (QNEW*QOLD+sqrt((QNEW*QNEW+AMNE*AMNE)*(QOLD*QOLD+AMNE*AMNE)));
488
489 // GCHAR=(QOLD**2-QNEW**2)/(
490 // & QNEW*SQRT(QOLD**2+AMCH2)+QOLD*SQRT(QNEW**2+AMCH2)
491 // & )
492 // GCHAR=SQRT(1D0+GCHAR**2)
493 // GNEU=(QOLD**2-QNEW**2)/(
494 // & QNEW*SQRT(QOLD**2+AMNE2)+QOLD*SQRT(QNEW**2+AMNE2)
495 // & )
496 // GNEU=SQRT(1D0+GNEU**2)
497 if(GNEU<1.||GCHAR<1.){
498 printf(" PHOTOS TRYPAR GBOOST LT 1., LIMIT OF PHASE SPACE %18.13f %18.13f %18.13f %18.13f %18.13f %18.13f %18.13f %18.13f \n"
499 ,GNEU,GCHAR,QNEW,QOLD,AMTO,AMTOST,AMNE,AMCH);
500 return;
501 }
502 PARCH =GCHAR+sqrt(GCHAR*GCHAR-1.0);
503 PARNEU=GNEU -sqrt(GNEU*GNEU -1.0);
504
505 // 2) REDUCTIEV BOOSTS
506 bostd3(PARNEU,VEC ,VEC );
507 bostd3(PARNEU,PNEUTR,PNEUTR);
508 bostd3(PARCH,PAA ,PAA );
509 spaj(2,PNEUTR,PAA,PP,PE);
510
511 // TIME FOR THE PHOTON that is electron pair
512 double PMOD=xlam(XMP*XMP,AMEL*AMEL,AMEL*AMEL)/XMP/2.0;
513 double S2=sqrt(1.0-C2*C2);
514 PP[4-j]=XMP/2.0;
515 PP[3-j]=PMOD*C2;
516 PP[2-j]=PMOD*S2*cos(FIX2);
517 PP[1-j]=PMOD*S2*sin(FIX2);
518 PE[4-j]= PP[4-j];
519 PE[3-j]=-PP[3-j];
520 PE[2-j]=-PP[2-j];
521 PE[1-j]=-PP[1-j];
522 // PHOTON ENERGY and momentum IN THE REDUCED SYSTEM
523 double PENE=(AMTO*AMTO-XMP*XMP-XMK2)/2.0/sqrt(XMK2);
524 double PPED=sqrt(PENE*PENE-XMP*XMP);
525 FI3=FIX1;
526 double COSTHG=C1;
527 double SINTHG=sqrt(1.0-C1*C1);
528 X1 = -COSTHG;
529 X2 = SINTHG;
530 TH3 =angxy(X1,X2);
531 PHOT[1-j]=PMOD*SINTHG*cos(FI3);
532 PHOT[2-j]=PMOD*SINTHG*sin(FI3);
533 // MINUS BECAUSE AXIS OPPOSITE TO PAA
534 PHOT[3-j]=-PMOD*COSTHG;
535 PHOT[4-j]=PENE;
536 // ROTATE TO PUT PHOTON ALONG THIRD AXIS
537 X1 = PHOT[1-j];
538 X2 = PHOT[2-j];
539 lortra(3,-FI3,PNEUTR,VEC,PAA,PP,PE);
540 rotod3(-FI3,PHOT,PHOT);
541 lortra(2,-TH3,PNEUTR,VEC,PAA,PP,PE);
542 rotod2(-TH3,PHOT,PHOT);
543 spaj(21,PNEUTR,PAA,PP,PE);
544 // ... now get the pair !
545 double PAIRB=PENE/XMP+PPED/XMP;
546 bostd3(PAIRB,PE,PE);
547 bostd3(PAIRB,PP,PP);
548 spaj(3,PNEUTR,PAA,PP,PE);
549 double GAMM=(PNEUTR[4-j]+PAA[4-j]+PP[4-j]+PE[4-j])/AMTO;
550
551 // TP and ZW: 25.II.2015: fix for cases when mother is very close
552 // to its rest frame and pair is generated after photon emission.
553 // Then GAMM can be slightly less than 1.0 due to rounding error
554 if( GAMM < 1.0 ) {
555 if( GAMM > 0.9999999 ) GAMM = 1.0;
556 else {
557 printf("Photos::partra: GAMM = %20.18e\n",GAMM);
558 printf(" BELOW 0.9999999 THRESHOLD!\n");
559 GAMM = 1.0;
560 }
561 }
562
563 BPAR=GAMM-sqrt(GAMM*GAMM-1.0);
564 lortra(1, BPAR,PNEUTR,VEC,PAA,PP,PE);
565 bostd3( BPAR,PHOT,PHOT);
566 spaj(4,PNEUTR,PAA,PP,PE);
567 // BACK IN THE TAU REST FRAME BUT PNEUTR NOT YET ORIENTED.
568 X1 = PNEUTR[1-j];
569 X2 = PNEUTR[2-j];
570 FI4 =angfi(X1,X2);
571 X1 = PNEUTR[3-j];
572 X2 = sqrt(PNEUTR[1-j]*PNEUTR[1-j]+PNEUTR[2-j]*PNEUTR[2-j]);
573 TH4 =angxy(X1,X2);
574 lortra(3, FI4,PNEUTR,VEC,PAA,PP,PE);
575 rotod3( FI4,PHOT,PHOT);
576 lortra(2,-TH4,PNEUTR,VEC,PAA,PP,PE);
577 rotod2(-TH4,PHOT,PHOT);
578 X1 = VEC[1-j];
579 X2 = VEC[2-j];
580 FI5=angfi(X1,X2);
581 lortra(3,-FI5,PNEUTR,VEC,PAA,PP,PE);
582 rotod3(-FI5,PHOT,PHOT);
583 // PAA RESTORES ORIGINAL DIRECTION
584 lortra(3, FI2,PNEUTR,VEC,PAA,PP,PE);
585 rotod3( FI2,PHOT,PHOT);
586 lortra(2, TH1,PNEUTR,VEC,PAA,PP,PE);
587 rotod2( TH1,PHOT,PHOT);
588 lortra(3, FI1,PNEUTR,VEC,PAA,PP,PE);
589 rotod3( FI1,PHOT,PHOT);
590 spaj(10,PNEUTR,PAA,PP,PE);
591 lortra(1,BSTB,PNEUTR,VEC,PAA,PP,PE);
592 lortra(2,TH0,PNEUTR,VEC,PAA,PP,PE);
593 lortra(3,FI0,PNEUTR,VEC,PAA,PP,PE);
594 spaj(11,PNEUTR,PAA,PP,PE);
595
596
597 // matrix element as formula 1 from journals.aps.org/prd/pdf/10.1103/PhysRevD.49.1178
598
599 double pq= PAA[3]*PP[3]-PAA[2]*PP[2]-PAA[1]*PP[1]-PAA[0]*PP[0];
600 pq=pq +PAA[3]*PE[3]-PAA[2]*PE[2]-PAA[1]*PE[1]-PAA[0]*PE[0];
601
602 double ppq= PNEUTR[3]*PP[3]-PNEUTR[2]*PP[2]-PNEUTR[1]*PP[1]-PNEUTR[0]*PP[0];
603 ppq=ppq+ PNEUTR[3]*PE[3]-PNEUTR[2]*PE[2]-PNEUTR[1]*PE[1]-PNEUTR[0]*PE[0];
604 double pq1 =PAA[3]*PP[3]-PAA[2]*PP[2]-PAA[1]*PP[1]-PAA[0]*PP[0];
605 double pq2 =PAA[3]*PE[3]-PAA[2]*PE[2]-PAA[1]*PE[1]-PAA[0]*PE[0];
606
607 double ppq1=PNEUTR[3]*PP[3]-PNEUTR[2]*PP[2]-PNEUTR[1]*PP[1]-PNEUTR[0]*PP[0];
608 double ppq2=PNEUTR[3]*PE[3]-PNEUTR[2]*PE[2]-PNEUTR[1]*PE[1]-PNEUTR[0]*PE[0];
609
610 double ppp=PNEUTR[3]*PAA[3]-PNEUTR[2]*PAA[2]-PNEUTR[1]*PAA[1]-PNEUTR[0]*PAA[0];
611
612 double YOT1=1./2./XMP/XMP/XMP/XMP*
613 ( 4*(pq1/pq-ppq1/ppq)*(pq2/pq-ppq2/ppq)
614 -XMP*XMP*(AMCH2/pq/pq+AMNE2/ppq/ppq-ppp/pq/ppq-ppp/pq/ppq) );
615 // *(1-XP/XPMAX+0.5*(XP/XPMAX)*(XP/XPMAX)); // A-P kernel divide by (1-XP/XPMAX)?
616 // if (*sameflav){
617 //printf ("samef= %d\n",*sameflav);
618 //exit(0);
619 //}
620 if(*sameflav){
621 // we interchange: PAA <--> pp
622 for( int k=0;k<=3;k++){
623 double stored=PAA[k];
624 PAA[k]=PE[k];
625 PE[k]=stored;
626 }
627
628 pq= PAA[3]*PP[3]-PAA[2]*PP[2]-PAA[1]*PP[1]-PAA[0]*PP[0];
629 pq=pq +PAA[3]*PE[3]-PAA[2]*PE[2]-PAA[1]*PE[1]-PAA[0]*PE[0];
630
631 ppq= PNEUTR[3]*PP[3]-PNEUTR[2]*PP[2]-PNEUTR[1]*PP[1]-PNEUTR[0]*PP[0];
632 ppq=ppq+ PNEUTR[3]*PE[3]-PNEUTR[2]*PE[2]-PNEUTR[1]*PE[1]-PNEUTR[0]*PE[0];
633 pq1 =PAA[3]*PP[3]-PAA[2]*PP[2]-PAA[1]*PP[1]-PAA[0]*PP[0];
634 pq2 =PAA[3]*PE[3]-PAA[2]*PE[2]-PAA[1]*PE[1]-PAA[0]*PE[0];
635
636 ppq1=PNEUTR[3]*PP[3]-PNEUTR[2]*PP[2]-PNEUTR[1]*PP[1]-PNEUTR[0]*PP[0];
637 ppq2=PNEUTR[3]*PE[3]-PNEUTR[2]*PE[2]-PNEUTR[1]*PE[1]-PNEUTR[0]*PE[0];
638
639 ppp=PNEUTR[3]*PAA[3]-PNEUTR[2]*PAA[2]-PNEUTR[1]*PAA[1]-PNEUTR[0]*PAA[0];
640
641 XMP=-(PP[0]+PE[0])*(PP[0]+PE[0])-(PP[1]+PE[1])*(PP[1]+PE[1])
642 -(PP[2]+PE[2])*(PP[2]+PE[2])+(PP[3]+PE[3])*(PP[3]+PE[3]);
643 XMP=sqrt(fabs(XMP));
644
645
646double YOT1p=1./2./XMP/XMP/XMP/XMP*
647 ( 4*(pq1/pq-ppq1/ppq)*(pq2/pq-ppq2/ppq)
648 -XMP*XMP*(AMCH2/pq/pq+AMNE2/ppq/ppq-ppp/pq/ppq-ppp/pq/ppq) );
649 // *(1-XP/XPMAX+0.5*(XP/XPMAX)*(XP/XPMAX)); // A-P kernel divide by (1-XP/XPMAX)?
650 double wtint=0.; // not yet installed
651 wtint=1;//(YOT1+YOT1p+wtint)/(YOT1+YOT1p);
652 YOT1=YOT1*wtint;
653
654 // we interchange: PAA <--> pp back into place
655 for( int k=0;k<=3;k++){
656 double stored=PAA[k];
657 PAA[k]=PE[k];
658 PE[k]=stored;
659 }
660 } // end sameflav
661
662 double WT=YOT1*YOT2*YOT3;
663
664 WT=WT/8/FREJECT; // origin must be understood
665
666 if(WT>1.0){
667 printf (" from Photos pairs.cxx WT= %15.8f \n",WT);
668
669}
670
671 if (RRR[8-j]>WT){
672 *JESLI=false;
673 return;
674 }
675
676 }
677
678} // namespace Photospp
679
static double(* randomDouble)()