Making the directory structure of AliFemtoUser flat. All files go into one common...
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoModelGausRinvFreezeOutGenerator.cxx
1 ////////////////////////////////////////////////////////////////////////////////
2 ///                                                                          ///
3 /// AliFemtoModelGausRinvFreezeOutGenerator - freeze-out                     ///
4 /// coordinates generator, generating a 3D gaussian ellipsoid in LCMS        ///
5 /// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
6 ///                                                                          ///
7 ////////////////////////////////////////////////////////////////////////////////
8 #ifdef __ROOT__
9   ClassImp(AliFemtoModelGausRinvFreezeOutGenerator, 1)
10 #endif
11
12 #include "math.h"
13 #include "AliFemtoModelGausRinvFreezeOutGenerator.h"
14 #include "AliFemtoModelHiddenInfo.h"
15 #include "AliFemtoLorentzVector.h"
16
17 //_______________________
18 AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator() :
19   fSizeInv(0)
20 {
21   // Default constructor
22   fRandom = new TRandom2();
23 }
24
25 //_______________________
26 AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator(const AliFemtoModelGausRinvFreezeOutGenerator &aModel):
27   fSizeInv(0)
28 {
29   // Copy constructor
30   fRandom = new TRandom2();
31   SetSizeInv(aModel.GetSizeInv());
32 }
33 //_______________________
34 AliFemtoModelGausRinvFreezeOutGenerator::~AliFemtoModelGausRinvFreezeOutGenerator()
35 {
36   if (fRandom) delete fRandom;
37 }
38 //_______________________
39 void AliFemtoModelGausRinvFreezeOutGenerator::GenerateFreezeOut(AliFemtoPair *aPair)
40 {
41   // Generate two particle emission points with respect
42   // to their pair momentum 
43   // The source is the 3D Gaussian ellipsoid in the LCMS frame
44   AliFemtoModelHiddenInfo *inf1 = (AliFemtoModelHiddenInfo *) aPair->Track1()->HiddenInfo();
45   AliFemtoModelHiddenInfo *inf2 = (AliFemtoModelHiddenInfo *) aPair->Track2()->HiddenInfo();
46
47   if ((!inf1) || (!inf2)) { cout << "Hidden info not created! "  << endl; exit(kFALSE); }
48
49   // Calculate sum momenta
50   Double_t tPx = inf1->GetTrueMomentum()->x() + inf2->GetTrueMomentum()->x();
51   Double_t tPy = inf1->GetTrueMomentum()->y() + inf2->GetTrueMomentum()->y();
52   Double_t tPz = inf1->GetTrueMomentum()->z() + inf2->GetTrueMomentum()->z();
53   Double_t tM1 = inf1->GetMass();
54   Double_t tM2 = inf2->GetMass();
55   Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->mag2());
56   Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->mag2());
57   Double_t tEs = tE1 + tE2;
58
59   Double_t tPt = sqrt(tPx*tPx + tPy*tPy);
60   Double_t tMt = sqrt(tEs*tEs - tPz*tPz);
61
62   // Generate positions in PRF from a Gaussian
63   Double_t tROutS =  fRandom->Gaus(0,fSizeInv); // reuse of long
64   Double_t tRSideS = fRandom->Gaus(0,fSizeInv);
65   Double_t tRLongS = fRandom->Gaus(0,fSizeInv);
66   Double_t tRTimeS = 0;
67       
68   Double_t tBetat  = tPt/tMt;
69   Double_t tGammat = 1.0/sqrt(1.0-tBetat*tBetat);
70
71   Double_t tBetaz  = tPz/tEs;
72   Double_t tGammaz = 1.0/sqrt(1.0-tBetaz*tBetaz);
73
74   Double_t tROut = tGammat * (tROutS + tBetat * tRTimeS);
75   Double_t tDtL  = tGammat * (tRTimeS + tBetat * tROutS);
76   Double_t tRSide = tRSideS;
77
78   Double_t tRLong = tGammaz * (tRLongS + tBetaz * tDtL);
79   Double_t tDt    = tGammaz * (tDtL + tBetaz * tRLongS);
80           
81   tPx /= tPt;
82   tPy /= tPt;
83           
84   Double_t tXout  = tROut*tPx-tRSide*tPy;
85   Double_t tXside = tROut*tPy+tRSide*tPx;
86   Double_t tXlong = tRLong;
87   Double_t tXtime = tDt;
88   
89   if (!(inf1->GetEmissionPoint())) {
90     AliFemtoLorentzVector tPos(0,0,0,0);
91     inf1->SetEmissionPoint(&tPos);
92   }
93   else
94     inf1->SetEmissionPoint(0,0,0,0);
95   if (!(inf2->GetEmissionPoint())) {
96     AliFemtoLorentzVector tPos(tXout,tXside,tXlong,tXtime);
97     inf2->SetEmissionPoint(&tPos);
98   }
99   else
100     inf2->SetEmissionPoint(tXout, tXside, tXlong, tXtime);
101 }
102
103 //_______________________
104 void AliFemtoModelGausRinvFreezeOutGenerator::SetSizeInv(Double_t aSizeInv)
105 {
106   fSizeInv = aSizeInv;
107 }
108 //_______________________
109 Double_t AliFemtoModelGausRinvFreezeOutGenerator::GetSizeInv() const
110 {
111   return fSizeInv;
112 }
113 //_______________________
114 AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::Clone() const
115
116   return GetGenerator(); 
117 }
118 //_______________________
119 inline AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::GetGenerator() const 
120
121   AliFemtoModelFreezeOutGenerator* tModel = new AliFemtoModelGausRinvFreezeOutGenerator(*this); 
122   return tModel; 
123 }