StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StHbtFsiLednicky.cxx
1 /***************************************************************************
2  *
3  *
4  *
5  * Author: Laurent Conin, Fabrice Retiere, Subatech, France
6  ***************************************************************************
7  *
8  * Description : Calculate correlation weight using R.Lednicky's code
9  * Use the fortran files : FsiLednickyWeight.F and FsiTools.F
10  *
11  ***************************************************************************
12  *
13  *
14  *
15  ***************************************************************************/
16 
17 #include "StHbtMaker/ThCorrFctn/StHbtFsiLednicky.h"
18 #include "StarCallf77.h"
19 #include <Stsstream.h>
20 
21 
22 #ifdef SOLARIS
23 # ifndef false
24 typedef int bool;
25 #define false 0
26 #define true 1
27 # endif
28 #endif
29 
30 // --- Prototype of the function used in the weight calculator
31 // (in FsiWeightLedinicky.F)
32 #define fsiin F77_NAME(fsiin,FSIIN)
33 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);}
34 #define llini F77_NAME(llini,LLINI)
35 extern "C" {void type_of_call F77_NAME(llini,LLINI)(const int &lll,const int &ns, const int &itest);}
36 
37 #define fsinucl F77_NAME(fsinucl,FSINUCL)
38 extern "C" {void type_of_call F77_NAME(fsinucl,FSINUCL)(const double &mn,const double &cn);}
39 #define fsimomentum F77_NAME(fsimomentum,FSIMOMENTUM)
40 extern "C" {void type_of_call F77_NAME(fsimomentum,FSIMOMENTUM)(double &p1,double &p2);}
41 #define fsiposition F77_NAME(fsiposition,FSIPOSITION)
42 extern "C" {void type_of_call F77_NAME(fsiposition,FSIPOSITION)(double &x1,double &x2);}
43 #define fsiw F77_NAME(fsiw,FSIW)
44 extern "C" {void type_of_call F77_NAME(fsiw,FSIW)(const int &i,double &weif,
45  double &wei,double &wein);}
46 #define ltran12 F77_NAME(ltran12,LTRAN12)
47 extern "C" {void type_of_call ltran12_();}
48 
49 // Test function for Lambda potential
50 #define printlam F77_NAME(printlam,PRINTLAM)
51 extern "C" {void type_of_call printlam_();}
52 
53 // --- Additional prototyping of some CERN functions (in FsiTool.F)
54 typedef float REAL;
55 typedef struct { REAL re; REAL im; } COMPLEX;
56 #define cgamma F77_NAME(cgamma,CGAMMA)
57 extern "C" {COMPLEX type_of_call cgamma_(COMPLEX*);}
58 
59 #ifdef __ROOT__
60 ClassImp(StHbtFsiLednicky)
61 #endif
62 
63 StHbtFsiLednicky::StHbtFsiLednicky() : StHbtFsiWeight(),
64  mItest(0),mIch(1),mIqs(1),mIsi(1),mI3c(0),
65  mNuclMass(1.),mNuclCharge(0.),
66  mSphereApp(false),mT0App(false),mNuclChargeSign(1), mLLMax(30),mNumbNonId(0) {
67 
68  mLLName=new char*[mLLMax+1];
69  mNumProcessPair=new int[mLLMax+1];
70  int i;
71  for (i=1;i<=mLLMax;i++) {mLLName[i]=new char[40];mNumProcessPair[i]=0;}
72  strcpy( mLLName[1],"neutron neutron");
73  strcpy( mLLName[2],"proton proton");
74  strcpy( mLLName[3],"neutron proton");
75  strcpy( mLLName[4],"alpha alpha");
76  strcpy( mLLName[5],"pi+ pi-");
77  strcpy( mLLName[6],"pi0 pi0");
78  strcpy( mLLName[7],"pi+ pi+");
79  strcpy( mLLName[8],"neutron deuteron");
80  strcpy( mLLName[9],"proton deuteron");
81  strcpy( mLLName[10],"pi+ K-");
82  strcpy( mLLName[11],"pi+ K+");
83  strcpy( mLLName[12],"pi+ proton");
84  strcpy( mLLName[13],"pi- proton");
85  strcpy( mLLName[14],"K+ K-");
86  strcpy( mLLName[15],"K+ K+");
87  strcpy( mLLName[16],"K+ proton");
88  strcpy( mLLName[17],"K- proton");
89  strcpy( mLLName[18],"deuteron deuteron");
90  strcpy( mLLName[19],"deuton alpha");
91  strcpy( mLLName[20],"triton triton");
92  strcpy( mLLName[21],"triton alpha");
93  strcpy( mLLName[22],"K0 K0");
94  strcpy( mLLName[23],"K0 K0b");
95  strcpy( mLLName[24],"deuteron triton");
96  strcpy( mLLName[25],"proton triton");
97  strcpy( mLLName[26],"proton alpha");
98  strcpy( mLLName[27],"proton lambda");
99  strcpy( mLLName[28],"neutron lambda");
100  strcpy( mLLName[29],"Lambda lambda");// gael 21May02
101  strcpy( mLLName[30],"Proton Anti-proton");// gael 21May02
102  FsiInit();
103  FsiNucl();
104 };
105 
106 
107 double StHbtFsiLednicky::GetWeight(const StHbtThPair* aThPair){
108 
109  if (!SetPid(aThPair->GetPid1(),aThPair->GetPid2())) {
110  mWeightDen=1.;
111  return 1;
112  } else { // Good Pid
113  const StHbtLorentzVector* p;
114  p=(aThPair->GetRealMomentum1());
115  double p1[]={p->px(),p->py(),p->pz()};
116  p=(aThPair->GetRealMomentum2());
117  double p2[]={p->px(),p->py(),p->pz()};
118  if ((p1[0]==p2[0])&&(p1[1]==p2[1])&&(p1[2]==p2[2])) {
119  mWeightDen=0.;
120  return 0;
121  }
122  if (mSwap) {
123  fsimomentum(*p2,*p1);
124  } else {
125  fsimomentum(*p1,*p2);
126  }
127  p=(aThPair->GetEmPoint1());
128 // cout << "Pid1:dans GetWeight = " << aThPair->GetPid1() << endl;
129 // cout << "Pid2:dans GetWeight = " << aThPair->GetPid2() << endl;
130 // cout << "LL:in GetWeight = " << mLL << endl;
131  double x1[]={p->x(),p->y(),p->z(),p->t()};
132  p=(aThPair->GetEmPoint2());
133  double x2[]={p->x(),p->y(),p->z(),p->t()};
134  if ((x1[0]==x2[0])&&(x1[1]==x2[1])&&(x1[2]==x2[2])&&(x1[3]==x2[3])) {
135  mWeightDen=0.;
136  return 0;
137  }
138  if (mSwap) {
139  fsiposition(*x1,*x2);
140  } else {
141  fsiposition(*x2,*x1);
142  }
143  FsiSetLL();
144  ltran12();
145  fsiw(1,mWeif,mWei,mWein);
146  if (mI3c==0) return mWein;
147  mWeightDen=mWeif;
148  return mWei;
149  };
150 };
151 
152 StHbtString StHbtFsiLednicky::Report() {
153  std::ostringstream tStr;
154  tStr << "Lednicky afterburner calculation for Correlation - Report" << endl;
155  tStr << " Setting : Quantum : " << ((mIqs) ? "On" : "Off");
156  tStr << " - Coulbomb : " << ((mIch) ? "On" : "Off") ;
157  tStr << " - Strong : " << ((mIsi) ? "On" : "Off");
158  tStr << endl;
159  tStr << " 3-Body : " << ((mI3c) ? "On" : "Off") ;
160  if (mI3c) tStr << " Mass=" << mNuclMass << " - Charge= " << mNuclCharge ;
161  tStr << endl;
162  tStr << " " << mNumProcessPair[0] << " Pairs have been Processed :" << endl;
163  int i;
164  for(i=1;i<=mLLMax;i++) {
165  if (mNumProcessPair[i])
166  tStr << " "<< setw(8) << mNumProcessPair[i] << " " << mLLName[i] << endl;
167  }
168  if (mNumbNonId)
169  tStr << " "<< setw(8) << mNumbNonId << " Non Identified" << endl;
170  StHbtString returnThis = tStr.str();
171  return returnThis;
172 }
173 
174 void StHbtFsiLednicky::FsiInit(){
175  cout << "*******************StHbtFsiLednicky check FsiInit ************" << endl;
176  cout <<"mItest dans FsiInit() = " << mItest << endl;
177  cout <<"mIch dans FsiInit() = " << mIch << endl;
178  cout <<"mIqs dans FsiInit() = " << mIqs << endl;
179  cout <<"mIsi dans FsiInit() = " << mIsi << endl;
180  cout <<"mI3c dans FsiInit() = " << mI3c << endl;
181  fsiin(mItest,mIch,mIqs,mIsi,mI3c);
182 };
183 
184 void StHbtFsiLednicky::FsiNucl(){
185  cout << "*******************StHbtFsiLednicky check FsiNucl ************" << endl;
186  cout <<"mNuclMass dans FsiNucl() = " << mNuclMass << endl;
187  cout <<"mNuclCharge dans FsiNucl() = " << mNuclCharge << endl;
188  cout <<"mNuclChargeSign dans FsiNucl() = " << mNuclChargeSign << endl;
189  fsinucl(mNuclMass,mNuclCharge*mNuclChargeSign);
190 };
191 
192 void StHbtFsiLednicky::FsiSetLL(){
193  int tNS;
194  if (mSphereApp||(mLL>5)) {
195  if (mT0App) { tNS=4;}
196  else {tNS=2;};
197  } else { tNS=1;};
198 // cout <<"mLL dans FsiSetLL() = "<< mLL << endl;
199 // cout <<"tNS dans FsiSetLL() = "<< tNS << endl;
200 // cout <<"mItest dans FsiSetLL() = "<< mItest << endl;
201  llini(mLL,tNS,mItest);
202 }
203 
204 bool StHbtFsiLednicky::SetPid(const int aPid1,const int aPid2) {
205 
206  static const int sPi0Pid=111;
207  static const int sPionPid=211;
208  static const int sK0Pid=311;
209  static const int sKPid=321;
210  static const int sNeutPid=2112;
211  static const int sProtPid=2212;
212  static const int sLamPid=3122;
213  // static const int sLamLamPid=3122;
214 
215  // cout << "Setting PID to " << aPid1 << " " << aPid2 << endl;
216 
217  int tPidl,tPidh;
218  int tChargeFactor=1;
219 
220  if (abs(aPid1)<abs(aPid2)) {
221  if (aPid1<0) tChargeFactor=-1;
222  tPidl=aPid1*tChargeFactor;
223  tPidh=aPid2*tChargeFactor;
224  mSwap=false;
225  } else {
226  if (aPid2<0) tChargeFactor=-1;
227  tPidl=aPid2*tChargeFactor;
228  tPidh=aPid1*tChargeFactor;
229  mSwap=true;
230  }
231  switch (tPidl) {
232  case sPionPid:
233  switch (tPidh) {
234  case -sPionPid: mLL=5; tChargeFactor*=1 ;break;
235  case sPionPid: mLL=7; tChargeFactor*=1 ;break;
236  case -sKPid: mLL=10;tChargeFactor*=1 ;break;
237  case sKPid: mLL=11;tChargeFactor*=1 ;break;
238  case sProtPid: mLL=12;tChargeFactor*=1 ;break;
239  case -sProtPid: mLL=13;tChargeFactor*=-1;break;
240  default: mLL=0;
241  }
242  break;
243  case sProtPid:
244  switch (tPidh) {
245  case sProtPid: mLL=2; tChargeFactor*=1 ;break;
246  case sLamPid: mLL=27;tChargeFactor*=1 ;break;
247  case -sProtPid: mLL=30;tChargeFactor*=1 ;break;
248  default: mLL=0;
249  };
250  break;
251  case sKPid:
252  switch (tPidh) {
253  case -sKPid: mLL=14;tChargeFactor*=1 ;break;
254  case sKPid: mLL=15;tChargeFactor*=1 ;break;
255  case sProtPid: mLL=16;tChargeFactor*=1 ;break;
256  case -sProtPid: mLL=17;tChargeFactor*=-1 ;break;
257  default: mLL=0;
258  };
259  break;
260  case sK0Pid:
261  switch (tPidh) {
262  case sK0Pid: mLL=22;tChargeFactor*=1 ;break;
263  case -sK0Pid: mLL=23;tChargeFactor*=1 ;break;
264  default: mLL=0;
265  };
266  break;
267  case sPi0Pid:
268  switch (tPidh) {
269  case sPi0Pid: mLL=6; tChargeFactor*=1 ;break;
270  default: mLL=0;
271  };
272  break;
273  case sNeutPid:
274  switch (tPidh) {
275  case sNeutPid: mLL=1; tChargeFactor*=1 ;break;
276  case sProtPid: mLL=3; tChargeFactor*=1 ;break;
277  case sLamPid: mLL=28;tChargeFactor*=1 ;break;
278  default: mLL=0;
279  };
280  break; //Gael 21May02
281  case sLamPid: //Gael 21May02
282  switch (tPidh) { //Gael 21May02
283  case sLamPid: mLL=29;tChargeFactor*=1 ;break;//Gael 21May02
284  default: mLL=0; //Gael 21May02
285  }; //Gael 21May02
286  break; //Gael 21May02
287  default: mLL=0;
288  }
289  if (tChargeFactor!=mNuclChargeSign) {
290  mNuclChargeSign=tChargeFactor;
291  FsiNucl();
292  }
293  (mNumProcessPair[0])++;
294  if (mLL) {
295  (mNumProcessPair[mLL])++;
296  return true;
297  } else {
298  mNumbNonId++;
299  return false;
300  }
301  cout << "*******************StHbtFsiLednicky check SetPid ************" << endl;
302  cout << "mLL=="<< mLL << endl;
303  cout << "mNuclCharge=="<< mNuclCharge << endl;
304 
305 }
306 StHbtFsiLednicky::~StHbtFsiLednicky()
307 { /* no-op */ };
308 
309 
310 inline void StHbtFsiLednicky::SetNuclCharge(const double aNuclCharge) {mNuclCharge=aNuclCharge;FsiNucl();};
311 inline void StHbtFsiLednicky::SetNuclMass(const double aNuclMass){mNuclMass=aNuclMass;FsiNucl();};
312 
313 inline void StHbtFsiLednicky::SetSphere(){mSphereApp=true;};
314 inline void StHbtFsiLednicky::SetSquare(){mSphereApp=false;};
315 inline void StHbtFsiLednicky::SetT0ApproxOn(){ mT0App=true;};
316 inline void StHbtFsiLednicky::SetT0ApproxOff(){ mT0App=false;};
317 inline void StHbtFsiLednicky::SetDefaultCalcPar(){
318  mItest=1;mIqs=1;mIsi=1;mI3c=0;mIch=1;FsiInit();
319  mSphereApp=false;mT0App=false;};
320 
321 inline void StHbtFsiLednicky::SetCoulOn(){mItest=1;mIch=1;FsiInit();};
322 inline void StHbtFsiLednicky::SetCoulOff(){mItest=1;mIch=0;FsiInit();};
323 inline void StHbtFsiLednicky::SetQuantumOn(){mItest=1;mIqs=1;FsiInit();};
324 inline void StHbtFsiLednicky::SetQuantumOff(){mItest=1;mIqs=0;FsiInit();};
325 inline void StHbtFsiLednicky::SetStrongOn(){mItest=1;mIsi=1;FsiInit();};
326 inline void StHbtFsiLednicky::SetStrongOff(){mItest=1;mIsi=0;FsiInit();};
327 inline void StHbtFsiLednicky::Set3BodyOn(){mItest=1;mI3c=1;FsiInit();FsiNucl();};
328 inline void StHbtFsiLednicky::Set3BodyOff(){mItest=1;mI3c=0;FsiInit();mWeightDen=1.;FsiNucl();};
329 
330 void StHbtFsiLednicky::PrintLambdas() {printlam();};
331 
332 
333 
Generated by   doxygen 1.8.5