]>
Commit | Line | Data |
---|---|---|
4fdf4eb3 | 1 | #include "AliHBTLLWeights.h" |
88cb7938 | 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) | |
c81f9591 | 86 | // and ializes various parameters. |
88cb7938 | 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 | |
1a1e58ac | 93 | // SPI=TMath::Sqrt(PI) |
88cb7938 | 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 | ||
4fdf4eb3 | 138 | |
7f92929e | 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_ | |
1a1e58ac | 147 | # define boosttoprf boosttoprf_ |
7f92929e | 148 | # define fsiw fsiw_ |
88cb7938 | 149 | # define setpdist setpdist_ |
7f92929e | 150 | # define type_of_call |
151 | #else | |
152 | # define led_bldata LED_BLDATA | |
153 | # define fsiini FSIINI | |
154 | # define ltran12 LTRAN12 | |
1a1e58ac | 155 | # define boosttoprf BOOSTTOPRF |
7f92929e | 156 | # define fsiw FSIW |
88cb7938 | 157 | # define setpdist SETPDIST |
7f92929e | 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(); | |
1a1e58ac | 164 | extern "C" void type_of_call boosttoprf(); |
7f92929e | 165 | extern "C" void type_of_call fsiw(); |
88cb7938 | 166 | extern "C" void type_of_call setpdist(Double_t& r); |
7f92929e | 167 | /**************************************************************/ |
168 | ||
88cb7938 | 169 | #include "AliHBTPair.h" |
78d7c6d3 | 170 | #include "AliVAODParticle.h" |
88cb7938 | 171 | #include "WLedCOMMONS.h" |
88cb7938 | 172 | #include <TRandom.h> |
173 | #include <TMath.h> | |
174 | #include <TPDGCode.h> | |
78d7c6d3 | 175 | #include <TParticlePDG.h> |
176 | #include <TDatabasePDG.h> | |
88cb7938 | 177 | |
178 | ||
7f92929e | 179 | ClassImp(AliHBTLLWeights) |
7f92929e | 180 | |
88cb7938 | 181 | AliHBTLLWeights* AliHBTLLWeights::fgLLWeights = 0x0; |
182 | const Double_t AliHBTLLWeights::fgkWcons = 1./0.1973; | |
7f92929e | 183 | |
88cb7938 | 184 | AliHBTLLWeights::AliHBTLLWeights(): |
185 | fTest(kTRUE), | |
186 | fColoumbSwitch(kTRUE), | |
187 | fQuantStatSwitch(kTRUE), | |
188 | fStrongInterSwitch(kTRUE), | |
111e505b | 189 | fColWithResidNuclSwitch(kFALSE), |
88cb7938 | 190 | fNuclMass(0.0), |
191 | fNuclCharge(0.0), | |
1a1e58ac | 192 | fRandomPosition(kNone), |
88cb7938 | 193 | fRadius(0.0), |
526c2bd5 | 194 | fOneMinusLambda(0.0), |
88cb7938 | 195 | fPID1(0), |
196 | fPID2(0), | |
197 | fSigma(0.0) | |
7f92929e | 198 | { |
88cb7938 | 199 | // Default Constructor |
ae4a4473 | 200 | if (fgLLWeights) |
201 | Fatal("AliHBTLLWeights","LLWeights already instatiated. Use AliHBTLLWeights::Instance()"); | |
7f92929e | 202 | } |
88cb7938 | 203 | /**************************************************************/ |
7f92929e | 204 | |
111e505b | 205 | AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/): |
dd82cadc | 206 | AliHBTWeights(), |
111e505b | 207 | fTest(kTRUE), |
208 | fColoumbSwitch(kTRUE), | |
209 | fQuantStatSwitch(kTRUE), | |
210 | fStrongInterSwitch(kTRUE), | |
211 | fColWithResidNuclSwitch(kFALSE), | |
212 | fNuclMass(0.0), | |
213 | fNuclCharge(0.0), | |
1a1e58ac | 214 | fRandomPosition(kNone), |
111e505b | 215 | fRadius(0.0), |
526c2bd5 | 216 | fOneMinusLambda(0.0), |
111e505b | 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 | ||
4fdf4eb3 | 234 | AliHBTLLWeights* AliHBTLLWeights::Instance() |
88cb7938 | 235 | { |
236 | // returns instance of class | |
237 | if (fgLLWeights) | |
238 | { | |
4fdf4eb3 | 239 | return fgLLWeights; |
88cb7938 | 240 | } |
241 | else | |
242 | { | |
243 | fgLLWeights = new AliHBTLLWeights(); | |
244 | return fgLLWeights; | |
245 | } | |
246 | } | |
526c2bd5 | 247 | /************************************************************/ |
7f92929e | 248 | |
dd82cadc | 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 | ||
1a1e58ac | 265 | Double_t AliHBTLLWeights::GetWeight(AliHBTPair* partpair) |
7f92929e | 266 | { |
88cb7938 | 267 | // calculates weight for a pair |
37e71815 | 268 | static const Double_t kcmtofm = 1.e13; |
269 | static const Double_t kcmtoOneOverGeV = kcmtofm*fgkWcons; | |
88cb7938 | 270 | |
78d7c6d3 | 271 | AliVAODParticle *part1 = partpair->Particle1(); |
272 | AliVAODParticle *part2 = partpair->Particle2(); | |
2f8eea63 | 273 | |
4fdf4eb3 | 274 | if ( (part1 == 0x0) || (part2 == 0x0)) |
88cb7938 | 275 | { |
276 | Error("GetWeight","Null particle pointer"); | |
277 | return 0.0; | |
278 | } | |
279 | ||
72e72d00 | 280 | if ( fPID1 != part1->GetPdgCode() ) return 1.0; |
281 | if ( fPID2 != part2->GetPdgCode() ) return 1.0; | |
88cb7938 | 282 | |
526c2bd5 | 283 | //takes a lot of time |
4fdf4eb3 | 284 | if ( (part1->Px() == part2->Px()) && |
285 | (part1->Py() == part2->Py()) && | |
286 | (part1->Pz() == part2->Pz()) ) | |
88cb7938 | 287 | { |
288 | return 0.0; | |
289 | } | |
290 | ||
4fdf4eb3 | 291 | if ((!fRandomPosition) && |
88cb7938 | 292 | (part1->Vx() == part2->Vx()) && |
293 | (part1->Vy() == part2->Vy()) && | |
294 | (part1->Vz() == part2->Vz()) ) | |
2f8eea63 | 295 | { |
4fdf4eb3 | 296 | return 0.0; |
2f8eea63 | 297 | } |
88cb7938 | 298 | |
526c2bd5 | 299 | if(fOneMinusLambda)//implemetation of non-zero intetcept parameter |
300 | { | |
301 | if( gRandom->Rndm() < fOneMinusLambda ) return 1.0; | |
302 | } | |
303 | ||
7f92929e | 304 | |
88cb7938 | 305 | //this must be after ltran12 because it would overwrite what we set below |
1a1e58ac | 306 | switch (fRandomPosition) |
88cb7938 | 307 | { |
1a1e58ac | 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; | |
88cb7938 | 370 | } |
1a1e58ac | 371 | |
88cb7938 | 372 | fsiw(); |
4fdf4eb3 | 373 | return LEDWEIGHT.WEIN; |
374 | } | |
7f92929e | 375 | /************************************************************/ |
88cb7938 | 376 | |
7f92929e | 377 | void AliHBTLLWeights::Init() |
4fdf4eb3 | 378 | { |
88cb7938 | 379 | //initial parameters of model |
380 | ||
4fdf4eb3 | 381 | FSI_NS.NS = fApproximationModel; |
382 | ||
88cb7938 | 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 | ||
4fdf4eb3 | 393 | if ( (fPID1 == 0) || (fPID2 == 0) ) |
88cb7938 | 394 | { |
395 | Fatal("Init","Particles types are not set"); | |
396 | return;//pro forma | |
397 | } | |
c81f9591 | 398 | |
399 | ||
4fdf4eb3 | 400 | FSI_NS.LL = GetPairCode(fPID1,fPID2); |
88cb7938 | 401 | |
4fdf4eb3 | 402 | if (FSI_NS.LL == 0) |
88cb7938 | 403 | { |
404 | Fatal("Init","Particles types are not supported"); | |
405 | return;//pro forma | |
406 | } | |
407 | ||
c81f9591 | 408 | Info("Init","Setting PIDs %d %d. LL Code is %d",fPID1,fPID2,FSI_NS.LL); |
409 | ||
88cb7938 | 410 | |
4fdf4eb3 | 411 | TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1); |
412 | if (tpart1 == 0x0) | |
88cb7938 | 413 | { |
414 | Fatal("init","We can not find particle with ID=%d in PDG DataBase",fPID1); | |
415 | return; | |
416 | } | |
417 | ||
4fdf4eb3 | 418 | FSI_POC.AM1=tpart1->Mass(); |
419 | FSI_POC.C1=tpart1->Charge(); | |
88cb7938 | 420 | |
4fdf4eb3 | 421 | TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2); |
88cb7938 | 422 | //lv |
4fdf4eb3 | 423 | if (tpart2 == 0x0) |
88cb7938 | 424 | { |
425 | Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2); | |
426 | return; | |
427 | } | |
428 | ||
4fdf4eb3 | 429 | FSI_POC.AM2=tpart2->Mass(); |
430 | FSI_POC.C1=tpart2->Charge(); | |
88cb7938 | 431 | |
4fdf4eb3 | 432 | led_bldata(); |
433 | fsiini(); | |
88cb7938 | 434 | |
435 | ||
436 | //constants for radii simulation | |
437 | ||
438 | if(fRandomPosition) | |
439 | { | |
440 | fSigma =TMath::Sqrt(2.)*fRadius; | |
441 | } | |
7f92929e | 442 | } |
88cb7938 | 443 | /************************************************************/ |
7f92929e | 444 | |
445 | Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair) | |
446 | { | |
88cb7938 | 447 | //returns Code corresponding to that pair |
448 | return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode()); | |
7f92929e | 449 | } |
88cb7938 | 450 | /************************************************************/ |
7f92929e | 451 | |
452 | Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2) | |
453 | { | |
88cb7938 | 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 | ||
7f92929e | 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 | ||
88cb7938 | 469 | //determine the order of selcetion in switch |
7f92929e | 470 | if (TMath::Abs(pid1) < TMath::Abs(pid2) ) |
88cb7938 | 471 | { |
472 | if (pid1<0) chargefactor=-1; | |
473 | hpid=pid2*chargefactor; | |
474 | lpid=pid1*chargefactor; | |
475 | swap = kFALSE; | |
476 | } | |
7f92929e | 477 | else |
88cb7938 | 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 | |
7f92929e | 491 | switch (hpid) //switch on first particle id |
88cb7938 | 492 | { |
493 | case kNeutron: | |
7f92929e | 494 | switch (lpid) |
88cb7938 | 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 | } | |
7f92929e | 512 | break; |
88cb7938 | 513 | |
514 | case kProton: | |
7f92929e | 515 | switch (lpid) |
88cb7938 | 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 | } | |
7f92929e | 530 | break; |
88cb7938 | 531 | |
532 | case kPiPlus: | |
533 | ||
7f92929e | 534 | switch (lpid) |
88cb7938 | 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 | } | |
7f92929e | 561 | break; |
88cb7938 | 562 | case kPi0: |
7f92929e | 563 | switch (lpid) |
88cb7938 | 564 | { |
565 | case kPi0: | |
566 | code = 6; | |
567 | break; | |
568 | ||
569 | default: | |
570 | return 0; //given pair not supported | |
571 | break; | |
572 | } | |
7f92929e | 573 | break; |
574 | ||
88cb7938 | 575 | case kKPlus: |
7f92929e | 576 | switch (lpid) |
88cb7938 | 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 | } | |
7f92929e | 594 | break; |
595 | ||
88cb7938 | 596 | case kKMinus: |
7f92929e | 597 | switch (lpid) |
88cb7938 | 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 | } | |
7f92929e | 608 | break; |
609 | ||
88cb7938 | 610 | case kK0: |
7f92929e | 611 | switch (lpid) |
88cb7938 | 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 | } | |
7f92929e | 625 | break; |
88cb7938 | 626 | |
627 | default: return 0; | |
628 | } | |
7f92929e | 629 | return code; |
630 | } | |
88cb7938 | 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 | /************************************************************/ | |
7f92929e | 661 | |
88cb7938 | 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 | ||
1a1e58ac | 676 | void AliHBTLLWeights::SetRandomPosition(ERandomizationWay rw) |
88cb7938 | 677 | { |
1a1e58ac | 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; | |
88cb7938 | 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 | } | |
1a1e58ac | 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 | } |