Update by Ludmila Malinina
[u/mrichter/AliRoot.git] / HBTAN / AliHBTLLWeights.cxx
1 #include "AliHBTLLWeights.h"
2 #include "AliPDG.h"
3 #include "AliHBTPair.h"
4 #include "AliHBTParticle.h"
5 #include <TList.h>
6 #include <TRandom.h>                                                                     
7 #include <TMath.h>                                                                       
8
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         return LEDWEIGHT.WEIN;
125
126  }
127  
128 /************************************************************/
129 void AliHBTLLWeights::Init()
130  {
131 //---------------------------------------------------------------------
132
133 //initial parameters of model
134
135       FSI_NS.NS = approximationModel;      
136       
137       if(!ftest){LEDWEIGHT.ITEST=0;}
138     
139       if(ftest){
140       LEDWEIGHT.ITEST=1;
141       if(fColoumbSwitch){FSI_NS.ICH =1;}
142       else{FSI_NS.ICH=0;}
143       if(fStrongInterSwitch){FSI_NS.ISI=1;}
144       else{FSI_NS.ISI=0;}
145       if(fQuantStatSwitch){FSI_NS.IQS=1;}
146       else{FSI_NS.IQS=0;}
147       if(fColWithResidNuclSwitch){FSI_NS.I3C=1;}
148       else{FSI_NS.I3C=0;}
149       }
150       
151       if(fRandomPosition){LEDWEIGHT.IRANPOS=1;}
152       else{LEDWEIGHT.IRANPOS=0;}
153
154
155      if ( (fPID1 == 0) || (fPID2 == 0) )
156       {
157         Fatal("Init","Particles types are not set");
158         return;//pro forma
159       }
160
161       FSI_NS.LL = GetPairCode(fPID1,fPID2);
162        
163       if (FSI_NS.LL == 0) 
164        {
165          Fatal("Init","Particles types are not supported");
166          return;//pro forma
167        }
168
169
170       TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
171       if (tpart1 == 0x0)
172        {
173          Fatal("init","We can not find particle with ID=%d in our DataBase",fPID1);
174          return;
175        }
176       
177       FSI_POC.AM1=tpart1->Mass();
178       FSI_POC.C1=tpart1->Charge(); 
179
180       TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
181 //mlv
182       
183       
184
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        fsigma =TMath::Sqrt(2.)*fRadius;     
202        fwcons =FSI_CONS.W;
203       } 
204
205
206 Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
207 {
208  return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
209 }
210
211 Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
212 {
213 //   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
214 //   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
215 //   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
216 //   NS=1 y/n:  +  +  +  +   +    -   -  -  -  -   -   -   -  -  -  -  -   -  -    -    -    -  -   - -  -      -      -
217
218 //alphas, deuterons and tyts are NOT supported here
219
220   Int_t chargefactor = 1;
221   Int_t hpid; //pid in higher row
222   Int_t lpid; //pid in lower row
223   Int_t code; //pairCode
224   
225   Bool_t swap;
226   
227 //determine the order of selcetion in switch  
228   if (TMath::Abs(pid1) < TMath::Abs(pid2) ) 
229    {
230     if (pid1<0) chargefactor=-1;
231     hpid=pid2*chargefactor;
232     lpid=pid1*chargefactor;
233     swap = kFALSE;
234    } 
235   else 
236    {
237     if (pid2<0) chargefactor=-1;
238     hpid=pid1*chargefactor;
239     lpid=pid2*chargefactor;
240     swap = kTRUE;
241    }
242
243 //mlv
244    hpid=pid1;
245    lpid=pid2;
246
247
248 //Determine the pair code
249   switch (hpid) //switch on first  particle id
250    {
251      case kNeutron:
252       switch (lpid)
253        {
254          case kNeutron: 
255            code = 1;  //neutron neutron
256            break;
257         
258          case kProton: 
259            code = 3;  //neutron proton
260            break;
261            
262          case kLambda0: 
263            code = 28;  //neutron lambda
264            break;
265            
266          default: 
267            return 0; //given pair not supported
268            break;
269        }
270       break;
271
272      case kProton:
273       switch (lpid)
274        {
275          case kProton:
276            code = 2; //proton proton
277            break;
278            
279          case kLambda0: 
280            code = 27;//proton lambda
281            break;
282            
283          default: 
284            return 0; //given pair not supported
285            break;
286            
287        }
288       break;
289
290      case kPiPlus:
291      
292       switch (lpid)
293        {
294          case kPiPlus:
295            code = 7; //piplus piplus
296            break;
297
298          case kPiMinus:
299            code = 5; //piplus piminus
300            break;
301         
302          case kKMinus:
303            code = 10; //piplus Kminus
304            break;
305
306          case kKPlus:
307            code = 11; //piplus Kplus
308            break;
309
310          case kProton:
311            code = 12; //piplus proton
312            chargefactor*=-1;
313            break;
314
315          default: 
316            return 0; //given pair not supported
317            break;
318        }
319       break;
320      case kPi0:
321       switch (lpid)
322        {
323          case kPi0:
324            code = 6;
325            break;
326            
327          default: 
328            return 0; //given pair not supported
329            break;
330        }
331       break;
332       
333      case kKPlus:
334       switch (lpid)
335        {
336          case kKMinus:
337            code = 14; //Kplus Kminus
338            break;
339
340          case kKPlus:
341            code = 15; //Kplus Kplus
342            break;
343
344          case kProton:
345            code = 16; //Kplus proton
346            break;
347            
348          default: 
349            return 0; //given pair not supported
350            break;
351        }
352       break;
353       
354      case kKMinus:
355       switch (lpid)
356        {
357          case kProton:
358            code = 17; //Kminus proton
359            chargefactor*=1;
360            break;
361            
362          default: 
363            return 0; //given pair not supported
364            break;
365        }
366       break;
367       
368      case kK0:
369       switch (lpid)
370        {
371          case kK0:
372            code = 2; //Kzero Kzero
373            break;
374          
375          case kK0Bar:
376            code = 17; //Kzero KzeroBar
377            break;
378
379          default: 
380            return 0; //given pair not supported
381            break;
382        }
383       break;
384
385      default: return 0;
386    }
387   return code;
388 }
389