HBTAN/AliHBTLLWeights.cxx

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