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