Initial check-in of the model classes
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / Model / AliFemtoModelWeightGeneratorBasic.cxx
CommitLineData
75c432a7 1////////////////////////////////////////////////////////////////////////////////
2/// ///
3/// AliFemtoModelWeightGeneratorBasic - basic femtoscopic weight generator ///
4/// only return a simple ///
5/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu ///
6/// ///
7////////////////////////////////////////////////////////////////////////////////
8#ifdef __ROOT__
9 ClassImp(AliFemtoModelWeightGeneratorBasic, 1)
10#endif
11
12#include "AliFemtoModelWeightGeneratorBasic.h"
13#include "AliFemtoModelHiddenInfo.h"
14
15//________________________
16AliFemtoModelWeightGeneratorBasic::AliFemtoModelWeightGeneratorBasic():
17 AliFemtoModelWeightGenerator()
18{
19 /* no-op */
20}
21//________________________
22AliFemtoModelWeightGeneratorBasic::AliFemtoModelWeightGeneratorBasic(const AliFemtoModelWeightGeneratorBasic &aModel) :
23 AliFemtoModelWeightGenerator(aModel)
24{
25 /* no-op */
26}
27
28//________________________
29AliFemtoModelWeightGeneratorBasic::~AliFemtoModelWeightGeneratorBasic()
30{
31 /* no-op */
32}
33
34//________________________
35Double_t AliFemtoModelWeightGeneratorBasic::GenerateWeight(AliFemtoPair *aPair)
36{
37 // Generate a simple femtoscopic weight coming only from
38 // quantum statistics - symmetrization or anti-symmetrization
39 // of the pair wave function
40
41 // Get hidden information pointers
42 AliFemtoModelHiddenInfo *inf1 = (AliFemtoModelHiddenInfo *) aPair->track1()->HiddenInfo();
43 AliFemtoModelHiddenInfo *inf2 = (AliFemtoModelHiddenInfo *) aPair->track2()->HiddenInfo();
44
45 // Calculate pair variables
46 Double_t tPx = inf1->GetTrueMomentum()->x()+inf2->GetTrueMomentum()->x();
47 Double_t tPy = inf1->GetTrueMomentum()->y()+inf2->GetTrueMomentum()->y();
48 Double_t tPz = inf1->GetTrueMomentum()->z()+inf2->GetTrueMomentum()->z();
49 // double tE = inf1->GetTrueMomentum()->e +inf2->GetTrueMomentum()->.e;
50 Double_t tM1 = inf1->GetMass();
51 Double_t tM2 = inf2->GetMass();
52 Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->mag2());
53 Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->mag2());
54 Double_t tE = tE1 + tE2;
55 Double_t tPt = tPx*tPx + tPy*tPy;
56 Double_t tMt = tE*tE - tPz*tPz;//mCVK;
57 Double_t tM = sqrt(tMt - tPt);
58 tMt = sqrt(tMt);
59 tPt = sqrt(tPt);
60 Double_t tBetat = tPt/tMt;
61
62 // Boost to LCMS
63 Double_t tBeta = tPz/tE;
64 Double_t tGamma = tE/tMt;
65 fKStarLong = tGamma * (inf1->GetTrueMomentum()->z() - tBeta * tE1);
66 Double_t tE1L = tGamma * (tE1 - tBeta * inf1->GetTrueMomentum()->z());
67
68 // Transform positions to LCMS
69// Double_t tP1zl = tGamma * (inf1->GetEmissionPoint()->z() - tBeta * inf1->GetEmissionPoint()->t());
70// Double_t tP1tl = tGamma * (inf1->GetEmissionPoint()->t() - tBeta * inf1->GetEmissionPoint()->z());
71
72// Double_t tP2zl = tGamma * (inf2->GetEmissionPoint()->z() - tBeta * inf2->GetEmissionPoint()->t());
73// Double_t tP2tl = tGamma * (inf2->GetEmissionPoint()->t() - tBeta * inf2->GetEmissionPoint()->z());
74
75// Double_t tP1pzl = tGamma * (inf1->GetTrueMomentum()->z() - tBeta * tE1);
76// Double_t tP1el = tGamma * (tE1 - tBeta * inf1->GetTrueMomentum()->z());
77
78// Double_t tP2pzl = tGamma * (inf2->GetTrueMomentum()->z() - tBeta * tE2);
79// Double_t tP2el = tGamma * (tE2 - tBeta * inf2->GetTrueMomentum()->z());
80
81 // Rotate in transverse plane
82 fKStarOut = ( inf1->GetTrueMomentum()->x()*tPx + inf1->GetTrueMomentum()->y()*tPy)/tPt;
83 fKStarSide = (-inf1->GetTrueMomentum()->x()*tPy + inf1->GetTrueMomentum()->y()*tPx)/tPt;
84
85// Double_t tP1pxl = fKStarOut;
86// Double_t tP1pyl = fKStarSide;
87
88// Double_t tP2pxl = (inf2->GetTrueMomentum()->x()*tPx + inf2->GetTrueMomentum()->y()*tPy)/tPt;
89// Double_t tP2pyl = (inf2->GetTrueMomentum()->y()*tPx - inf2->GetTrueMomentum()->x()*tPy)/tPt;;
90
91// Double_t tKO = tP1pxl - tP2pxl;
92// Double_t tKS = tP1pyl - tP2pyl;
93// Double_t tKL = tP1pzl - tP2pzl;
94// Double_t tDE = tP1el - tP2el;
95
96 // save the rotated coordinates in LCMS variables
97// Double_t tP1xl = ( inf1->GetEmissionPoint()->x()*tPx + inf1->GetEmissionPoint()->y()*tPy)/tPt;
98// Double_t tP1yl = (-inf1->GetEmissionPoint()->x()*tPy + inf1->GetEmissionPoint()->y()*tPx)/tPt;
99
100// Double_t tP2xl = ( inf2->GetEmissionPoint()->x()*tPx + inf2->GetEmissionPoint()->y()*tPy)/tPt;
101// Double_t tP2yl = (-inf2->GetEmissionPoint()->x()*tPy + inf2->GetEmissionPoint()->y()*tPx)/tPt;
102
103 // Boost to pair cms
104 fKStarOut = tMt/tM * (fKStarOut - tPt/tMt * tE1L);
105
106 tBetat = tPt/tMt;
107// Double_t tGammat = 1.0/sqrt(1.0-tBetat*tBetat);
108
109// Double_t tP1xp = tGammat*(tP1xl - tBetat*tP1tl);
110// Double_t tP1tp = tGammat*(tP1tl - tBetat*tP1xl);
111
112// Double_t tP2xp = tGammat*(tP2xl - tBetat*tP2tl);
113// Double_t tP2tp = tGammat*(tP2tl - tBetat*tP2xl);
114
115// Double_t tRO = (tP1xl - tP2xl)/0.197327;
116// Double_t tRS = (tP1yl - tP2yl)/0.197327;
117// Double_t tRL = (tP1zl - tP2zl)/0.197327;
118// Double_t tDT = (tP1tl - tP2tl)/0.197327;
119
120 Double_t tDX = inf1->GetEmissionPoint()->x()-inf2->GetEmissionPoint()->x();
121 Double_t tDY = inf1->GetEmissionPoint()->y()-inf2->GetEmissionPoint()->y();
122 Double_t tRLong = inf1->GetEmissionPoint()->z()-inf2->GetEmissionPoint()->z();
123 Double_t tDTime = inf1->GetEmissionPoint()->t()-inf2->GetEmissionPoint()->t();
124
125 Double_t tROut = (tDX*tPx + tDY*tPy)/tPt;
126 Double_t tRSide = (-tDX*tPy + tDY*tPx)/tPt;
127
128 fRStarSide = tRSide;
129 Double_t tRSS = fRStarSide/0.197327;
130
131 fRStarLong = tGamma*(tRLong - tBeta* tDTime);
132 Double_t tDTimePairLCMS = tGamma*(tDTime - tBeta* tRLong);
133
134 Double_t tRLS = fRStarLong/0.197327;
135 tBeta = tPt/tMt;
136 tGamma = tMt/tM;
137
138 fRStarOut = tGamma*(tROut - tBeta* tDTimePairLCMS);
139 Double_t tROS = fRStarOut/0.197327;
140// Double_t tDTimePairCMS = tGamma*(tDTimePairLCMS - tBeta* tROut);
141 fRStar = ::sqrt(fRStarOut*fRStarOut + fRStarSide*fRStarSide +
142 fRStarLong*fRStarLong);
143 fKStar = ::sqrt(fKStarOut*fKStarOut + fKStarSide*fKStarSide + fKStarLong*fKStarLong);
144// Double_t tRSt = fRStar/0.197327;
145
146 if ((fPairType == kPionPlusPionPlus) || (fPairType == kKaonPlusKaonPlus))
147 return 1.0 + cos (2*(fKStarOut * tROS + fKStarSide * tRSS + fKStarLong * tRLS));
148 else if (fPairType == kProtonProton)
149 return 1.0 - 0.5 * cos (2*(fKStarOut * tROS + fKStarSide * tRSS + fKStarLong * tRLS));
150 else
151 return 1.0;
152}
153
154//________________________
155void AliFemtoModelWeightGeneratorBasic::SetPairType(Int_t aPairType)
156{
157 AliFemtoModelWeightGenerator::SetPairType(aPairType);
158}
159//________________________
160void AliFemtoModelWeightGeneratorBasic::SetPairTypeFromPair(AliFemtoPair *aPair)
161{
162 AliFemtoModelWeightGenerator::SetPairTypeFromPair(aPair);
163}
164//________________________
165Int_t AliFemtoModelWeightGeneratorBasic::GetPairType()
166{
167 return AliFemtoModelWeightGenerator::GetPairType();
168}
169//________________________
170AliFemtoModelWeightGenerator* AliFemtoModelWeightGeneratorBasic::Clone() const
171{
172 return GetGenerator();
173}
174//________________________
175AliFemtoModelWeightGenerator* AliFemtoModelWeightGeneratorBasic::GetGenerator() const
176{
177 AliFemtoModelWeightGeneratorBasic *tGen = new AliFemtoModelWeightGeneratorBasic(*this);
178 return tGen;
179}