PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelGausRinvFreezeOutGenerator.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 ///                                                                          ///
   3 /// AliFemtoModelGausRinvFreezeOutGenerator - freeze-out                     ///
   4 /// coordinates generator, generating a 3D gaussian ellipsoid in LCMS        ///
   5 /// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
   6 ///                                                                          ///
   7 ////////////////////////////////////////////////////////////////////////////////
   8 #ifdef __ROOT__
   9   ClassImp(AliFemtoModelGausRinvFreezeOutGenerator, 1)
  10 #endif
  11
  12 #include "math.h"
  13 #include "AliFemtoModelGausRinvFreezeOutGenerator.h"
  14 #include "AliFemtoModelHiddenInfo.h"
  15 #include "AliFemtoModelGlobalHiddenInfo.h"
  16 #include "AliFemtoLorentzVector.h"
  17
  18 //_______________________
  19 AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator() :
  20   fSizeInv(0),
  21   fSelectPrimary(false)
  22 {
  23   // Default constructor
  24   fRandom = new TRandom2();
  25 }
  26
  27 //_______________________
  28 AliFemtoModelGausRinvFreezeOutGenerator::AliFemtoModelGausRinvFreezeOutGenerator(const AliFemtoModelGausRinvFreezeOutGenerator &aModel):
  29   AliFemtoModelFreezeOutGenerator(),
  30   fSizeInv(0),
  31   fSelectPrimary(false)
  32 {
  33   // Copy constructor
  34   fRandom = new TRandom2();
  35   SetSizeInv(aModel.GetSizeInv());
  36 }
  37 //_______________________
  38 AliFemtoModelGausRinvFreezeOutGenerator::~AliFemtoModelGausRinvFreezeOutGenerator()
  39 {
  40   if (fRandom) delete fRandom;
  41 }
  42 //_______________________
  43 void AliFemtoModelGausRinvFreezeOutGenerator::GenerateFreezeOut(AliFemtoPair *aPair)
  44 {
  45   AliFemtoModelHiddenInfo *inf1 = (AliFemtoModelHiddenInfo *) aPair->Track1()->HiddenInfo();
  46   AliFemtoModelHiddenInfo *inf2 = (AliFemtoModelHiddenInfo *) aPair->Track2()->HiddenInfo();
  47
  48   if ((!inf1) || (!inf2)) { cout << "Hidden info not created! "  << endl; exit(kFALSE); }
  49
  50   if (fSelectPrimary) {
  51     const AliFemtoModelGlobalHiddenInfo *infg1 = dynamic_cast<const AliFemtoModelGlobalHiddenInfo *>(aPair->Track1()->HiddenInfo());
  52     const AliFemtoModelGlobalHiddenInfo *infg2 = dynamic_cast<const AliFemtoModelGlobalHiddenInfo *>(aPair->Track2()->HiddenInfo());
  53
  54     if ((infg1) && (infg2)) {
  55       // assume the emission point is in [cm] and try to judge if
  56       // both particles are primary
  57       Double_t dist1 = infg1->GetGlobalEmissionPoint()->perp();
  58       Double_t dist2 = infg2->GetGlobalEmissionPoint()->perp();
  59
  60       if ((dist1 > 0.05) && (dist2 > 0.05)) {
  61         // At least one particle is non primary
  62         if (!(inf1->GetEmissionPoint())) {
  63           AliFemtoLorentzVector tPos(-1000,1000,-500,0);
  64           inf1->SetEmissionPoint(&tPos);
  65         }
  66         else
  67           inf1->SetEmissionPoint(-1000,1000,-500,0);
  68         if (!(inf2->GetEmissionPoint())) {
  69           AliFemtoLorentzVector tPos(fRandom->Gaus(0,1000.0),fRandom->Gaus(0,1000),fRandom->Gaus(0,1000),0.0);
  70           inf2->SetEmissionPoint(&tPos);
  71         }
  72         else
  73           inf2->SetEmissionPoint(fRandom->Gaus(0,1000), fRandom->Gaus(0,1000), fRandom->Gaus(0,1000), 0.0);
  74
  75         return;
  76       }
  77     }
  78   }
  79
  80   // Generate two particle emission points with respect
  81   // to their pair momentum
  82   // The source is the 3D Gaussian ellipsoid in the LCMS frame
  83
  84   // Calculate sum momenta
  85   Double_t tPx = inf1->GetTrueMomentum()->x() + inf2->GetTrueMomentum()->x();
  86   Double_t tPy = inf1->GetTrueMomentum()->y() + inf2->GetTrueMomentum()->y();
  87   Double_t tPz = inf1->GetTrueMomentum()->z() + inf2->GetTrueMomentum()->z();
  88   Double_t tM1 = inf1->GetMass();
  89   Double_t tM2 = inf2->GetMass();
  90   Double_t tE1 = sqrt(tM1*tM1 + inf1->GetTrueMomentum()->mag2());
  91   Double_t tE2 = sqrt(tM2*tM2 + inf2->GetTrueMomentum()->mag2());
  92   Double_t tEs = tE1 + tE2;
  93
  94   Double_t tPt = sqrt(tPx*tPx + tPy*tPy);
  95   Double_t tMt = sqrt(tEs*tEs - tPz*tPz);
  96
  97   // Generate positions in PRF from a Gaussian
  98   Double_t tROutS =  fRandom->Gaus(0,fSizeInv); // reuse of long
  99   Double_t tRSideS = fRandom->Gaus(0,fSizeInv);
 100   Double_t tRLongS = fRandom->Gaus(0,fSizeInv);
 101   Double_t tRTimeS = 0;
 102
 103   Double_t tBetat  = tPt/tMt;
 104   Double_t tGammat = 1.0/sqrt(1.0-tBetat*tBetat);
 105
 106   Double_t tBetaz  = tPz/tEs;
 107   Double_t tGammaz = 1.0/sqrt(1.0-tBetaz*tBetaz);
 108
 109   Double_t tROut = tGammat * (tROutS + tBetat * tRTimeS);
 110   Double_t tDtL  = tGammat * (tRTimeS + tBetat * tROutS);
 111   Double_t tRSide = tRSideS;
 112
 113   Double_t tRLong = tGammaz * (tRLongS + tBetaz * tDtL);
 114   Double_t tDt    = tGammaz * (tDtL + tBetaz * tRLongS);
 115
 116   tPx /= tPt;
 117   tPy /= tPt;
 118
 119   Double_t tXout  = tROut*tPx-tRSide*tPy;
 120   Double_t tXside = tROut*tPy+tRSide*tPx;
 121   Double_t tXlong = tRLong;
 122   Double_t tXtime = tDt;
 123
 124   if (!(inf1->GetEmissionPoint())) {
 125     AliFemtoLorentzVector tPos(0,0,0,0);
 126     inf1->SetEmissionPoint(&tPos);
 127   }
 128   else
 129     inf1->SetEmissionPoint(0,0,0,0);
 130   if (!(inf2->GetEmissionPoint())) {
 131     AliFemtoLorentzVector tPos(tXout,tXside,tXlong,tXtime);
 132     inf2->SetEmissionPoint(&tPos);
 133   }
 134   else
 135     inf2->SetEmissionPoint(tXout, tXside, tXlong, tXtime);
 136 }
 137
 138 //_______________________
 139 void AliFemtoModelGausRinvFreezeOutGenerator::SetSizeInv(Double_t aSizeInv)
 140 {
 141   fSizeInv = aSizeInv;
 142 }
 143 //_______________________
 144 Double_t AliFemtoModelGausRinvFreezeOutGenerator::GetSizeInv() const
 145 {
 146   return fSizeInv;
 147 }
 148 //_______________________
 149 AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::Clone() const
 150 {
 151   return GetGenerator();
 152 }
 153 //_______________________
 154 inline AliFemtoModelFreezeOutGenerator* AliFemtoModelGausRinvFreezeOutGenerator::GetGenerator() const
 155 {
 156   AliFemtoModelFreezeOutGenerator* tModel = new AliFemtoModelGausRinvFreezeOutGenerator(*this);
 157   return tModel;
 158 }
 159 //_______________________
 160 void AliFemtoModelGausRinvFreezeOutGenerator::SetSelectPrimaryFromHidden(bool aUse)
 161 {
 162   fSelectPrimary = aUse;
 163 }
 164 Bool_t AliFemtoModelGausRinvFreezeOutGenerator::GetSelectPrimaryFromHidden()
 165 {
 166   return fSelectPrimary;
 167 }
 168