]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HBTAN/AliHBTLLWeights.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
example analysis macro extended for monitor functions
[u/mrichter/AliRoot.git] / HBTAN / AliHBTLLWeights.cxx
CommitLineData
7f92929e 1#include <iostream.h>
2#include "AliHBTLLWeights.h"
3#include "AliPDG.h"
4#include "AliHBTPair.h"
5#include "AliHBTParticle.h"
6#include <TList.h>
7#include <TRandom.h>
8#include <TMath.h>
9
10/*******************************************************************/
11/****** ROUTINES USED FOR COMMUNUCATION ********/
12/******************** WITH FORTRAN ********************/
13/*******************************************************************/
14#ifndef WIN32
15# define led_bldata led_bldata_
16# define fsiini fsiini_
17# define ltran12 ltran12_
18# define fsiw fsiw_
19# define type_of_call
20#else
21# define led_bldata LED_BLDATA
22# define fsiini FSIINI
23# define ltran12 LTRAN12
24# define fsiw FSIW
25# define type_of_call _stdcall
26#endif
27/****************************************************************/
28extern "C" void type_of_call led_bldata();
29extern "C" void type_of_call fsiini();
30extern "C" void type_of_call ltran12();
31extern "C" void type_of_call fsiw();
32/**************************************************************/
33
34ClassImp(AliHBTLLWeights)
35
36//Interface to Richard Lednicky's weight alghorithm (Fortran implementation)
37//Ludmila Malinina JINR
38//Piotr Krzysztof Skowronski CERN/WUT
39
40AliHBTLLWeights* AliHBTLLWeights::fgLLWeights=NULL;
41
42AliHBTLLWeights::AliHBTLLWeights()
43{
44// Default Constructor
45 fPID1 = 0;
46 fPID2 = 0;
47 SetRandomPosition();
48 SetColWithResidNuclSwitch();
49 SetStrongInterSwitch();
50 SetQuantumStatistics();
51 SetColoumb();
52 SetTest();
53}
54
55
56 AliHBTLLWeights* AliHBTLLWeights::Instance()
57{
58 if (fgLLWeights)
59 {
60 return fgLLWeights;
61 }
62 else
63 {
64 fgLLWeights = new AliHBTLLWeights();
65 return fgLLWeights;
66 }
67}
68
69
70Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair)
71{
72 AliHBTParticle *part1 = partpair->Particle1();
73 AliHBTParticle *part2 = partpair->Particle2();
74
75 if ( (part1 == 0x0) || (part2 == 0x0))
76 {
77 Error("GetWeight","Null particle pointer");
78 return 0.0;
79 }
80
81 Double_t part1Momentum[]={part1->Px(),part1->Py(),part1->Pz()};
82 Double_t part2Momentum[]={part2->Px(),part2->Py(),part2->Pz()};
83
84
85 if ( (part1->Px() == part2->Px()) && (part1->Py() == part2->Py())
86 && (part1->Pz() == part2->Pz()) )
87 {//if particles have the same momentum
88 return 0.0;
89 }
90
91
92 if ((!fRandomPosition) && (part1->Vx() == part2->Vx()) && (part1->Vy() == part2->Vy())
93 && (part1->Vz() == part2->Vz()) )
94 { //if particles have the same position
95 return 0.0;
96 }
97
98//put momenta of particles in LAB into Fortran commons
99
100 FSI_MOM.P1X=part1Momentum[0];
101 FSI_MOM.P1Y=part1Momentum[1];
102 FSI_MOM.P1Z=part1Momentum[2];
103
104 FSI_MOM.P2X=part2Momentum[0];
105 FSI_MOM.P2Y=part2Momentum[1];
106 FSI_MOM.P2Z=part2Momentum[2];
107
108 if (fRandomPosition)
109 {
110 Double_t rxcm = fsigma*gRandom->Gaus();
111 Double_t rycm = fsigma*gRandom->Gaus();
112 Double_t rzcm = fsigma*gRandom->Gaus();
113
114 FSI_PRF.X=rxcm*fwcons;
115 FSI_PRF.Y=rycm*fwcons;
116 FSI_PRF.Z=rzcm*fwcons;
117 FSI_PRF.T=0.;
118
119 Double_t rps=rxcm*rxcm+rycm*rycm+rzcm*rzcm;
120 Double_t rp=TMath::Sqrt(rps);
121 FSI_PRF.RP=rp;
122 FSI_PRF.RPS=rps;
123 }
124
125 ltran12();
126 fsiw();
127 return LEDWEIGHT.WEIN;
128 }
129
130/************************************************************/
131void AliHBTLLWeights::Init()
132 {
133//---------------------------------------------------------------------
134
135//initial parameters of model
136
137 FSI_NS.NS = approximationModel;
138
139 if(!ftest) LEDWEIGHT.ITEST=0;
140
141 if(ftest)
142 {
143 LEDWEIGHT.ITEST=1;
144 if(fColoumbSwitch) FSI_NS.ICH =1;
145 else FSI_NS.ICH=0;
146
147 if(fStrongInterSwitch) FSI_NS.ISI=1;
148 else FSI_NS.ISI=0;
149
150 if(fQuantStatSwitch) FSI_NS.IQS=1;
151 else FSI_NS.IQS=0;
152
153 if(fColWithResidNuclSwitch) FSI_NS.I3C=1;
154 else FSI_NS.I3C=0;
155 }
156
157 if(fRandomPosition) LEDWEIGHT.IRANPOS=1;
158 else LEDWEIGHT.IRANPOS=0;
159
160 if ( (fPID1 == 0) || (fPID2 == 0) )
161 {
162 Fatal("Init","Particles types are not set");
163 return;//pro forma
164 }
165
166 FSI_NS.LL = GetPairCode(fPID1,fPID2);
167
168 if (FSI_NS.LL == 0)
169 {
170 Fatal("Init","Particles types are not supported");
171 return;//pro forma
172 }
173
174 TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
175 if (tpart1 == 0x0)
176 {
177 Fatal("init","We can not find particle with ID=%d in our DataBase",fPID1);
178 return;
179 }
180
181 FSI_POC.AM1=tpart1->Mass();
182 FSI_POC.C1=tpart1->Charge();
183
184 TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
185 if (tpart2 == 0x0)
186 {
187 Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2);
188 return;
189 }
190
191 FSI_POC.AM2=tpart2->Mass();
192 FSI_POC.C1=tpart2->Charge();
193
194 led_bldata();
195 fsiini();
196
197
198//constants for radii simulation
199
200 if(fRandomPosition)
201 {
202 fsigma =TMath::Sqrt(2.)*fRadius;
203 fwcons =FSI_CONS.W;
204 }
205}
206
207Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
208{
209 return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
210}
211
212Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
213{
214// 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
215// 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
216// 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
217// NS=1 y/n: + + + + + - - - - - - - - - - - - - - - - - - - - - - -
218
219//alphas, deuterons and tyts are NOT supported here
220
221 Int_t chargefactor = 1;
222 Int_t hpid; //pid in higher row
223 Int_t lpid; //pid in lower row
224 Int_t code; //pairCode
225
226 Bool_t swap;
227
228//determine the order of selcetion in switch
229 if (TMath::Abs(pid1) < TMath::Abs(pid2) )
230 {
231 if (pid1<0) chargefactor=-1;
232 hpid=pid2*chargefactor;
233 lpid=pid1*chargefactor;
234 swap = kFALSE;
235 }
236 else
237 {
238 if (pid2<0) chargefactor=-1;
239 hpid=pid1*chargefactor;
240 lpid=pid2*chargefactor;
241 swap = kTRUE;
242 }
243
244//Determine the pair code
245 switch (hpid) //switch on first particle id
246 {
247 case kNeutron:
248 switch (lpid)
249 {
250 case kNeutron:
251 code = 1; //neutron neutron
252 break;
253
254 case kProton:
255 code = 3; //neutron proton
256 break;
257
258 case kLambda0:
259 code = 28; //neutron lambda
260 break;
261
262 default:
263 return 0; //given pair not supported
264 break;
265 }
266 break;
267
268 case kProton:
269 switch (lpid)
270 {
271 case kProton:
272 code = 2; //proton proton
273 break;
274
275 case kLambda0:
276 code = 27;//proton lambda
277 break;
278
279 default:
280 return 0; //given pair not supported
281 break;
282
283 }
284 break;
285
286 case kPiPlus:
287 switch (lpid)
288 {
289 case kPiPlus:
290 code = 7; //piplus piplus
291 break;
292
293 case kPiMinus:
294 code = 5; //piplus piminus
295 break;
296
297 case kKMinus:
298 code = 10; //piplus Kminus
299 break;
300
301 case kKPlus:
302 code = 11; //piplus Kplus
303 break;
304
305 case kProton:
306 code = 12; //piplus proton
307 chargefactor*=-1;
308 break;
309
310 default:
311 return 0; //given pair not supported
312 break;
313 }
314 break;
315 case kPi0:
316 switch (lpid)
317 {
318 case kPi0:
319 code = 6;
320 break;
321
322 default:
323 return 0; //given pair not supported
324 break;
325 }
326 break;
327
328 case kKPlus:
329 switch (lpid)
330 {
331 case kKMinus:
332 code = 14; //Kplus Kminus
333 break;
334
335 case kKPlus:
336 code = 15; //Kplus Kplus
337 break;
338
339 case kProton:
340 code = 16; //Kplus proton
341 break;
342
343 default:
344 return 0; //given pair not supported
345 break;
346 }
347 break;
348
349 case kKMinus:
350 switch (lpid)
351 {
352 case kProton:
353 code = 17; //Kminus proton
354 chargefactor*=1;
355 break;
356
357 default:
358 return 0; //given pair not supported
359 break;
360 }
361 break;
362
363 case kK0:
364 switch (lpid)
365 {
366 case kK0:
367 code = 2; //Kzero Kzero
368 break;
369
370 case kK0Bar:
371 code = 17; //Kzero KzeroBar
372 break;
373
374 default:
375 return 0; //given pair not supported
376 break;
377 }
378 break;
379
380 default: return 0;
381 }
382 return code;
383}
384