Changes for #91126: Request to commit a patch for AMPT

diff --git a/TAmpt/AMPT/amptsetdef.f b/TAmpt/AMPT/amptsetdef.f
index 75a5c33bd09ba1c48ec4cc4ecae1b2fe3747f91d..7313444227cda14e4c745adcbf2223a5a384e6f2 100644 (file)
--- a/TAmpt/AMPT/amptsetdef.f
+++ b/TAmpt/AMPT/amptsetdef.f
@@ -68,7 +68,7 @@ c     quark coalescence radii in momentum and space for string melting:
       dpcoal   = 1d6           ! dpcoal in GeV
       drcoal   = 1d6           ! drcoal in fm
 c     flag: read in HIJING random # seed at runtime(1) or from input.ampt(D=0):
-      ihjsed   = 0             ! ihjsed: take HIJING seed from below (D=0)or at runtime(11)
+      ihjsed   = 11            ! ihjsed: take HIJING seed from below (D=0)or at runtime(11)
 c     2 seeds for random number generators in HIJING/hadron cascade and ZPC:
       nseed    = 53153511      ! random seed for HIJING
       iseedp   = 8             ! random seed for parton cascade

diff --git a/TAmpt/AMPT/hijing1.383_ampt.f b/TAmpt/AMPT/hijing1.383_ampt.f
index f7dd2ff1be46da2d711c2f46f6eab4599d2e295f..55c7bbe061e4247cd0b5e062227c8f2b309fcca4 100644 (file)
--- a/TAmpt/AMPT/hijing1.383_ampt.f
+++ b/TAmpt/AMPT/hijing1.383_ampt.f
@@ -507,8 +507,8 @@ C       BY THE ANGLE PHI  FOR EACH COLLISION.******************
 C
         BB=SQRT(BMIN**2+RANART(NSEED)*(BMAX**2-BMIN**2))
 cbz6/28/99 flow1
-        PHI=2.0*HIPR1(40)*RANART(NSEED)
-c        PHI=0.
+c        PHI=2.0*HIPR1(40)*RANART(NSEED)
+        PHI=0.
 cbz6/28/99 flow1 end
         BBX=BB*COS(PHI)
         BBY=BB*SIN(PHI)

diff --git a/TAmpt/AliGenAmpt.cxx b/TAmpt/AliGenAmpt.cxx
index 05a2a62cec6d2e99f5df70dbac06ed7e4c6cc7c6..b166057e372b1a6333580b980c9d20e6ed4eb35a 100644 (file)
--- a/TAmpt/AliGenAmpt.cxx
+++ b/TAmpt/AliGenAmpt.cxx
@@ -83,7 +83,8 @@ AliGenAmpt::AliGenAmpt()
     fStringB(0.9),
     fEventTime(0.),
     fHeader(new AliGenAmptEventHeader("Ampt")),
