Added preprocessor conditionals to support ROOT > 5.11.2.
[u/mrichter/AliRoot.git] / HBTAN / AliHBTLLWeights.cxx
1 #include "AliHBTLLWeights.h"
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 //_________________________________________________________________________
18 ///////////////////////////////////////////////////////////////////////////
19 //
20 //  class AliHBTLLWeights
21 //
22 //  This class introduces the weight's calculation 
23 //  according to the Lednicky's algorithm.
24 //  
25 //  
26 //  fsiw.f, fsiini.f  
27 //
28 //  Description from fortran code by author R. Lednicky
29 //
30 //  Calculates final state interaction (FSI) weights                     
31 //  WEIF = weight due to particle - (effective) nucleus FSI (p-N)
32 //  WEI  = weight due to p-p-N FSI
33 //  WEIN = weight due to p-p FSI; note that WEIN=WEI if I3C=0;
34 //                                note that if I3C=1 the calculation of
35 //                                WEIN can be skipped by putting J=0
36 //.......................................................................
37 //  Correlation Functions:
38 //  CF(p-p-N)   = sum(WEI)/sum(WEIF)
39 //  CF(p-p)     = sum(WEIN)/sum(1); here the nucleus is completely
40 //                                  inactive
41 //  CF(p-p-"N") = sum(WEIN*WEIF')/sum(WEIF'), where WEIN and WEIF'
42 //                are not correlated (calculated at different emission
43 //                points, e.g., for different events);
44 //                thus here the nucleus affects one-particle
45 //                spectra but not the correlation
46 //.......................................................................
47 //  User must supply data file <fn> on unit NUNIT (e.g. =11) specifying
48 //  LL   : particle pair
49 //  NS   : approximation used to calculate Bethe-Salpeter amplitude
50 //  ITEST: test switch
51 //         If ITEST=1 then also following parameters are required
52 //  ICH  : 1(0) Coulomb interaction between the two particles ON (OFF)
53 //  IQS  : 1(0) quantum statistics for the two particles ON (OFF)
54 //  ISI  : 1(0) strong interaction between the two particles ON (OFF)
55 //  I3C  : 1(0) Coulomb interaction with residual nucleus ON (OFF)
56 //  This data file can contain other information useful for the user.
57 //  It is read by subroutines READINT4 and READREA8(4) (or READ_FILE).
58 //  -------------------------------------------------------------------
59 //-   LL       1  2  3  4   5    6   7  8  9 10  11  12  13  14 15 16 17
60 //-   part. 1: n  p  n alfa pi+ pi0 pi+ n  p pi+ pi+ pi+ pi- K+ K+ K+ K-
61 //-   part. 2: n  p  p alfa pi- pi0 pi+ d  d  K-  K+  p   p  K- K+ p  p
62 //   NS=1 y/n: +  +  +  +   +    -   -  -  -  -   -   -   -  -  -  -  -
63 //  -------------------------------------------------------------------
64 //-   LL       18  19 20  21  22 23 24 25 26    27     28
65 //-   part. 1:  d  d   t  t   K0 K0  d p  p      p      n
66 //-   part. 2:  d alfa t alfa K0 K0b t t alfa lambda lambda
67 //   NS=1 y/n:  -  -   -  -   -  -   - -  -      +      +
68 //  -------------------------------------------------------------------
69 //   NS=1  Square well potential,
70 //   NS=3  not used
71 //   NS=4  scattered wave approximated by the spherical wave,
72 //   NS=2  same as NS=4 but the approx. of equal emission times in PRF
73 //       not required (t=0 approx. used in all other cases).
74 //   Note: if NS=2,4, the B-S amplitude diverges at zero distance r* in
75 //       the two-particle c.m.s.; user can specify a cutoff AA in
76 //       SUBROUTINE FSIINI, for example:
77 //       IF(NS.EQ.2.OR.NS.EQ.4)AA=5.D0 !! in 1/GeV --> AA=1. fm
78 //  ------------------------------------------------------------------
79 //  ITEST=1 any values of parameters ICH, IQS, ISI, I3C are allowed
80 //          and should be given in data file <fn>
81 //  ITEST=0 physical values of these parameters are put automatically
82 //          in FSIINI (their values are not required in data file)
83 //=====================================================================
84 //  At the beginning of calculation user should call FSIINI,
85 //  which reads LL, NS, ITEST (and eventually ICH, IQS, ISI, I3C)
86 //  and ializes various parameters.
87 //  In particular the constants in
88 //    COMMON/FSI_CONS/PI,PI2,SPI,DR,W
89 //  may be useful for the user:
90 //   W=1/.1973D0    ! from fm to 1/GeV
91 //   PI=4*DATAN(1.D0)
92 //   PI2=2*PI
93 //   SPI=TMath::Sqrt(PI)
94 //   DR=180.D0/PI   ! from radian to degree
95 //    _______________________________________________________
96 //  !! |Important note: all real quantities are assumed REAL*8 | !!
97 //    -------------------------------------------------------
98 //  For each event user should fill in the following information
99 //  in COMMONs (all COMMONs in FSI calculation start with FSI_):
100 //  ...................................................................
101 //   COMMON/FSI_POC/AMN,AM1,AM2,CN,C1,C2,AC1,AC2
102 //  Only
103 //       AMN  = mass of the effective nucleus   [GeV/c**2]
104 //       CN   = charge of the effective nucleus [elem. charge units]
105 //  are required
106 //  ...................................................................
107 //   COMMON/FSI_MOM/P1X,P1Y,P1Z,E1,P1, !part. momenta in the rest frame 
108 //  1               P2X,P2Y,P2Z,E2,P2  !of effective nucleus (NRF)
109 //  Only the components
110 //                      PiX,PiY,PiZ  [GeV/c]
111 //  in NRF are required.
112 //  To make the corresponding Lorentz transformation user can use the
113 //  subroutines LTRAN and LTRANB
114 //  ...................................................................
115 //  COMMON/FSI_COOR/X1,Y1,Z1,T1,R1,     ! 4-coord. of emission
116 //  1               X2,Y2,Z2,T2,R2      ! points in NRF
117 //  The componets
118 //                     Xi,Yi,Zi  [fm]
119 //  and emission times
120 //                        Ti   [fm/c]
121 //  should be given in NRF with the origin assumed at the center
122 //  of the effective nucleus. If the effect of residual nucleus is
123 //  not calculated within FSIW, the NRF can be any fixed frame.
124 //  --------------------------------------------------------------------
125 //  Before calling FSIW the user must call
126 //   CALL LTRAN12
127 //  Besides Lorentz transformation to pair rest frame:
128 //  (p1-p2)/2 --> k* it also transforms 4-coordinates of
129 //  emission points from fm to 1/GeV and calculates Ei,Pi and Ri.
130 //  Note that |k*|=AK in COMMON/FSI_PRF/
131 //  --------------------------------------------------------------------
132 //  After making some additional filtering using k* (say k* < k*max)
133 //  or direction of vector k*,
134 //  user can finally call FSIW to calculate the FSI weights
135 //  to be used to construct the correlation function
136 //======================================================================
137
138
139 /*******************************************************************/
140 /******      ROUTINES    USED    FOR     COMMUNUCATION      ********/
141 /********************     WITH      FORTRAN     ********************/
142 /*******************************************************************/
143 #ifndef WIN32
144 # define led_bldata led_bldata_
145 # define fsiini fsiini_
146 # define ltran12 ltran12_
147 # define boosttoprf boosttoprf_
148 # define fsiw fsiw_
149 # define setpdist setpdist_
150 # define type_of_call
151 #else
152 # define led_bldata LED_BLDATA
153 # define fsiini FSIINI
154 # define ltran12 LTRAN12
155 # define boosttoprf BOOSTTOPRF
156 # define fsiw FSIW
157 # define setpdist SETPDIST
158 # define type_of_call _stdcall
159 #endif
160 /****************************************************************/
161 extern "C" void type_of_call led_bldata(); 
162 extern "C" void type_of_call fsiini();
163 extern "C" void type_of_call ltran12();
164 extern "C" void type_of_call boosttoprf();
165 extern "C" void type_of_call fsiw();
166 extern "C" void type_of_call setpdist(Double_t& r);
167 /**************************************************************/
168
169 #include "AliHBTPair.h"
170 #include "AliVAODParticle.h"
171 #include "WLedCOMMONS.h"
172 #include <TRandom.h>   
173 #include <TMath.h>     
174 #include <TPDGCode.h>
175 #include <TParticlePDG.h>
176 #include <TDatabasePDG.h>
177
178
179 ClassImp(AliHBTLLWeights)  
180  
181 AliHBTLLWeights* AliHBTLLWeights::fgLLWeights = 0x0; 
182 const Double_t AliHBTLLWeights::fgkWcons = 1./0.1973;
183
184 AliHBTLLWeights::AliHBTLLWeights():
185  fTest(kTRUE),
186  fColoumbSwitch(kTRUE),
187  fQuantStatSwitch(kTRUE),
188  fStrongInterSwitch(kTRUE),
189  fColWithResidNuclSwitch(kFALSE),
190  fNuclMass(0.0),
191  fNuclCharge(0.0),
192  fRandomPosition(kNone),
193  fRadius(0.0),
194  fOneMinusLambda(0.0),
195  fPID1(0),
196  fPID2(0),
197  fSigma(0.0)
198 {
199 // Default Constructor 
200   if (fgLLWeights)
201    Fatal("AliHBTLLWeights","LLWeights already instatiated. Use AliHBTLLWeights::Instance()");
202 }
203 /**************************************************************/
204
205 AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/):
206  AliHBTWeights(),
207  fTest(kTRUE),
208  fColoumbSwitch(kTRUE),
209  fQuantStatSwitch(kTRUE),
210  fStrongInterSwitch(kTRUE),
211  fColWithResidNuclSwitch(kFALSE),
212  fNuclMass(0.0),
213  fNuclCharge(0.0),
214  fRandomPosition(kNone),
215  fRadius(0.0),
216  fOneMinusLambda(0.0),
217  fPID1(0),
218  fPID2(0),
219  fSigma(0.0)
220 {
221   //Copy ctor needed by the coding conventions but not used
222   Fatal("AliHBTLLWeights","copy ctor not implemented");
223 }
224 /************************************************************/
225
226 AliHBTLLWeights& AliHBTLLWeights::operator=(const AliHBTLLWeights& /*source*/)
227 {
228   //Assignment operator needed by the coding conventions but not used
229   Fatal("AliHBTLLWeights","assignment operator not implemented");
230   return * this;
231 }
232 /************************************************************/
233
234 AliHBTLLWeights* AliHBTLLWeights::Instance()
235 {     
236 // returns instance of class 
237  if (fgLLWeights) 
238   {
239     return fgLLWeights;
240   } 
241  else 
242   {
243    fgLLWeights = new AliHBTLLWeights();            
244    return fgLLWeights; 
245   } 
246 }     
247 /************************************************************/
248
249 void AliHBTLLWeights::Set()
250 {
251  //sets this as weighitng class
252  Info("Set","Setting Lednicky-Lyuboshitz as Weighing Class");
253  
254  if ( fgWeights == 0x0 )  
255   {
256     fgWeights = AliHBTLLWeights::Instance();
257     return;
258   }  
259  if ( fgWeights == AliHBTLLWeights::Instance() ) return;
260  delete fgWeights;
261  fgWeights = AliHBTLLWeights::Instance();
262 }
263 /************************************************************/
264
265 Double_t AliHBTLLWeights::GetWeight(AliHBTPair* partpair)
266 {
267 // calculates weight for a pair
268   static const Double_t kcmtofm = 1.e13;
269   static const Double_t kcmtoOneOverGeV = kcmtofm*fgkWcons;  
270   
271   AliVAODParticle *part1 = partpair->Particle1();
272   AliVAODParticle *part2 = partpair->Particle2();
273
274   if ( (part1 == 0x0) || (part2 == 0x0))
275    {
276      Error("GetWeight","Null particle pointer");
277      return 0.0;
278    }
279
280   if ( fPID1 != part1->GetPdgCode() ) return 1.0; 
281   if ( fPID2 != part2->GetPdgCode() ) return 1.0; 
282
283 //takes a lot of time
284   if ( (part1->Px() == part2->Px()) && 
285        (part1->Py() == part2->Py()) && 
286        (part1->Pz() == part2->Pz()) )
287    {
288      return 0.0;
289    }
290
291   if ((!fRandomPosition) && 
292       (part1->Vx()  == part2->Vx()) && 
293       (part1->Vy()  == part2->Vy()) && 
294       (part1->Vz()  == part2->Vz()) )
295     {        
296       return 0.0;
297     }
298
299   if(fOneMinusLambda)//implemetation of non-zero intetcept parameter 
300    {
301      if( gRandom->Rndm() < fOneMinusLambda ) return 1.0;
302    }
303     
304
305   //this must  be after ltran12 because it would overwrite what we set below
306   switch (fRandomPosition)
307    {
308     case kGaussianQInv:
309      {
310        // Set Momenta in the common block
311        FSI_MOM.P1X = part1->Px();
312        FSI_MOM.P1Y = part1->Py();
313        FSI_MOM.P1Z = part1->Pz();
314
315        FSI_MOM.P2X = part2->Px();
316        FSI_MOM.P2Y = part2->Py();
317        FSI_MOM.P2Z = part2->Pz();
318       
319        //boost it to PRF
320        ltran12();
321        boosttoprf();
322        
323        //Set particle positions now so they are Gaussian in PRF
324        RandomPairDistances();
325        
326        FSI_PRF.RPS= FSI_COOR.X2*FSI_COOR.X2 + FSI_COOR.Y2*FSI_COOR.Y2 +FSI_COOR.Z2*FSI_COOR.Z2;
327        FSI_PRF.RP=TMath::Sqrt(FSI_PRF.RPS);
328        
329        break;
330       } 
331     case kGaussianOSL:
332       { 
333        //boost pair to LCMS and set such a momenta in the common block
334        Int_t retv = SetMomentaInLCMS(part1,part2);
335        if (retv) return 1;
336        //random particle positions/distance so they are Gaussian in LCMS
337        RandomPairDistances();
338        
339        //Boost to PRF
340        ltran12();
341        boosttoprf();
342     
343        break;
344       }
345     case kNone:
346     default:
347        //set momenta and paricle positions as they are
348        FSI_MOM.P1X = part1->Px();
349        FSI_MOM.P1Y = part1->Py();
350        FSI_MOM.P1Z = part1->Pz();
351
352        FSI_MOM.P2X = part2->Px();
353        FSI_MOM.P2Y = part2->Py();
354        FSI_MOM.P2Z = part2->Pz();
355
356        FSI_COOR.X1 = part1->Vx()*kcmtoOneOverGeV;
357        FSI_COOR.Y1 = part1->Vy()*kcmtoOneOverGeV;
358        FSI_COOR.Z1 = part1->Vz()*kcmtoOneOverGeV;
359        FSI_COOR.T1 = part1->T()*kcmtoOneOverGeV;
360
361        FSI_COOR.X2 = part2->Vx()*kcmtoOneOverGeV;
362        FSI_COOR.Y2 = part2->Vy()*kcmtoOneOverGeV;
363        FSI_COOR.Z2 = part2->Vz()*kcmtoOneOverGeV;
364        FSI_COOR.T2 = part2->T()*kcmtoOneOverGeV;
365
366        ltran12();
367        boosttoprf();
368        
369        break;
370    }
371
372   fsiw();
373   return LEDWEIGHT.WEIN;
374 }
375 /************************************************************/
376
377 void AliHBTLLWeights::Init()
378 {
379 //initial parameters of model
380
381   FSI_NS.NS = fApproximationModel;      
382   
383   LEDWEIGHT.ITEST = fTest;  
384   if(fTest)
385    {
386      FSI_NS.ICH = fColoumbSwitch;
387      FSI_NS.ISI = fStrongInterSwitch;
388      FSI_NS.IQS = fQuantStatSwitch;
389      FSI_NS.I3C = fColWithResidNuclSwitch;
390      LEDWEIGHT.IRANPOS = fRandomPosition;
391    }
392  
393   if ( (fPID1 == 0) || (fPID2 == 0) )
394    {
395      Fatal("Init","Particles types are not set");
396      return;//pro forma
397    }
398   
399   
400   FSI_NS.LL = GetPairCode(fPID1,fPID2);
401        
402   if (FSI_NS.LL == 0) 
403    {
404      Fatal("Init","Particles types are not supported");
405      return;//pro forma
406    }
407
408   Info("Init","Setting PIDs %d %d. LL Code is %d",fPID1,fPID2,FSI_NS.LL);
409
410
411   TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
412   if (tpart1 == 0x0)
413    {
414      Fatal("init","We can not find particle with ID=%d in PDG DataBase",fPID1);
415      return;
416    }
417       
418   FSI_POC.AM1=tpart1->Mass();
419   FSI_POC.C1=tpart1->Charge(); 
420
421   TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
422 //lv
423   if (tpart2 == 0x0)
424    {
425      Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2);
426      return;
427    }
428
429   FSI_POC.AM2=tpart2->Mass();
430   FSI_POC.C1=tpart2->Charge();
431
432   led_bldata();
433   fsiini();
434
435
436 //constants for radii simulation 
437
438   if(fRandomPosition)
439    {
440      fSigma =TMath::Sqrt(2.)*fRadius;     
441    } 
442
443 /************************************************************/
444
445 Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
446 {
447 //returns Code corresponding to that pair
448  return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
449 }
450 /************************************************************/
451
452 Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
453 {
454 // returns code corresponding to the pair of PIDs
455 //   pairCode   1  2  3  4   5    6   7  8  9 10  11  12  13  14 15 16 17 18  19  20   21   22  23 24 25 26    27     28
456 //   hpid:      n  p  n alfa pi+ pi0 pi+ n  p pi+ pi+ pi+ pi- K+ K+ K+ K-  d  d    t   t    K0  K0  d p  p      p      n
457 //   lpid:      n  p  p alfa pi- pi0 pi+ d  d  K-  K+  p   p  K- K+ p  p   d alfa  t  alfa  K0  K0b t t alfa lambda lambda
458 //   NS=1 y/n:  +  +  +  +   +    -   -  -  -  -   -   -   -  -  -  -  -   -  -    -    -    -  -   - -  -      -      -
459
460 //alphas, deuterons and tyts are NOT supported here
461
462   Int_t chargefactor = 1;
463   Int_t hpid; //pid in higher row
464   Int_t lpid; //pid in lower row
465   Int_t code; //pairCode
466   
467   Bool_t swap;
468   
469 //determine the order of selcetion in switch  
470   if (TMath::Abs(pid1) < TMath::Abs(pid2) ) 
471    {
472     if (pid1<0) chargefactor=-1;
473     hpid=pid2*chargefactor;
474     lpid=pid1*chargefactor;
475     swap = kFALSE;
476    } 
477   else 
478    {
479     if (pid2<0) chargefactor=-1;
480     hpid=pid1*chargefactor;
481     lpid=pid2*chargefactor;
482     swap = kTRUE;
483    }
484
485 //mlv
486    hpid=pid1;
487    lpid=pid2;
488
489
490 //Determine the pair code
491   switch (hpid) //switch on first  particle id
492    {
493      case kNeutron:
494       switch (lpid)
495        {
496          case kNeutron: 
497            code = 1;  //neutron neutron
498            break;
499         
500          case kProton: 
501            code = 3;  //neutron proton
502            break;
503            
504          case kLambda0: 
505            code = 28;  //neutron lambda
506            break;
507            
508          default: 
509            return 0; //given pair not supported
510            break;
511        }
512       break;
513
514      case kProton:
515       switch (lpid)
516        {
517          case kProton:
518            code = 2; //proton proton
519            break;
520            
521          case kLambda0: 
522            code = 27;//proton lambda
523            break;
524            
525          default: 
526            return 0; //given pair not supported
527            break;
528            
529        }
530       break;
531
532      case kPiPlus:
533      
534       switch (lpid)
535        {
536          case kPiPlus:
537            code = 7; //piplus piplus
538            break;
539
540          case kPiMinus:
541            code = 5; //piplus piminus
542            break;
543         
544          case kKMinus:
545            code = 10; //piplus Kminus
546            break;
547
548          case kKPlus:
549            code = 11; //piplus Kplus
550            break;
551
552          case kProton:
553            code = 12; //piplus proton
554            chargefactor*=-1;
555            break;
556
557          default: 
558            return 0; //given pair not supported
559            break;
560        }
561       break;
562      case kPi0:
563       switch (lpid)
564        {
565          case kPi0:
566            code = 6;
567            break;
568            
569          default: 
570            return 0; //given pair not supported
571            break;
572        }
573       break;
574       
575      case kKPlus:
576       switch (lpid)
577        {
578          case kKMinus:
579            code = 14; //Kplus Kminus
580            break;
581
582          case kKPlus:
583            code = 15; //Kplus Kplus
584            break;
585
586          case kProton:
587            code = 16; //Kplus proton
588            break;
589            
590          default: 
591            return 0; //given pair not supported
592            break;
593        }
594       break;
595       
596      case kKMinus:
597       switch (lpid)
598        {
599          case kProton:
600            code = 17; //Kminus proton
601            chargefactor*=1;
602            break;
603            
604          default: 
605            return 0; //given pair not supported
606            break;
607        }
608       break;
609       
610      case kK0:
611       switch (lpid)
612        {
613          case kK0:
614            code = 2; //Kzero Kzero
615            break;
616          
617          case kK0Bar:
618            code = 17; //Kzero KzeroBar
619            break;
620
621          default: 
622            return 0; //given pair not supported
623            break;
624        }
625       break;
626
627      default: return 0;
628    }
629   return code;
630 }
631 /************************************************************/
632
633 void AliHBTLLWeights::SetTest(Bool_t rtest)
634 {
635   //Sets fTest member
636   fTest = rtest;
637
638 /************************************************************/
639
640 void AliHBTLLWeights::SetColoumb(Bool_t col)
641 {
642   // (ICH in fortran code) Coulomb interaction between the two particles ON (OFF)
643   fColoumbSwitch = col;
644 }
645 /************************************************************/
646
647 void AliHBTLLWeights::SetQuantumStatistics(Bool_t qss)
648 {
649   //IQS: quantum statistics for the two particles ON (OFF) 
650   //if non-identical particles automatically off
651   fQuantStatSwitch = qss;
652 }
653 /************************************************************/
654
655 void AliHBTLLWeights::SetStrongInterSwitch(Bool_t sis)
656 {
657   //ISI: strong interaction between the two particles ON (OFF)
658   fStrongInterSwitch = sis;
659 }
660 /************************************************************/
661
662 void AliHBTLLWeights::SetColWithResidNuclSwitch(Bool_t crn)
663 {
664   //I3C: Coulomb interaction with residual nucleus ON (OFF)  
665   fColWithResidNuclSwitch = crn;
666 }
667 /************************************************************/
668
669 void AliHBTLLWeights::SetApproxModel(Int_t ap)
670 {
671   //sets  Model of Approximation (NS in Fortran code)
672   fApproximationModel=ap;
673 }
674 /************************************************************/
675      
676 void AliHBTLLWeights::SetRandomPosition(ERandomizationWay rw)
677
678 // rw can be: kGaussianQInv - so 1D Qinv correlation function has a Gaussian shape 
679 //                            (and also Qlong and Qside, but not Qout)
680 //            kGaussianOSL - so 3D Qout-Qside-Qlong correlation function has a Gaussian shape
681 //            kNone  - no randomization performed (DEFAULT)
682  fRandomPosition = rw;
683 }
684 /************************************************************/
685
686 void AliHBTLLWeights::SetR1dw(Double_t R)
687 {
688   //spherical source model radii
689   fRadius=R;
690 }
691 /************************************************************/
692
693 void AliHBTLLWeights::SetParticlesTypes(Int_t pid1, Int_t pid2)
694 {
695   //set AliRoot particles types   
696   fPID1 = pid1; 
697   fPID2 = pid2;
698 }
699 /************************************************************/
700     
701 void AliHBTLLWeights::SetNucleusCharge(Double_t ch)
702 {
703   // not used now  (see comments in fortran code)
704   fNuclCharge=ch;
705 }
706 /************************************************************/
707
708 void AliHBTLLWeights::SetNucleusMass(Double_t mass)
709 {
710   // (see comments in fortran code)
711   fNuclMass=mass;
712 }
713 /************************************************************/
714
715 Int_t AliHBTLLWeights::SetMomentaInLCMS(AliVAODParticle* part1, AliVAODParticle* part2)
716 {
717 //sets paricle momenta in the common block in LCMS frame
718 //---> Particle energies ---------
719
720   Double_t p1x = part1->Px();
721   Double_t p1y = part1->Py();
722   Double_t p1z = part1->Pz();
723       
724   Double_t p2x = part2->Px();
725   Double_t p2y = part2->Py();
726   Double_t p2z = part2->Pz();
727
728   Double_t am1 = part1->Mass();
729   Double_t am2 = part2->Mass();
730   
731   Double_t p1s=p1x*p1x+p1y*p1y+p1z*p1z;
732   Double_t p2s=p2x*p2x+p2y*p2y+p2z*p2z;
733   Double_t e1=TMath::Sqrt(am1*am1+p1s);
734   Double_t e2=TMath::Sqrt(am2*am2+p2s);
735 //---> pair parameters -----------
736   Double_t e12=e1+e2;       // energy
737   Double_t p12x=p1x+p2x;    // px
738   Double_t p12y=p1y+p2y;    // py
739   Double_t p12z=p1z+p2z;    // pz
740   Double_t p12s=p12x*p12x+p12y*p12y+p12z*p12z;
741   Double_t p12 =TMath::Sqrt(p12s);// momentum
742   Double_t v12 =p12/e12;    // velocity
743
744   Double_t cth =p12z/p12;   // cos(theta)
745 //  Double_t sth =TMath::Sqrt(1. - cth*cth); //sin
746   Double_t v12z=v12*cth;    // longit. v
747   Double_t gamz=1./TMath::Sqrt(1.-v12z*v12z);
748 //--      v12t=v12*sth    // transv. v in cms (not needed)
749
750
751   Double_t p12ts=p12x*p12x+p12y*p12y;
752   Double_t p12t=TMath::Sqrt(p12ts); //pt
753 //===> azimuthal rotation (pt||x) ============
754    if(p12t != 0.0) 
755     {
756      Double_t cphi=p12x/p12t;  // cos(phi)
757      Double_t sphi=p12y/p12t;  // sin(phi)
758      Rotate(p1x,p1y,sphi,cphi,p1x,p1y);       
759      Rotate(p2x,p2y,sphi,cphi,p2x,p2y);
760 //     Rotate(x1,y1,sphi,cphi,x1,y1);
761 //     Rotate(x2,y2,sphi,cphi,x2,y2);
762      }    
763    else // rotation impossible 
764     {
765       return 1;
766     }
767
768 //===> co-moving ref. frame       ============
769    
770    Double_t nothing;
771    Boost(p1z,e1,v12z,gamz,p1z,nothing);
772    Boost(p2z,e2,v12z,gamz,p2z,nothing);
773
774 //   p1s=p1x*p1x+p1y*p1y+p1z*p1z;
775 //   p2s=p2x*p2x+p2y*p2y+p2z*p2z;
776 //   e1=TMath::Sqrt(am1*am1+p1s);
777 //   e2=TMath::Sqrt(am2*am2+p2s);
778 //   Boost(z1,t1,v12z,gamz,z1,t1);
779 //   Boost(z2,t2,v12z,gamz,z2,t2);
780
781
782      FSI_MOM.P1X = p1x;
783      FSI_MOM.P1Y = p1y;
784      FSI_MOM.P1Z = p1z;
785
786      FSI_MOM.P2X = p2x;
787      FSI_MOM.P2Y = p2y;
788      FSI_MOM.P2Z = p2z;
789      
790 //     Info("Simulate"," %f %f %f ",p1x + p2x, p1y + p2y, p1z + p2z);
791      return 0;
792 }
793
794 /**********************************************************/  
795
796 void AliHBTLLWeights::RandomPairDistances()
797 {
798  //randomizes pair distances
799   Double_t rxcm = fSigma*gRandom->Gaus();
800   Double_t rycm = fSigma*gRandom->Gaus();
801   Double_t rzcm = fSigma*gRandom->Gaus();
802
803
804   FSI_COOR.X1 = 0;
805   FSI_COOR.Y1 = 0;
806   FSI_COOR.Z1 = 0;
807   FSI_COOR.T1 = 0;
808
809   FSI_COOR.X2 = rxcm*fgkWcons;
810   FSI_COOR.Y2 = rycm*fgkWcons;
811   FSI_COOR.Z2 = rzcm*fgkWcons;
812   FSI_COOR.T2 = 0;
813
814 }
815
816
817 void AliHBTLLWeights::Rotate(Double_t x, Double_t y, Double_t sphi, Double_t cphi, Double_t& xr, Double_t& yr)
818 {
819 //rotates 
820  yr=y*cphi-x*sphi;
821  xr=x*cphi+y*sphi;
822 }
823
824 void AliHBTLLWeights::Boost(Double_t z, Double_t t, Double_t beta, Double_t gamma, Double_t& zt, Double_t& yt)
825 {
826 //boosts
827   zt=gamma*(z-beta*t);
828   yt=gamma*(t-beta*z);
829 }