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