]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenHIJINGpara.cxx
J/Psi from B (kJpsiFromB) added
[u/mrichter/AliRoot.git] / EVGEN / AliGenHIJINGpara.cxx
index 5488b103f8bc7007dcccd8fa698a4dd96d9ad419..272a1868d5e9997d747ca0bad5273ccfb67cd9a3 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.13  2002/10/14 14:55:35  hristov
-Merging the VirtualMC branch to the main development branch (HEAD)
-
-Revision 1.11.4.1  2002/07/24 08:56:28  alibrary
-Updating EVGEN on TVirtulaMC
-
-Revision 1.12  2002/06/19 06:56:34  hristov
-Memory leak corrected
-
-Revision 1.11  2002/03/20 10:21:13  hristov
-Set fPtMax to 15 GeV in order to avoid some numerical problems
-
-Revision 1.10  2001/10/15 16:44:46  morsch
-- Possibility for vertex distribution truncation.
-- Write mc header with vertex position.
-
-Revision 1.9  2001/07/27 17:09:36  morsch
-Use local SetTrack, KeepTrack and SetHighWaterMark methods
-to delegate either to local stack or to stack owned by AliRun.
-(Piotr Skowronski, A.M.)
-
-Revision 1.8  2001/07/20 11:03:58  morsch
-Issue warning message if used outside allowed eta range (-8 to 8).
-
-Revision 1.7  2001/07/17 12:41:01  morsch
-- Calculation of fraction of event corresponding to selected pt-range corrected
-(R. Turrisi)
-- Parent weight corrected.
-
-Revision 1.6  2001/05/16 14:57:10  alibrary
-New files for folders and Stack
-
-Revision 1.5  2000/12/21 16:24:06  morsch
-Coding convention clean-up
-
-Revision 1.4  2000/11/30 07:12:50  alibrary
-Introducing new Rndm and QA classes
-
-Revision 1.3  2000/10/02 21:28:06  fca
-Removal of useless dependecies via forward declarations
-
-Revision 1.2  2000/07/11 18:24:55  fca
-Coding convention corrections + few minor bug fixes
-
-Revision 1.1  2000/06/09 20:20:30  morsch
-Same class as previously in AliSimpleGen.cxx
-All coding rule violations except RS3 corrected (AM)
-
-*/
+/* $Id$ */
 
 // Parameterisation of pi and K, eta and pt distributions
 // used for the ALICE TDRs.
@@ -87,25 +37,28 @@ All coding rule violations except RS3 corrected (AM)
 //                                                               //
 ///////////////////////////////////////////////////////////////////
 
-#include "AliGenHIJINGpara.h"
-#include "AliGenEventHeader.h"
-#include "AliRun.h"
-#include "AliConst.h"
-#include "AliDecayer.h"
-#include "AliDecayerPythia.h"
-#include "AliPDG.h"
-
-#include <TF1.h>
 #include <TArrayF.h>
-#include <TDatabasePDG.h>
 #include <TClonesArray.h>
+#include <TDatabasePDG.h>
+#include <TF1.h>
 #include <TParticle.h>
+#include <TPDGCode.h>
+#include <TVirtualMC.h>
+
+#include "AliConst.h"
+#include "AliDecayer.h"
+#include "AliGenEventHeader.h"
+#include "AliGenHIJINGpara.h"
+#include "AliRun.h"
 
 ClassImp(AliGenHIJINGpara)
 
-AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para)
+
+AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para):
+    AliGenerator(para)
 {
-// copy constructor
+// Copy constructor
+    para.Copy(*this);
 }
 
 //_____________________________________________________________________________
@@ -211,11 +164,12 @@ AliGenHIJINGpara::AliGenHIJINGpara()
     //
     // Default constructor
     //
-    fPtpi    = 0;
-    fPtka    = 0;
-    fETApic  = 0;
-    fETAkac  = 0;
-    fDecayer = 0;
+    fPtpi    =  0;
+    fPtka    =  0;
+    fETApic  =  0;
+    fETAkac  =  0;
+    fDecayer =  0;
+    fNt      = -1;
     SetCutVertexZ();
     SetPtRange();
     SetPi0Decays();
