relative momentum between the p-n pairs should be computed in their CM frame
authormorsch <andreas.morsch@cern.ch>
Sat, 5 Apr 2014 20:07:08 +0000 (22:07 +0200)
committermorsch <andreas.morsch@cern.ch>
Sat, 5 Apr 2014 20:07:08 +0000 (22:07 +0200)
Eulogio Serradilla <eulogio.serradilla@cern.ch>

EVGEN/AliGenDeuteron.cxx
EVGEN/AliGenDeuteron.h

index 6bccb91..03b9794 100644 (file)
@@ -38,8 +38,8 @@
 
 #include "Riostream.h"
 #include "TParticle.h"
-#include "AliRun.h"
 #include "AliStack.h"
+#include "AliRun.h"
 #include "TMath.h"
 #include "TMCProcess.h"
 #include "TList.h"
@@ -59,8 +59,7 @@ AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t c
  :fSign(1)
  ,fPmax(pmax)
  ,fRmax(rmax)
- ,fSpinProb(0.75)
- ,fMaxRadius(1000.)
+ ,fSpinProb(1)
  ,fClusterType(cluster)
  ,fModel(0)
  ,fTimeLength(2.5)
@@ -127,12 +126,6 @@ void AliGenDeuteron::Generate()
                        if(fB >= 2.*fR) continue; // no collision
                }
                
-               // primary vertex
-               TArrayF primVtx;
-               gener->GetActiveEventHeader()->PrimaryVertex(primVtx);
-               TVector3 r0(primVtx[0],primVtx[1],primVtx[2]);
-               
-               // emerging nucleons from the collision
                fCurStack = gener->GetStack(ns);
                if(!fCurStack)
                {
@@ -145,15 +138,12 @@ void AliGenDeuteron::Generate()
                TList* neutrons = new TList();
                neutrons->SetOwner(kFALSE);
                
-               for (Int_t i=0; i < fCurStack->GetNtrack(); ++i)
+               // particles produced by the generator
+               for (Int_t i=0; i < fCurStack->GetNprimary(); ++i)
                {
                        TParticle* iParticle = fCurStack->Particle(i);
                        
-                       if(iParticle->TestBit(kDoneBit)) continue;
-                       
-                       // select only nucleons within the freeze-out volume
-                       TVector3 r(iParticle->Vx(),iParticle->Vy(),iParticle->Vz());
-                       if((r-r0).Mag() > fMaxRadius*1.e-13) continue;
+                       if(iParticle->GetStatusCode() != 1) continue;
                        
                        Int_t pdgCode = iParticle->GetPdgCode();
                        if(pdgCode == fSign*2212)// (anti)proton
@@ -197,14 +187,17 @@ Double_t AliGenDeuteron::GetCoalescenceProbability(const TParticle* nucleon1, co
 // returns < 0 if coalescence is not possible
 // otherwise returns a coalescence probability
 //
+       const Double_t kProtonMass  = 0.938272013;
+       const Double_t kNeutronMass = 0.939565378;
+       
        TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz());
        TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz());
        
        TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz());
        TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz());
        
-       Double_t deltaP = (p2-p1).Mag();       // relative momentum
-       if( deltaP >= fPmax)       return -1.;
+       Double_t deltaP = 2.*this->GetPcm(p1, kProtonMass, p2, kNeutronMass); // relative momentum in CM frame
+       if( deltaP >= fPmax) return -1.;
        
        Double_t deltaR = (v2-v1).Mag();       // relative distance (cm)
        if(deltaR >= fRmax*1.e-13) return -1.;
@@ -288,7 +281,7 @@ void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons)
        Double_t* iWeight = new Double_t[nMaxIntPair];
        
        Int_t iIdx = -1;
-       while(true)
+       while(kTRUE)
        {
                Int_t j = -1;
                Double_t wMax = 0;
@@ -414,3 +407,32 @@ void AliGenDeuteron::FixProductionVertex(TParticle* i)
        i->SetProductionVertex(i->Vx() + 1.e-13*x, i->Vy() + 1.e-13*y, i->Vz() + 1.e-13*z, i->T());
 }
 
