StRoot: StHbtMaker/ThCorrFctn/StHbtFsiLednicky.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: 
00004  *
00005  * Author: Laurent Conin, Fabrice Retiere, Subatech, France
00006  ***************************************************************************
00007  *
00008  * Description : Calculate correlation weight using R.Lednicky's code 
00009  *  Use the fortran files : FsiLednickyWeight.F and FsiTools.F
00010  *
00011  ***************************************************************************
00012  *
00013  * $Log: 
00014  *
00015  ***************************************************************************/
00016 
00017 #include "StHbtMaker/ThCorrFctn/StHbtFsiLednicky.h"
00018 #include "StarCallf77.h"
00019 #include <Stsstream.h>
00020 
00021 
00022 #ifdef SOLARIS
00023 # ifndef false
00024 typedef int bool;
00025 #define false 0
00026 #define true 1
00027 # endif
00028 #endif
00029 
00030 // --- Prototype of the function used in the weight calculator 
00031 //     (in FsiWeightLedinicky.F)
00032 #define fsiin F77_NAME(fsiin,FSIIN)
00033 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);}
00034 #define llini F77_NAME(llini,LLINI)
00035 extern "C" {void type_of_call F77_NAME(llini,LLINI)(const int &lll,const int &ns, const int &itest);}
00036 
00037 #define fsinucl F77_NAME(fsinucl,FSINUCL)
00038 extern "C" {void type_of_call  F77_NAME(fsinucl,FSINUCL)(const double &mn,const double &cn);}
00039 #define fsimomentum F77_NAME(fsimomentum,FSIMOMENTUM)
00040 extern "C" {void type_of_call F77_NAME(fsimomentum,FSIMOMENTUM)(double &p1,double &p2);}
00041 #define fsiposition F77_NAME(fsiposition,FSIPOSITION)
00042 extern "C" {void type_of_call F77_NAME(fsiposition,FSIPOSITION)(double &x1,double &x2);}
00043 #define fsiw F77_NAME(fsiw,FSIW)
00044 extern "C" {void type_of_call F77_NAME(fsiw,FSIW)(const int &i,double &weif,
00045                                                   double &wei,double &wein);}
00046 #define ltran12 F77_NAME(ltran12,LTRAN12)
00047 extern "C" {void type_of_call ltran12_();}
00048 
00049 // Test function for Lambda potential
00050 #define printlam F77_NAME(printlam,PRINTLAM)
00051 extern "C" {void type_of_call printlam_();}
00052 
00053 // --- Additional prototyping of some CERN functions (in FsiTool.F)
00054 typedef float   REAL;
00055 typedef struct { REAL re; REAL im; } COMPLEX;
00056 #define cgamma F77_NAME(cgamma,CGAMMA)
00057 extern "C" {COMPLEX type_of_call cgamma_(COMPLEX*);}
00058 
00059 #ifdef __ROOT__
00060 ClassImp(StHbtFsiLednicky)
00061 #endif
00062 
00063 StHbtFsiLednicky::StHbtFsiLednicky() : StHbtFsiWeight(),
00064   mItest(0),mIch(1),mIqs(1),mIsi(1),mI3c(0),
00065   mNuclMass(1.),mNuclCharge(0.),
00066   mSphereApp(false),mT0App(false),mNuclChargeSign(1), mLLMax(30),mNumbNonId(0) {
00067 
00068   mLLName=new char*[mLLMax+1];
00069   mNumProcessPair=new int[mLLMax+1];
00070   int i;
00071   for (i=1;i<=mLLMax;i++) {mLLName[i]=new char[40];mNumProcessPair[i]=0;}
00072   strcpy( mLLName[1],"neutron neutron");
00073   strcpy( mLLName[2],"proton proton");
00074   strcpy( mLLName[3],"neutron proton");
00075   strcpy( mLLName[4],"alpha alpha");
00076   strcpy( mLLName[5],"pi+ pi-");
00077   strcpy( mLLName[6],"pi0 pi0");
00078   strcpy( mLLName[7],"pi+ pi+");
00079   strcpy( mLLName[8],"neutron deuteron");
00080   strcpy( mLLName[9],"proton deuteron");
00081   strcpy( mLLName[10],"pi+ K-");
00082   strcpy( mLLName[11],"pi+ K+");
00083   strcpy( mLLName[12],"pi+ proton");
00084   strcpy( mLLName[13],"pi- proton");
00085   strcpy( mLLName[14],"K+ K-");
00086   strcpy( mLLName[15],"K+ K+");
00087   strcpy( mLLName[16],"K+ proton");
00088   strcpy( mLLName[17],"K- proton");
00089   strcpy( mLLName[18],"deuteron deuteron");
00090   strcpy( mLLName[19],"deuton alpha");
00091   strcpy( mLLName[20],"triton triton");
00092   strcpy( mLLName[21],"triton alpha");
00093   strcpy( mLLName[22],"K0 K0");
00094   strcpy( mLLName[23],"K0 K0b");
00095   strcpy( mLLName[24],"deuteron triton");
00096   strcpy( mLLName[25],"proton triton");
00097   strcpy( mLLName[26],"proton alpha");
00098   strcpy( mLLName[27],"proton lambda");
00099   strcpy( mLLName[28],"neutron lambda");
00100   strcpy( mLLName[29],"Lambda lambda");// gael 21May02
00101   strcpy( mLLName[30],"Proton Anti-proton");// gael 21May02
00102   FsiInit();
00103   FsiNucl();
00104 };
00105 
00106 
00107 double StHbtFsiLednicky::GetWeight(const StHbtThPair* aThPair){
00108 
00109   if (!SetPid(aThPair->GetPid1(),aThPair->GetPid2())) {
00110     mWeightDen=1.;
00111     return 1;    
00112   } else { // Good Pid
00113     const StHbtLorentzVector*  p;
00114     p=(aThPair->GetRealMomentum1());
00115     double p1[]={p->px(),p->py(),p->pz()};
00116     p=(aThPair->GetRealMomentum2());
00117     double p2[]={p->px(),p->py(),p->pz()};
00118     if ((p1[0]==p2[0])&&(p1[1]==p2[1])&&(p1[2]==p2[2])) {
00119       mWeightDen=0.;
00120       return 0;  
00121     } 
00122     if (mSwap) {
00123       fsimomentum(*p2,*p1);
00124     } else {
00125       fsimomentum(*p1,*p2);
00126     }
00127     p=(aThPair->GetEmPoint1());
00128 //     cout << "Pid1:dans GetWeight = " << aThPair->GetPid1() << endl;
00129 //     cout << "Pid2:dans GetWeight = " << aThPair->GetPid2() << endl;
00130 //     cout << "LL:in GetWeight = " << mLL << endl;
00131     double x1[]={p->x(),p->y(),p->z(),p->t()};
00132     p=(aThPair->GetEmPoint2());
00133     double x2[]={p->x(),p->y(),p->z(),p->t()};
00134     if ((x1[0]==x2[0])&&(x1[1]==x2[1])&&(x1[2]==x2[2])&&(x1[3]==x2[3])) {
00135       mWeightDen=0.;
00136       return 0;  
00137     } 
00138     if (mSwap) {
00139       fsiposition(*x1,*x2);
00140     } else {
00141       fsiposition(*x2,*x1);
00142     }
00143     FsiSetLL();
00144     ltran12();
00145     fsiw(1,mWeif,mWei,mWein);
00146     if (mI3c==0) return mWein;
00147     mWeightDen=mWeif;
00148     return mWei;
00149   };
00150 };
00151 
00152 StHbtString StHbtFsiLednicky::Report() {
00153   ostrstream tStr; 
00154   tStr << "Lednicky afterburner calculation for  Correlation -  Report" << endl;
00155   tStr << "    Setting : Quantum : " << ((mIqs) ? "On" : "Off"); 
00156   tStr << " - Coulbomb : " << ((mIch) ? "On" : "Off") ;
00157   tStr << " - Strong : " << ((mIsi) ? "On" : "Off");
00158   tStr << endl;
00159   tStr << "              3-Body : " << ((mI3c) ? "On"  : "Off") ;
00160   if (mI3c) tStr << " Mass=" <<  mNuclMass << " - Charge= " << mNuclCharge ;
00161   tStr << endl;
00162   tStr << "    " << mNumProcessPair[0] << " Pairs have been Processed :" << endl;
00163   int i;
00164   for(i=1;i<=mLLMax;i++) { 
00165     if (mNumProcessPair[i])
00166       tStr << "         "<<  setw(8) << mNumProcessPair[i] << " " << mLLName[i] << endl;
00167   }
00168   if (mNumbNonId)
00169     tStr << "         "<< setw(8) << mNumbNonId << " Non Identified" << endl;
00170   StHbtString returnThis = tStr.str();
00171   return returnThis;
00172 }
00173 
00174 void StHbtFsiLednicky::FsiInit(){
00175   cout << "*******************StHbtFsiLednicky check FsiInit ************" << endl;
00176   cout <<"mItest dans FsiInit() = " << mItest << endl;
00177   cout <<"mIch dans FsiInit() = " << mIch << endl;
00178   cout <<"mIqs dans FsiInit() = " << mIqs << endl;
00179   cout <<"mIsi dans FsiInit() = " << mIsi << endl;
00180   cout <<"mI3c dans FsiInit() = " << mI3c << endl;
00181   fsiin(mItest,mIch,mIqs,mIsi,mI3c);
00182 };
00183 
00184 void StHbtFsiLednicky::FsiNucl(){
00185   cout << "*******************StHbtFsiLednicky check FsiNucl ************" << endl;
00186   cout <<"mNuclMass dans FsiNucl() = " << mNuclMass << endl;
00187   cout <<"mNuclCharge dans FsiNucl() = " << mNuclCharge << endl;
00188   cout <<"mNuclChargeSign dans FsiNucl() = " << mNuclChargeSign << endl;
00189   fsinucl(mNuclMass,mNuclCharge*mNuclChargeSign);
00190 };
00191 
00192 void StHbtFsiLednicky::FsiSetLL(){
00193   int tNS;
00194   if (mSphereApp||(mLL>5)) {
00195     if (mT0App) { tNS=4;} 
00196     else {tNS=2;};
00197   } else { tNS=1;};
00198 //   cout <<"mLL dans FsiSetLL() = "<< mLL << endl;
00199 //   cout <<"tNS dans FsiSetLL() = "<< tNS << endl;
00200 //   cout <<"mItest dans FsiSetLL() = "<< mItest << endl;
00201   llini(mLL,tNS,mItest);
00202 }
00203          
00204 bool StHbtFsiLednicky::SetPid(const int aPid1,const int aPid2) {
00205 
00206   static const int sPi0Pid=111;
00207   static const int sPionPid=211; 
00208   static const int sK0Pid=311;
00209   static const int sKPid=321;
00210   static const int sNeutPid=2112;
00211   static const int sProtPid=2212;
00212   static const int sLamPid=3122;
00213   //  static const int sLamLamPid=3122;
00214 
00215   //  cout << "Setting PID to " << aPid1 << " " << aPid2 << endl;
00216 
00217   int tPidl,tPidh;
00218   int tChargeFactor=1;
00219   
00220   if (abs(aPid1)<abs(aPid2)) {
00221     if (aPid1<0) tChargeFactor=-1;
00222     tPidl=aPid1*tChargeFactor;
00223     tPidh=aPid2*tChargeFactor;
00224     mSwap=false;
00225   } else {
00226     if (aPid2<0) tChargeFactor=-1;
00227     tPidl=aPid2*tChargeFactor;
00228     tPidh=aPid1*tChargeFactor;
00229     mSwap=true;
00230   }
00231   switch (tPidl) {
00232   case sPionPid:
00233     switch (tPidh) {
00234     case -sPionPid:   mLL=5; tChargeFactor*=1 ;break;
00235     case sPionPid:    mLL=7; tChargeFactor*=1 ;break;
00236     case -sKPid:      mLL=10;tChargeFactor*=1 ;break;  
00237     case sKPid:       mLL=11;tChargeFactor*=1 ;break;  
00238     case sProtPid:    mLL=12;tChargeFactor*=1 ;break;
00239     case -sProtPid:   mLL=13;tChargeFactor*=-1;break;
00240     default: mLL=0;
00241     }
00242     break;
00243   case sProtPid:
00244     switch (tPidh) {
00245     case sProtPid:    mLL=2; tChargeFactor*=1 ;break;
00246     case sLamPid:     mLL=27;tChargeFactor*=1 ;break;
00247     case -sProtPid:   mLL=30;tChargeFactor*=1 ;break;
00248     default: mLL=0;
00249     };
00250     break;
00251   case sKPid:
00252     switch (tPidh) {
00253     case -sKPid:      mLL=14;tChargeFactor*=1 ;break;
00254     case sKPid:       mLL=15;tChargeFactor*=1 ;break;
00255     case sProtPid:    mLL=16;tChargeFactor*=1 ;break;
00256     case -sProtPid:   mLL=17;tChargeFactor*=-1 ;break;
00257     default: mLL=0;
00258     };
00259     break;    
00260   case sK0Pid:
00261     switch (tPidh) {
00262     case sK0Pid:         mLL=22;tChargeFactor*=1 ;break;
00263     case -sK0Pid:        mLL=23;tChargeFactor*=1 ;break;
00264     default: mLL=0;
00265     };
00266     break;   
00267   case sPi0Pid:
00268     switch (tPidh) {
00269     case sPi0Pid:        mLL=6; tChargeFactor*=1 ;break;
00270     default: mLL=0;
00271     };
00272     break;
00273   case sNeutPid:
00274     switch (tPidh) {
00275     case sNeutPid:      mLL=1; tChargeFactor*=1 ;break;
00276     case sProtPid:      mLL=3; tChargeFactor*=1 ;break;
00277     case sLamPid:       mLL=28;tChargeFactor*=1 ;break;
00278     default: mLL=0;
00279     };
00280     break;                                             //Gael 21May02 
00281   case sLamPid:                                        //Gael 21May02 
00282     switch (tPidh) {                                   //Gael 21May02 
00283     case sLamPid:       mLL=29;tChargeFactor*=1 ;break;//Gael 21May02  
00284     default: mLL=0;                                    //Gael 21May02 
00285     };                                                 //Gael 21May02 
00286     break;                                             //Gael 21May02 
00287   default: mLL=0;
00288   }
00289   if (tChargeFactor!=mNuclChargeSign) {
00290     mNuclChargeSign=tChargeFactor;
00291     FsiNucl();
00292   }
00293   (mNumProcessPair[0])++;
00294   if (mLL) {
00295     (mNumProcessPair[mLL])++;
00296     return true;
00297   } else {
00298     mNumbNonId++;
00299     return false;
00300   }
00301   cout << "*******************StHbtFsiLednicky check SetPid ************" << endl;
00302   cout << "mLL=="<< mLL << endl;
00303   cout << "mNuclCharge=="<< mNuclCharge << endl;
00304 
00305 }    
00306 StHbtFsiLednicky::~StHbtFsiLednicky() 
00307 { /* no-op */ };
00308 
00309 
00310 inline  void StHbtFsiLednicky::SetNuclCharge(const double aNuclCharge) {mNuclCharge=aNuclCharge;FsiNucl();};
00311 inline  void StHbtFsiLednicky::SetNuclMass(const double aNuclMass){mNuclMass=aNuclMass;FsiNucl();};
00312 
00313 inline  void StHbtFsiLednicky::SetSphere(){mSphereApp=true;};
00314 inline  void StHbtFsiLednicky::SetSquare(){mSphereApp=false;};
00315 inline  void StHbtFsiLednicky::SetT0ApproxOn(){ mT0App=true;};
00316 inline  void StHbtFsiLednicky::SetT0ApproxOff(){ mT0App=false;};
00317 inline  void StHbtFsiLednicky::SetDefaultCalcPar(){
00318   mItest=1;mIqs=1;mIsi=1;mI3c=0;mIch=1;FsiInit();
00319   mSphereApp=false;mT0App=false;};
00320 
00321 inline  void StHbtFsiLednicky::SetCoulOn(){mItest=1;mIch=1;FsiInit();};
00322 inline  void StHbtFsiLednicky::SetCoulOff(){mItest=1;mIch=0;FsiInit();};
00323 inline  void StHbtFsiLednicky::SetQuantumOn(){mItest=1;mIqs=1;FsiInit();};
00324 inline  void StHbtFsiLednicky::SetQuantumOff(){mItest=1;mIqs=0;FsiInit();};
00325 inline  void StHbtFsiLednicky::SetStrongOn(){mItest=1;mIsi=1;FsiInit();};
00326 inline  void StHbtFsiLednicky::SetStrongOff(){mItest=1;mIsi=0;FsiInit();};
00327 inline  void StHbtFsiLednicky::Set3BodyOn(){mItest=1;mI3c=1;FsiInit();FsiNucl();};
00328 inline  void StHbtFsiLednicky::Set3BodyOff(){mItest=1;mI3c=0;FsiInit();mWeightDen=1.;FsiNucl();};
00329 
00330 void StHbtFsiLednicky::PrintLambdas() {printlam();};
00331 
00332   
00333