HBTAN/AliHBTLLWeights.cxx

   1 #include "AliHBTLLWeights.h"
   2 #include "AliPDG.h"
   3 #include "AliHBTPair.h"
   4 #include "AliHBTParticle.h"
   5 #include <TList.h>
   6 #include <TRandom.h>
   7 #include <TMath.h>
   8
   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         return LEDWEIGHT.WEIN;
 125
 126  }
 127
 128 /************************************************************/
 129 void AliHBTLLWeights::Init()
 130  {
 131 //---------------------------------------------------------------------
 132
 133 //initial parameters of model
 134
 135       FSI_NS.NS = approximationModel;
 136
 137       if(!ftest){LEDWEIGHT.ITEST=0;}
 138
 139       if(ftest){
 140       LEDWEIGHT.ITEST=1;
 141       if(fColoumbSwitch){FSI_NS.ICH =1;}
 142       else{FSI_NS.ICH=0;}
 143       if(fStrongInterSwitch){FSI_NS.ISI=1;}
 144       else{FSI_NS.ISI=0;}
 145       if(fQuantStatSwitch){FSI_NS.IQS=1;}
 146       else{FSI_NS.IQS=0;}
 147       if(fColWithResidNuclSwitch){FSI_NS.I3C=1;}
 148       else{FSI_NS.I3C=0;}
 149       }
 150
 151       if(fRandomPosition){LEDWEIGHT.IRANPOS=1;}
 152       else{LEDWEIGHT.IRANPOS=0;}
 153
 154
 155      if ( (fPID1 == 0) || (fPID2 == 0) )
 156       {
 157         Fatal("Init","Particles types are not set");
 158         return;//pro forma
 159       }
 160
 161       FSI_NS.LL = GetPairCode(fPID1,fPID2);
 162
 163       if (FSI_NS.LL == 0)
 164        {
 165          Fatal("Init","Particles types are not supported");
 166          return;//pro forma
 167        }
 168
 169
 170       TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1);
 171       if (tpart1 == 0x0)
 172        {
 173          Fatal("init","We can not find particle with ID=%d in our DataBase",fPID1);
 174          return;
 175        }
 176
 177       FSI_POC.AM1=tpart1->Mass();
 178       FSI_POC.C1=tpart1->Charge();
 179
 180       TParticlePDG* tpart2 = TDatabasePDG::Instance()->GetParticle(fPID2);
 181 //mlv
 182
 183
 184
 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        fsigma =TMath::Sqrt(2.)*fRadius;
 202        fwcons =FSI_CONS.W;
 203       }
 204 }
 205
 206 Int_t AliHBTLLWeights::GetPairCode(const AliHBTPair* partpair)
 207 {
 208  return GetPairCode(partpair->Particle1()->GetPdgCode(),partpair->Particle2()->GetPdgCode());
 209 }
 210
 211 Int_t AliHBTLLWeights::GetPairCode(Int_t pid1,Int_t pid2)
 212 {
 213 //   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
 214 //   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
 215 //   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
 216 //   NS=1 y/n:  +  +  +  +   +    -   -  -  -  -   -   -   -  -  -  -  -   -  -    -    -    -  -   - -  -      -      -
 217
 218 //alphas, deuterons and tyts are NOT supported here
 219
 220   Int_t chargefactor = 1;
 221   Int_t hpid; //pid in higher row
 222   Int_t lpid; //pid in lower row
 223   Int_t code; //pairCode
 224
 225   Bool_t swap;
 226
 227 //determine the order of selcetion in switch
 228   if (TMath::Abs(pid1) < TMath::Abs(pid2) )
 229    {
 230     if (pid1<0) chargefactor=-1;
 231     hpid=pid2*chargefactor;
 232     lpid=pid1*chargefactor;
 233     swap = kFALSE;
 234    }
 235   else
 236    {
 237     if (pid2<0) chargefactor=-1;
 238     hpid=pid1*chargefactor;
 239     lpid=pid2*chargefactor;
 240     swap = kTRUE;
 241    }
 242
 243 //mlv
 244    hpid=pid1;
 245    lpid=pid2;
 246
 247
 248 //Determine the pair code
 249   switch (hpid) //switch on first  particle id
 250    {
 251      case kNeutron:
 252       switch (lpid)
 253        {
 254          case kNeutron:
 255            code = 1;  //neutron neutron
 256            break;
 257
 258          case kProton:
 259            code = 3;  //neutron proton
 260            break;
 261
 262          case kLambda0:
 263            code = 28;  //neutron lambda
 264            break;
 265
 266          default:
 267            return 0; //given pair not supported
 268            break;
 269        }
 270       break;
 271
 272      case kProton:
 273       switch (lpid)
 274        {
 275          case kProton:
 276            code = 2; //proton proton
 277            break;
 278
 279          case kLambda0:
 280            code = 27;//proton lambda
 281            break;
 282
 283          default:
 284            return 0; //given pair not supported
 285            break;
 286
 287        }
 288       break;
 289
 290      case kPiPlus:
 291
 292       switch (lpid)
 293        {
 294          case kPiPlus:
 295            code = 7; //piplus piplus
 296            break;
 297
 298          case kPiMinus:
 299            code = 5; //piplus piminus
 300            break;
 301
 302          case kKMinus:
 303            code = 10; //piplus Kminus
 304            break;
 305
 306          case kKPlus:
 307            code = 11; //piplus Kplus
 308            break;
 309
 310          case kProton:
 311            code = 12; //piplus proton
 312            chargefactor*=-1;
 313            break;
 314
 315          default:
 316            return 0; //given pair not supported
 317            break;
 318        }
 319       break;
 320      case kPi0:
 321       switch (lpid)
 322        {
 323          case kPi0:
 324            code = 6;
 325            break;
 326
 327          default:
 328            return 0; //given pair not supported
 329            break;
 330        }
 331       break;
 332
 333      case kKPlus:
 334       switch (lpid)
 335        {
 336          case kKMinus:
 337            code = 14; //Kplus Kminus
 338            break;
 339
 340          case kKPlus:
 341            code = 15; //Kplus Kplus
 342            break;
 343
 344          case kProton:
 345            code = 16; //Kplus proton
 346            break;
 347
 348          default:
 349            return 0; //given pair not supported
 350            break;
 351        }
 352       break;
 353
 354      case kKMinus:
 355       switch (lpid)
 356        {
 357          case kProton:
 358            code = 17; //Kminus proton
 359            chargefactor*=1;
 360            break;
 361
 362          default:
 363            return 0; //given pair not supported
 364            break;
 365        }
 366       break;
 367
 368      case kK0:
 369       switch (lpid)
 370        {
 371          case kK0:
 372            code = 2; //Kzero Kzero
 373            break;
 374
 375          case kK0Bar:
 376            code = 17; //Kzero KzeroBar
 377            break;
 378
 379          default:
 380            return 0; //given pair not supported
 381            break;
 382        }
 383       break;
 384
 385      default: return 0;
 386    }
 387   return code;
 388 }
 389