HBTAN/AliHBTLLWeights.cxx

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