StRoot: StHbtMaker/ThCorrFctn/StHbtThPairDoubleGauss.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: 
00004  *
00005  * Author: Laurent Conin, Fabrice Retiere, Subatech, France
00006  *         Adam Kisiel, Warsaw University of Technology
00007  ***************************************************************************
00008  *
00009  * Description : implementation of StHbtThPAirGAuss
00010  *
00011  ***************************************************************************
00012  *
00013  * $Log: 
00014  *
00015  ***************************************************************************/
00016 
00017 #include "StHbtMaker/ThCorrFctn/StHbtThPairDoubleGauss.h"
00018 #include "StHbtMaker/Infrastructure/StHbtParticle.hh"
00019 #include "StHbtMaker/Base/StHbtFsiWeight.hh"
00020 
00021 #ifdef __ROOT__
00022 ClassImp(StHbtThPairDoubleGauss)
00023 #endif
00024 
00025 StHbtThPairDoubleGauss::StHbtThPairDoubleGauss() : StHbtThPair(){
00026   mEmPoint1=&mPos1;
00027   mEmPoint2=&mPos2;
00028   mUseHidMom=false;
00029   mUseHidPid=false;
00030   mMomentum1=&mMom1;
00031   mMomentum2=&mMom2;
00032   mRef=RCMSDG;
00033   mSizeX1=mSizeY1=mSizeZ1=1.;
00034   mTime1=0;
00035   mSizeX2=mSizeY2=mSizeZ2=1.;
00036   mTime2=0;
00037   mProb1=0.5;
00038   mBetaRCMS=0.;
00039   mGammaRCMS=1.;
00040   mCoreHalo=1;
00041   mXShift = 0.0;
00042   mYShift = 0.0;
00043   mZShift = 0.0;
00044   mTShift = 0.0;
00045   mRand = *(new TRandom2(31811));
00046 
00047   Int_t nbins = 100;
00048   Double_t posMin = -10.0;
00049   Double_t posMax = 10.0;
00050   Double_t ptMin = 0.0;
00051   Double_t ptMax = 1.0;
00052   
00053   mPosDist1 = new StHbt3DHisto("PosDist1","PosDist1",nbins,posMin,posMax,nbins,posMin,posMax,nbins,posMin,posMax);
00054   mPosDist2 = new StHbt3DHisto("PosDist2","PosDist2",nbins,posMin,posMax,nbins,posMin,posMax,nbins,posMin,posMax);
00055   mTDist1 = new StHbt1DHisto("TDist1","PosDist1",nbins,posMin,posMax);
00056   mTDist2 = new StHbt1DHisto("TDist2","PosDist2",nbins,posMin,posMax);
00057   mPosPtDist1 = new StHbt2DHisto("PosPtDist1","PosPtDist1",nbins,posMin,posMax,nbins,ptMin,ptMax);
00058   mPosPtDist2 = new StHbt2DHisto("PosPtDist2","PosPtDist2",nbins,posMin,posMax,nbins,ptMin,ptMax);
00059 }; 
00060 
00061 StHbtThPairDoubleGauss::~StHbtThPairDoubleGauss() 
00062 {
00063   delete mPosDist1;
00064   delete mPosDist2;
00065   delete mTDist1;
00066   delete mTDist2;
00067   delete mPosPtDist1;
00068   delete mPosPtDist2;
00069 };
00070 
00071 void StHbtThPairDoubleGauss::Set(const StHbtPair* aPair){
00072   SetMomentum_PID (aPair);
00073   SetPosition(aPair);
00074   BoostPosition();
00075 
00076   mMomParCalculated=0;
00077   mPosParCalculated=0;
00078 
00079   mMeasPair=aPair;
00080   mWeightOk=false;
00081 
00082 }
00083 
00084 void StHbtThPairDoubleGauss::SetMomentum_PID( const StHbtPair* aPair ){
00085 
00086   StHbtSmearedHiddenInfo* tSmearedHidInf1=0;
00087   StHbtSmearedHiddenInfo* tSmearedHidInf2=0;
00088   const StHbtEvtGenHiddenInfo* tEvtGenHidInf1=0;
00089   const StHbtEvtGenHiddenInfo* tEvtGenHidInf2=0;
00090   StHbtShiftedHiddenInfo* tShiftedHidInf1=0;
00091   StHbtShiftedHiddenInfo* tShiftedHidInf2=0;
00092 
00093   if (mUseHidMom||mUseHidPid) {
00094     if (mHiddenInfoType == SMEAR) {
00095       tSmearedHidInf1=dynamic_cast< StHbtSmearedHiddenInfo*>
00096         (aPair->track1()->getHiddenInfo()); 
00097       tSmearedHidInf2=dynamic_cast< StHbtSmearedHiddenInfo*>
00098         (aPair->track2()->getHiddenInfo());
00099       if ((tSmearedHidInf1==0)||(tSmearedHidInf2==0)) {
00100         if (tSmearedHidInf1==0) {
00101           StHbtMomRes *mMomRes1 = new StHbtMomRes(mPid1);
00102           mMomRes1->setMult(mResMult);
00103            StHbtLorentzVector *temp;
00104            temp = GenerateFreezeOut(1);
00105            aPair->track1()->SetHiddenInfo(new StHbtSmearedHiddenInfo(aPair->track1()->FourMomentum(), *temp, mPid1, &mRand, mMomRes1));
00106           tSmearedHidInf1 = dynamic_cast< StHbtSmearedHiddenInfo*>(aPair->track1()->getHiddenInfo());
00107           delete mMomRes1;
00108            delete temp;
00109         } 
00110         if (tSmearedHidInf2==0) {
00111           StHbtMomRes *mMomRes2 = new StHbtMomRes(mPid2);
00112           mMomRes2->setMult(mResMult);
00113            StHbtLorentzVector *temp;
00114            temp = GenerateFreezeOut(2);
00115           aPair->track2()->SetHiddenInfo(new StHbtSmearedHiddenInfo(aPair->track2()->FourMomentum(), *temp, mPid2, &mRand, mMomRes2));
00116           tSmearedHidInf2 = dynamic_cast< StHbtSmearedHiddenInfo*>(aPair->track2()->getHiddenInfo());
00117           delete mMomRes2;
00118            delete temp;
00119         } 
00120       }
00121     }
00122     else if (mHiddenInfoType == SHIFT) {
00123       tShiftedHidInf1=dynamic_cast< StHbtShiftedHiddenInfo*>
00124         (aPair->track1()->getHiddenInfo()); 
00125       tShiftedHidInf2=dynamic_cast< StHbtShiftedHiddenInfo*>
00126         (aPair->track2()->getHiddenInfo());
00127       if ((tShiftedHidInf1==0)||(tShiftedHidInf2==0)) {
00128         if (tShiftedHidInf1==0) {
00129           StHbtMomRes *mMomRes1 = new StHbtMomRes(mPid1);
00130           mMomRes1->setMult(mResMult);
00131           aPair->track1()->SetHiddenInfo(new StHbtShiftedHiddenInfo(aPair->track1()->FourMomentum(), mPid1, &mRand, mMomRes1, mShift, PHISHIFT));
00132           tShiftedHidInf1 = dynamic_cast< StHbtShiftedHiddenInfo*>(aPair->track1()->getHiddenInfo());
00133           delete mMomRes1;
00134         } 
00135         if (tShiftedHidInf2==0) {
00136           StHbtMomRes *mMomRes2 = new StHbtMomRes(mPid2);
00137           mMomRes2->setMult(mResMult);
00138           aPair->track2()->SetHiddenInfo(new StHbtShiftedHiddenInfo(aPair->track2()->FourMomentum(), mPid2, &mRand, mMomRes2, mShift, PHISHIFT));
00139           tShiftedHidInf2 = dynamic_cast< StHbtShiftedHiddenInfo*>(aPair->track2()->getHiddenInfo());
00140           delete mMomRes2;
00141         } 
00142       }
00143     }
00144     else {
00145       tEvtGenHidInf1=dynamic_cast<const StHbtEvtGenHiddenInfo*>
00146         (aPair->track1()->HiddenInfo()); 
00147       tEvtGenHidInf2=dynamic_cast<const StHbtEvtGenHiddenInfo*>
00148         (aPair->track2()->HiddenInfo());
00149     }
00150   }
00151   
00152   if (mUseHidMom&&tSmearedHidInf2&&tSmearedHidInf1) {
00153     //    cout << "Setting Smeared HiddenMomentum " << endl;
00154     mMomentum1=&(tSmearedHidInf1->mSmearedMom);
00155     mMomentum2=&(tSmearedHidInf2->mSmearedMom);
00156   }
00157   else if (mUseHidMom&&tShiftedHidInf2&&tShiftedHidInf1) {
00158     mMomentum1=&(tShiftedHidInf1->mShiftedMom);
00159     mMomentum2=&(tShiftedHidInf2->mShiftedMom);
00160   }
00161   else if (mUseHidMom&&tEvtGenHidInf2&&tEvtGenHidInf1) {
00162     //    mMomentum1=&(tEvtGenHidInf1->getFreezeOutMomEn());
00163     //    mMomentum2=&(tEvtGenHidInf2->getFreezeOutMomEn());
00164     mMomentum1= new StHbtLorentzVector(*tEvtGenHidInf1->getFreezeOutMomEn());
00165     mMomentum2= new StHbtLorentzVector(*tEvtGenHidInf2->getFreezeOutMomEn());
00166   }
00167   else{
00168     mMom1=aPair->track1()->FourMomentum();
00169     mMom1.setE(::sqrt(mMassSq1+mMom1.vect().mag2()));
00170     mMom2=aPair->track2()->FourMomentum();
00171     mMom2.setE(::sqrt(mMassSq2+mMom2.vect().mag2()));
00172     if ((!mUseHidMom)&&(mBetaRCMS>0)) {
00173       double tE=mMom1.e();
00174       mMom1.setE(mGammaRCMS*(tE-mBetaRCMS*mMom1.pz()));
00175       mMom1.setPz(mGammaRCMS*(mMom1.pz()-mBetaRCMS*tE));
00176       tE=mMom2.e();
00177       mMom2.setE(mGammaRCMS*(tE-mBetaRCMS*mMom2.pz()));
00178       mMom2.setPz(mGammaRCMS*(mMom2.pz()-mBetaRCMS*tE));                 
00179     }
00180   }
00181 }
00182     
00183 void StHbtThPairDoubleGauss::SetPosition( const StHbtPair* aPair) {
00184   /*
00185    * Set the random freeze-out space coordinates
00186    * according to either of two models:
00187    * Core-Halo - the coordinates are randomly generated
00188    * either with core (1st) or halo (2nd) parametrization
00189    * TwoSources - The first particle coordinates are generated 
00190    * from the 1st set of parameters, coordinates for te second - 
00191    * from the second set.
00192    */
00193 
00194   const StHbtSmearedHiddenInfo* tSmearedHidInf1=0;
00195   const StHbtSmearedHiddenInfo* tSmearedHidInf2=0;
00196 
00197   if (mHiddenInfoType == SMEAR) {
00198     tSmearedHidInf1 = dynamic_cast<const StHbtSmearedHiddenInfo*>(aPair->track1()->HiddenInfo());
00199     tSmearedHidInf2 = dynamic_cast<const StHbtSmearedHiddenInfo*>(aPair->track2()->HiddenInfo());
00200   }
00201 
00202   if ((mHiddenInfoType == SMEAR) && (tSmearedHidInf1 != 0) && (tSmearedHidInf2 != 0))
00203     {
00204       //      cout << "Setting position for pair:  " << this << endl;
00205       mPos1 = tSmearedHidInf1->getFreezeOut();
00206       mPos2 = tSmearedHidInf2->getFreezeOut();
00207       //      cout << "Set the position 1 to: X: " << mPos1.x() << "   Y: " << mPos1.y() << "   Z: " << mPos1.z() << "   T: " << mPos1.t() << endl;
00208       //      cout << "Set the position 2 to: X: " << mPos2.x() << "   Y: " << mPos2.y() << "   Z: " << mPos2.z() << "   T: " << mPos2.t() << endl;
00209 
00210     }
00211   else
00212     {
00213        StHbtLorentzVector* temp;
00214        temp = GenerateFreezeOut(1);
00215        mPos1 = *temp;
00216        delete temp;
00217        temp = GenerateFreezeOut(2);
00218        mPos2 = *temp;
00219        delete temp;
00220       //      cout << mPos1.x() << " " << mPos1.y() << endl;
00221     }
00222 
00223 };
00224 
00225 StHbtLorentzVector* StHbtThPairDoubleGauss::GenerateFreezeOut(int partno) {
00226 
00227   StHbtLorentzVector* mPos;
00228   mPos = new StHbtLorentzVector();
00229   float x,y,z,t;
00230 
00231   //  cout << "Generating Freeze-out! " << partno  << endl;
00232   mRand.Rannor(x,y);
00233   mRand.Rannor(z,t);
00234   if (mCoreHalo == 1) {
00235      if (mRand.Uniform() < mProb1) {
00236        mPos->setX(x*mSizeX1);
00237        mPos->setY(y*mSizeY1);
00238        mPos->setZ(z*mSizeZ1);
00239        mPos->setT(t*mTime1);
00240      }
00241      else {
00242        mPos->setX(x*mSizeX2);
00243        mPos->setY(y*mSizeY2);
00244        mPos->setZ(z*mSizeZ2);
00245        mPos->setT(t*mTime2);
00246      }
00247   }
00248   else if (mCoreHalo ==2) 
00249     {
00250       Double_t phi;
00251       phi = mRand.Rndm();
00252       
00253       mPos->setX(x*mSizeX1*cos(phi*2*TMath::Pi()));
00254       mPos->setY(x*mSizeX1*sin(phi*2*TMath::Pi()));
00255       mPos->setZ(z*mSizeZ2);
00256       mPos->setT(t*mTime2);      
00257     }
00258   else if (mCoreHalo ==3) 
00259     {
00260       if (partno == 1)
00261         {
00262           mPos->setX(0.0);
00263           mPos->setY(0.0);
00264           mPos->setZ(0.0);
00265           mPos->setT(0.0);      
00266         }
00267       else
00268         {
00269           double tROut =  mRand.Gaus(mTime1,mSizeX1); // reuse of long
00270           double tRSide = mRand.Gaus(0,mSizeX1);
00271           double tRLong = mRand.Gaus(0,mSizeX1);
00272           double tDTime = tROut;
00273       
00274           double tPx = mMomentum1->x()+mMomentum2->x();
00275           double tPy = mMomentum1->y()+mMomentum2->y();
00276           double tPz = mMomentum1->z()+mMomentum2->z();
00277           double tE  = mMomentum1->e()+mMomentum2->e();
00278 
00279           double tPt = tPx*tPx+tPy*tPy;
00280           double tMt = tE*tE-tPz*tPz; 
00281           double tM = ::sqrt(tMt - tPt);
00282           tPt = ::sqrt(tPt);
00283           tMt = ::sqrt(tMt);
00284 
00285           tROut *= (tMt/tM); // Rout*gammaT
00286           tDTime *= (tPt/tM); // Rout*betaT*gammaT
00287           double ttDTime = tDTime; 
00288           tDTime += (tPz/tE*tRLong);
00289           tDTime *= (tE/tMt);
00290           tRLong += (tPz/tE*ttDTime); 
00291           tRLong *= (tE/tMt);
00292           
00293           tPx /= tPt;
00294           tPy /= tPt;
00295           
00296           mPos->setX(tROut*tPx-tRSide*tPy);
00297           mPos->setY(tROut*tPy+tRSide*tPx);
00298           mPos->setZ(tRLong);
00299           mPos->setT(tDTime);
00300 
00301           //       mPos->setX(x*mSizeX1*cos(phi*2*TMath::Pi()));
00302           //       mPos->setY(x*mSizeX1*sin(phi*2*TMath::Pi()));
00303           //       mPos->setZ(z*mSizeZ2);
00304           //       mPos->setT(t*mTime2);      
00305         }
00306     }
00307   else
00308   {
00309     if (partno == 1)
00310       {
00311         mPos->setX(x*mSizeX1);
00312         mPos->setY(y*mSizeY1);
00313         mPos->setZ(z*mSizeZ1);
00314         mPos->setT(t*mTime1);
00315         //      mPos->setT(mTime1);
00316       }
00317     else
00318       {
00319         mPos->setX(mXShift+x*mSizeX2);
00320         mPos->setY(mYShift+y*mSizeY2);
00321         mPos->setZ(mZShift+z*mSizeZ2);
00322         mPos->setT(mTShift+t*mTime2);
00323         //      mPos->setT(mTime2);
00324       }
00325   }
00326   if (partno == 1)
00327     {
00328       mTDist1->Fill(mPos->t());
00329       mPosDist1->Fill(mPos->x(),mPos->y(),mPos->z());
00330       //      mPosPtDist1->Fill(hypot(mPos->x(),mPos->y()), hypot(partMom->x(),partMom->y()));
00331     }
00332   else
00333     {
00334       mTDist2->Fill(mPos->t());
00335       mPosDist2->Fill(mPos->x(),mPos->y(),mPos->z());
00336       //      mPosPtDist2->Fill(hypot(mPos->x(),mPos->y()), hypot(partMom->y(),partMom->x()));
00337     }
00338   return mPos;
00339 }
00340 
00341 void StHbtThPairDoubleGauss::BoostPosition(){
00342   double tBeta;
00343   double tGamma;
00344   double tT;
00345   switch (mRef) {
00346   case RCMSDG: break;
00347   case LCMSDG:
00348     tBeta=(mMomentum1->pz()+mMomentum2->pz())/(mMomentum1->e()+mMomentum2->e());
00349      tGamma=::sqrt(1/1-tBeta*tBeta);
00350      tT=mPos1.t();
00351     mPos1.setT(tGamma*(tT-tBeta*mPos1.z()));
00352     mPos1.setZ(tGamma*(tBeta*mPos1.z()-tBeta*tT));
00353     tT=mPos2.t();
00354     mPos2.setT(tGamma*(tT-tBeta*mPos2.z()));
00355     mPos2.setZ(tGamma*(mPos2.z()-tBeta*tT)); 
00356     break;
00357   case PRFDG:
00358     StHbtLorentzVector tBoost=*mMomentum1+*mMomentum2;
00359     mPos1=mPos1.boost(tBoost);
00360     mPos2=mPos2.boost(tBoost);  
00361     break;
00362   }
00363 };
00364 
00365 inline void StHbtThPairDoubleGauss::SetResolutionMult(const double mult) {
00366   mResMult = mult;
00367 }
00368 
00369 inline void StHbtThPairDoubleGauss::SetMomentumShift(const double shift)
00370 {
00371   mShift = shift;
00372 }
00373 
00374 inline void StHbtThPairDoubleGauss::UseSmearedHiddenInfo()
00375 {
00376   mHiddenInfoType = SMEAR;
00377 }
00378 
00379 inline void StHbtThPairDoubleGauss::UseShiftedHiddenInfo()
00380 {
00381   mHiddenInfoType = SHIFT;
00382 }
00383 
00384 inline void StHbtThPairDoubleGauss::UseEvtGenHiddenInfo()
00385 {
00386   mHiddenInfoType = EVTGEN;
00387 }
00388 
00389 inline  void          StHbtThPairDoubleGauss::SetSizes(double aXYZ1, double aT1, double aXYZ2, double aT2) 
00390 {mSizeX1=mSizeY1=mSizeZ1=aXYZ1; mTime1=aT1; mSizeX2=mSizeY2=mSizeZ2=aXYZ2; mTime2=aT2; }
00391 inline  void          StHbtThPairDoubleGauss::SetSizes(double aX1,double aY1,double aZ1, double aT1,double aX2,double aY2,double aZ2, double aT2)
00392 {mSizeX1=aX1; mSizeY1=aY1; mSizeZ1=aZ1; mTime1=aT1; mSizeX2=aX2; mSizeY2=aY2; mSizeZ2=aZ2; mTime2=aT2; }
00393 inline  void          StHbtThPairDoubleGauss::SetSize1(double aSize,double aTime) {mSizeX1=aSize;mSizeY1=aSize;mSizeZ1=aSize;mTime1=aTime;};
00394 inline  void          StHbtThPairDoubleGauss::SetSize1(double aSizeX,double aSizeY, double aSizeZ,double aTime)
00395 {mSizeX1=aSizeX;mSizeY1=aSizeY;mSizeZ1=aSizeZ;mTime1=aTime;};
00396 inline  void          StHbtThPairDoubleGauss::SetSize2(double aSize,double aTime) {mSizeX2=aSize;mSizeY2=aSize;mSizeZ2=aSize;mTime2=aTime;};
00397 inline  void          StHbtThPairDoubleGauss::SetSize2(double aSizeX,double aSizeY, double aSizeZ,double aTime)
00398 {mSizeX2=aSizeX;mSizeY2=aSizeY;mSizeZ2=aSizeZ;mTime2=aTime;};
00399 
00400 inline  void          StHbtThPairDoubleGauss::UseHiddenMomentum(){mUseHidMom=1;};
00401 inline  void          StHbtThPairDoubleGauss::UseParticleMomentum(){
00402   mUseHidMom=0;
00403   if (mUseHidPid){
00404     mMomentum1=&mMom1;
00405     mMomentum2=&mMom2;
00406   }
00407 };
00408 
00409 inline  void          StHbtThPairDoubleGauss::UseHiddenPid() {
00410   mUseHidPid=true;
00411   if (mUseHidMom){
00412     mMomentum1=&mMom1;
00413     mMomentum2=&mMom2;
00414   }
00415 };
00416 
00417 inline  void          StHbtThPairDoubleGauss::UseFixedPid( int const tPid1,double const tMass1,int const tPid2, double const tMass2) {
00418   mUseHidPid=false;mPid1=tPid1;mPid2=tPid2;mMassSq1=tMass1*tMass1;mMassSq2=tMass2*tMass2; }; 
00419 inline  void          StHbtThPairDoubleGauss::UseFixedPid( int const tPid1,double const tMass1) {
00420   mUseHidPid=false;mPid1=tPid1;mPid2=tPid1;mMassSq1=tMass1*tMass1;mMassSq2=tMass1*tMass1; }; 
00421 
00422 inline  void          StHbtThPairDoubleGauss::SetRCMS() {mRef=RCMSDG;};
00423 inline  void          StHbtThPairDoubleGauss::SetLCMS(){mRef=LCMSDG;};
00424 inline  void          StHbtThPairDoubleGauss::SetPRF(){mRef=PRFDG;};
00425 
00426 inline  void          StHbtThPairDoubleGauss::SetBoostRCMS(double aPlab,double aMBeam, double aMTarget){
00427   double tEBeamLab=::sqrt(aPlab*aPlab+aMBeam*aMBeam);
00428   mGammaRCMS=(tEBeamLab+aMTarget)/::sqrt(aMBeam*aMBeam+aMTarget*aMTarget+2*tEBeamLab*aMTarget);
00429   mBetaRCMS=::sqrt(1.-1/(mGammaRCMS*mGammaRCMS));
00430 }
00431 
00432 inline  void          StHbtThPairDoubleGauss::SetFirstProb(double aProb1) { 
00433   if (aProb1<0.0) mProb1 = 0.0;
00434   else if (aProb1>1.0) mProb1 = 1.0;
00435   else mProb1 = aProb1;
00436 }
00437 
00438 inline  void          StHbtThPairDoubleGauss::SetCoreHalo()   { mCoreHalo=1; }
00439 inline  void          StHbtThPairDoubleGauss::SetTwoSources() { mCoreHalo=0; }
00440 inline  void          StHbtThPairDoubleGauss::SetRadialGaus() { mCoreHalo=2; }
00441 inline  void          StHbtThPairDoubleGauss::SetPRFGaus()    { mCoreHalo=3; }
00442 inline  void          StHbtThPairDoubleGauss::SetPositionShift(double aX, double aY, double aZ, double aT)
00443 {
00444   mXShift = aX;
00445   mYShift = aY;
00446   mZShift = aZ;
00447   mTShift = aT;
00448 }
00449 
00450 void StHbtThPairDoubleGauss::Write()
00451 {
00452   mPosDist1->Write();
00453   mPosDist2->Write();
00454   mTDist1->Write();
00455   mTDist2->Write();
00456   mPosPtDist1->Write();
00457   mPosPtDist2->Write();
00458 
00459 }