+Double_t AliGenDeuteron::GetS(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const
+{
+//
+// square of the center of mass energy
+//
+       Double_t E1 = TMath::Sqrt( p1x*p1x + p1y*p1y + p1z*p1z + m1*m1);
+       Double_t E2 = TMath::Sqrt( p2x*p2x + p2y*p2y + p2z*p2z + m2*m2);
+       
+       return (E1+E2)*(E1+E2) - ((p1x+p2x)*(p1x+p2x) + (p1y+p2y)*(p1y+p2y) + (p1z+p2z)*(p1z+p2z));
+}
+
+Double_t AliGenDeuteron::GetPcm(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const
+{
+//
+// momentum in the CM frame
+//
+       Double_t s = this->GetS(p1x, p1y, p1z, m1, p2x, p2y, p2z, m2);
+       
+       return TMath::Sqrt( (s-(m1-m2)*(m1-m2))*(s-(m1+m2)*(m1+m2)) )/(2.*TMath::Sqrt(s));
+}
+
+Double_t AliGenDeuteron::GetPcm(const TVector3& p1, Double_t m1, const TVector3& p2, Double_t m2) const
+{
+//
+// momentum in the CM frame
+//
+       return this->GetPcm(p1.Px(),p1.Py(),p1.Pz(),m1,p2.Px(),p2.Py(),p2.Pz(),m2);
+}
+
index 9d974bf..ff2f21a 100644 (file)
@@ -17,7 +17,7 @@ class AliGenDeuteron: public AliGenerator
 
  public:
 
-       AliGenDeuteron(Int_t sign=1, Double_t pmax=0.3, Double_t rmax=2.1, Int_t cluster=0 );
+       AliGenDeuteron(Int_t sign=1, Double_t pmax=0.2, Double_t rmax=2.1, Int_t cluster=0 );
        virtual ~AliGenDeuteron();
 
        virtual void Init();
@@ -27,7 +27,6 @@ class AliGenDeuteron: public AliGenerator
        Double_t GetCoalescenceMomentum() const { return fPmax; }
        Double_t GetCoalescenceDistance() const { return fRmax; }
        Double_t GetSpinProbability() const { return fSpinProb; }
-       Double_t GetFreezeOutRadius() const { return fMaxRadius; }
        Double_t GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const;
        Int_t GetClusterType() const { return fClusterType; }
        Int_t GetFreezeOutModel() const { return fModel; }
@@ -37,14 +36,13 @@ class AliGenDeuteron: public AliGenerator
        void SetCoalescenceMomentum(Double_t p) { fPmax = p; }
        void SetCoalescenceDistance(Double_t r) { fRmax = r; }
        void SetSpinProbability(Double_t s) { fSpinProb = s; }
-       void SetFreezeOutRadius(Double_t r) { fMaxRadius = r; }
        void SetClusterType(Int_t cluster) { fClusterType = cluster; }
        void SetFreezeOutModel(Int_t model, Double_t timeLength=2.5) { fModel = model; fTimeLength=timeLength;}
        
  public:
 
-       enum { kFirstPartner, kLowestMomentum, kLowestDistance, kBoth };
-       enum { kNone, kThermal, kExpansion };
+       enum { kFirstPartner=0, kLowestMomentum, kLowestDistance, kBoth };
+       enum { kNone=0, kThermal, kExpansion };
        
  private:
  
@@ -56,13 +54,16 @@ class AliGenDeuteron: public AliGenerator
        void WeightMatrix(const TList* protons, const TList* neutrons);
        void PushDeuteron(TParticle* parent1, TParticle* parent2);
        
+       Double_t GetS(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const;
+       Double_t GetPcm(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const;
+       Double_t GetPcm(const TVector3& p1, Double_t m1, const TVector3& p2, Double_t m2) const;
+       
  private:
        
        Int_t fSign;          // +1 for deuterons, -1 for antideuterons
        Double_t fPmax;       // Maximum p-n momentum difference (GeV/c)
        Double_t fRmax;       // Maximum p-n distance (fm)
        Double_t fSpinProb;   // cluster formation probability due to spin
-       Double_t fMaxRadius;  // Maximum freeze-out radius (fm)
        Int_t fClusterType;   // Probability criteria to find clusters
        Int_t fModel;         // Model to override generator source
        Double_t fTimeLength; // Thermal and chemical freeze-out time (fm/c)
@@ -72,7 +73,7 @@ class AliGenDeuteron: public AliGenerator
        AliStack* fCurStack;  //! current event stack
        Int_t fNtrk;          //! number of the stored track
        
-       ClassDef(AliGenDeuteron,1)
+       ClassDef(AliGenDeuteron,2)
 };
 
 #endif // ALIGENDEUTERON_H