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 }