PWGCF/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   AliFemtoTrack *inf1 = (AliFemtoTrack *) aPair->Track1()->Track();
 294   AliFemtoTrack *inf2 = (AliFemtoTrack *) aPair->Track2()->Track();
 295
 296   // Calculate pair variables
 297   Double_t tPx = inf1->GetTrueMomentum()->x()+inf2->GetTrueMomentum()->x();
 298   Double_t tPy = inf1->GetTrueMomentum()->y()+inf2->GetTrueMomentum()->y();
 299   Double_t tPz = inf1->GetTrueMomentum()->z()+inf2->GetTrueMomentum()->z();
 300   Double_t tM1 = inf1->GetMass();
 301   Double_t tM2 = inf2->GetMass();
 302   Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->Mag2());
 303   Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->Mag2());
 304   Double_t tE  = tE1 + tE2;
 305   Double_t tPt = tPx*tPx + tPy*tPy;
 306   Double_t tMt = tE*tE - tPz*tPz;//mCVK;
 307   Double_t tM  = sqrt(tMt - tPt);
 308   tMt = sqrt(tMt);
 309   tPt = sqrt(tPt);
 310   Double_t tBetat = tPt/tMt;
 311
 312   // Boost to LCMS
 313   Double_t tBeta = tPz/tE;
 314   Double_t tGamma = tE/tMt;
 315   fKStarLong = tGamma * (inf1->GetTrueMomentum()->z() - tBeta * tE1);
 316   Double_t tE1L = tGamma * (tE1  - tBeta * inf1->GetTrueMomentum()->z());
 317
 318   // Rotate in transverse plane
 319   fKStarOut  = ( inf1->GetTrueMomentum()->x()*tPx + inf1->GetTrueMomentum()->y()*tPy)/tPt;
 320   fKStarSide = (-inf1->GetTrueMomentum()->x()*tPy + inf1->GetTrueMomentum()->y()*tPx)/tPt;
 321
 322   // Boost to pair cms
 323   fKStarOut = tMt/tM * (fKStarOut - tPt/tMt * tE1L);
 324
 325   tBetat = tPt/tMt;
 326
 327   Double_t tDX = inf1->GetEmissionPoint()->x()-inf2->GetEmissionPoint()->x();
 328   Double_t tDY = inf1->GetEmissionPoint()->y()-inf2->GetEmissionPoint()->y();
 329   Double_t tRLong = inf1->GetEmissionPoint()->z()-inf2->GetEmissionPoint()->z();
 330   Double_t tDTime = inf1->GetEmissionPoint()->t()-inf2->GetEmissionPoint()->t();
 331
 332   Double_t tROut = (tDX*tPx + tDY*tPy)/tPt;
 333   Double_t tRSide = (-tDX*tPy + tDY*tPx)/tPt;
 334
 335 //   cout << "Got points 1 " << inf1->GetEmissionPoint()->x() << "  " <<  inf1->GetEmissionPoint()->y() << " "  << inf1->GetEmissionPoint()->z() << "  " << inf1->GetEmissionPoint()->t() << endl;
 336
 337 //   cout << "Got points 2 " << inf2->GetEmissionPoint()->x() << "  " << inf2->GetEmissionPoint()->y() << " " << inf2->GetEmissionPoint()->z() << "  " << inf2->GetEmissionPoint()->t() << endl;
 338
 339   fRStarSide = tRSide;
 340
 341   fRStarLong = tGamma*(tRLong - tBeta* tDTime);
 342   Double_t tDTimePairLCMS = tGamma*(tDTime - tBeta* tRLong);
 343
 344   tBeta = tPt/tMt;
 345   tGamma = tMt/tM;
 346
 347   fRStarOut = tGamma*(tROut - tBeta* tDTimePairLCMS);
 348   fRStar = ::sqrt(fRStarOut*fRStarOut + fRStarSide*fRStarSide +
 349                            fRStarLong*fRStarLong);
 350   fKStar = ::sqrt(fKStarOut*fKStarOut + fKStarSide*fKStarSide + fKStarLong*fKStarLong);
 351
 352 //   cout << "Got out side " << fRStarOut << " " << fRStarSide << endl;
 353
 354   if (!SetPid(inf1->GetPDGPid(),inf2->GetPDGPid())) {
 355     fWeightDen=1.;
 356     return 1;
 357   }
 358   else { // Good Pid
 359     AliFemtoThreeVector*  p;
 360     p=(inf1->GetTrueMomentum());
 361     double p1[]={p->x(),p->y(),p->z()};
 362     p=(inf2->GetTrueMomentum());
 363     double p2[]={p->x(),p->y(),p->z()};
 364     if ((p1[0]==p2[0])&&(p1[1]==p2[1])&&(p1[2]==p2[2])) {
 365       fWeightDen=0.;
 366       return 0;
 367     }
 368     if (fSwap) {
 369       fsimomentum(*p2,*p1);
 370     } else {
 371       fsimomentum(*p1,*p2);
 372     }
 373     AliFemtoLorentzVector* tPoint;
 374     tPoint=(inf1->GetEmissionPoint());
 375 //     cout << "Pid1:dans GetWeight = " << aThPair->GetPid1() << endl;
 376 //     cout << "Pid2:dans GetWeight = " << aThPair->GetPid2() << endl;
 377 //     cout << "LL:in GetWeight = " << mLL << endl;
 378
 379     double x1[]={tPoint->x(),tPoint->y(),tPoint->z(),tPoint->t()};
 380     tPoint=(inf2->GetEmissionPoint());
 381     double x2[]={tPoint->x(),tPoint->y(),tPoint->z(),tPoint->t()};
 382     if ((x1[0]==x2[0])&&(x1[1]==x2[1])&&(x1[2]==x2[2])&&(x1[3]==x2[3])) {
 383       fWeightDen=0.;
 384       return 0;
 385     }
 386     if (fSwap) {
 387       fsiposition(*x2,*x1);
 388     } else {
 389       fsiposition(*x1,*x2);
 390     }
 391     FsiSetLL();
 392     ltran12();
 393     fsiw(1,fWeif,fWei,fWein);
 394
 395     if (fI3c==0) return fWein;
 396     fWeightDen=fWeif;
 397     return fWei;
 398   }
 399 }
 400
 401 AliFemtoString AliFemtoModelWeightGeneratorLednicky::Report() {
 402   // create report
 403   ostringstream tStr;
 404   tStr << "Lednicky afterburner calculation for  Correlation -  Report" << endl;
 405   tStr << "    Setting : Quantum : " << ((fIqs) ? "On" : "Off");
 406   tStr << " - Coulbomb : " << ((fIch) ? "On" : "Off") ;
 407   tStr << " - Strong : " << ((fIsi) ? "On" : "Off");
 408   tStr << endl;
 409   tStr << "              3-Body : " << ((fI3c) ? "On"  : "Off") ;
 410   if (fI3c) tStr << " Mass=" <<  fNuclMass << " - Charge= " << fNuclCharge ;
 411   tStr << endl;
 412   tStr << "    " << fNumProcessPair[0] << " Pairs have been Processed :" << endl;
 413   int i;
 414   for(i=1;i<=fLLMax;i++) {
 415     if (fNumProcessPair[i])
 416       tStr << "         " << fNumProcessPair[i] << " " << fLLName[i] << endl;
 417   }
 418   if (fNumbNonId)
 419     tStr << "         "<< fNumbNonId << " Non Identified" << endl;
 420   AliFemtoString returnThis = tStr.str();
 421   return returnThis;
 422 }
 423
 424 void AliFemtoModelWeightGeneratorLednicky::FsiInit(){
 425   // Initialize weight generation module
 426 //   cout << "*******************AliFemtoModelWeightGeneratorLednicky check FsiInit ************" << endl;
 427 //   cout <<"mItest dans FsiInit() = " << fItest << endl;
 428 //   cout <<"mIch dans FsiInit() = " << fIch << endl;
 429 //   cout <<"mIqs dans FsiInit() = " << fIqs << endl;
 430 //   cout <<"mIsi dans FsiInit() = " << fIsi << endl;
 431 //   cout <<"mI3c dans FsiInit() = " << fI3c << endl;
 432   fsiin(fItest,fIch,fIqs,fIsi,fI3c);
 433 }
 434
 435 void AliFemtoModelWeightGeneratorLednicky::FsiNucl(){
 436   // initialize weight generation taking into account the residual charge
 437 //   cout << "*******************AliFemtoModelWeightGeneratorLednicky check FsiNucl ************" << endl;
 438 //   cout <<"fNuclMass dans FsiNucl() = " << fNuclMass << endl;
 439 //   cout <<"fNuclCharge dans FsiNucl() = " << fNuclCharge << endl;
 440 //   cout <<"fNuclChargeSign dans FsiNucl() = " << fNuclChargeSign << endl;
 441   fsinucl(fNuclMass,fNuclCharge*fNuclChargeSign);
 442 }
 443
 444 void AliFemtoModelWeightGeneratorLednicky::FsiSetLL(){
 445   // set internal pair type for the module
 446   int tNS;
 447   if (fSphereApp||(fLL>5)) {
 448     if (fT0App) { tNS=4;}
 449     else {tNS=2;}
 450   } else { tNS=1;}
 451    //cout <<"fLL dans FsiSetLL() = "<< fLL << endl;
 452    //cout <<"tNS dans FsiSetLL() = "<< tNS << endl;
 453    //cout <<"fItest dans FsiSetLL() = "<< fItest << endl;
 454   llini(fLL,tNS,fItest);
 455   //cout<<" end of FsiSetLL"<<endl;
 456 }
 457
 458 bool AliFemtoModelWeightGeneratorLednicky::SetPid(const int aPid1,const int aPid2) {
 459   // set calculated system for basing on particles' pids
 460   static const int ksPi0Pid=111;
 461   static const int ksPionPid=211;
 462   static const int ksK0Pid=311;
 463   static const int ksKPid=321;
 464   static const int ksNeutPid=2112;
 465   static const int ksProtPid=2212;
 466   static const int ksLamPid=3122;
 467   //  static const int sLamLamPid=3122;
 468
 469    // cout << "Setting PID to " << aPid1 << " " << aPid2 << endl;
 470
 471   int tPidl,tPidh;
 472   int tChargeFactor=1;
 473
 474   if (abs(aPid1)<abs(aPid2)) {
 475     if (aPid1<0) tChargeFactor=-1;
 476     tPidl=aPid1*tChargeFactor;
 477     tPidh=aPid2*tChargeFactor;
 478     fSwap=false;
 479   } else {
 480     if (aPid2<0) tChargeFactor=-1;
 481     tPidl=aPid2*tChargeFactor;
 482     tPidh=aPid1*tChargeFactor;
 483     fSwap=true;
 484   }
 485   switch (tPidl) {
 486   case ksPionPid:
 487     switch (tPidh) {
 488     case -ksPionPid:   fLL=5; tChargeFactor*=1 ;break;
 489     case ksPionPid:    fLL=7; tChargeFactor*=1 ;break;
 490     case -ksKPid:      fLL=10;tChargeFactor*=1 ;break;
 491     case ksKPid:       fLL=11;tChargeFactor*=1 ;break;
 492     case ksProtPid:    fLL=12;tChargeFactor*=1 ;break;
 493     case -ksProtPid:   fLL=13;tChargeFactor*=-1;break;
 494     default: fLL=0;
 495     }
 496     break;
 497   case ksProtPid:
 498     switch (tPidh) {
 499     case ksProtPid:    fLL=2; tChargeFactor*=1 ;break;
 500     case ksLamPid:     fLL=27;tChargeFactor*=1 ;break;
 501     case -ksProtPid:   fLL=30;tChargeFactor*=1 ;break;
 502     default: fLL=0;
 503     }
 504     break;
 505   case ksKPid:
 506     switch (tPidh) {
 507     case -ksKPid:      fLL=14;tChargeFactor*=1 ;break;
 508     case ksKPid:       fLL=15;tChargeFactor*=1 ;break;
 509     case ksProtPid:    fLL=16;tChargeFactor*=1 ;break;
 510     case -ksProtPid:   fLL=17;tChargeFactor*=-1 ;break;
 511     default: fLL=0;
 512     }
 513     break;
 514   case ksK0Pid:
 515     switch (tPidh) {
 516     case ksK0Pid:         fLL=22;tChargeFactor*=1 ;break;
 517     case -ksK0Pid:        fLL=23;tChargeFactor*=1 ;break;
 518     default: fLL=0;
 519     }
 520     break;
 521   case ksPi0Pid:
 522     switch (tPidh) {
 523     case ksPi0Pid:        fLL=6; tChargeFactor*=1 ;break;
 524     default: fLL=0;
 525     }
 526     break;
 527   case ksNeutPid:
 528     switch (tPidh) {
 529     case ksNeutPid:      fLL=1; tChargeFactor*=1 ;break;
 530     case ksProtPid:      fLL=3; tChargeFactor*=1 ;break;
 531     case ksLamPid:       fLL=28;tChargeFactor*=1 ;break;
 532     default: fLL=0;
 533     }
 534     break;                                             //Gael 21May02
 535   case ksLamPid:                                        //Gael 21May02
 536     switch (tPidh) {                                   //Gael 21May02
 537     case ksLamPid:       fLL=29;tChargeFactor*=1 ;break;//Gael 21May02
 538     default: fLL=0;                                    //Gael 21May02
 539     }                                                 //Gael 21May02
 540     break;                                             //Gael 21May02
 541   default: fLL=0;
 542   }
 543   if (tChargeFactor!=fNuclChargeSign) {
 544     fNuclChargeSign=tChargeFactor;
 545     FsiNucl();
 546   }
 547   (fNumProcessPair[0])++;
 548   if (fLL) {
 549     (fNumProcessPair[fLL])++;
 550     return true;
 551   } else {
 552     fNumbNonId++;
 553     return false;
 554   }
 555 //   cout << "*******************AliFemtoModelWeightGeneratorLednicky check SetPid ************" << endl;
 556 //   cout << "fLL=="<< fLL << endl;
 557 //   cout << "fNuclCharge=="<< fNuclCharge << endl;
 558
 559 }
 560 AliFemtoModelWeightGeneratorLednicky::~AliFemtoModelWeightGeneratorLednicky()
 561 {
 562   if (fLLName) delete [] fLLName;
 563   if (fNumProcessPair) delete [] fNumProcessPair;
 564 /* no-op */
 565 }
 566
 567 //_____________________________________________
 568 void     AliFemtoModelWeightGeneratorLednicky::SetPairType(Int_t aPairType)
 569 {
 570   // set calculated system basing on the pair type
 571   fPairType = aPairType;
 572   if (fPairType == fgkPionPlusPionPlus) SetPid(211,211);
 573   if (fPairType == fgkPionPlusPionMinus ) SetPid(211, -211);
 574   if (fPairType == fgkKaonPlusKaonPlus ) SetPid(321, 321);
 575   if (fPairType == fgkKaonPlusKaonMinus ) SetPid(321, -321);
 576   if (fPairType == fgkProtonProton ) SetPid(2212, 2212);
 577   if (fPairType == fgkProtonAntiproton ) SetPid(2212, -2212);
 578   if (fPairType == fgkPionPlusKaonPlus ) SetPid(211, 321);
 579   if (fPairType == fgkPionPlusKaonMinus ) SetPid(211, -321);
 580   if (fPairType == fgkPionPlusProton ) SetPid(211, 2212);
 581   if (fPairType == fgkPionPlusAntiproton ) SetPid(211, -2212);
 582   if (fPairType == fgkKaonPlusProton ) SetPid(321, 2212);
 583   if (fPairType == fgkKaonPlusAntiproton ) SetPid(321, -2212);
 584 }
 585
 586 //_____________________________________________
 587 Int_t    AliFemtoModelWeightGeneratorLednicky::GetPairType() const
 588 {
 589   // return pair type
 590   return fPairType;
 591 }
 592
 593 //_____________________________________________
 594 void     AliFemtoModelWeightGeneratorLednicky::SetPairTypeFromPair(AliFemtoPair *aPair)
 595 {
 596   // set calculated system based on the hidden info in the pair
 597   AliFemtoModelHiddenInfo *inf1 = ( AliFemtoModelHiddenInfo *) aPair->Track1()->HiddenInfo();
 598   AliFemtoModelHiddenInfo *inf2 = ( AliFemtoModelHiddenInfo *) aPair->Track2()->HiddenInfo();
 599
 600   const Int_t ktPid1 = inf1->GetPDGPid();
 601   const Int_t ktPid2 = inf2->GetPDGPid();
 602
 603   if      (((ktPid1 ==   211) && (ktPid2 ==   211)) ||
 604            ((ktPid1 ==  -211) && (ktPid2 ==  -211)))
 605     fPairType = fgkPionPlusPionPlus;
 606   else if (((ktPid1 ==  -211) && (ktPid2 ==   211)) ||
 607            ((ktPid1 ==   211) && (ktPid2 ==  -211)))
 608     fPairType = fgkPionPlusPionMinus;
 609   else if (((ktPid1 ==   321) && (ktPid2 ==   321)) ||
 610            ((ktPid1 ==  -321) && (ktPid2 ==  -321)))
 611     fPairType = fgkKaonPlusKaonPlus;
 612   else if (((ktPid1 ==  -321) && (ktPid2 ==   321)) ||
 613            ((ktPid1 ==   321) && (ktPid2 ==  -321)))
 614     fPairType = fgkKaonPlusKaonMinus;
 615   else if (((ktPid1 ==  2212) && (ktPid2 ==  2212)) ||
 616            ((ktPid1 == -2212) && (ktPid2 == -2212)))
 617     fPairType = fgkProtonProton;
 618   else if (((ktPid1 == -2212) && (ktPid2 ==  2212)) ||
 619            ((ktPid1 ==  2212) && (ktPid2 == -2212)))
 620     fPairType = fgkProtonAntiproton;
 621   else if (((ktPid1 ==   211) && (ktPid2 ==   321)) ||
 622            ((ktPid1 ==  -211) && (ktPid2 ==  -321)))
 623     fPairType = fgkPionPlusKaonPlus;
 624   else if (((ktPid1 ==  -211) && (ktPid2 ==   321)) ||
 625            ((ktPid1 ==   211) && (ktPid2 ==  -321)))
 626     fPairType = fgkPionPlusKaonMinus;
 627   else if (((ktPid1 ==   211) && (ktPid2 ==  2212)) ||
 628            ((ktPid1 ==  -211) && (ktPid2 == -2212)))
 629     fPairType = fgkPionPlusProton;
 630   else if (((ktPid1 ==  -211) && (ktPid2 ==  2212)) ||
 631            ((ktPid1 ==   211) && (ktPid2 == -2212)))
 632     fPairType = fgkPionPlusAntiproton;
 633   else if (((ktPid1 ==   321) && (ktPid2 ==  2212)) ||
 634            ((ktPid1 ==  -321) && (ktPid2 == -2212)))
 635     fPairType = fgkKaonPlusProton;
 636   else if (((ktPid1 ==  -321) && (ktPid2 ==  2212)) ||
 637            ((ktPid1 ==   321) && (ktPid2 == -2212)))
 638     fPairType = fgkKaonPlusAntiproton;
 639   SetPid(ktPid1, ktPid2);
 640 }
 641
 642 void AliFemtoModelWeightGeneratorLednicky::SetNuclCharge(const double aNuclCharge) {fNuclCharge=aNuclCharge;FsiNucl();}
 643 void AliFemtoModelWeightGeneratorLednicky::SetNuclMass(const double aNuclMass){fNuclMass=aNuclMass;FsiNucl();}
 644
 645 void AliFemtoModelWeightGeneratorLednicky::SetSphere(){fSphereApp=true;}
 646 void AliFemtoModelWeightGeneratorLednicky::SetSquare(){fSphereApp=false;}
 647 void AliFemtoModelWeightGeneratorLednicky::SetT0ApproxOn(){ fT0App=true;}
 648 void AliFemtoModelWeightGeneratorLednicky::SetT0ApproxOff(){ fT0App=false;}
 649 void AliFemtoModelWeightGeneratorLednicky::SetDefaultCalcPar(){
 650   fItest=1;fIqs=1;fIsi=1;fI3c=0;fIch=1;FsiInit();
 651   fSphereApp=false;fT0App=false;}
 652
 653 void AliFemtoModelWeightGeneratorLednicky::SetCoulOn()    {fItest=1;fIch=1;FsiInit();}
 654 void AliFemtoModelWeightGeneratorLednicky::SetCoulOff()   {fItest=1;fIch=0;FsiInit();}
 655 void AliFemtoModelWeightGeneratorLednicky::SetQuantumOn() {fItest=1;fIqs=1;FsiInit();}
 656 void AliFemtoModelWeightGeneratorLednicky::SetQuantumOff(){fItest=1;fIqs=0;FsiInit();}
 657 void AliFemtoModelWeightGeneratorLednicky::SetStrongOn()  {fItest=1;fIsi=1;FsiInit();}
 658 void AliFemtoModelWeightGeneratorLednicky::SetStrongOff() {fItest=1;fIsi=0;FsiInit();}
 659 void AliFemtoModelWeightGeneratorLednicky::Set3BodyOn()   {fItest=1;fI3c=1;FsiInit();FsiNucl();}
 660 void AliFemtoModelWeightGeneratorLednicky::Set3BodyOff()  {fItest=1;fI3c=0;FsiInit();fWeightDen=1.;FsiNucl();}
 661
 662 Double_t AliFemtoModelWeightGeneratorLednicky::GetKStar() const {return AliFemtoModelWeightGenerator::GetKStar();}
 663 Double_t AliFemtoModelWeightGeneratorLednicky::GetKStarOut() const { return AliFemtoModelWeightGenerator::GetKStarOut(); }
 664 Double_t AliFemtoModelWeightGeneratorLednicky::GetKStarSide() const { return AliFemtoModelWeightGenerator::GetKStarSide(); }
 665 Double_t AliFemtoModelWeightGeneratorLednicky::GetKStarLong() const { return AliFemtoModelWeightGenerator::GetKStarLong(); }
 666 Double_t AliFemtoModelWeightGeneratorLednicky::GetRStar() const { return AliFemtoModelWeightGenerator::GetRStar(); }
 667 Double_t AliFemtoModelWeightGeneratorLednicky::GetRStarOut() const { return AliFemtoModelWeightGenerator::GetRStarOut(); }
 668 Double_t AliFemtoModelWeightGeneratorLednicky::GetRStarSide() const { return AliFemtoModelWeightGenerator::GetRStarSide(); }
 669 Double_t AliFemtoModelWeightGeneratorLednicky::GetRStarLong() const { return AliFemtoModelWeightGenerator::GetRStarLong(); }
 670
 671 AliFemtoModelWeightGenerator* AliFemtoModelWeightGeneratorLednicky::Clone() const {
 672   AliFemtoModelWeightGenerator* tmp = new AliFemtoModelWeightGeneratorLednicky(*this);
 673   return tmp;
 674 }