@@ -228,13 +182,14 @@ AliGenHIJINGpara::AliGenHIJINGpara(Int_t npart)
   // 
   // Standard constructor
   //
-    fName="HIGINGpara";
+    fName="HIJINGpara";
     fTitle="HIJING Parametrisation Particle Generator";
-    fPtpi    = 0;
-    fPtka    = 0;
-    fETApic  = 0;
-    fETAkac  = 0;
-    fDecayer = 0;
+    fPtpi    =  0;
+    fPtka    =  0;
+    fETApic  =  0;
+    fETAkac  =  0;
+    fDecayer =  0;
+    fNt      = -1;
     SetCutVertexZ();
     SetPtRange();
     SetPi0Decays();
@@ -267,11 +222,11 @@ void AliGenHIJINGpara::Init()
     fETApic = new TF1("etapic",&etapic,etaMin,etaMax,0);
     fETAkac = new TF1("etakac",&etakac,etaMin,etaMax,0);
 
-    TF1 etaPic0("etapic",&etapic,-7,7,0);
-    TF1 etaKac0("etakac",&etakac,-7,7,0);
+    TF1 etaPic0("etaPic0",&etapic,-7,7,0);
+    TF1 etaKac0("etaKac0",&etakac,-7,7,0);
 
-    TF1 ptPic0("ptpi",&ptpi,0.,15.,0);
-    TF1 ptKac0("ptka",&ptka,0.,15.,0);
+    TF1 ptPic0("ptPic0",&ptpi,0.,15.,0);
+    TF1 ptKac0("ptKac0",&ptka,0.,15.,0);
 
     Float_t intETApi  = etaPic0.Integral(-0.5, 0.5);
     Float_t intETAka  = etaKac0.Integral(-0.5, 0.5);
@@ -291,8 +246,16 @@ void AliGenHIJINGpara::Init()
 //  Fraction of events corresponding to the selected phi-range    
     Float_t phiFrac    = (fPhiMax-fPhiMin)/2/TMath::Pi();
 
+    
     fParentWeight = Float_t(fNpart)/(intETASel*ptFrac*phiFrac);
     
+    if (fAnalog != 0) {
+       fPtWgtPi = (fPtMax - fPtMin) / fPtpi->Integral(0., 20.);
+       fPtWgtKa = (fPtMax - fPtMin) / fPtka->Integral(0., 20.);
+       fParentWeight = Float_t(fNpart)/(intETASel*phiFrac);
+    }
+    
+    
     printf("%s: The number of particles in the selected kinematic region corresponds to %f percent of a full event\n ", 
           ClassName(),100.*fParentWeight);
 
@@ -306,7 +269,8 @@ void AliGenHIJINGpara::Init()
     }
 //
 //
-    if (fPi0Decays) fDecayer = new AliDecayerPythia();
+    if (fPi0Decays && gMC)
+       fDecayer = gMC->GetDecayer();
 }
 
 
@@ -326,10 +290,10 @@ void AliGenHIJINGpara::Generate()
     const Int_t kKaons[4] = {kK0Long, kK0Short, kKPlus, kKMinus};
     //
     Float_t origin[3];
-    Float_t pt, pl, ptot;
+    Float_t pt, pl, ptot, wgt;
     Float_t phi, theta;
     Float_t p[3];
-    Int_t i, part, nt, j;
+    Int_t i, part, j;
     //
     TF1 *ptf;
     TF1 *etaf;
@@ -339,39 +303,40 @@ void AliGenHIJINGpara::Generate()
     for (j=0;j<3;j++) origin[j]=fOrigin[j];
 
     if(fVertexSmear == kPerEvent) {
-       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++) origin[j] += dv[j];
+       Vertex();
+       for (j=0; j < 3; j++) origin[j] = fVertex[j];
     } // if kPerEvent
     TArrayF eventVertex;
     eventVertex.Set(3);
     eventVertex[0] = origin[0];
     eventVertex[1] = origin[1];
     eventVertex[2] = origin[2];
