Merging correctly the changes with the previous code clean-up
[u/mrichter/AliRoot.git] / HBTAN / AliHBTLLWeights.cxx
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 //------------------------------------------------------------------
11
12 #include <TMath.h>
13 #include <TPDGCode.h>
14 #include <TRandom.h>
15
16 #include "AliHBTLLWeights.h"
17 #include "AliHBTPair.h"
18 #include "AliHBTParticle.h"
19 #include "WLedCOMMONS.h"
20
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 /****************************************************************/
39 extern "C" void type_of_call led_bldata(); 
40 extern "C" void type_of_call fsiini();
41 extern "C" void type_of_call ltran12();
42 extern "C" void type_of_call fsiw();
43 /**************************************************************/
44
45 ClassImp(AliHBTLLWeights)  
46  
47 AliHBTLLWeights* AliHBTLLWeights::fgLLWeights=NULL; 
48
49 AliHBTLLWeights::AliHBTLLWeights()
50 {
51   // Default Constructor 
52   fPID1 = 0;
53   fPID2 = 0;
54   SetRandomPosition();
55   SetColWithResidNuclSwitch();
56   SetStrongInterSwitch();
57   SetQuantumStatistics();
58   SetColoumb();
59   SetTest();
60   
61 }
62
63
64 AliHBTLLWeights* 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 }
74                                             
75
76 Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair)
77 {
78   // Returns the weignt of the pair "partpair"
79   AliHBTParticle *part1 = partpair->Particle1();
80   AliHBTParticle *part2 = partpair->Particle2();
81
82   if ( (part1 == 0x0) || (part2 == 0x0))
83     {
84       Error("GetWeight","Null particle pointer");
85       return 0.0;
86     }
87   
88    
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()) )
95     {
96       return 0.0;
97     }
98   
99              
100   if ((!fRandomPosition) && 
101       (part1->Vx()  == part2->Vx()) && (part1->Vy()  == part2->Vy())
102       && (part1->Vz()  == part2->Vz()) )
103     {        
104       return 0.0;
105     }
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();
137
138   if(flambda<1){
139     if(gRandom->Rndm()<(1-flambda))LEDWEIGHT.WEIN=1.;}
140   
141   return LEDWEIGHT.WEIN;
142 }
143
144 /************************************************************/
145 void AliHBTLLWeights::Init()
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   } 
219
220
221 Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
222 {
223   // Return the code of the pair "partpair"
224   return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
225 }
226
227 Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
228 {
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   
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   
243   //determine the order of selcetion in switch  
244   if (TMath::Abs(pid1) < TMath::Abs(pid2) ) 
245     {
246       if (pid1<0) chargefactor=-1;
247       hpid=pid2*chargefactor;
248       lpid=pid1*chargefactor;
249       swap = kFALSE;
250     } 
251   else 
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
265   switch (hpid) //switch on first  particle id
266     {
267     case kNeutron:
268       switch (lpid)
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         }
286       break;
287       
288     case kProton:
289       switch (lpid)
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         }
304       break;
305       
306     case kPiPlus:
307       
308       switch (lpid)
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         }
335       break;
336     case kPi0:
337       switch (lpid)
338         {
339         case kPi0:
340           code = 6;
341           break;
342           
343         default: 
344           return 0; //given pair not supported
345           break;
346         }
347       break;
348       
349     case kKPlus:
350       switch (lpid)
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         }
368       break;
369       
370     case kKMinus:
371       switch (lpid)
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         }
382       break;
383       
384     case kK0:
385       switch (lpid)
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         }
399       break;
400       
401     default: return 0;
402     }
403   return code;
404 }
405