Add Boost() method to boost all particles to LHC lab frame. Needed for asymmetric...
[u/mrichter/AliRoot.git] / EVGEN / AliGenHijing.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 Revision 1.30  2001/10/21 18:35:56  hristov
19 Several pointers were set to zero in the default constructors to avoid memory management problems
20
21 Revision 1.29  2001/10/15 08:12:24  morsch
22 - Vertex smearing with truncated gaussian.
23 - Store triggered jet info before and after final state radiation into mc-heade
24
25 Revision 1.28  2001/10/08 11:55:25  morsch
26 Store 4-momenta of trigegred jets in event header.
27 Possibility to switch of initial and final state radiation.
28
29 Revision 1.27  2001/10/08 07:13:14  morsch
30 Add setter for minimum transverse momentum of triggered jet.
31
32 Revision 1.26  2001/10/04 08:12:24  morsch
33 Redefinition of stable condition.
34
35 Revision 1.25  2001/07/27 17:09:36  morsch
36 Use local SetTrack, KeepTrack and SetHighWaterMark methods
37 to delegate either to local stack or to stack owned by AliRun.
38 (Piotr Skowronski, A.M.)
39
40 Revision 1.24  2001/07/20 09:34:56  morsch
41 Count the number of spectator neutrons and protons and add information
42 to the event header. (Chiara Oppedisano)
43
44 Revision 1.23  2001/07/13 17:30:22  morsch
45 Derive from AliGenMC.
46
47 Revision 1.22  2001/06/11 13:09:23  morsch
48 - Store cross-Section and number of binary collisions as a function of impact parameter
49 - Pass AliGenHijingEventHeader to gAlice.
50
51 Revision 1.21  2001/02/28 17:35:24  morsch
52 Consider elastic interactions (ks = 1 and ks = 11) as spectator (Chiara Oppedisano)
53
54 Revision 1.20  2001/02/14 15:50:40  hristov
55 The last particle in event marked using SetHighWaterMark
56
57 Revision 1.19  2000/12/21 16:24:06  morsch
58 Coding convention clean-up
59
60 Revision 1.18  2000/12/06 17:46:30  morsch
61 Avoid random numbers 1 and 0.
62
63 Revision 1.17  2000/12/04 11:22:03  morsch
64 Init of sRandom as in 1.15
65
66 Revision 1.16  2000/12/02 11:41:39  morsch
67 Use SetRandom() to initialize random number generator in constructor.
68
69 Revision 1.15  2000/11/30 20:29:02  morsch
70 Initialise static variable sRandom in constructor: sRandom = fRandom;
71
72 Revision 1.14  2000/11/30 07:12:50  alibrary
73 Introducing new Rndm and QA classes
74
75 Revision 1.13  2000/11/09 17:40:27  morsch
76 Possibility to select/unselect spectator protons and neutrons.
77 Method SetSpectators(Int_t spect) added. (FCA, Ch. Oppedisano)
78
79 Revision 1.12  2000/10/20 13:38:38  morsch
80 Debug printouts commented.
81
82 Revision 1.11  2000/10/20 13:22:26  morsch
83 - skip particle type 92 (string)
84 - Charmed and beauty baryions (5122, 4122) are considered as stable consistent with
85   mesons.
86
87 Revision 1.10  2000/10/17 15:10:20  morsch
88 Write first all the parent particles to the stack and then the final state particles.
89
90 Revision 1.9  2000/10/17 13:38:59  morsch
91 Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..)     (FCA)
92
93 Revision 1.8  2000/10/17 12:46:31  morsch
94 Protect EvaluateCrossSections() against division by zero.
95
96 Revision 1.7  2000/10/02 21:28:06  fca
97 Removal of useless dependecies via forward declarations
98
99 Revision 1.6  2000/09/11 13:23:37  morsch
100 Write last seed to file (fortran lun 50) and reed back from same lun using calls to
101 luget_hijing and luset_hijing.
102
103 Revision 1.5  2000/09/07 16:55:40  morsch
104 fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
105
106 Revision 1.4  2000/07/11 18:24:56  fca
107 Coding convention corrections + few minor bug fixes
108
109 Revision 1.3  2000/06/30 12:08:36  morsch
110 In member data: char* replaced by TString, Init takes care of resizing the strings to
111 8 characters required by Hijing.
112
113 Revision 1.2  2000/06/15 14:15:05  morsch
114 Add possibility for heavy flavor selection: charm and beauty.
115
116 Revision 1.1  2000/06/09 20:47:27  morsch
117 AliGenerator interface class to HIJING using THijing (test version)
118
119 */
120
121
122
123 // Generator using HIJING as an external generator
124 // The main HIJING options are accessable for the user through this interface.
125 // Uses the THijing implementation of TGenerator.
126 //
127 // andreas.morsch@cern.ch
128
129 #include "AliGenHijing.h"
130 #include "AliGenHijingEventHeader.h"
131 #include "AliRun.h"
132 #include "AliPDG.h"
133
134 #include <TArrayI.h>
135 #include <TParticle.h>
136 #include <THijing.h>
137 #include <TGraph.h>
138 #include <TLorentzVector.h>
139
140
141  ClassImp(AliGenHijing)
142
143 AliGenHijing::AliGenHijing()
144                  :AliGenMC()
145 {
146 // Constructor
147 }
148
149 AliGenHijing::AliGenHijing(Int_t npart)
150     :AliGenMC(npart)
151 {
152 // Default PbPb collisions at 5. 5 TeV
153 //
154     SetEnergyCMS();
155     SetImpactParameterRange();
156     SetTarget();
157     SetProjectile();
158     SetBoostLHC();
159     fKeep       =  0;
160     fQuench     =  1;
161     fShadowing  =  1;
162     fTrigger    =  0;
163     fDecaysOff  =  1;
164     fEvaluate   =  0;
165     fSelectAll  =  0;
166     fFlavor     =  0;
167     fSpectators =  1;
168     fDsigmaDb   =  0;
169     fDnDb       =  0;
170     fPtMinJet   = -2.5;         
171     fRadiation  =  1;
172     fEventVertex.Set(3);
173 //
174 // Set random number generator   
175     sRandom = fRandom;
176     fHijing = 0;
177
178 }
179
180 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
181 {
182 // copy constructor
183 }
184
185
186 AliGenHijing::~AliGenHijing()
187 {
188 // Destructor
189     if ( fDsigmaDb) delete  fDsigmaDb;  
190     if ( fDnDb)     delete  fDnDb;  
191 }
192
193 void AliGenHijing::Init()
194 {
195 // Initialisation
196     fFrame.Resize(8);
197     fTarget.Resize(8);
198     fProjectile.Resize(8);
199     
200     SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
201                       fAProjectile, fZProjectile, fATarget, fZTarget, 
202                       fMinImpactParam, fMaxImpactParam));
203
204     fHijing=(THijing*) fgMCEvGen;
205     fHijing->SetIHPR2(2,  fRadiation);
206     fHijing->SetIHPR2(3,  fTrigger);
207     fHijing->SetIHPR2(4,  fQuench);
208     fHijing->SetIHPR2(6,  fShadowing);
209     fHijing->SetIHPR2(12, fDecaysOff);    
210     fHijing->SetIHPR2(21, fKeep);
211     fHijing->SetHIPR1(10, fPtMinJet);   
212     fHijing->Initialize();
213
214     
215 //
216     if (fEvaluate) EvaluateCrossSections();
217 //
218 //
219 //  Initialize random generator
220 }
221
222 void AliGenHijing::Generate()
223 {
224 // Generate one event
225
226   Float_t polar[3]    =   {0,0,0};
227   Float_t origin[3]   =   {0,0,0};
228   Float_t origin0[3]  =   {0,0,0};
229   Float_t p[3], random[6];
230   Float_t tof;
231
232   static TClonesArray *particles;
233 //  converts from mm/c to s
234   const Float_t kconv = 0.001/2.999792458e8;
235 //
236   Int_t nt  = 0;
237   Int_t jev = 0;
238   Int_t j, kf, ks, imo;
239   kf = 0;
240     
241   if(!particles) particles = new TClonesArray("TParticle",10000);
242     
243   fTrials = 0;
244   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
245   if(fVertexSmear == kPerEvent) {
246       Float_t dv[3];
247       dv[2] = 1.e10;
248       while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
249           Rndm(random,6);
250           for (j=0; j < 3; j++) {
251               dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
252                   TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
253           }
254       }
255       for (j=0; j < 3; j++) origin0[j] += dv[j];
256   } else if (fVertexSmear == kPerTrack) {
257 //          fHijing->SetMSTP(151,0);
258       for (j = 0; j < 3; j++) {
259 //            fHijing->SetPARP(151+j, fOsigma[j]*10.);
260       }
261   }
262   while(1)
263   {
264 //    Generate one event
265 // --------------------------------------------------------------------------
266       fSpecn    = 0;  
267       fSpecp    = 0;
268 // --------------------------------------------------------------------------
269       fHijing->GenerateEvent();
270       fTrials++;
271       fHijing->ImportParticles(particles,"All");
272       if (fLHC) Boost(particles);
273       
274       Int_t np = particles->GetEntriesFast();
275       printf("\n **************************************************%d\n",np);
276       Int_t nc = 0;
277       if (np == 0 ) continue;
278       Int_t i;
279       Int_t * newPos = new Int_t[np];
280
281       for (i = 0; i < np; i++) *(newPos+i) = i;
282 //      Get event vertex
283 //
284       TParticle *  iparticle = (TParticle *) particles->At(0);
285       fEventVertex[0] = origin0[0];
286       fEventVertex[1] = origin0[1];     
287       fEventVertex[2] = origin0[2];
288       
289 //
290 //      First write parent particles
291 //
292
293       for (i = 0; i < np; i++) {
294           iparticle       = (TParticle *) particles->At(i);
295 // Is this a parent particle ?
296           if (Stable(iparticle)) continue;
297 //
298         Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
299         Bool_t  selected             =  kTRUE;
300         Bool_t  hasSelectedDaughters =  kFALSE;
301             
302             
303         kf        = iparticle->GetPdgCode();
304         ks        = iparticle->GetStatusCode();
305         if (kf == 92) continue;
306             
307         if (!fSelectAll) selected = KinematicSelection(iparticle, 0)&&SelectFlavor(kf);
308         hasSelectedDaughters = DaughtersSelection(iparticle, particles);
309 //
310 // Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
311 //
312         if (selected || hasSelectedDaughters) {
313             nc++;
314             p[0] = iparticle->Px();
315             p[1] = iparticle->Py();
316             p[2] = iparticle->Pz();
317             origin[0] = origin0[0]+iparticle->Vx()/10;
318             origin[1] = origin0[1]+iparticle->Vy()/10;
319             origin[2] = origin0[2]+iparticle->Vz()/10;
320             tof = kconv*iparticle->T();
321             imo = -1;
322             TParticle* mother = 0;
323             if (hasMother) {
324                 imo = iparticle->GetFirstMother();
325                 mother = (TParticle *) particles->At(imo);
326                 imo = (mother->GetPdgCode() != 92) ? imo =* (newPos+imo) : -1;
327                 
328             }
329 // Put particle on the stack ... 
330 //              printf("\n set track mother: %d %d %d %d %d %d ",i,imo, kf, nt+1, selected, hasSelectedDaughters);
331             
332             SetTrack(0,imo,kf,p,origin,polar, tof,kPPrimary,nt);
333 // ... and keep it there
334             KeepTrack(nt);
335 //
336             *(newPos+i)=nt;
337         } // selected
338       } // particle loop parents
339 //
340 // Now write the final state particles
341 //
342
343       for (i = 0; i<np; i++) {
344         TParticle *  iparticle       = (TParticle *) particles->At(i);
345 // Is this a final state particle ?
346         if (!Stable(iparticle)) continue;
347
348         Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
349         Bool_t  selected             =  kTRUE;
350         kf        = iparticle->GetPdgCode();
351         ks        = iparticle->GetStatusCode();
352 // --------------------------------------------------------------------------
353 // Count spectator neutrons and protons
354         if(ks == 0 || ks == 1 || ks == 10 || ks == 11){
355           if(kf == kNeutron) fSpecn += 1;
356           if(kf == kProton)  fSpecp += 1;
357         }
358 // --------------------------------------------------------------------------
359 //          
360         if (!fSelectAll) {
361             selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
362             if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10
363                                                       && ks != 11);
364         }
365 //
366 // Put particle on the stack if selected
367 //
368         if (selected) {
369             nc++;
370             p[0] = iparticle->Px();
371             p[1] = iparticle->Py();
372             p[2] = iparticle->Pz();
373             origin[0] = origin0[0]+iparticle->Vx()/10;
374             origin[1] = origin0[1]+iparticle->Vy()/10;
375             origin[2] = origin0[2]+iparticle->Vz()/10;
376             tof = kconv*iparticle->T();
377             imo = -1;
378             TParticle* mother = 0;
379             if (hasMother) {
380                 imo = iparticle->GetFirstMother();
381                 mother = (TParticle *) particles->At(imo);
382                 imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1;
383             }   
384 // Put particle on the stack
385             SetTrack(fTrackIt,imo,kf,p,origin,polar,
386                                                      tof,kPNoProcess,nt);
387             KeepTrack(nt);
388             *(newPos+i)=nt;
389         } // selected
390       } // particle loop final state
391  
392       delete[] newPos;
393
394       printf("\n I've put %i particles on the stack \n",nc);
395       if (nc > 0) {
396           jev += nc;
397           if (jev >= fNpart || fNpart == -1) {
398               fKineBias = Float_t(fNpart)/Float_t(fTrials);
399               printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
400               break;
401           }
402       }
403   } // event loop
404   MakeHeader();
405   
406   SetHighWaterMark(nt);
407 }
408
409 void AliGenHijing::KeepFullEvent()
410 {
411     fKeep=1;
412 }
413
414 void AliGenHijing::EvaluateCrossSections()
415 {
416 //     Glauber Calculation of geometrical x-section
417 //
418     Float_t xTot       = 0.;          // barn
419     Float_t xTotHard   = 0.;          // barn 
420     Float_t xPart      = 0.;          // barn
421     Float_t xPartHard  = 0.;          // barn 
422     Float_t sigmaHard  = 0.1;         // mbarn
423     Float_t bMin       = 0.;
424     Float_t bMax       = fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
425     const Float_t kdib = 0.2;
426     Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
427
428
429     printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
430     printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
431     Int_t i;
432     Float_t oldvalue= 0.;
433
434     Float_t* b   = new Float_t[kMax];
435     Float_t* si1 = new Float_t[kMax];    
436     Float_t* si2 = new Float_t[kMax];    
437     
438     for (i = 0; i < kMax; i++)
439     {
440         Float_t xb  = bMin+i*kdib;
441         Float_t ov;
442         ov=fHijing->Profile(xb);
443         Float_t gb  =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
444         Float_t gbh =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
445         xTot+=gb;
446         xTotHard += gbh;
447         if (xb > fMinImpactParam && xb < fMaxImpactParam)
448         {
449             xPart += gb;
450             xPartHard += gbh;
451         }
452         
453         if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
454         oldvalue = xTot;
455         printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
456         printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
457         if (i>0) {
458             si1[i] = gb/kdib;
459             si2[i] = gbh/gb;
460             b[i]  = xb;
461         }
462     }
463
464     printf("\n Total cross section (barn): %f \n",xTot);
465     printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
466     printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
467     printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
468
469 //  Store result as a graph
470     b[0] = 0;
471     si1[0] = 0;
472     si2[0]=si2[1];
473     
474     fDsigmaDb  = new TGraph(i, b, si1);
475     fDnDb      = new TGraph(i, b, si2);
476 }
477
478 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
479 {
480 //
481 // Looks recursively if one of the daughters has been selected
482 //
483 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
484     Int_t imin = -1;
485     Int_t imax = -1;
486     Int_t i;
487     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
488     Bool_t selected = kFALSE;
489     if (hasDaughters) {
490         imin = iparticle->GetFirstDaughter();
491         imax = iparticle->GetLastDaughter();       
492         for (i = imin; i <= imax; i++){
493             TParticle *  jparticle = (TParticle *) particles->At(i);    
494             Int_t ip = jparticle->GetPdgCode();
495             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
496                 selected=kTRUE; break;
497             }
498             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
499         }
500     } else {
501         return kFALSE;
502     }
503     return selected;
504 }
505
506
507 Bool_t AliGenHijing::SelectFlavor(Int_t pid)
508 {
509 // Select flavor of particle
510 // 0: all
511 // 4: charm and beauty
512 // 5: beauty
513     if (fFlavor == 0) return kTRUE;
514     
515     Int_t ifl = TMath::Abs(pid/100);
516     if (ifl > 10) ifl/=10;
517     return (fFlavor == ifl);
518 }
519
520 Bool_t AliGenHijing::Stable(TParticle*  particle)
521 {
522 // Return true for a stable particle
523 //
524     if (particle->GetFirstDaughter() < 0 )
525     {
526         return kTRUE;
527     } else {
528         return kFALSE;
529     }
530 }
531
532
533 void AliGenHijing::Boost(TClonesArray* particles)
534 {
535 //
536 // Boost cms into LHC lab frame
537 //
538     Double_t dy    = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) / 
539                                       (Double_t(fZTarget)    * Double_t(fAProjectile)));
540     Double_t beta  = TMath::TanH(dy);
541     Double_t gamma = 1./TMath::Sqrt(1.-beta*beta);
542     Double_t gb    = gamma * beta;
543
544     printf("\n Boosting particles to lab frame %f %f %f", dy, beta, gamma);
545     
546     Int_t i;
547     Int_t np = particles->GetEntriesFast();
548     for (i = 0; i < np; i++) 
549     {
550         TParticle* iparticle = (TParticle*) particles->At(i);
551
552         Double_t e   = iparticle->Energy();
553         Double_t px  = iparticle->Px();
554         Double_t py  = iparticle->Py();
555         Double_t pz  = iparticle->Pz();
556
557         Double_t eb  = gamma * e -      gb * pz;
558         Double_t pzb =   -gb * e +   gamma * pz;
559
560         iparticle->SetMomentum(px, py, pzb, eb);
561     }
562 }
563
564
565 void AliGenHijing::MakeHeader()
566 {
567 // Builds the event header, to be called after each event
568     AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing");
569     ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT());
570     ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
571     ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
572     ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
573     ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
574     ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(),
575                                                        fHijing->GetN01(),
576                                                        fHijing->GetN10(),
577                                                        fHijing->GetN11());
578     ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
579
580 // 4-momentum vectors of the triggered jets.
581 //
582 // Before final state gluon radiation.
583     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
584                                               fHijing->GetHINT1(22),
585                                               fHijing->GetHINT1(23),
586                                               fHijing->GetHINT1(24));
587
588     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
589                                               fHijing->GetHINT1(32),
590                                               fHijing->GetHINT1(33),
591                                               fHijing->GetHINT1(34));
592 // After final state gluon radiation.
593     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
594                                               fHijing->GetHINT1(27),
595                                               fHijing->GetHINT1(28),
596                                               fHijing->GetHINT1(29));
597
598     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
599                                               fHijing->GetHINT1(37),
600                                               fHijing->GetHINT1(38),
601                                               fHijing->GetHINT1(39));
602     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
603 // Event Vertex
604     header->SetPrimaryVertex(fEventVertex);
605     gAlice->SetGenEventHeader(header);    
606 }
607
608
609 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
610 {
611 // Assignment operator
612     return *this;
613 }
614
615 #ifndef WIN32
616 # define rluget_hijing rluget_hijing_
617 # define rluset_hijing rluset_hijing_
618 # define rlu_hijing    rlu_hijing_
619 # define type_of_call
620 #else
621 # define rluget_hijing RLUGET_HIJING
622 # define rluset_hijing RLUSET_HIJING
623 # define rlu_hijing    RLU_HIJING
624 # define type_of_call _stdcall
625 #endif
626
627
628 extern "C" {
629   void type_of_call rluget_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
630   {printf("Dummy version of rluget_hijing reached\n");}
631
632   void type_of_call rluset_hijing(Int_t & /*lfn*/, Int_t & /*move*/)
633   {printf("Dummy version of rluset_hijing reached\n");}
634
635   Double_t type_of_call rlu_hijing(Int_t & /*idum*/) 
636   {
637       Float_t r;
638       do r=sRandom->Rndm(); while(0 >= r || r >= 1);
639       return r;
640   }
641 }