-    fDecay(kTRUE)
+    fDecay(kTRUE),
+    fRotating(kFALSE)
 {
   // Constructor
   fEnergyCMS = 2760.;
@@ -137,7 +138,8 @@ AliGenAmpt::AliGenAmpt(Int_t npart)
     fStringB(0.9),
     fEventTime(0.),
     fHeader(new AliGenAmptEventHeader("Ampt")),
-    fDecay(kTRUE)
+    fDecay(kTRUE),
+    fRotating(kFALSE)
 {
   // Default PbPb collisions at 2.76 TeV
 
@@ -228,6 +230,9 @@ void AliGenAmpt::Init()
   fAmpt->Initialize();
   if (fEvaluate) 
     EvaluateCrossSections();
+
+  fAmpt->SetReactionPlaneAngle(0.0);
+  fRotating=kFALSE;
 }
 
 void AliGenAmpt::Generate()
@@ -264,6 +269,13 @@ void AliGenAmpt::Generate()
   Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
 
   while(1) {
+
+    // Generate random reaction plane angle if requested
+    if( fRotating ) {     
+      TRandom *r=AliAmptRndm::GetAmptRandom();
+      fAmpt->SetReactionPlaneAngle(TMath::TwoPi()*r->Rndm());
+    }
+
     // Generate one event
     Int_t fpemask = gSystem->GetFPEMask();
     gSystem->SetFPEMask(0);
@@ -271,6 +283,8 @@ void AliGenAmpt::Generate()
     gSystem->SetFPEMask(fpemask);
     fTrials++;
     fNprimaries = 0;
+
+
     fAmpt->ImportParticles(&fParticles,"All");
     Int_t np = fParticles.GetEntriesFast();
     if (np == 0 ) 
@@ -670,7 +684,8 @@ void AliGenAmpt::MakeHeader()
                         fAmpt->GetN11());
   fHeader->SetSpectators(fProjectileSpecn, fProjectileSpecp,
                         fTargetSpecn,fTargetSpecp);
-  fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
+  //fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
+  fHeader->SetReactionPlaneAngle(fAmpt->GetReactionPlaneAngle());
   //printf("Impact Parameter %13.3f \n", fAmpt->GetHINT1(19));
 
   // 4-momentum vectors of the triggered jets.

diff --git a/TAmpt/AliGenAmpt.h b/TAmpt/AliGenAmpt.h
index bee50ea6c52e3d2892178b0a253149f819c8152d..2ef9b9610ad823a0a1889c34eb5ecf2e723f9a02 100644 (file)
--- a/TAmpt/AliGenAmpt.h
+++ b/TAmpt/AliGenAmpt.h
@@ -88,6 +88,7 @@ class AliGenAmpt : public AliGenMC
     virtual TGraph* BinaryCollisions() {return fDnDb;}
     virtual Bool_t  CheckTrigger();
     virtual Bool_t  IsThisAKnownParticle(TParticle *thisGuy);
+    virtual void    SetRandomReactionPlane(Bool_t on=kTRUE){fRotating=on;};
 
   protected:
     Bool_t      SelectFlavor(Int_t pid);
@@ -139,6 +140,8 @@ class AliGenAmpt : public AliGenMC
     Float_t                  fEventTime;       // The event time
     AliGenHijingEventHeader *fHeader;          // header
     Bool_t                   fDecay;           // decay "long-lived" particles
+    Bool_t                   fRotating;        // Random reaction plane on/off
+ 
 
   private:
     AliGenAmpt(const AliGenAmpt &Ampt);

diff --git a/TAmpt/TAmpt.cxx b/TAmpt/TAmpt.cxx
index ad487c05886569ed6ef754c9b0d0cc6c73d60269..28aa16708ace654a05b3cb8d6a80d596ce1c0d01 100644 (file)
--- a/TAmpt/TAmpt.cxx
+++ b/TAmpt/TAmpt.cxx
@@ -77,7 +77,8 @@ TAmpt::TAmpt()
     fIat(208),
     fIzt(82),
     fBmin(0.),
-    fBmax(5.)
+    fBmax(5.),
+    fPsi(0.)
 {
   // Default constructor 
   amptsetdef();
@@ -98,7 +99,8 @@ TAmpt::TAmpt(Double_t efrm, const char *frame="CMS",
     fIat(iat),
     fIzt(izt),
     fBmin(bmin),
-    fBmax(bmax)
+    fBmax(bmax),
+    fPsi(0.)
 {
   // TAmpt constructor: 
   // Note that there may be only one functional TAmpt object
@@ -136,6 +138,13 @@ TObjArray* TAmpt::ImportParticles(Option_t */*option*/)
     //Double_t vy = HBT.xlast[i][1]*1e-12;//mm
     //Double_t vz = HBT.xlast[i][2]*1e-12;//mm
     //Double_t vt = HBT.xlast[i][3]*1e-12;//mm/c
+
+    // Rotate xy components of vectors to new reaction plane:
+    Double_t pxp=TMath::Cos(fPsi)*px - TMath::Sin(fPsi)*py;
+    Double_t pyp=TMath::Sin(fPsi)*px + TMath::Cos(fPsi)*py;
+    Double_t vxp=TMath::Cos(fPsi)*vx - TMath::Sin(fPsi)*vy;
+    Double_t vyp=TMath::Sin(fPsi)*vx + TMath::Cos(fPsi)*vy;
+
     Int_t pdg   = invflv(HBT.lblast[i]);
     TParticle *p = new TParticle(pdg,
                                  status,
@@ -143,12 +152,12 @@ TObjArray* TAmpt::ImportParticles(Option_t */*option*/)
                                  -1,
                                  -1,
                                  -1,
-                                 px,
-                                 py,
+                                 pxp,
+                                 pyp,
                                  pz,
                                  TMath::Sqrt(ma*ma+px*px+py*py+pz*pz),
-                                 vx,
-                                 vy,
+                                 vxp,
+                                 vyp,
                                  vz,
                                  vt);
     if((px==0)&&(py==0)) {
@@ -184,14 +193,21 @@ Int_t TAmpt::ImportParticles(TClonesArray *particles, Option_t */*option*/)
     Double_t py = HBT.plast[i][1];//GeV/c
     Double_t pz = HBT.plast[i][2];//GeV/c
     Double_t ma = HBT.plast[i][3];//GeV/c/c
-//    Double_t vx = 0;//HBT.xlast[i][0]*1e-12;//mm
-//    Double_t vy = 0;//HBT.xlast[i][1]*1e-12;//mm
-//    Double_t vz = 0;//HBT.xlast[i][2]*1e-12;//mm
-//    Double_t vt = 0;//HBT.xlast[i][3]*1e-12;//mm/c
-    Double_t vx = HBT.xlast[i][0]*1e-12;//mm
-    Double_t vy = HBT.xlast[i][1]*1e-12;//mm
-    Double_t vz = HBT.xlast[i][2]*1e-12;//mm
-    Double_t vt = HBT.xlast[i][3]*1e-12;//mm/c
+    Double_t vx = 0;//HBT.xlast[i][0]*1e-12;//mm
+    Double_t vy = 0;//HBT.xlast[i][1]*1e-12;//mm
+    Double_t vz = 0;//HBT.xlast[i][2]*1e-12;//mm
+    Double_t vt = 0;//HBT.xlast[i][3]*1e-12;//mm/c
+    //Double_t vx = HBT.xlast[i][0]*1e-12;//mm
+    //Double_t vy = HBT.xlast[i][1]*1e-12;//mm
+    //Double_t vz = HBT.xlast[i][2]*1e-12;//mm
+    //Double_t vt = HBT.xlast[i][3]*1e-12;//mm/c
+
+    // Rotate xy components of vectors to new reaction plane:
+    Double_t pxp=TMath::Cos(fPsi)*px - TMath::Sin(fPsi)*py;
+    Double_t pyp=TMath::Sin(fPsi)*px + TMath::Cos(fPsi)*py;
+    Double_t vxp=TMath::Cos(fPsi)*vx - TMath::Sin(fPsi)*vy;
+    Double_t vyp=TMath::Sin(fPsi)*vx + TMath::Cos(fPsi)*vy;
+
     Int_t pdg  = invflv(HBT.lblast[i]);
     //printf("i %d pdg %d px %f py %f pz %f vx %f vy %f vz %f vt %f\n", i, pdg, px, py, pz, vx, vy, vz, vt);
     new(particlesR[i]) TParticle(pdg,
@@ -200,12 +216,12 @@ Int_t TAmpt::ImportParticles(TClonesArray *particles, Option_t */*option*/)
                                  -1,
                                  -1,
                                  -1,
-                                 px,
-                                 py,
+                                 pxp,
+                                 pyp,
                                  pz,
                                  TMath::Sqrt(ma*ma+px*px+py*py+pz*pz),
-                                 vx,
-                                 vy,
+                                 vxp,
+                                 vyp,
                                  vz,
                                  vt);
     if((px==0)&&(py==0)){
@@ -235,6 +251,11 @@ Int_t TAmpt::ImportNucleons(TClonesArray *nucleons, Option_t */*option*/)
     Double_t x = HJCRDN.yp[i][0] + 0.5*GetBB();
     Double_t y = HJCRDN.yp[i][1];
     Double_t z = HJCRDN.yp[i][2];
+
+    // Rotate xy components of vectors to new reaction plane:
+    Double_t xp=TMath::Cos(fPsi)*x - TMath::Sin(fPsi)*y;
+    Double_t yp=TMath::Sin(fPsi)*x + TMath::Cos(fPsi)*y;
+
     Int_t    p = HSTRNG.nfp[3][i];
     Int_t    s = HSTRNG.nfp[4][i];
     new(nucleonsR[i]) TParticle(p,
@@ -247,8 +268,8 @@ Int_t TAmpt::ImportNucleons(TClonesArray *nucleons, Option_t */*option*/)
                                  0,
                                  0,
                                  0,
-                                 x,
-                                 y,
+                                 xp,
+                                 yp,
                                  z,
                                  0);
     nucleonsR[i]->SetUniqueID(1);
@@ -258,6 +279,11 @@ Int_t TAmpt::ImportNucleons(TClonesArray *nucleons, Option_t */*option*/)
     Double_t x = HJCRDN.yt[i][0] - 0.5*HPARNT.hint1[18];
     Double_t y = HJCRDN.yt[i][1];
     Double_t z = HJCRDN.yt[i][2];
+
+    // Rotate xy components of vectors to new reaction plane:
+    Double_t xp=TMath::Cos(fPsi)*x - TMath::Sin(fPsi)*y;
+    Double_t yp=TMath::Sin(fPsi)*x + TMath::Cos(fPsi)*y;
+
     Int_t    p = HSTRNG.nft[3][i];
     Int_t    s = HSTRNG.nft[4][i];
     new(nucleonsR[nA+i]) TParticle(p,
@@ -270,8 +296,8 @@ Int_t TAmpt::ImportNucleons(TClonesArray *nucleons, Option_t */*option*/)
                                    0,
                                    0,
                                    0,
-                                   x,
-                                   y,
+                                   xp,
+                                   yp,
                                    z,
                                    0);
     nucleonsR[nA+i]->SetUniqueID(-1);

diff --git a/TAmpt/TAmpt.h b/TAmpt/TAmpt.h
index b6c6c685359df8de8cc57d49dc00267ca860b319..84a09f79376f0107762a0a824ad865c263b42bc4 100644 (file)
--- a/TAmpt/TAmpt.h
+++ b/TAmpt/TAmpt.h
@@ -138,6 +138,8 @@ public:
   virtual void        SetXmu(Float_t m);
   virtual void        SetAlpha(Float_t alpha);
   virtual void        SetStringFrag(Float_t a, Float_t b);
+  virtual void        SetReactionPlaneAngle(Double_t psi){fPsi=psi;};  // Set reaction plane angle
+  virtual Double_t    GetReactionPlaneAngle(){return fPsi;};           // Get reaction plane angle
 
 protected:
   Double_t     fEfrm;     // Energy in the centre of mass (CMS) or lab-frame (LAB)
@@ -150,7 +152,7 @@ protected:
   Int_t        fIzt;      // Charge number of target
   Double_t     fBmin;     // Minimum impact parameter
   Double_t     fBmax;     // Maximum impact parameter
- 
-  ClassDef(TAmpt,2)  //Interface to Ampt Event Generator
+  Double_t     fPsi;      // Reaction plane angle
+  ClassDef(TAmpt,3)  //Interface to Ampt Event Generator
 };
 #endif