HBTAN/AliHBTLLWeights.cxx

   1 #include <Riostream.h>
   2 #include "AliHBTLLWeights.h"
   3 #include "AliPDG.h"
   4 #include "AliHBTPair.h"
   5 #include "AliHBTParticle.h"
   6 #include <TList.h>
   7 #include <TRandom.h>
   8 #include <TMath.h>
   9
  10 /*******************************************************************/
  11 /******      ROUTINES    USED    FOR     COMMUNUCATION      ********/
  12 /********************     WITH      FORTRAN     ********************/
  13 /*******************************************************************/
  14 #ifndef WIN32
  15 # define led_bldata led_bldata_
  16 # define fsiini fsiini_
  17 # define ltran12 ltran12_
  18 # define fsiw fsiw_
  19 # define type_of_call
  20 #else
  21 # define led_bldata LED_BLDATA
  22 # define fsiini FSIINI
  23 # define ltran12 LTRAN12
  24 # define fsiw FSIW
  25 # define type_of_call _stdcall
  26 #endif
  27 /****************************************************************/
  28 extern "C" void type_of_call led_bldata();
  29 extern "C" void type_of_call fsiini();
  30 extern "C" void type_of_call ltran12();
  31 extern "C" void type_of_call fsiw();
  32 /**************************************************************/
  33
  34 ClassImp(AliHBTLLWeights)
  35
  36 //Interface to Richard Lednicky's weight alghorithm (Fortran implementation)
  37 //Ludmila Malinina JINR
  38 //Piotr Krzysztof Skowronski CERN/WUT
  39
  40 AliHBTLLWeights* AliHBTLLWeights::fgLLWeights=NULL;
  41
  42 AliHBTLLWeights::AliHBTLLWeights()
  43 {
  44 // Default Constructor
  45     fPID1 = 0;
  46     fPID2 = 0;
  47     SetRandomPosition();
  48     SetColWithResidNuclSwitch();
  49     SetStrongInterSwitch();
  50     SetQuantumStatistics();
  51     SetColoumb();
  52     SetTest();
  53 }
  54
  55
  56  AliHBTLLWeights* AliHBTLLWeights::Instance()
  57 {
  58   if (fgLLWeights)
  59    {
  60      return fgLLWeights;
  61    }
  62   else
  63    {
  64     fgLLWeights = new AliHBTLLWeights();
  65     return fgLLWeights;
  66   }
  67 }
  68
  69
  70 Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair)
  71 {
  72     AliHBTParticle *part1 = partpair->Particle1();
  73     AliHBTParticle *part2 = partpair->Particle2();
  74
  75     if ( (part1 == 0x0) || (part2 == 0x0))
  76      {
  77       Error("GetWeight","Null particle pointer");
  78       return 0.0;
  79      }
  80
  81     Double_t part1Momentum[]={part1->Px(),part1->Py(),part1->Pz()};
  82     Double_t part2Momentum[]={part2->Px(),part2->Py(),part2->Pz()};
  83
  84
  85     if ( (part1->Px() == part2->Px()) && (part1->Py() == part2->Py())
  86        && (part1->Pz() == part2->Pz()) )
  87       {//if particles have the same momentum
  88         return 0.0;
  89       }
  90
  91
  92     if ((!fRandomPosition) && (part1->Vx()  == part2->Vx()) && (part1->Vy()  == part2->Vy())
  93           && (part1->Vz()  == part2->Vz()) )
  94       { //if particles have the same position
  95        return 0.0;
  96       }
  97
  98 //put momenta of particles in LAB into Fortran commons
  99
 100     FSI_MOM.P1X=part1Momentum[0];
 101     FSI_MOM.P1Y=part1Momentum[1];
 102     FSI_MOM.P1Z=part1Momentum[2];
 103
 104     FSI_MOM.P2X=part2Momentum[0];
 105     FSI_MOM.P2Y=part2Momentum[1];
 106     FSI_MOM.P2Z=part2Momentum[2];
 107
 108     if (fRandomPosition)
 109      {
 110        Double_t rxcm = fsigma*gRandom->Gaus();
 111        Double_t rycm = fsigma*gRandom->Gaus();
 112        Double_t rzcm = fsigma*gRandom->Gaus();
 113
 114        FSI_PRF.X=rxcm*fwcons;
 115        FSI_PRF.Y=rycm*fwcons;
 116        FSI_PRF.Z=rzcm*fwcons;
 117        FSI_PRF.T=0.;
 118
 119        Double_t rps=rxcm*rxcm+rycm*rycm+rzcm*rzcm;
 120        Double_t rp=TMath::Sqrt(rps);
 121        FSI_PRF.RP=rp;
 122        FSI_PRF.RPS=rps;
 123      }
 124
 125     ltran12();
 126     fsiw();
 127     return LEDWEIGHT.WEIN;
 128  }
 129
 130 /************************************************************/
 131 void AliHBTLLWeights::Init()
 132  {
 133 //---------------------------------------------------------------------
 134
 135 //initial parameters of model
 136
 137   FSI_NS.NS = approximationModel;
 138
 139   if(!ftest) LEDWEIGHT.ITEST=0;
 140
 141   if(ftest)
 142    {
 143      LEDWEIGHT.ITEST=1;
 144      if(fColoumbSwitch) FSI_NS.ICH =1;
 145      else FSI_NS.ICH=0;
 146
 147      if(fStrongInterSwitch) FSI_NS.ISI=1;
 148      else FSI_NS.ISI=0;
 149
 150      if(fQuantStatSwitch) FSI_NS.IQS=1;
 151      else FSI_NS.IQS=0;
 152
 153      if(fColWithResidNuclSwitch) FSI_NS.I3C=1;
 154      else FSI_NS.I3C=0;
 155    }
 156
 157   if(fRandomPosition) LEDWEIGHT.IRANPOS=1;
 158   else LEDWEIGHT.IRANPOS=0;
 159
 160   if ( (fPID1 == 0) || (fPID2 == 0) )
 161    {
 162     Fatal("Init","Particles types are not set");
 163     return;//pro forma
 164    }
 165
 166   FSI_NS.LL = GetPairCode(fPID1,fPID2);
 167
 168   if (FSI_NS.LL == 0)
 169    {
 170      Fatal("Init","Particles types are not supported");
 171      return;//pro forma
 172    }
 173
 174   TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
 175   if (tpart1 == 0x0)
 176    {
 177      Fatal("init","We can not find particle with ID=%d in our DataBase",fPID1);
 178      return;
 179    }
 180
 181   FSI_POC.AM1=tpart1->Mass();
 182   FSI_POC.C1=tpart1->Charge();
 183
 184   TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
 185   if (tpart2 == 0x0)
 186    {
 187      Fatal("init","We can not find particle with ID=%d in our DataBase",fPID2);
 188      return;
 189    }
 190
 191   FSI_POC.AM2=tpart2->Mass();
 192   FSI_POC.C1=tpart2->Charge();
 193
 194   led_bldata();
 195   fsiini();
 196
 197
 198 //constants for radii simulation
 199
 200   if(fRandomPosition)
 201    {
 202     fsigma =TMath::Sqrt(2.)*fRadius;
 203     fwcons =FSI_CONS.W;
 204    }
 205 }
 206
 207 Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
 208 {
 209  return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
 210 }
 211
 212 Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
 213 {
 214 //   pairCode   1  2  3  4   5    6   7  8  9 10  11  12  13  14 15 16 17 18  19  20   21   22  23 24 25 26    27     28
 215 //   hpid:      n  p  n alfa pi+ pi0 pi+ n  p pi+ pi+ pi+ pi- K+ K+ K+ K-  d  d    t   t    K0  K0  d p  p      p      n
 216 //   lpid:      n  p  p alfa pi- pi0 pi+ d  d  K-  K+  p   p  K- K+ p  p   d alfa  t  alfa  K0  K0b t t alfa lambda lambda
 217 //   NS=1 y/n:  +  +  +  +   +    -   -  -  -  -   -   -   -  -  -  -  -   -  -    -    -    -  -   - -  -      -      -
 218
 219 //alphas, deuterons and tyts are NOT supported here
 220
 221   Int_t chargefactor = 1;
 222   Int_t hpid; //pid in higher row
 223   Int_t lpid; //pid in lower row
 224   Int_t code; //pairCode
 225
 226   Bool_t swap;
 227
 228 //determine the order of selcetion in switch
 229   if (TMath::Abs(pid1) < TMath::Abs(pid2) )
 230    {
 231     if (pid1<0) chargefactor=-1;
 232     hpid=pid2*chargefactor;
 233     lpid=pid1*chargefactor;
 234     swap = kFALSE;
 235    }
 236   else
 237    {
 238     if (pid2<0) chargefactor=-1;
 239     hpid=pid1*chargefactor;
 240     lpid=pid2*chargefactor;
 241     swap = kTRUE;
 242    }
 243
 244 //Determine the pair code
 245   switch (hpid) //switch on first  particle id
 246    {
 247      case kNeutron:
 248       switch (lpid)
 249        {
 250          case kNeutron:
 251            code = 1;  //neutron neutron
 252            break;
 253
 254          case kProton:
 255            code = 3;  //neutron proton
 256            break;
 257
 258          case kLambda0:
 259            code = 28;  //neutron lambda
 260            break;
 261
 262          default:
 263            return 0; //given pair not supported
 264            break;
 265        }
 266       break;
 267
 268      case kProton:
 269       switch (lpid)
 270        {
 271          case kProton:
 272            code = 2; //proton proton
 273            break;
 274
 275          case kLambda0:
 276            code = 27;//proton lambda
 277            break;
 278
 279          default:
 280            return 0; //given pair not supported
 281            break;
 282
 283        }
 284       break;
 285
 286      case kPiPlus:
 287       switch (lpid)
 288        {
 289          case kPiPlus:
 290            code = 7; //piplus piplus
 291            break;
 292
 293          case kPiMinus:
 294            code = 5; //piplus piminus
 295            break;
 296
 297          case kKMinus:
 298            code = 10; //piplus Kminus
 299            break;
 300
 301          case kKPlus:
 302            code = 11; //piplus Kplus
 303            break;
 304
 305          case kProton:
 306            code = 12; //piplus proton
 307            chargefactor*=-1;
 308            break;
 309
 310          default:
 311            return 0; //given pair not supported
 312            break;
 313        }
 314       break;
 315      case kPi0:
 316       switch (lpid)
 317        {
 318          case kPi0:
 319            code = 6;
 320            break;
 321
 322          default:
 323            return 0; //given pair not supported
 324            break;
 325        }
 326       break;
 327
 328      case kKPlus:
 329       switch (lpid)
 330        {
 331          case kKMinus:
 332            code = 14; //Kplus Kminus
 333            break;
 334
 335          case kKPlus:
 336            code = 15; //Kplus Kplus
 337            break;
 338
 339          case kProton:
 340            code = 16; //Kplus proton
 341            break;
 342
 343          default:
 344            return 0; //given pair not supported
 345            break;
 346        }
 347       break;
 348
 349      case kKMinus:
 350       switch (lpid)
 351        {
 352          case kProton:
 353            code = 17; //Kminus proton
 354            chargefactor*=1;
 355            break;
 356
 357          default:
 358            return 0; //given pair not supported
 359            break;
 360        }
 361       break;
 362
 363      case kK0:
 364       switch (lpid)
 365        {
 366          case kK0:
 367            code = 2; //Kzero Kzero
 368            break;
 369
 370          case kK0Bar:
 371            code = 17; //Kzero KzeroBar
 372            break;
 373
 374          default:
 375            return 0; //given pair not supported
 376            break;
 377        }
 378       break;
 379
 380      default: return 0;
 381    }
 382   return code;
 383 }
 384