Lines getting the matched track moved to a method in AliCalorimeterUtils. Lines copie...
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoModelWeightGeneratorLednicky.cxx
1 ///////////////////////////////////////////////////////////////////////////
2 //                                                                       //
3 // AliFemtoModelWeightGeneratorLednicky : the most advanced weight       //
4 // generator available. Supports a large number of different pair types  //
5 // and interaction types. Can calculate pair weights coming from         //
6 // quantum statistics, coulomb interation and strong interaction ot any  //
7 // combination of the three, as applicable.                              //
8 // This class is a wrapper for the fortran code provided by Richard      //
9 // Lednicky.                                                             //
10 //                                                                       //
11 ///////////////////////////////////////////////////////////////////////////
12
13 //#include "StHbtMaker/ThCorrFctn/AliFemtoModelWeightGeneratorLednicky.h"
14 #include "AliFemtoModelWeightGeneratorLednicky.h"
15 #include "AliFemtoModelHiddenInfo.h"
16 #include "AliFemtoPair.h"
17 //#include "StarCallf77.h"
18 //#include <strstream.h>
19 //#include <iomanip.h>
20 //#include <stream>
21 //#include <iomanip>
22 #include <sstream>
23
24 #ifdef SOLARIS
25 # ifndef false
26 typedef int bool;
27 #define false 0
28 #define true 1
29 # endif
30 #endif
31
32 #ifdef WIN32
33 # ifdef CERNLIB_MSSTDCALL
34 #  define F77_UCASE
35 #  define type_of_call _stdcall
36 #  ifndef CERNLIB_QXCAPT
37 #    define CERNLIB_QXCAPT
38 #  endif
39 # else
40 #  define F77_LCASE
41 #  ifndef CERNLIB_QXNO_SC
42 #    define CERNLIB_QXNO_SC
43 #  endif
44 # endif
45 # define type_of_call  _stdcall
46 # define DEFCHARD   const char* , const int        
47 # define DEFCHARL          
48 # define PASSCHARD(string) string, strlen(string) 
49 # define PASSCHARL(string) 
50 #else
51 # define DEFCHARD     const char* 
52 # define DEFCHARL   , const int 
53 # define PASSCHARD(string) string 
54 # define PASSCHARL(string) , strlen(string) 
55 #endif
56 #ifdef CERNLIB_QXCAPT
57 #  define F77_NAME(name,NAME) NAME
58 #else
59 #  if defined(CERNLIB_QXNO_SC)
60 #    define F77_NAME(name,NAME) name
61 #  else
62 #    define F77_NAME(name,NAME) name##_
63 #  endif
64 #endif
65 #ifndef type_of_call
66 # define type_of_call
67 #endif
68
69 // --- Prototype of the function used in the weight calculator 
70 //     (in FsiWeightLedinicky.F)
71 #define fsiin F77_NAME(fsiin,FSIIN)
72 extern "C" {void type_of_call F77_NAME(fsiin,FSIIN)(const int &itest,const int &ich, const int &iqs, const int &isi,const int &i3c);}
73 #define llini F77_NAME(llini,LLINI)
74 extern "C" {void type_of_call F77_NAME(llini,LLINI)(const int &lll,const int &ns, const int &itest);}
75
76 #define fsinucl F77_NAME(fsinucl,FSINUCL)
77 extern "C" {void type_of_call  F77_NAME(fsinucl,FSINUCL)(const double &mn,const double &cn);}
78 #define fsimomentum F77_NAME(fsimomentum,FSIMOMENTUM)
79 extern "C" {void type_of_call F77_NAME(fsimomentum,FSIMOMENTUM)(double &p1,double &p2);}
80 #define fsiposition F77_NAME(fsiposition,FSIPOSITION)
81 extern "C" {void type_of_call F77_NAME(fsiposition,FSIPOSITION)(double &x1,double &x2);}
82 #define fsiw F77_NAME(fsiw,FSIW)
83 extern "C" {void type_of_call F77_NAME(fsiw,FSIW)(const int &i,double &weif,
84                                                   double &wei,double &wein);}
85 #define ltran12 F77_NAME(ltran12,LTRAN12)
86 extern "C" {void type_of_call ltran12_();}
87
88 // Test function for Lambda potential
89 //#define printlam F77_NAME(printlam,PRINTLAM)
90 //extern "C" {void type_of_call printlam_();}
91 //there is not PRINTLAM in *.F file
92
93 // --- Additional prototyping of some CERN functions (in FsiTool.F)
94 typedef float   REAL;
95 typedef struct { REAL re; REAL im; } COMPLEX;
96 #define cgamma F77_NAME(cgamma,CGAMMA)
97 extern "C" {COMPLEX type_of_call cgamma_(COMPLEX*);}
98
99 #ifdef __ROOT__
100 ClassImp(AliFemtoModelWeightGeneratorLednicky)
101 #endif
102
103 AliFemtoModelWeightGeneratorLednicky::AliFemtoModelWeightGeneratorLednicky() : 
104   AliFemtoModelWeightGenerator(),
105   fWei(0), fWein(0), fWeif(0), fWeightDen(0), 
106   fItest(0),fIch(1),fIqs(1),fIsi(1),fI3c(0),
107   fNuclMass(1.),fNuclCharge(0.),
108   fSphereApp(false),fT0App(false) ,
109   fLL(0), fNuclChargeSign(1), fSwap(0), fLLMax(30), fLLName(0), 
110   fNumProcessPair(0), fNumbNonId(0)
111 {
112   // default constructor
113   fLLName=new char*[fLLMax+1];
114   fNumProcessPair=new int[fLLMax+1];
115   int i;
116   for (i=1;i<=fLLMax;i++) {fLLName[i]=new char[40];fNumProcessPair[i]=0;}
117   strncpy( fLLName[1],"neutron neutron",40);
118   strncpy( fLLName[2],"proton proton",40);
119   strncpy( fLLName[3],"neutron proton",40);
120   strncpy( fLLName[4],"alpha alpha",40);
121   strncpy( fLLName[5],"pi+ pi-",40);
122   strncpy( fLLName[6],"pi0 pi0",40);
123   strncpy( fLLName[7],"pi+ pi+",40);
124   strncpy( fLLName[8],"neutron deuteron",40);
125   strncpy( fLLName[9],"proton deuteron",40);
126   strncpy( fLLName[10],"pi+ K-",40);
127   strncpy( fLLName[11],"pi+ K+",40);
128   strncpy( fLLName[12],"pi+ proton",40);
129   strncpy( fLLName[13],"pi- proton",40);
130   strncpy( fLLName[14],"K+ K-",40);
131   strncpy( fLLName[15],"K+ K+",40);
132   strncpy( fLLName[16],"K+ proton",40);
133   strncpy( fLLName[17],"K- proton",40);
134   strncpy( fLLName[18],"deuteron deuteron",40);
135   strncpy( fLLName[19],"deuton alpha",40);
136   strncpy( fLLName[20],"triton triton",40);
137   strncpy( fLLName[21],"triton alpha",40);
138   strncpy( fLLName[22],"K0 K0",40);
139   strncpy( fLLName[23],"K0 K0b",40);
140   strncpy( fLLName[24],"deuteron triton",40);
141   strncpy( fLLName[25],"proton triton",40);
142   strncpy( fLLName[26],"proton alpha",40);
143   strncpy( fLLName[27],"proton lambda",40);
144   strncpy( fLLName[28],"neutron lambda",40);
145   strncpy( fLLName[29],"Lambda lambda",40);// gael 21May02
146   strncpy( fLLName[30],"Proton Anti-proton",40);// gael 21May02
147   FsiInit();
148   FsiNucl();
149 }
150 //______________________
151 AliFemtoModelWeightGeneratorLednicky::AliFemtoModelWeightGeneratorLednicky(const AliFemtoModelWeightGeneratorLednicky &aWeight):
152   AliFemtoModelWeightGenerator(),
153   fWei(0), fWein(0), fWeif(0), fWeightDen(0), 
154   fItest(0),fIch(1),fIqs(1),fIsi(1),fI3c(0),
155   fNuclMass(1.),fNuclCharge(0.),
156   fSphereApp(false),fT0App(false) ,
157   fLL(0), fNuclChargeSign(1), fSwap(0), fLLMax(30), fLLName(0), 
158   fNumProcessPair(0), fNumbNonId(0)
159 {
160   // copy constructor
161   fWei = aWeight.fWei; 
162   fWein = aWeight.  fWein;
163   fWeif = aWeight. fWeif;
164   fWeightDen = aWeight.fWeightDen;
165   
166   fItest = aWeight.fItest;
167   fIch = aWeight.fIch;
168   fIqs = aWeight.fIqs;
169   fIsi = aWeight.fIsi;
170   fI3c = aWeight.fI3c;
171   fNuclMass = aWeight.fNuclMass;
172   fNuclCharge = aWeight.fNuclCharge;
173   fSphereApp = aWeight.fSphereApp;
174   fT0App = aWeight.fT0App; 
175   fLL = aWeight.fLL;
176   fNuclChargeSign = aWeight.fNuclChargeSign;
177   fSwap = aWeight.fSwap;
178   fLLName = aWeight.fLLName; 
179   fNumProcessPair = aWeight.fNumProcessPair;
180   fNumbNonId = aWeight.fNumbNonId;
181   fLLName=new char*[fLLMax+1];
182   fNumProcessPair=new int[fLLMax+1];
183   int i;
184   for (i=1;i<=fLLMax;i++) {fLLName[i]=new char[40];fNumProcessPair[i]=0;}
185   strncpy( fLLName[1],"neutron neutron",40);
186   strncpy( fLLName[2],"proton proton",40);
187   strncpy( fLLName[3],"neutron proton",40);
188   strncpy( fLLName[4],"alpha alpha",40);
189   strncpy( fLLName[5],"pi+ pi-",40);
190   strncpy( fLLName[6],"pi0 pi0",40);
191   strncpy( fLLName[7],"pi+ pi+",40);
192   strncpy( fLLName[8],"neutron deuteron",40);
193   strncpy( fLLName[9],"proton deuteron",40);
194   strncpy( fLLName[10],"pi+ K-",40);
195   strncpy( fLLName[11],"pi+ K+",40);
196   strncpy( fLLName[12],"pi+ proton",40);
197   strncpy( fLLName[13],"pi- proton",40);
198   strncpy( fLLName[14],"K+ K-",40);
199   strncpy( fLLName[15],"K+ K+",40);
200   strncpy( fLLName[16],"K+ proton",40);
201   strncpy( fLLName[17],"K- proton",40);
202   strncpy( fLLName[18],"deuteron deuteron",40);
203   strncpy( fLLName[19],"deuton alpha",40);
204   strncpy( fLLName[20],"triton triton",40);
205   strncpy( fLLName[21],"triton alpha",40);
206   strncpy( fLLName[22],"K0 K0",40);
207   strncpy( fLLName[23],"K0 K0b",40);
208   strncpy( fLLName[24],"deuteron triton",40);
209   strncpy( fLLName[25],"proton triton",40);
210   strncpy( fLLName[26],"proton alpha",40);
211   strncpy( fLLName[27],"proton lambda",40);
212   strncpy( fLLName[28],"neutron lambda",40);
213   strncpy( fLLName[29],"Lambda lambda",40);// gael 21May02
214   strncpy( fLLName[30],"Proton Anti-proton",40);// gael 21May02
215   FsiInit();
216   FsiNucl();
217 }
218
219 AliFemtoModelWeightGeneratorLednicky& AliFemtoModelWeightGeneratorLednicky::operator=(const AliFemtoModelWeightGeneratorLednicky& aWeight)
220 {
221   // assignment operator
222   if (this == &aWeight)
223     return *this;
224
225   fWei = aWeight.fWei; 
226   fWein = aWeight.  fWein;
227   fWeif = aWeight. fWeif;
228   fWeightDen = aWeight.fWeightDen;
229   
230   fItest = aWeight.fItest;
231   fIch = aWeight.fIch;
232   fIqs = aWeight.fIqs;
233   fIsi = aWeight.fIsi;
234   fI3c = aWeight.fI3c;
235   fNuclMass = aWeight.fNuclMass;
236   fNuclCharge = aWeight.fNuclCharge;
237   fSphereApp = aWeight.fSphereApp;
238   fT0App = aWeight.fT0App; 
239   fLL = aWeight.fLL;
240   fNuclChargeSign = aWeight.fNuclChargeSign;
241   fSwap = aWeight.fSwap;
242   //  fLLName = aWeight.fLLName; 
243   fNumProcessPair = aWeight.fNumProcessPair;
244   fNumbNonId = aWeight.fNumbNonId;
245   if (fLLName) free(fLLName);
246   fLLName=new char*[fLLMax+1];
247   if (fNumProcessPair) free(fNumProcessPair);
248   fNumProcessPair=new int[fLLMax+1];
249   int i;
250   for (i=1;i<=fLLMax;i++) {fLLName[i]=new char[40];fNumProcessPair[i]=0;}
251   strncpy( fLLName[1],"neutron neutron",40);
252   strncpy( fLLName[2],"proton proton",40);
253   strncpy( fLLName[3],"neutron proton",40);
254   strncpy( fLLName[4],"alpha alpha",40);
255   strncpy( fLLName[5],"pi+ pi-",40);
256   strncpy( fLLName[6],"pi0 pi0",40);
257   strncpy( fLLName[7],"pi+ pi+",40);
258   strncpy( fLLName[8],"neutron deuteron",40);
259   strncpy( fLLName[9],"proton deuteron",40);
260   strncpy( fLLName[10],"pi+ K-",40);
261   strncpy( fLLName[11],"pi+ K+",40);
262   strncpy( fLLName[12],"pi+ proton",40);
263   strncpy( fLLName[13],"pi- proton",40);
264   strncpy( fLLName[14],"K+ K-",40);
265   strncpy( fLLName[15],"K+ K+",40);
266   strncpy( fLLName[16],"K+ proton",40);
267   strncpy( fLLName[17],"K- proton",40);
268   strncpy( fLLName[18],"deuteron deuteron",40);
269   strncpy( fLLName[19],"deuton alpha",40);
270   strncpy( fLLName[20],"triton triton",40);
271   strncpy( fLLName[21],"triton alpha",40);
272   strncpy( fLLName[22],"K0 K0",40);
273   strncpy( fLLName[23],"K0 K0b",40);
274   strncpy( fLLName[24],"deuteron triton",40);
275   strncpy( fLLName[25],"proton triton",40);
276   strncpy( fLLName[26],"proton alpha",40);
277   strncpy( fLLName[27],"proton lambda",40);
278   strncpy( fLLName[28],"neutron lambda",40);
279   strncpy( fLLName[29],"Lambda lambda",40);// gael 21May02
280   strncpy( fLLName[30],"Proton Anti-proton",40);// gael 21May02
281   FsiInit();
282   FsiNucl();
283   
284   return *this;
285 }
286
287
288 double AliFemtoModelWeightGeneratorLednicky::GenerateWeight(AliFemtoPair* aPair)
289 {
290   // Get hidden information pointers
291   AliFemtoModelHiddenInfo *inf1 = (AliFemtoModelHiddenInfo *) aPair->Track1()->HiddenInfo();
292   AliFemtoModelHiddenInfo *inf2 = (AliFemtoModelHiddenInfo *) aPair->Track2()->HiddenInfo();
293
294   // Calculate pair variables
295   Double_t tPx = inf1->GetTrueMomentum()->x()+inf2->GetTrueMomentum()->x();
296   Double_t tPy = inf1->GetTrueMomentum()->y()+inf2->GetTrueMomentum()->y();
297   Double_t tPz = inf1->GetTrueMomentum()->z()+inf2->GetTrueMomentum()->z();
298   Double_t tM1 = inf1->GetMass();
299   Double_t tM2 = inf2->GetMass();
300   Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->Mag2());
301   Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->Mag2());
302   Double_t tE  = tE1 + tE2;
303   Double_t tPt = tPx*tPx + tPy*tPy;
304   Double_t tMt = tE*tE - tPz*tPz;//mCVK;
305   Double_t tM  = sqrt(tMt - tPt);
306   tMt = sqrt(tMt);
307   tPt = sqrt(tPt);
308   Double_t tBetat = tPt/tMt;
309
310   // Boost to LCMS
311   Double_t tBeta = tPz/tE;
312   Double_t tGamma = tE/tMt;         
313   fKStarLong = tGamma * (inf1->GetTrueMomentum()->z() - tBeta * tE1);
314   Double_t tE1L = tGamma * (tE1  - tBeta * inf1->GetTrueMomentum()->z());
315     
316   // Rotate in transverse plane
317   fKStarOut  = ( inf1->GetTrueMomentum()->x()*tPx + inf1->GetTrueMomentum()->y()*tPy)/tPt;
318   fKStarSide = (-inf1->GetTrueMomentum()->x()*tPy + inf1->GetTrueMomentum()->y()*tPx)/tPt;
319       
320   // Boost to pair cms
321   fKStarOut = tMt/tM * (fKStarOut - tPt/tMt * tE1L);
322   
323   tBetat = tPt/tMt;
324   
325   Double_t tDX = inf1->GetEmissionPoint()->x()-inf2->GetEmissionPoint()->x();
326   Double_t tDY = inf1->GetEmissionPoint()->y()-inf2->GetEmissionPoint()->y();
327   Double_t tRLong = inf1->GetEmissionPoint()->z()-inf2->GetEmissionPoint()->z();
328   Double_t tDTime = inf1->GetEmissionPoint()->t()-inf2->GetEmissionPoint()->t();
329
330   Double_t tROut = (tDX*tPx + tDY*tPy)/tPt;
331   Double_t tRSide = (-tDX*tPy + tDY*tPx)/tPt;
332
333 //   cout << "Got points 1 " << inf1->GetEmissionPoint()->x() << "  " <<  inf1->GetEmissionPoint()->y() << " "  << inf1->GetEmissionPoint()->z() << "  " << inf1->GetEmissionPoint()->t() << endl;
334
335 //   cout << "Got points 2 " << inf2->GetEmissionPoint()->x() << "  " << inf2->GetEmissionPoint()->y() << " " << inf2->GetEmissionPoint()->z() << "  " << inf2->GetEmissionPoint()->t() << endl;
336
337   fRStarSide = tRSide;
338
339   fRStarLong = tGamma*(tRLong - tBeta* tDTime);
340   Double_t tDTimePairLCMS = tGamma*(tDTime - tBeta* tRLong);
341
342   tBeta = tPt/tMt;
343   tGamma = tMt/tM;
344
345   fRStarOut = tGamma*(tROut - tBeta* tDTimePairLCMS);
346   fRStar = ::sqrt(fRStarOut*fRStarOut + fRStarSide*fRStarSide +
347                            fRStarLong*fRStarLong);
348   fKStar = ::sqrt(fKStarOut*fKStarOut + fKStarSide*fKStarSide + fKStarLong*fKStarLong);
349
350 //   cout << "Got out side " << fRStarOut << " " << fRStarSide << endl;
351
352   if (!SetPid(inf1->GetPDGPid(),inf2->GetPDGPid())) {
353     fWeightDen=1.;
354     return 1;    
355   } 
356   else { // Good Pid
357     AliFemtoThreeVector*  p;
358     p=(inf1->GetTrueMomentum());
359     double p1[]={p->x(),p->y(),p->z()};
360     p=(inf2->GetTrueMomentum());
361     double p2[]={p->x(),p->y(),p->z()};
362     if ((p1[0]==p2[0])&&(p1[1]==p2[1])&&(p1[2]==p2[2])) {
363       fWeightDen=0.;
364       return 0;  
365     } 
366     if (fSwap) {
367       fsimomentum(*p2,*p1);
368     } else {
369       fsimomentum(*p1,*p2);
370     }
371     AliFemtoLorentzVector* tPoint;
372     tPoint=(inf1->GetEmissionPoint());
373 //     cout << "Pid1:dans GetWeight = " << aThPair->GetPid1() << endl;
374 //     cout << "Pid2:dans GetWeight = " << aThPair->GetPid2() << endl;
375 //     cout << "LL:in GetWeight = " << mLL << endl;
376
377     double x1[]={tPoint->x(),tPoint->y(),tPoint->z(),tPoint->t()};
378     tPoint=(inf2->GetEmissionPoint());
379     double x2[]={tPoint->x(),tPoint->y(),tPoint->z(),tPoint->t()};
380     if ((x1[0]==x2[0])&&(x1[1]==x2[1])&&(x1[2]==x2[2])&&(x1[3]==x2[3])) {
381       fWeightDen=0.;
382       return 0;  
383     } 
384     if (fSwap) {
385       fsiposition(*x2,*x1);
386     } else {
387       fsiposition(*x1,*x2);
388     }
389     FsiSetLL();
390     ltran12();
391     fsiw(1,fWeif,fWei,fWein);
392
393     if (fI3c==0) return fWein;
394     fWeightDen=fWeif;
395     return fWei;
396   }
397 }
398
399 AliFemtoString AliFemtoModelWeightGeneratorLednicky::Report() {
400   // create report
401   ostringstream tStr; 
402   tStr << "Lednicky afterburner calculation for  Correlation -  Report" << endl;
403   tStr << "    Setting : Quantum : " << ((fIqs) ? "On" : "Off"); 
404   tStr << " - Coulbomb : " << ((fIch) ? "On" : "Off") ;
405   tStr << " - Strong : " << ((fIsi) ? "On" : "Off");
406   tStr << endl;
407   tStr << "              3-Body : " << ((fI3c) ? "On"  : "Off") ;
408   if (fI3c) tStr << " Mass=" <<  fNuclMass << " - Charge= " << fNuclCharge ;
409   tStr << endl;
410   tStr << "    " << fNumProcessPair[0] << " Pairs have been Processed :" << endl;
411   int i;
412   for(i=1;i<=fLLMax;i++) { 
413     if (fNumProcessPair[i])
414       tStr << "         " << fNumProcessPair[i] << " " << fLLName[i] << endl;
415   }
416   if (fNumbNonId)
417     tStr << "         "<< fNumbNonId << " Non Identified" << endl;
418   AliFemtoString returnThis = tStr.str();
419   return returnThis;
420 }
421
422 void AliFemtoModelWeightGeneratorLednicky::FsiInit(){
423   // Initialize weight generation module
424 //   cout << "*******************AliFemtoModelWeightGeneratorLednicky check FsiInit ************" << endl;
425 //   cout <<"mItest dans FsiInit() = " << fItest << endl;
426 //   cout <<"mIch dans FsiInit() = " << fIch << endl;
427 //   cout <<"mIqs dans FsiInit() = " << fIqs << endl;
428 //   cout <<"mIsi dans FsiInit() = " << fIsi << endl;
429 //   cout <<"mI3c dans FsiInit() = " << fI3c << endl;
430   fsiin(fItest,fIch,fIqs,fIsi,fI3c);
431 }
432
433 void AliFemtoModelWeightGeneratorLednicky::FsiNucl(){
434   // initialize weight generation taking into account the residual charge
435 //   cout << "*******************AliFemtoModelWeightGeneratorLednicky check FsiNucl ************" << endl;
436 //   cout <<"fNuclMass dans FsiNucl() = " << fNuclMass << endl;
437 //   cout <<"fNuclCharge dans FsiNucl() = " << fNuclCharge << endl;
438 //   cout <<"fNuclChargeSign dans FsiNucl() = " << fNuclChargeSign << endl;
439   fsinucl(fNuclMass,fNuclCharge*fNuclChargeSign);
440 }
441
442 void AliFemtoModelWeightGeneratorLednicky::FsiSetLL(){
443   // set internal pair type for the module
444   int tNS;
445   if (fSphereApp||(fLL>5)) {
446     if (fT0App) { tNS=4;} 
447     else {tNS=2;}
448   } else { tNS=1;}
449    //cout <<"fLL dans FsiSetLL() = "<< fLL << endl;
450    //cout <<"tNS dans FsiSetLL() = "<< tNS << endl;
451    //cout <<"fItest dans FsiSetLL() = "<< fItest << endl;
452   llini(fLL,tNS,fItest);
453   //cout<<" end of FsiSetLL"<<endl;
454 }
455          
456 bool AliFemtoModelWeightGeneratorLednicky::SetPid(const int aPid1,const int aPid2) {
457   // set calculated system for basing on particles' pids
458   static const int ksPi0Pid=111;
459   static const int ksPionPid=211; 
460   static const int ksK0Pid=311;
461   static const int ksKPid=321;
462   static const int ksNeutPid=2112;
463   static const int ksProtPid=2212;
464   static const int ksLamPid=3122;
465   //  static const int sLamLamPid=3122;
466
467    // cout << "Setting PID to " << aPid1 << " " << aPid2 << endl;
468
469   int tPidl,tPidh;
470   int tChargeFactor=1;
471   
472   if (abs(aPid1)<abs(aPid2)) {
473     if (aPid1<0) tChargeFactor=-1;
474     tPidl=aPid1*tChargeFactor;
475     tPidh=aPid2*tChargeFactor;
476     fSwap=false;
477   } else {
478     if (aPid2<0) tChargeFactor=-1;
479     tPidl=aPid2*tChargeFactor;
480     tPidh=aPid1*tChargeFactor;
481     fSwap=true;
482   }
483   switch (tPidl) {
484   case ksPionPid:
485     switch (tPidh) {
486     case -ksPionPid:   fLL=5; tChargeFactor*=1 ;break;
487     case ksPionPid:    fLL=7; tChargeFactor*=1 ;break;
488     case -ksKPid:      fLL=10;tChargeFactor*=1 ;break;  
489     case ksKPid:       fLL=11;tChargeFactor*=1 ;break;  
490     case ksProtPid:    fLL=12;tChargeFactor*=1 ;break;
491     case -ksProtPid:   fLL=13;tChargeFactor*=-1;break;
492     default: fLL=0;
493     }
494     break;
495   case ksProtPid:
496     switch (tPidh) {
497     case ksProtPid:    fLL=2; tChargeFactor*=1 ;break;
498     case ksLamPid:     fLL=27;tChargeFactor*=1 ;break;
499     case -ksProtPid:   fLL=30;tChargeFactor*=1 ;break;
500     default: fLL=0;
501     }
502     break;
503   case ksKPid:
504     switch (tPidh) {
505     case -ksKPid:      fLL=14;tChargeFactor*=1 ;break;
506     case ksKPid:       fLL=15;tChargeFactor*=1 ;break;
507     case ksProtPid:    fLL=16;tChargeFactor*=1 ;break;
508     case -ksProtPid:   fLL=17;tChargeFactor*=-1 ;break;
509     default: fLL=0;
510     }
511     break;    
512   case ksK0Pid:
513     switch (tPidh) {
514     case ksK0Pid:         fLL=22;tChargeFactor*=1 ;break;
515     case -ksK0Pid:        fLL=23;tChargeFactor*=1 ;break;
516     default: fLL=0;
517     }
518     break;   
519   case ksPi0Pid:
520     switch (tPidh) {
521     case ksPi0Pid:        fLL=6; tChargeFactor*=1 ;break;
522     default: fLL=0;
523     }
524     break;
525   case ksNeutPid:
526     switch (tPidh) {
527     case ksNeutPid:      fLL=1; tChargeFactor*=1 ;break;
528     case ksProtPid:      fLL=3; tChargeFactor*=1 ;break;
529     case ksLamPid:       fLL=28;tChargeFactor*=1 ;break;
530     default: fLL=0;
531     }
532     break;                                             //Gael 21May02 
533   case ksLamPid:                                        //Gael 21May02 
534     switch (tPidh) {                                   //Gael 21May02 
535     case ksLamPid:       fLL=29;tChargeFactor*=1 ;break;//Gael 21May02  
536     default: fLL=0;                                    //Gael 21May02 
537     }                                                 //Gael 21May02 
538     break;                                             //Gael 21May02 
539   default: fLL=0;
540   }
541   if (tChargeFactor!=fNuclChargeSign) {
542     fNuclChargeSign=tChargeFactor;
543     FsiNucl();
544   }
545   (fNumProcessPair[0])++;
546   if (fLL) {
547     (fNumProcessPair[fLL])++;
548     return true;
549   } else {
550     fNumbNonId++;
551     return false;
552   }
553 //   cout << "*******************AliFemtoModelWeightGeneratorLednicky check SetPid ************" << endl;
554 //   cout << "fLL=="<< fLL << endl;
555 //   cout << "fNuclCharge=="<< fNuclCharge << endl;
556
557 }    
558 AliFemtoModelWeightGeneratorLednicky::~AliFemtoModelWeightGeneratorLednicky() 
559
560   if (fLLName) delete [] fLLName;
561   if (fNumProcessPair) delete [] fNumProcessPair;
562 /* no-op */ 
563 }
564
565 //_____________________________________________
566 void     AliFemtoModelWeightGeneratorLednicky::SetPairType(Int_t aPairType)
567 {
568   // set calculated system basing on the pair type
569   fPairType = aPairType;
570   if (fPairType == fgkPionPlusPionPlus) SetPid(211,211);
571   if (fPairType == fgkPionPlusPionMinus ) SetPid(211, -211);
572   if (fPairType == fgkKaonPlusKaonPlus ) SetPid(321, 321);
573   if (fPairType == fgkKaonPlusKaonMinus ) SetPid(321, -321);
574   if (fPairType == fgkProtonProton ) SetPid(2212, 2212);
575   if (fPairType == fgkProtonAntiproton ) SetPid(2212, -2212);
576   if (fPairType == fgkPionPlusKaonPlus ) SetPid(211, 321);
577   if (fPairType == fgkPionPlusKaonMinus ) SetPid(211, -321);
578   if (fPairType == fgkPionPlusProton ) SetPid(211, 2212);
579   if (fPairType == fgkPionPlusAntiproton ) SetPid(211, -2212);
580   if (fPairType == fgkKaonPlusProton ) SetPid(321, 2212);
581   if (fPairType == fgkKaonPlusAntiproton ) SetPid(321, -2212);
582 }
583
584 //_____________________________________________
585 Int_t    AliFemtoModelWeightGeneratorLednicky::GetPairType() const
586 {
587   // return pair type
588   return fPairType;
589 }
590
591 //_____________________________________________
592 void     AliFemtoModelWeightGeneratorLednicky::SetPairTypeFromPair(AliFemtoPair *aPair)
593 {
594   // set calculated system based on the hidden info in the pair
595   AliFemtoModelHiddenInfo *inf1 = ( AliFemtoModelHiddenInfo *) aPair->Track1()->HiddenInfo();
596   AliFemtoModelHiddenInfo *inf2 = ( AliFemtoModelHiddenInfo *) aPair->Track2()->HiddenInfo();
597
598   const Int_t ktPid1 = inf1->GetPDGPid();
599   const Int_t ktPid2 = inf2->GetPDGPid();
600
601   if      (((ktPid1 ==   211) && (ktPid2 ==   211)) ||
602            ((ktPid1 ==  -211) && (ktPid2 ==  -211)))
603     fPairType = fgkPionPlusPionPlus;
604   else if (((ktPid1 ==  -211) && (ktPid2 ==   211)) ||
605            ((ktPid1 ==   211) && (ktPid2 ==  -211)))
606     fPairType = fgkPionPlusPionMinus;
607   else if (((ktPid1 ==   321) && (ktPid2 ==   321)) ||
608            ((ktPid1 ==  -321) && (ktPid2 ==  -321)))
609     fPairType = fgkKaonPlusKaonPlus;
610   else if (((ktPid1 ==  -321) && (ktPid2 ==   321)) ||
611            ((ktPid1 ==   321) && (ktPid2 ==  -321)))
612     fPairType = fgkKaonPlusKaonMinus;
613   else if (((ktPid1 ==  2212) && (ktPid2 ==  2212)) ||
614            ((ktPid1 == -2212) && (ktPid2 == -2212)))
615     fPairType = fgkProtonProton;
616   else if (((ktPid1 == -2212) && (ktPid2 ==  2212)) ||
617            ((ktPid1 ==  2212) && (ktPid2 == -2212)))
618     fPairType = fgkProtonAntiproton;
619   else if (((ktPid1 ==   211) && (ktPid2 ==   321)) ||
620            ((ktPid1 ==  -211) && (ktPid2 ==  -321)))
621     fPairType = fgkPionPlusKaonPlus;
622   else if (((ktPid1 ==  -211) && (ktPid2 ==   321)) ||
623            ((ktPid1 ==   211) && (ktPid2 ==  -321)))
624     fPairType = fgkPionPlusKaonMinus;
625   else if (((ktPid1 ==   211) && (ktPid2 ==  2212)) ||
626            ((ktPid1 ==  -211) && (ktPid2 == -2212)))
627     fPairType = fgkPionPlusProton;
628   else if (((ktPid1 ==  -211) && (ktPid2 ==  2212)) ||
629            ((ktPid1 ==   211) && (ktPid2 == -2212)))
630     fPairType = fgkPionPlusAntiproton;
631   else if (((ktPid1 ==   321) && (ktPid2 ==  2212)) ||
632            ((ktPid1 ==  -321) && (ktPid2 == -2212)))
633     fPairType = fgkKaonPlusProton;
634   else if (((ktPid1 ==  -321) && (ktPid2 ==  2212)) ||
635            ((ktPid1 ==   321) && (ktPid2 == -2212)))
636     fPairType = fgkKaonPlusAntiproton;
637   SetPid(ktPid1, ktPid2);
638 }
639
640 void AliFemtoModelWeightGeneratorLednicky::SetNuclCharge(const double aNuclCharge) {fNuclCharge=aNuclCharge;FsiNucl();}
641 void AliFemtoModelWeightGeneratorLednicky::SetNuclMass(const double aNuclMass){fNuclMass=aNuclMass;FsiNucl();}
642
643 void AliFemtoModelWeightGeneratorLednicky::SetSphere(){fSphereApp=true;}
644 void AliFemtoModelWeightGeneratorLednicky::SetSquare(){fSphereApp=false;}
645 void AliFemtoModelWeightGeneratorLednicky::SetT0ApproxOn(){ fT0App=true;}
646 void AliFemtoModelWeightGeneratorLednicky::SetT0ApproxOff(){ fT0App=false;}
647 void AliFemtoModelWeightGeneratorLednicky::SetDefaultCalcPar(){
648   fItest=1;fIqs=1;fIsi=1;fI3c=0;fIch=1;FsiInit();
649   fSphereApp=false;fT0App=false;}
650
651 void AliFemtoModelWeightGeneratorLednicky::SetCoulOn()    {fItest=1;fIch=1;FsiInit();}
652 void AliFemtoModelWeightGeneratorLednicky::SetCoulOff()   {fItest=1;fIch=0;FsiInit();}
653 void AliFemtoModelWeightGeneratorLednicky::SetQuantumOn() {fItest=1;fIqs=1;FsiInit();}
654 void AliFemtoModelWeightGeneratorLednicky::SetQuantumOff(){fItest=1;fIqs=0;FsiInit();}
655 void AliFemtoModelWeightGeneratorLednicky::SetStrongOn()  {fItest=1;fIsi=1;FsiInit();}
656 void AliFemtoModelWeightGeneratorLednicky::SetStrongOff() {fItest=1;fIsi=0;FsiInit();}
657 void AliFemtoModelWeightGeneratorLednicky::Set3BodyOn()   {fItest=1;fI3c=1;FsiInit();FsiNucl();}
658 void AliFemtoModelWeightGeneratorLednicky::Set3BodyOff()  {fItest=1;fI3c=0;FsiInit();fWeightDen=1.;FsiNucl();}
659
660 Double_t AliFemtoModelWeightGeneratorLednicky::GetKStar() const {return AliFemtoModelWeightGenerator::GetKStar();}
661 Double_t AliFemtoModelWeightGeneratorLednicky::GetKStarOut() const { return AliFemtoModelWeightGenerator::GetKStarOut(); }
662 Double_t AliFemtoModelWeightGeneratorLednicky::GetKStarSide() const { return AliFemtoModelWeightGenerator::GetKStarSide(); }
663 Double_t AliFemtoModelWeightGeneratorLednicky::GetKStarLong() const { return AliFemtoModelWeightGenerator::GetKStarLong(); }
664 Double_t AliFemtoModelWeightGeneratorLednicky::GetRStar() const { return AliFemtoModelWeightGenerator::GetRStar(); }
665 Double_t AliFemtoModelWeightGeneratorLednicky::GetRStarOut() const { return AliFemtoModelWeightGenerator::GetRStarOut(); }
666 Double_t AliFemtoModelWeightGeneratorLednicky::GetRStarSide() const { return AliFemtoModelWeightGenerator::GetRStarSide(); }
667 Double_t AliFemtoModelWeightGeneratorLednicky::GetRStarLong() const { return AliFemtoModelWeightGenerator::GetRStarLong(); }
668
669 AliFemtoModelWeightGenerator* AliFemtoModelWeightGeneratorLednicky::Clone() const {
670   AliFemtoModelWeightGenerator* tmp = new AliFemtoModelWeightGeneratorLednicky(*this);
671   return tmp;
672 }