-
+     
     for(i=0;i<fNpart;i++) {
        while(1) {
-           Rndm(random,3);
+           Rndm(random,4);
            if(random[0]<kBorne) {
                part=kPions[Int_t (random[1]*3)];
                ptf=fPtpi;
                etaf=fETApic;
+               wgt = fPtWgtPi;
            } else {
                part=kKaons[Int_t (random[1]*4)];
                ptf=fPtka;
                etaf=fETAkac;
+               wgt = fPtWgtKa;
            }
            phi=fPhiMin+random[2]*(fPhiMax-fPhiMin);
            theta=2*TMath::ATan(TMath::Exp(-etaf->GetRandom()));
            if(theta<fThetaMin || theta>fThetaMax) continue;
-           pt=ptf->GetRandom();
+        
+           if (fAnalog == 0) { 
+               pt = ptf->GetRandom();
+           } else {
+               pt = fPtMin + random[3] * (fPtMax - fPtMin);
+           }
+           
+           
            pl=pt/TMath::Tan(theta);
            ptot=TMath::Sqrt(pt*pt+pl*pl);
            if(ptot<fPMin || ptot>fPMax) continue;
@@ -385,13 +350,31 @@ void AliGenHIJINGpara::Generate()
                        TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
                }
            }
-           SetTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,nt,fParentWeight);
+           
+           if (fAnalog == 0) { 
+               wgt = fParentWeight;
+           } else {
+               wgt *= (fParentWeight * ptf->Eval(pt));
+           }
+           
+           
+           if (part == kPi0 && fPi0Decays){
 //
 //          Decay pi0 if requested
-           if (part == kPi0 && fPi0Decays) DecayPi0(origin, p);
+               PushTrack(0,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight);
+               KeepTrack(fNt);
+               DecayPi0(origin, p);
+           } else {
+               PushTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight);
+               KeepTrack(fNt);
+           }
+
            break;
        }
+       SetHighWaterMark(fNt);
     }
+//
+
 // Header
     AliGenEventHeader* header = new AliGenEventHeader("HIJINGparam");
 // Event Vertex
@@ -399,12 +382,6 @@ void AliGenHIJINGpara::Generate()
     gAlice->SetGenEventHeader(header); 
 }
 
-AliGenHIJINGpara& AliGenHIJINGpara::operator=(const  AliGenHIJINGpara& rhs)
-{
-// Assignment operator
-    return *this;
-}
-
 void AliGenHIJINGpara::SetPtRange(Float_t ptmin, Float_t ptmax) {
     AliGenerator::SetPtRange(ptmin, ptmax);
 }
@@ -419,25 +396,33 @@ void AliGenHIJINGpara::DecayPi0(Float_t* orig, Float_t * p)
     if(!particles) particles = new TClonesArray("TParticle",1000);
 //    
     const Float_t kMass = TDatabasePDG::Instance()->GetParticle(kPi0)->Mass();
-    Float_t       e     = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] +p[2] * p[2] * kMass * kMass);
+    Float_t       e     = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]+ kMass * kMass);
 //
 //  Decay the pi0    
     TLorentzVector pmom(p[0], p[1], p[2], e);
     fDecayer->Decay(kPi0, &pmom);
+    
 //
 // Put decay particles on the stack
 //
-    Int_t nt = 0;
     Float_t polar[3] = {0., 0., 0.};
     Int_t np = fDecayer->ImportParticles(particles);
-    
+    Int_t nt;    
     for (Int_t i = 1; i < np; i++)
     {
-       TParticle* iParticle =  (TParticle *) particles->At(1);
-       p[0]=iParticle->Px();
-       p[1]=iParticle->Py();
-       p[2]=iParticle->Pz();
+       TParticle* iParticle =  (TParticle *) particles->At(i);
+       p[0] = iParticle->Px();
+       p[1] = iParticle->Py();
+       p[2] = iParticle->Pz();
        Int_t part = iParticle->GetPdgCode();
-       SetTrack(fTrackIt, 0, part, p, orig, polar, 0, kPDecay, nt, fParentWeight);
+
+       PushTrack(fTrackIt, fNt, part, p, orig, polar, 0, kPDecay, nt, fParentWeight);
+       KeepTrack(nt);
     }
+    fNt = nt;
+}
+
+void AliGenHIJINGpara::Copy(AliGenHIJINGpara &) const
+{
+  Fatal("Copy","Not implemented!\n");
 }