Cosmetic corrections
[u/mrichter/AliRoot.git] / HBTAN / AliHBTLLWeights.cxx
1 #include <Riostream.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 /****************************************************************/
28 extern "C" void type_of_call led_bldata(); 
29 extern "C" void type_of_call fsiini();
30 extern "C" void type_of_call ltran12();
31 extern "C" void type_of_call fsiw();
32 /**************************************************************/
33
34 ClassImp(AliHBTLLWeights)  
35
36 //Interface to Richard Lednicky's weight alghorithm (Fortran implementation)
37 //Ludmila Malinina JINR
38 //Piotr Krzysztof Skowronski CERN/WUT
39  
40 AliHBTLLWeights* AliHBTLLWeights::fgLLWeights=NULL; 
41
42 AliHBTLLWeights::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
70 Double_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 /************************************************************/
131 void 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
207 Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
208 {
209  return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
210 }
211
212 Int_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