X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=HBTAN%2FAliHBTLLWeights.cxx;h=484355726ea7f00ae395ca472bd2f86f2600abba;hb=9988704266cf8e4b13017cc4e53d86463471b873;hp=e53acfd388fe063cb864252153e6a21bf54ccdf7;hpb=88cb7938ca21d4a80991d4e7aa564008c29340f7;p=u%2Fmrichter%2FAliRoot.git diff --git a/HBTAN/AliHBTLLWeights.cxx b/HBTAN/AliHBTLLWeights.cxx index e53acfd388f..484355726ea 100644 --- a/HBTAN/AliHBTLLWeights.cxx +++ b/HBTAN/AliHBTLLWeights.cxx @@ -83,14 +83,14 @@ //===================================================================== // At the beginning of calculation user should call FSIINI, // which reads LL, NS, ITEST (and eventually ICH, IQS, ISI, I3C) -// and initializes various parameters. +// and ializes various parameters. // In particular the constants in // COMMON/FSI_CONS/PI,PI2,SPI,DR,W // may be useful for the user: // W=1/.1973D0 ! from fm to 1/GeV // PI=4*DATAN(1.D0) // PI2=2*PI -// SPI=DSQRT(PI) +// SPI=TMath::Sqrt(PI) // DR=180.D0/PI ! from radian to degree // _______________________________________________________ // !! |Important note: all real quantities are assumed REAL*8 | !! @@ -144,6 +144,7 @@ # define led_bldata led_bldata_ # define fsiini fsiini_ # define ltran12 ltran12_ +# define boosttoprf boosttoprf_ # define fsiw fsiw_ # define setpdist setpdist_ # define type_of_call @@ -151,6 +152,7 @@ # define led_bldata LED_BLDATA # define fsiini FSIINI # define ltran12 LTRAN12 +# define boosttoprf BOOSTTOPRF # define fsiw FSIW # define setpdist SETPDIST # define type_of_call _stdcall @@ -159,17 +161,19 @@ extern "C" void type_of_call led_bldata(); extern "C" void type_of_call fsiini(); extern "C" void type_of_call ltran12(); +extern "C" void type_of_call boosttoprf(); extern "C" void type_of_call fsiw(); extern "C" void type_of_call setpdist(Double_t& r); /**************************************************************/ #include "AliHBTPair.h" -#include "AliHBTParticle.h" +#include "AliVAODParticle.h" #include "WLedCOMMONS.h" -#include #include #include #include +#include +#include ClassImp(AliHBTLLWeights) @@ -182,19 +186,53 @@ AliHBTLLWeights::AliHBTLLWeights(): fColoumbSwitch(kTRUE), fQuantStatSwitch(kTRUE), fStrongInterSwitch(kTRUE), - fColWithResidNuclSwitch(kTRUE), + fColWithResidNuclSwitch(kFALSE), fNuclMass(0.0), fNuclCharge(0.0), - fRandomPosition(kFALSE), + fRandomPosition(kNone), fRadius(0.0), + fOneMinusLambda(0.0), + fApproximationModel(0), fPID1(0), fPID2(0), fSigma(0.0) { // Default Constructor + if (fgLLWeights) + Fatal("AliHBTLLWeights","LLWeights already instatiated. Use AliHBTLLWeights::Instance()"); } /**************************************************************/ +AliHBTLLWeights::AliHBTLLWeights(const AliHBTLLWeights &/*source*/): + AliHBTWeights(), + fTest(kTRUE), + fColoumbSwitch(kTRUE), + fQuantStatSwitch(kTRUE), + fStrongInterSwitch(kTRUE), + fColWithResidNuclSwitch(kFALSE), + fNuclMass(0.0), + fNuclCharge(0.0), + fRandomPosition(kNone), + fRadius(0.0), + fOneMinusLambda(0.0), + fApproximationModel(0), + fPID1(0), + fPID2(0), + fSigma(0.0) +{ + //Copy ctor needed by the coding conventions but not used + Fatal("AliHBTLLWeights","copy ctor not implemented"); +} +/************************************************************/ + +AliHBTLLWeights& AliHBTLLWeights::operator=(const AliHBTLLWeights& /*source*/) +{ + //Assignment operator needed by the coding conventions but not used + Fatal("AliHBTLLWeights","assignment operator not implemented"); + return * this; +} +/************************************************************/ + AliHBTLLWeights* AliHBTLLWeights::Instance() { // returns instance of class @@ -208,15 +246,32 @@ AliHBTLLWeights* AliHBTLLWeights::Instance() return fgLLWeights; } } - +/************************************************************/ -Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair) +void AliHBTLLWeights::Set() +{ + //sets this as weighitng class + Info("Set","Setting Lednicky-Lyuboshitz as Weighing Class"); + + if ( fgWeights == 0x0 ) + { + fgWeights = AliHBTLLWeights::Instance(); + return; + } + if ( fgWeights == AliHBTLLWeights::Instance() ) return; + delete fgWeights; + fgWeights = AliHBTLLWeights::Instance(); +} +/************************************************************/ + +Double_t AliHBTLLWeights::GetWeight(AliHBTPair* partpair) { // calculates weight for a pair - static const Double_t cmtofm = 1.e13; + static const Double_t kcmtofm = 1.e13; + static const Double_t kcmtoOneOverGeV = kcmtofm*fgkWcons; - AliHBTParticle *part1 = partpair->Particle1(); - AliHBTParticle *part2 = partpair->Particle2(); + AliVAODParticle *part1 = partpair->Particle1(); + AliVAODParticle *part2 = partpair->Particle2(); if ( (part1 == 0x0) || (part2 == 0x0)) { @@ -224,8 +279,10 @@ Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair) return 0.0; } + if ( fPID1 != part1->GetPdgCode() ) return 1.0; + if ( fPID2 != part2->GetPdgCode() ) return 1.0; -//eats a lot of time +//takes a lot of time if ( (part1->Px() == part2->Px()) && (part1->Py() == part2->Py()) && (part1->Pz() == part2->Pz()) ) @@ -241,43 +298,79 @@ Double_t AliHBTLLWeights::GetWeight(const AliHBTPair* partpair) return 0.0; } - FSI_MOM.P1X = part1->Px(); - FSI_MOM.P1Y = part1->Py(); - FSI_MOM.P1Z = part1->Pz(); - - FSI_MOM.P2X = part2->Px(); - FSI_MOM.P2Y = part2->Py(); - FSI_MOM.P2Z = part2->Pz(); - - FSI_COOR.X1 = part1->Vx()*cmtofm; - FSI_COOR.Y1 = part1->Vy()*cmtofm; - FSI_COOR.Z1 = part1->Vz()*cmtofm; - FSI_COOR.T1 = part1->T(); - - FSI_COOR.X2 = part2->Vx()*cmtofm; - FSI_COOR.Y2 = part2->Vy()*cmtofm; - FSI_COOR.Z2 = part2->Vz()*cmtofm; - FSI_COOR.T2 = part2->T(); - - ltran12(); + if(fOneMinusLambda)//implemetation of non-zero intetcept parameter + { + if( gRandom->Rndm() < fOneMinusLambda ) return 1.0; + } + //this must be after ltran12 because it would overwrite what we set below - if (fRandomPosition) + switch (fRandomPosition) { - Double_t rxcm = fSigma*gRandom->Gaus(); - Double_t rycm = fSigma*gRandom->Gaus(); - Double_t rzcm = fSigma*gRandom->Gaus(); - - FSI_PRF.X=rxcm*fgkWcons; - FSI_PRF.Y=rycm*fgkWcons; - FSI_PRF.Z=rzcm*fgkWcons; - FSI_PRF.T=0.; - - Double_t rps=rxcm*rxcm+rycm*rycm+rzcm*rzcm; - Double_t rp=TMath::Sqrt(rps); - setpdist(rp); + case kGaussianQInv: + { + // Set Momenta in the common block + FSI_MOM.P1X = part1->Px(); + FSI_MOM.P1Y = part1->Py(); + FSI_MOM.P1Z = part1->Pz(); + + FSI_MOM.P2X = part2->Px(); + FSI_MOM.P2Y = part2->Py(); + FSI_MOM.P2Z = part2->Pz(); + + //boost it to PRF + ltran12(); + boosttoprf(); + + //Set particle positions now so they are Gaussian in PRF + RandomPairDistances(); + + FSI_PRF.RPS= FSI_COOR.X2*FSI_COOR.X2 + FSI_COOR.Y2*FSI_COOR.Y2 +FSI_COOR.Z2*FSI_COOR.Z2; + FSI_PRF.RP=TMath::Sqrt(FSI_PRF.RPS); + + break; + } + case kGaussianOSL: + { + //boost pair to LCMS and set such a momenta in the common block + Int_t retv = SetMomentaInLCMS(part1,part2); + if (retv) return 1; + //random particle positions/distance so they are Gaussian in LCMS + RandomPairDistances(); + + //Boost to PRF + ltran12(); + boosttoprf(); + + break; + } + case kNone: + default: + //set momenta and paricle positions as they are + FSI_MOM.P1X = part1->Px(); + FSI_MOM.P1Y = part1->Py(); + FSI_MOM.P1Z = part1->Pz(); + + FSI_MOM.P2X = part2->Px(); + FSI_MOM.P2Y = part2->Py(); + FSI_MOM.P2Z = part2->Pz(); + + FSI_COOR.X1 = part1->Vx()*kcmtoOneOverGeV; + FSI_COOR.Y1 = part1->Vy()*kcmtoOneOverGeV; + FSI_COOR.Z1 = part1->Vz()*kcmtoOneOverGeV; + FSI_COOR.T1 = part1->T()*kcmtoOneOverGeV; + + FSI_COOR.X2 = part2->Vx()*kcmtoOneOverGeV; + FSI_COOR.Y2 = part2->Vy()*kcmtoOneOverGeV; + FSI_COOR.Z2 = part2->Vz()*kcmtoOneOverGeV; + FSI_COOR.T2 = part2->T()*kcmtoOneOverGeV; + + ltran12(); + boosttoprf(); + + break; } - + fsiw(); return LEDWEIGHT.WEIN; } @@ -304,6 +397,8 @@ void AliHBTLLWeights::Init() Fatal("Init","Particles types are not set"); return;//pro forma } + + FSI_NS.LL = GetPairCode(fPID1,fPID2); if (FSI_NS.LL == 0) @@ -312,6 +407,8 @@ void AliHBTLLWeights::Init() return;//pro forma } + Info("Init","Setting PIDs %d %d. LL Code is %d",fPID1,fPID2,FSI_NS.LL); + TParticlePDG* tpart1 = TDatabasePDG::Instance()->GetParticle(fPID1); if (tpart1 == 0x0) @@ -578,13 +675,13 @@ void AliHBTLLWeights::SetApproxModel(Int_t ap) } /************************************************************/ -void AliHBTLLWeights::SetRandomPosition(Bool_t rp) +void AliHBTLLWeights::SetRandomPosition(ERandomizationWay rw) { - //ON=kTRUE(OFF=kFALSE) - //ON -- calculation of the Gauss source radii - //if the generator don't allows the source generation (for example MeVSim) - //if ON the following parameters are requested: - fRandomPosition = rp; +// rw can be: kGaussianQInv - so 1D Qinv correlation function has a Gaussian shape +// (and also Qlong and Qside, but not Qout) +// kGaussianOSL - so 3D Qout-Qside-Qlong correlation function has a Gaussian shape +// kNone - no randomization performed (DEFAULT) + fRandomPosition = rw; } /************************************************************/ @@ -615,3 +712,120 @@ void AliHBTLLWeights::SetNucleusMass(Double_t mass) // (see comments in fortran code) fNuclMass=mass; } +/************************************************************/ + +Int_t AliHBTLLWeights::SetMomentaInLCMS(AliVAODParticle* part1, AliVAODParticle* part2) +{ +//sets paricle momenta in the common block in LCMS frame +//---> Particle energies --------- + + Double_t p1x = part1->Px(); + Double_t p1y = part1->Py(); + Double_t p1z = part1->Pz(); + + Double_t p2x = part2->Px(); + Double_t p2y = part2->Py(); + Double_t p2z = part2->Pz(); + + Double_t am1 = part1->Mass(); + Double_t am2 = part2->Mass(); + + Double_t p1s=p1x*p1x+p1y*p1y+p1z*p1z; + Double_t p2s=p2x*p2x+p2y*p2y+p2z*p2z; + Double_t e1=TMath::Sqrt(am1*am1+p1s); + Double_t e2=TMath::Sqrt(am2*am2+p2s); +//---> pair parameters ----------- + Double_t e12=e1+e2; // energy + Double_t p12x=p1x+p2x; // px + Double_t p12y=p1y+p2y; // py + Double_t p12z=p1z+p2z; // pz + Double_t p12s=p12x*p12x+p12y*p12y+p12z*p12z; + Double_t p12 =TMath::Sqrt(p12s);// momentum + Double_t v12 =p12/e12; // velocity + + Double_t cth =p12z/p12; // cos(theta) +// Double_t sth =TMath::Sqrt(1. - cth*cth); //sin + Double_t v12z=v12*cth; // longit. v + Double_t gamz=1./TMath::Sqrt(1.-v12z*v12z); +//-- v12t=v12*sth // transv. v in cms (not needed) + + + Double_t p12ts=p12x*p12x+p12y*p12y; + Double_t p12t=TMath::Sqrt(p12ts); //pt +//===> azimuthal rotation (pt||x) ============ + if(p12t != 0.0) + { + Double_t cphi=p12x/p12t; // cos(phi) + Double_t sphi=p12y/p12t; // sin(phi) + Rotate(p1x,p1y,sphi,cphi,p1x,p1y); + Rotate(p2x,p2y,sphi,cphi,p2x,p2y); +// Rotate(x1,y1,sphi,cphi,x1,y1); +// Rotate(x2,y2,sphi,cphi,x2,y2); + } + else // rotation impossible + { + return 1; + } + +//===> co-moving ref. frame ============ + + Double_t nothing; + Boost(p1z,e1,v12z,gamz,p1z,nothing); + Boost(p2z,e2,v12z,gamz,p2z,nothing); + +// p1s=p1x*p1x+p1y*p1y+p1z*p1z; +// p2s=p2x*p2x+p2y*p2y+p2z*p2z; +// e1=TMath::Sqrt(am1*am1+p1s); +// e2=TMath::Sqrt(am2*am2+p2s); +// Boost(z1,t1,v12z,gamz,z1,t1); +// Boost(z2,t2,v12z,gamz,z2,t2); + + + FSI_MOM.P1X = p1x; + FSI_MOM.P1Y = p1y; + FSI_MOM.P1Z = p1z; + + FSI_MOM.P2X = p2x; + FSI_MOM.P2Y = p2y; + FSI_MOM.P2Z = p2z; + +// Info("Simulate"," %f %f %f ",p1x + p2x, p1y + p2y, p1z + p2z); + return 0; +} + +/**********************************************************/ + +void AliHBTLLWeights::RandomPairDistances() +{ + //randomizes pair distances + Double_t rxcm = fSigma*gRandom->Gaus(); + Double_t rycm = fSigma*gRandom->Gaus(); + Double_t rzcm = fSigma*gRandom->Gaus(); + + + FSI_COOR.X1 = 0; + FSI_COOR.Y1 = 0; + FSI_COOR.Z1 = 0; + FSI_COOR.T1 = 0; + + FSI_COOR.X2 = rxcm*fgkWcons; + FSI_COOR.Y2 = rycm*fgkWcons; + FSI_COOR.Z2 = rzcm*fgkWcons; + FSI_COOR.T2 = 0; + +} + + +void AliHBTLLWeights::Rotate(Double_t x, Double_t y, Double_t sphi, Double_t cphi, Double_t& xr, Double_t& yr) +{ +//rotates + yr=y*cphi-x*sphi; + xr=x*cphi+y*sphi; +} + +void AliHBTLLWeights::Boost(Double_t z, Double_t t, Double_t beta, Double_t gamma, Double_t& zt, Double_t& yt) +{ +//boosts + zt=gamma*(z-beta*t); + yt=gamma*(t-beta*z); +}