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