Update for change in AliTRDsimple

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

diff --git a/EVGEN/AliGenHijing.cxx b/EVGEN/AliGenHijing.cxx
index fba6a3adf9903b21e06b183d6ba5d047cd128380..9bd097b9daf447a3d6a98261d599f92bc5e6dd68 100644 (file)
--- a/EVGEN/AliGenHijing.cxx
+++ b/EVGEN/AliGenHijing.cxx
@@ -15,6 +15,17 @@
 
 /*
 $Log$
+Revision 1.30  2001/10/21 18:35:56  hristov
+Several pointers were set to zero in the default constructors to avoid memory management problems
+
+Revision 1.29  2001/10/15 08:12:24  morsch
+- Vertex smearing with truncated gaussian.
+- Store triggered jet info before and after final state radiation into mc-heade
+
+Revision 1.28  2001/10/08 11:55:25  morsch
+Store 4-momenta of trigegred jets in event header.
+Possibility to switch of initial and final state radiation.
+
 Revision 1.27  2001/10/08 07:13:14  morsch
 Add setter for minimum transverse momentum of triggered jet.
 
@@ -144,6 +155,7 @@ AliGenHijing::AliGenHijing(Int_t npart)
     SetImpactParameterRange();
     SetTarget();
     SetProjectile();
+    SetBoostLHC();
     fKeep       =  0;
     fQuench     =  1;
     fShadowing  =  1;
@@ -157,9 +169,12 @@ AliGenHijing::AliGenHijing(Int_t npart)
     fDnDb       =  0;
     fPtMinJet   = -2.5;        
     fRadiation  =  1;
+    fEventVertex.Set(3);
 //
 // Set random number generator   
     sRandom = fRandom;
+    fHijing = 0;
+
 }
 
 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
@@ -228,11 +243,16 @@ void AliGenHijing::Generate()
   fTrials = 0;
   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
   if(fVertexSmear == kPerEvent) {
-       Rndm(random,6);
-       for (j=0; j < 3; j++) {
-           origin0[j] += fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
-               TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
-       }
+      Float_t dv[3];
+      dv[2] = 1.e10;
+      while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) {
+         Rndm(random,6);
+         for (j=0; j < 3; j++) {
+             dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
+                 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
+         }
+      }
+      for (j=0; j < 3; j++) origin0[j] += dv[j];
   } else if (fVertexSmear == kPerTrack) {
 //         fHijing->SetMSTP(151,0);
       for (j = 0; j < 3; j++) {
@@ -249,6 +269,8 @@ void AliGenHijing::Generate()
       fHijing->GenerateEvent();
       fTrials++;
       fHijing->ImportParticles(particles,"All");
+      if (fLHC) Boost(particles);
+      
       Int_t np = particles->GetEntriesFast();
       printf("\n **************************************************%d\n",np);
       Int_t nc = 0;
@@ -257,14 +279,21 @@ void AliGenHijing::Generate()
       Int_t * newPos = new Int_t[np];
 
       for (i = 0; i < np; i++) *(newPos+i) = i;
+//      Get event vertex
+//
+      TParticle *  iparticle = (TParticle *) particles->At(0);
+      fEventVertex[0] = origin0[0];
+      fEventVertex[1] = origin0[1];    
+      fEventVertex[2] = origin0[2];
+      
 //
 //      First write parent particles
 //
 
       for (i = 0; i < np; i++) {
-         TParticle *  iparticle       = (TParticle *) particles->At(i);
+         iparticle       = (TParticle *) particles->At(i);
 // Is this a parent particle ?
-        if (Stable(iparticle)) continue;
+         if (Stable(iparticle)) continue;
 //
         Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
         Bool_t  selected             =  kTRUE;
@@ -500,6 +529,39 @@ Bool_t AliGenHijing::Stable(TParticle*  particle)
     }
 }
 
+
+void AliGenHijing::Boost(TClonesArray* particles)
+{
+//
+// Boost cms into LHC lab frame
+//
+    Double_t dy    = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) / 
+                                     (Double_t(fZTarget)    * Double_t(fAProjectile)));
+    Double_t beta  = TMath::TanH(dy);
+    Double_t gamma = 1./TMath::Sqrt(1.-beta*beta);
+    Double_t gb    = gamma * beta;
+
+    printf("\n Boosting particles to lab frame %f %f %f", dy, beta, gamma);
+    
+    Int_t i;
+    Int_t np = particles->GetEntriesFast();
+    for (i = 0; i < np; i++) 
+    {
+       TParticle* iparticle = (TParticle*) particles->At(i);
+
+       Double_t e   = iparticle->Energy();
+       Double_t px  = iparticle->Px();
+       Double_t py  = iparticle->Py();
+       Double_t pz  = iparticle->Pz();
+
+       Double_t eb  = gamma * e -      gb * pz;
+       Double_t pzb =   -gb * e +   gamma * pz;
+
+       iparticle->SetMomentum(px, py, pzb, eb);
+    }
+}
+
+
 void AliGenHijing::MakeHeader()
 {
 // Builds the event header, to be called after each event
@@ -515,6 +577,9 @@ void AliGenHijing::MakeHeader()
                                                       fHijing->GetN11());
     ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp);
 
+// 4-momentum vectors of the triggered jets.
+//
+// Before final state gluon radiation.
     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
                                              fHijing->GetHINT1(22),
                                              fHijing->GetHINT1(23),
@@ -524,11 +589,23 @@ void AliGenHijing::MakeHeader()
                                              fHijing->GetHINT1(32),
                                              fHijing->GetHINT1(33),
                                              fHijing->GetHINT1(34));
-
-    ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2);
+// After final state gluon radiation.
+    TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
+                                             fHijing->GetHINT1(27),
+                                             fHijing->GetHINT1(28),
+                                             fHijing->GetHINT1(29));
+
+    TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
+                                             fHijing->GetHINT1(37),
+                                             fHijing->GetHINT1(38),
+                                             fHijing->GetHINT1(39));
+    ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
+// Event Vertex
+    header->SetPrimaryVertex(fEventVertex);
     gAlice->SetGenEventHeader(header);    
 }
 
+
 AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
 {
 // Assignment operator