primary particle criteria correction (ESD)

[u/mrichter/AliRoot.git] / EVGEN / AliGenDeuteron.cxx

diff --git a/EVGEN/AliGenDeuteron.cxx b/EVGEN/AliGenDeuteron.cxx
index a3ff77cd62d49922d23c18f6d8ffc23eb2da38b4..124ee13ab002daf710e0a9fd5cc10874b38bce48 100644 (file)
--- a/EVGEN/AliGenDeuteron.cxx
+++ b/EVGEN/AliGenDeuteron.cxx
@@ -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"
@@ -51,6 +51,7 @@
 #include "AliGenCocktailEntry.h"
 #include "AliGenCocktailAfterBurner.h"
 #include "AliGenDeuteron.h"
+#include "AliLog.h"
 
 ClassImp(AliGenDeuteron)
 
@@ -58,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)
@@ -67,7 +67,6 @@ AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t c
  ,fR(0)
  ,fPsiR(0)
  ,fCurStack(0)
- ,fNtrk(0)
 {
 //
 // constructor
@@ -126,12 +125,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)
                {
@@ -144,15 +137,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
@@ -196,14 +186,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.;
@@ -287,7 +280,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;
@@ -355,11 +348,19 @@ void AliGenDeuteron::PushDeuteron(TParticle* parent1, TParticle* parent2)
        // E^2 = p^2 + m^2
        Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass);
        
+       Int_t ntrk = 0;
+       Double_t weight = 1;
+       Int_t is = 1; // final state particle
+       
        // Add a new (anti)deuteron to current event stack
        fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
                         pN.X(), pN.Y(), pN.Z(), energy,
                         vN.X(), vN.Y(), vN.Z(), parent1->T(),
-                        0., 0., 0., kPNCapture, fNtrk, 1., 0);
+                        0., 0., 0., kPNCapture, ntrk, weight, is);
+       
+       // change the status code of the parents
+       parent1->SetStatusCode(kCluster);
+       parent2->SetStatusCode(kCluster);
        
        // Set kDoneBit for the parents
        parent1->SetBit(kDoneBit);
@@ -413,3 +414,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);
+}
+