]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenPythia.cxx
Correcting coding convention violations
[u/mrichter/AliRoot.git] / EVGEN / AliGenPythia.cxx
index 03dbfcdd286db6de331b1134babf28768b408727..8ab6626d0f418926ea99c72d4a6ec8c2b5d72840 100644 (file)
 
 /*
 $Log$
+Revision 1.29  2000/12/04 11:22:03  morsch
+Init of sRandom as in 1.15
+
+Revision 1.28  2000/12/02 11:41:39  morsch
+Use SetRandom() to initialize random number generator in constructor.
+
+Revision 1.27  2000/11/30 20:29:02  morsch
+Initialise static variable sRandom in constructor: sRandom = fRandom;
+
+Revision 1.26  2000/11/30 07:12:50  alibrary
+Introducing new Rndm and QA classes
+
+Revision 1.25  2000/10/18 19:11:27  hristov
+Division by zero fixed
+
+Revision 1.24  2000/09/18 10:41:35  morsch
+Add possibility to use nuclear structure functions from PDF library V8.
+
+Revision 1.23  2000/09/14 14:05:40  morsch
+dito
+
+Revision 1.22  2000/09/14 14:02:22  morsch
+- Correct conversion from mm to cm when passing particle vertex to MC.
+- Correct handling of fForceDecay == all.
+
+Revision 1.21  2000/09/12 14:14:55  morsch
+Call fDecayer->ForceDecay() at the beginning of Generate().
+
+Revision 1.20  2000/09/06 14:29:33  morsch
+Use AliPythia for event generation an AliDecayPythia for decays.
+Correct handling of "nodecay" option
+
+Revision 1.19  2000/07/11 18:24:56  fca
+Coding convention corrections + few minor bug fixes
+
+Revision 1.18  2000/06/30 12:40:34  morsch
+Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
+Geant units (cm).
+
 Revision 1.17  2000/06/09 20:34:07  morsch
 All coding rule violations except RS3 corrected
 
@@ -41,16 +80,20 @@ Introduction of the Copyright and cvs Log
 */
 
 #include "AliGenPythia.h"
+#include "AliDecayerPythia.h"
 #include "AliRun.h"
 #include "AliPythia.h"
+#include "AliPDG.h"
 #include <TParticle.h>
+#include <TSystem.h>
 
  ClassImp(AliGenPythia)
 
 AliGenPythia::AliGenPythia()
                  :AliGenerator()
 {
-// Constructor
+// Default Constructor
+  fDecayer = new AliDecayerPythia();
 }
 
 AliGenPythia::AliGenPythia(Int_t npart)
@@ -61,6 +104,8 @@ AliGenPythia::AliGenPythia(Int_t npart)
 // structure function GRVHO
 //
     fXsection  = 0.;
+    fNucA1=0;
+    fNucA2=0;
     fParentSelect.Set(5);
     fChildSelect.Set(5);
     for (Int_t i=0; i<5; i++) fParentSelect[i]=fChildSelect[i]=0;
@@ -69,6 +114,9 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetForceDecay();
     SetPtHard();
     SetEnergyCMS();
+    fDecayer = new AliDecayerPythia();
+    // Set random number generator 
+    sRandom=fRandom;
 }
 
 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
@@ -84,19 +132,27 @@ AliGenPythia::~AliGenPythia()
 void AliGenPythia::Init()
 {
 // Initialisation
-    SetMC(new AliPythia());
+  SetMC(AliPythia::Instance());
     fPythia=(AliPythia*) fgMCEvGen;
 //
     fParentWeight=1./Float_t(fNpart);
 //
 //  Forward Paramters to the AliPythia object
-    fPythia->DefineParticles();
+    //    gSystem->Exec("ln -s $ALICE_ROOT/data/Decay.table fort.1");
+    //    fPythia->Pyupda(2,1);    
+    //    gSystem->Exec("rm fort.1");
+
+    fDecayer->SetForceDecay(fForceDecay);    
+    fDecayer->Init();
+
+
     fPythia->SetCKIN(3,fPtHardMin);
     fPythia->SetCKIN(4,fPtHardMax);    
+    if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);  
     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
