- Update 11a10* (MC_PbPb) (also update for splines)

[u/mrichter/AliRoot.git] / TAmpt / AliGenAmpt.cxx

diff --git a/TAmpt/AliGenAmpt.cxx b/TAmpt/AliGenAmpt.cxx
index 05a2a62cec6d2e99f5df70dbac06ed7e4c6cc7c6..0424e0be49c057e59bc07c805adff84320c3a3a0 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()
@@ -241,9 +246,6 @@ void AliGenAmpt::Generate()
   Float_t p[3];
   Float_t tof;
 
-  //  converts from mm/c to s
-  const Float_t kconv = 0.001/2.99792458e8;
-
   Int_t nt  = 0;
   Int_t jev = 0;
   Int_t j, kf, ks, ksp, imo;
@@ -264,6 +266,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,11 +280,23 @@ void AliGenAmpt::Generate()
     gSystem->SetFPEMask(fpemask);
     fTrials++;
     fNprimaries = 0;
+
+
     fAmpt->ImportParticles(&fParticles,"All");
     Int_t np = fParticles.GetEntriesFast();
     if (np == 0 ) 
       continue;
-
+    //
+    //RS>>: Decayers now returns cm and sec. Since TAmpt returns mm and mm/c, convert its result to cm and sec here
+    const Float_t kconvT=0.001/2.999792458e8; // mm/c to seconds conversion
+    const Float_t kconvL=1./10; // mm to cm conversion
+    for (int ip=np;ip--;) {
+      TParticle* part = (TParticle*)fParticles[ip];
+      if (!part) continue;
+      part->SetProductionVertex(part->Vx()*kconvL,part->Vy()*kconvL,part->Vz()*kconvL,kconvT*part->T());
+    }
+    // RS<<
+    //
     if (fTrigger != kNoTrigger) {
       if (!CheckTrigger()) 
         continue;
@@ -338,19 +359,27 @@ void AliGenAmpt::Generate()
                  //arr.Print();
                  // iparticle->SetStatusCode(2);  to be compatible with Hijing
                  iparticle->SetFirstDaughter(np2);
-                 for (Int_t jj = 1; jj < ndecayed; jj++) {
+                 for (Int_t jj = 1; jj < ndecayed; jj++) {
                  TParticle *jp = (TParticle *)arr.At(jj);
                    if (jp->GetFirstMother()!=1)
                      continue;
-                 TParticle *newp = new(fParticles[np2]) TParticle(jp->GetPdgCode(),
+
+                   TParticle *newp = new(fParticles[np2]) TParticle(jp->GetPdgCode(),
                                                                   0, //1,  //to be compatible with Hijing
                                                                   i,
                                                                   -1,
                                                                   -1,
                                                                   -1,
                                                                   jp->Px(),jp->Py(),jp->Pz(),jp->Energy(),
-                                                                  jp->Vx(),jp->Vy(),jp->Vz(),jp->T());
-                 newp->SetUniqueID( jp->GetStatusCode() );
+                                                                    jp->Vx(),jp->Vy(),jp->Vz(),jp->T());
+                   //take care of the phi
+                   //if((kf == 333)||(kf == 313)) {
+                   if(IsThisAKnownParticle(iparticle)) {
+                     //Printf("=============PANOS===================");
+                     //Printf("Phi detected - daughet is: %d",jp->GetPdgCode());
+                     newp->SetUniqueID(4);
+                   }
+                   else newp->SetUniqueID( jp->GetStatusCode() );
                  np2++;
                } // end of jj->nDecayedParticles
                iparticle->SetLastDaughter(np2-1);
@@ -461,10 +490,10 @@ void AliGenAmpt::Generate()
         p[0] = iparticle->Px();
         p[1] = iparticle->Py();
         p[2] = iparticle->Pz() * sign;
-        origin[0] = origin0[0]+iparticle->Vx()/10;
-        origin[1] = origin0[1]+iparticle->Vy()/10;
-        origin[2] = origin0[2]+iparticle->Vz()/10;
-       tof = time0+kconv * iparticle->T();
+        origin[0] = origin0[0]+iparticle->Vx();
+        origin[1] = origin0[1]+iparticle->Vy();
+        origin[2] = origin0[2]+iparticle->Vz();
+       tof = time0 + iparticle->T();
 
         imo = -1;
         TParticle* mother = 0;
@@ -670,7 +699,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.