StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
pairs.cxx
1 #include "Photos.h"
2 #include "pairs.h"
3 
4 #include <cmath>
5 #include <stdio.h>
6 #include <stdlib.h>
7 using namespace Photospp;
8 
9 namespace 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 
50 void 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
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 
78 void 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 
95 void 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 
117 void 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 
154 double 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 
161 void 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 
203 void 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 ){
276 printf (" 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 
646 double 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 
Photospp::Photos::randomDouble
static double(* randomDouble)()
Definition: Photos.h:228
Generated by   doxygen 1.8.5