-    fPythia->ForceDecay(fForceDecay);
-    fPythia->Lulist(0);
-    fPythia->Pystat(2);
+
+    //    fPythia->Pylist(0);
+    //    fPythia->Pystat(2);
 //  Parent and Children Selection
     switch (fProcess) 
     {
@@ -140,7 +196,7 @@ void AliGenPythia::Init()
     case dielectron:
     case b_jpsi_dielectron:
     case b_psip_dielectron:
-       fChildSelect[0]=11;     
+       fChildSelect[0]=kElectron;      
        break;
     case semimuonic:
     case dimuon:
@@ -148,17 +204,22 @@ void AliGenPythia::Init()
     case b_psip_dimuon:
     case pitomu:
     case katomu:
-       fChildSelect[0]=13;
+       fChildSelect[0]=kMuonMinus;
        break;
+    case hadronicD:
+      fChildSelect[0]=kPiPlus;
+      fChildSelect[1]=kKPlus;
+      break;
     case all:
     case nodecay:
-       break;
+      break;
     }
 }
 
 void AliGenPythia::Generate()
 {
 // Generate one event
+    fDecayer->ForceDecay();
 
     Float_t polar[3] =   {0,0,0};
     Float_t origin[3]=   {0,0,0};
@@ -181,19 +242,19 @@ void AliGenPythia::Generate()
     
     fTrials=0;
     for (j=0;j<3;j++) origin0[j]=fOrigin[j];
-    if(fVertexSmear==perEvent) {
+    if(fVertexSmear==kPerEvent) {
        fPythia->SetMSTP(151,1);
        for (j=0;j<3;j++) {
            fPythia->SetPARP(151+j, fOsigma[j]/10.);
        }
-    } else if (fVertexSmear==perTrack) {
+    } else if (fVertexSmear==kPerTrack) {
        fPythia->SetMSTP(151,0);
     }
     
     while(1)
     {
        fPythia->Pyevnt();
-//     fPythia->Lulist(1);
+//     fPythia->Pylist(1);
        fTrials++;
        fPythia->ImportParticles(particles,"All");
        Int_t np = particles->GetEntriesFast();
@@ -207,7 +268,7 @@ void AliGenPythia::Generate()
                kf = CheckPDGCode(iparticle->GetPdgCode());
                if (ks==21) continue;
 
-               fChildWeight=(fPythia->GetBraPart(kf))*fParentWeight;     
+               fChildWeight=(fDecayer->GetPartialBranchingRatio(kf))*fParentWeight;      
 //
 // Parent
                if (ParentSelected(TMath::Abs(kf))) {
@@ -227,7 +288,8 @@ void AliGenPythia::Generate()
                            pP[3] = massP;
                            gAlice->SetTrack(0,-1,-1,
                                             pP,originP,polar,
-                                            0,"Hard Scat.",ntP,fParentWeight);
+                                            0,kPPrimary,ntP,fParentWeight);
+//                                          0,"Hard Scat.",ntP,fParentWeight);
                            gAlice->KeepTrack(ntP);
                        }
                        nc++;
@@ -236,42 +298,46 @@ void AliGenPythia::Generate()
                        p[0]=iparticle->Px();
                        p[1]=iparticle->Py();
                        p[2]=iparticle->Pz();
-                       origin[0]=origin0[0]+iparticle->Vx()*10.;
-                       origin[1]=origin0[1]+iparticle->Vy()*10.;
-                       origin[2]=origin0[2]+iparticle->Vz()*10.;
+                       origin[0]=origin0[0]+iparticle->Vx()/10.;
+                       origin[1]=origin0[1]+iparticle->Vy()/10.;
+                       origin[2]=origin0[2]+iparticle->Vz()/10.;
 
                        Int_t ifch=iparticle->GetFirstDaughter();
                        Int_t ilch=iparticle->GetLastDaughter();        
-                       if (ifch !=0 && ilch !=0) {
-                           gAlice->SetTrack(0,ntP,kf,
+
+                       if ((ifch !=0 && ilch !=0) || fForceDecay == nodecay) {
+                         Int_t trackit=0;
+                         if (fForceDecay == nodecay) trackit = 1;
+                           gAlice->SetTrack(trackit,ntP,kf,
                                             p,origin,polar,
-                                            0,"Primary",nt,fParentWeight);
+                                            0,kPPrimary,nt,fParentWeight);
                            gAlice->KeepTrack(nt);
                            Int_t iparent = nt;
 //
 // Children        
-
-                           for (j=ifch; j<=ilch; j++)
-                           {
-                               TParticle *  ichild = 
+                           if (fForceDecay != nodecay) {
+                             for (j=ifch; j<=ilch; j++)
+                               {
+                                 TParticle *  ichild = 
                                    (TParticle *) particles->At(j-1);
-                               kf = CheckPDGCode(ichild->GetPdgCode());
+                                 kf = CheckPDGCode(ichild->GetPdgCode());
 //
 // 
-                               if (ChildSelected(TMath::Abs(kf))) {
-                                   origin[0]=origin0[0]+ichild->Vx()*10.;
-                                   origin[1]=origin0[1]+ichild->Vy()*10.;
-                                   origin[2]=origin0[2]+ichild->Vz()*10.;              
+                                 if (ChildSelected(TMath::Abs(kf))) {
+                                   origin[0]=origin0[0]+ichild->Vx()/10.;
+                                   origin[1]=origin0[1]+ichild->Vy()/10.;
+                                   origin[2]=origin0[2]+ichild->Vz()/10.;              
                                    p[0]=ichild->Px();
                                    p[1]=ichild->Py();
                                    p[2]=ichild->Pz();
                                    Float_t tof=kconv*ichild->T();
                                    gAlice->SetTrack(fTrackIt, iparent, kf,
                                                     p,origin,polar,
-                                                    tof,"Decay",nt,fChildWeight);
+                                                    tof,kPDecay,nt,fChildWeight);
                                    gAlice->KeepTrack(nt);
-                               } // select child
-                           } // child loop
+                                 } // select child
+                               } // child loop
+                           } 
                        }
                    } // kinematic selection
                } // select particle
@@ -289,12 +355,12 @@ void AliGenPythia::Generate()
                        p[0]=iparticle->Px();
                        p[1]=iparticle->Py();
                        p[2]=iparticle->Pz();
-                       origin[0]=origin0[0]+iparticle->Vx()*10;
-                       origin[1]=origin0[1]+iparticle->Vy()*10;
-                       origin[2]=origin0[2]+iparticle->Vz()*10;
+                       origin[0]=origin0[0]+iparticle->Vx()/10.;
+                       origin[1]=origin0[1]+iparticle->Vy()/10.;
+                       origin[2]=origin0[2]+iparticle->Vz()/10.;
                        Float_t tof=kconv*iparticle->T();
                        gAlice->SetTrack(fTrackIt,-1,kf,p,origin,polar,
-                                        tof,"Primary",nt);
+                                        tof,kPPrimary,nt);
                        gAlice->KeepTrack(nt);
                } // select particle
            } // particle loop 
@@ -305,6 +371,8 @@ void AliGenPythia::Generate()
            if (jev >= fNpart || fNpart == -1) {
                fKineBias=Float_t(fNpart)/Float_t(fTrials);
                printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
+// Print x-section summary
+               fPythia->Pystat(1);
                break;
            }
        }
@@ -328,6 +396,8 @@ Bool_t AliGenPythia::ParentSelected(Int_t ip)
 Bool_t AliGenPythia::ChildSelected(Int_t ip)
 {
 // True if particle is in list of decay products to be selected
+    if (fForceDecay == all) return kTRUE;
+    
     for (Int_t i=0; i<5; i++)
     {
        if (fChildSelect[i]==ip) return kTRUE;
@@ -371,11 +441,16 @@ Bool_t AliGenPythia::KinematicSelection(TParticle *particle)
 
 //
 // rapidity cut
-    Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
-    if (y > fYMax || y < fYMin)
-    {
+    if (e==pz) {
+      return kFALSE;
+    }
+    else {
+      Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
+      if (y > fYMax || y < fYMin)
+        {
 //     printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
-       return kFALSE;
+          return kFALSE;
+        }
     }
 
 //
@@ -407,6 +482,8 @@ Int_t AliGenPythia::CheckPDGCode(Int_t pdgcode)
 //
 //  If the particle is in a diffractive state, then take action accordingly
   switch (pdgcode) {
+  case 91:
+    return 92;
   case 110:
     //rho_diff0 -- difficult to translate, return rho0
     return 113;
@@ -432,7 +509,16 @@ Int_t AliGenPythia::CheckPDGCode(Int_t pdgcode)
   //non diffractive state -- return code unchanged
   return pdgcode;
 }
-                 
+
+    
+void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
+{
+// Treat protons as inside nuclei with mass numbers a1 and a2  
+    fNucA1 = a1;
+    fNucA2 = a2;
+}
+       
+         
 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 {
 // Assignment operator
@@ -441,7 +527,43 @@ AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 
 
 
-
-
+void AliGenPythia::Streamer(TBuffer &R__b)
+{
+   // Stream an object of class AliGenPythia.
+
+   if (R__b.IsReading()) {
+      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
+      AliGenerator::Streamer(R__b);
+      R__b >> (Int_t&)fProcess;
+      R__b >> (Int_t&)fStrucFunc;
+      R__b >> (Int_t&)fForceDecay;
+      R__b >> fEnergyCMS;
+      R__b >> fKineBias;
+      R__b >> fTrials;
+      fParentSelect.Streamer(R__b);
+      fChildSelect.Streamer(R__b);
+      R__b >> fXsection;
+//      (AliPythia::Instance())->Streamer(R__b);
+      R__b >> fPtHardMin;
+      R__b >> fPtHardMax;
+//      if (fDecayer) fDecayer->Streamer(R__b);
+   } else {
+      R__b.WriteVersion(AliGenPythia::IsA());
+      AliGenerator::Streamer(R__b);
+      R__b << (Int_t)fProcess;
+      R__b << (Int_t)fStrucFunc;
+      R__b << (Int_t)fForceDecay;
+      R__b << fEnergyCMS;
+      R__b << fKineBias;
+      R__b << fTrials;
+      fParentSelect.Streamer(R__b);
+      fChildSelect.Streamer(R__b);
+      R__b << fXsection;
+//      R__b << fPythia;
+      R__b << fPtHardMin;
+      R__b << fPtHardMax;
+      //     fDecayer->Streamer(R__b);
+   }
+}