Use fVertex instead of fEventVertex

[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index 258b2e86d73571eb0b692aa70578c8219cffe5e0..5532a4fe50170462bf7254322f7450aeea2cf73d 100644 (file)
--- a/PYTHIA6/AliGenPythia.cxx
+++ b/PYTHIA6/AliGenPythia.cxx
@@ -13,209 +13,7 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 

-/*
-$Log$
-Revision 1.2  2003/03/27 09:28:38  morsch
-Charm production process configuration for pPb and pp added
-(A. Dainese, N. Carrer)
-
-Revision 1.1  2003/03/15 15:00:48  morsch
-Classed imported from EVGEN.
-
-Revision 1.69  2003/01/14 10:50:19  alibrary
-Cleanup of STEER coding conventions
-
-Revision 1.68  2002/12/11 09:16:16  morsch
-Use GetJets to fill header.
-
-Revision 1.67  2002/12/09 15:24:09  morsch
-Same trigger routine can use Pycell or Pyclus.
-
-Revision 1.66  2002/12/09 08:22:56  morsch
-UA1 jet finder (Pycell) for software triggering added.
-
-Revision 1.65  2002/11/19 08:57:10  morsch
-Configuration of pt-kick added.
-
-Revision 1.64  2002/11/15 00:43:06  morsch
-Changes for kPyJets
-- initial and final state g-radiation + pt-kick default
-- trigger based on parton clusters (using pyclus)
-- trigger jets are stored in header.
-
-Revision 1.63  2002/10/14 14:55:35  hristov
-Merging the VirtualMC branch to the main development branch (HEAD)
-
-Revision 1.52.4.4  2002/10/10 16:40:08  hristov
-Updating VirtualMC to v3-09-02
-
-Revision 1.62  2002/09/24 10:00:01  morsch
-CheckTrigger() corrected.
-
-Revision 1.61  2002/07/30 09:52:38  morsch
-Call SetGammaPhiRange() and SetGammaEtaRange() in the constructor.
-
-Revision 1.60  2002/07/19 14:49:03  morsch
-Typo corrected.
-
-Revision 1.59  2002/07/19 14:35:36  morsch
-Count total number of trials. Print mean Q, x1, x2.
-
-Revision 1.58  2002/07/17 10:04:09  morsch
-SetYHard method added.
-
-Revision 1.57  2002/05/22 13:22:53  morsch
-Process kPyMbNonDiffr added.
-
-Revision 1.56  2002/04/26 10:30:01  morsch
-Option kPyBeautyPbMNR added. (N. Carrer).
-
-Revision 1.55  2002/04/17 10:23:56  morsch
-Coding Rule violations corrected.
-
-Revision 1.54  2002/03/28 11:49:10  morsch
-Pass status code in SetTrack.
-
-Revision 1.53  2002/03/25 14:51:13  morsch
-New stack-fill and count options introduced (N. Carrer).
-
-Revision 1.51  2002/03/06 08:46:57  morsch
-- Loop until np-1
-- delete dyn. alloc. arrays (N. Carrer)
-
-Revision 1.50  2002/03/03 13:48:50  morsch
-Option  kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
-NLO calculations (Nicola Carrer).
-
-Revision 1.49  2002/02/08 16:50:50  morsch
-Add name and title in constructor.
-
-Revision 1.48  2001/12/20 11:44:28  morsch
-Add kinematic bias for direct gamma production.
-
-Revision 1.47  2001/12/19 14:45:00  morsch
-Store number of trials in header.
-
-Revision 1.46  2001/12/19 10:36:19  morsch
-Add possibility if jet kinematic biasing.
-
-Revision 1.45  2001/11/28 08:06:52  morsch
-Use fMaxLifeTime parameter.
-
-Revision 1.44  2001/11/27 13:13:07  morsch
-Maximum lifetime for long-lived particles to be put on the stack is parameter.
-It can be set via SetMaximumLifetime(..).
-
-Revision 1.43  2001/10/21 18:35:56  hristov
-Several pointers were set to zero in the default constructors to avoid memory management problems
-
-Revision 1.42  2001/10/15 08:21:55  morsch
-Vertex truncation settings moved to AliGenMC.
-
-Revision 1.41  2001/10/08 08:45:42  morsch
-Possibility of vertex cut added.
-
-Revision 1.40  2001/09/25 11:30:23  morsch
-Pass event vertex to header.
-
-Revision 1.39  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.38  2001/07/13 10:58:54  morsch
-- Some coded moved to AliGenMC
-- Improved handling of secondary vertices.
-
-Revision 1.37  2001/06/28 11:17:28  morsch
-SetEventListRange setter added. Events in specified range are listed for
-debugging. (Yuri Kharlov)
-
-Revision 1.36  2001/03/30 07:05:49  morsch
-Final print-out in finish run.
-Write parton system for jet-production (preliminary solution).
-
-Revision 1.35  2001/03/09 13:03:40  morsch
-Process_t and Struc_Func_t moved to AliPythia.h
-
-Revision 1.34  2001/02/14 15:50:40  hristov
-The last particle in event marked using SetHighWaterMark
-
-Revision 1.33  2001/01/30 09:23:12  hristov
-Streamers removed (R.Brun)
-
-Revision 1.32  2001/01/26 19:55:51  hristov
-Major upgrade of AliRoot code
-
-Revision 1.31  2001/01/17 10:54:31  hristov
-Better protection against FPE
-
-Revision 1.30  2000/12/18 08:55:35  morsch
-Make AliPythia dependent generartors work with new scheme of random number generation
-
-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
-
-Revision 1.16  2000/05/15 15:04:20  morsch
-The full event is written for fNtrack = -1
-Coding rule violations corrected.
-
-Revision 1.15  2000/04/26 10:14:24  morsch
-Particles array has one entry more than pythia particle list. Upper bound of
-particle loop changed to np-1 (R. Guernane, AM)
-
-Revision 1.14  2000/04/05 08:36:13  morsch
-Check status code of particles in Pythia event
-to avoid double counting as partonic state and final state particle.
-
-Revision 1.13  1999/11/09 07:38:48  fca
-Changes for compatibility with version 2.23 of ROOT
-
-Revision 1.12  1999/11/03 17:43:20  fca
-New version from G.Martinez & A.Morsch
-
-Revision 1.11  1999/09/29 09:24:14  fca
-Introduction of the Copyright and cvs Log
-*/
+/* $Id$ */

 
 //
 // Generator using the TPythia interface (via AliPythia)
 
 //
 // Generator using the TPythia interface (via AliPythia)


@@ -238,6 +36,7 @@ Introduction of the Copyright and cvs Log
 #include "AliGenPythia.h"
 #include "AliGenPythiaEventHeader.h"
 #include "AliPythia.h"
 #include "AliGenPythia.h"
 #include "AliGenPythiaEventHeader.h"
 #include "AliPythia.h"

+#include "AliPythiaRndm.h"

 #include "AliRun.h"
 
  ClassImp(AliGenPythia)
 #include "AliRun.h"
 
  ClassImp(AliGenPythia)


@@ -248,6 +47,7 @@ AliGenPythia::AliGenPythia()
 // Default Constructor
   fParticles = 0;
   fPythia    = 0;
 // Default Constructor
   fParticles = 0;
   fPythia    = 0;

+  fHeader = 0;

   fDecayer   = new AliDecayerPythia();
   SetEventListRange();
   SetJetPhiRange();
   fDecayer   = new AliDecayerPythia();
   SetEventListRange();
   SetJetPhiRange();


@@ -256,6 +56,9 @@ AliGenPythia::AliGenPythia()
   SetGammaPhiRange();
   SetGammaEtaRange();
   SetPtKick();
   SetGammaPhiRange();
   SetGammaEtaRange();
   SetPtKick();

+  fSetNuclei = kFALSE;
+  if (!AliPythiaRndm::GetPythiaRandom()) 
+    AliPythiaRndm::SetPythiaRandom(GetRandom());

 }
 
 AliGenPythia::AliGenPythia(Int_t npart)
 }
 
 AliGenPythia::AliGenPythia(Int_t npart)


@@ -268,8 +71,6 @@ AliGenPythia::AliGenPythia(Int_t npart)
     fName = "Pythia";
     fTitle= "Particle Generator using PYTHIA";
     fXsection  = 0.;
     fName = "Pythia";
     fTitle= "Particle Generator using PYTHIA";
     fXsection  = 0.;

-    fNucA1=0;
-    fNucA2=0;

     SetProcess();
     SetStrucFunc();
     SetForceDecay();
     SetProcess();
     SetStrucFunc();
     SetForceDecay();


@@ -278,11 +79,12 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetEnergyCMS();
     fDecayer = new AliDecayerPythia();
     // Set random number generator 
     SetEnergyCMS();
     fDecayer = new AliDecayerPythia();
     // Set random number generator 

-    sRandom=fRandom;
+    if (!AliPythiaRndm::GetPythiaRandom()) 
+      AliPythiaRndm::SetPythiaRandom(GetRandom());

     fFlavorSelect   = 0;
     // Produced particles  
     fParticles = new TClonesArray("TParticle",1000);
     fFlavorSelect   = 0;
     // Produced particles  
     fParticles = new TClonesArray("TParticle",1000);

-    fEventVertex.Set(3);
+    fHeader = 0;

     SetEventListRange();
     SetJetPhiRange();
     SetJetEtaRange();
     SetEventListRange();
     SetJetPhiRange();
     SetJetEtaRange();


@@ -298,6 +100,9 @@ AliGenPythia::AliGenPythia(Int_t npart)
     fFragmentation = kTRUE;
     // Default counting mode
     fCountMode = kCountAll;
     fFragmentation = kTRUE;
     // Default counting mode
     fCountMode = kCountAll;

+    // Pycel
+    SetPycellParameters();
+    fSetNuclei = kFALSE;

 }
 
 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 }
 
 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)


@@ -311,6 +116,21 @@ AliGenPythia::~AliGenPythia()
 // Destructor
 }
 
 // Destructor
 }
 

+void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
+                                      Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
+{
+// Set pycell parameters
+    fPycellEtaMax    =  etamax;
+    fPycellNEta      =  neta;
+    fPycellNPhi      =  nphi;
+    fPycellThreshold =  thresh;
+    fPycellEtSeed    =  etseed;
+    fPycellMinEtJet  =  minet;
+    fPycellMaxRadius =  r;
+}
+
+
+

 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 {
   // Set a range of event numbers, for which a table
 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 {
   // Set a range of event numbers, for which a table


@@ -339,7 +159,7 @@ void AliGenPythia::Init()
     fPythia->SetCKIN(7,fYHardMin);
     fPythia->SetCKIN(8,fYHardMax);
     
     fPythia->SetCKIN(7,fYHardMin);
     fPythia->SetCKIN(8,fYHardMax);
     

-    if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);  
+    if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);  

     // Fragmentation?
     if (fFragmentation) {
       fPythia->SetMSTP(111,1);
     // Fragmentation?
     if (fFragmentation) {
       fPythia->SetMSTP(111,1);


@@ -418,6 +238,23 @@ void AliGenPythia::Init()
        break;
     }
 //
        break;
     }
 //

+//
+//  JetFinder for Trigger
+//
+//  Configure detector (EMCAL like)
+//
+       fPythia->SetPARU(51, fPycellEtaMax);
+       fPythia->SetMSTU(51, fPycellNEta);
+       fPythia->SetMSTU(52, fPycellNPhi);
+//
+//  Configure Jet Finder
+//  
+       fPythia->SetPARU(58,  fPycellThreshold);
+       fPythia->SetPARU(52,  fPycellEtSeed);
+       fPythia->SetPARU(53,  fPycellMinEtJet);
+       fPythia->SetPARU(54,  fPycellMaxRadius);
+       fPythia->SetMSTU(54,  2);
+//

 //  This counts the total number of calls to Pyevnt() per run.
     fTrialsRun = 0;
     fQ         = 0.;
 //  This counts the total number of calls to Pyevnt() per run.
     fTrialsRun = 0;
     fQ         = 0.;


@@ -425,12 +262,22 @@ void AliGenPythia::Init()
     fX2        = 0.;    
     fNev       = 0 ;
 //    
     fX2        = 0.;    
     fNev       = 0 ;
 //    

+//
+//

     AliGenMC::Init();
     AliGenMC::Init();

+//
+//
+//  
+    if (fSetNuclei) {
+       fDyBoost = 0;
+       Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
+    }

 }
 
 void AliGenPythia::Generate()
 {
 // Generate one event
 }
 
 void AliGenPythia::Generate()
 {
 // Generate one event

+    

     fDecayer->ForceDecay();
 
     Float_t polar[3]   =   {0,0,0};
     fDecayer->ForceDecay();
 
     Float_t polar[3]   =   {0,0,0};


@@ -463,6 +310,7 @@ void AliGenPythia::Generate()
        
        fPythia->ImportParticles(fParticles,"All");
 
        
        fPythia->ImportParticles(fParticles,"All");
 

+       Boost();

 //
 //
 //
 //
 //
 //


@@ -475,9 +323,9 @@ void AliGenPythia::Generate()
        Float_t distz = iparticle->Vz()/10.;
        if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
 //
        Float_t distz = iparticle->Vz()/10.;
        if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
 //

-       fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
-       fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
-       fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
+       fVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
+       fVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
+       fVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);

 //
        Int_t* pParent   = new Int_t[np];
        Int_t* pSelected = new Int_t[np];
 //
        Int_t* pParent   = new Int_t[np];
        Int_t* pSelected = new Int_t[np];


@@ -594,7 +442,7 @@ void AliGenPythia::Generate()
 // Track final state particle
                if (ks == 1) trackIt[i] = 1;
 // Track semi-stable particles
 // Track final state particle
                if (ks == 1) trackIt[i] = 1;
 // Track semi-stable particles

-               if ((ks ==1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
+               if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;

 // Track particles selected by process if undecayed. 
                if (fForceDecay == kNoDecay) {
                    if (ParentSelected(kf)) trackIt[i] = 1;
 // Track particles selected by process if undecayed. 
                if (fForceDecay == kNoDecay) {
                    if (ParentSelected(kf)) trackIt[i] = 1;


@@ -759,24 +607,27 @@ void AliGenPythia::AdjustWeights()
 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 {
 // Treat protons as inside nuclei with mass numbers a1 and a2  
 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 {
 // Treat protons as inside nuclei with mass numbers a1 and a2  

-    fNucA1 = a1;
-    fNucA2 = a2;
+
+    fAProjectile = a1;
+    fATarget     = a2;
+    fSetNuclei   = kTRUE;

 }
 
 
 void AliGenPythia::MakeHeader()
 {
 // Builds the event header, to be called after each event
 }
 
 
 void AliGenPythia::MakeHeader()
 {
 // Builds the event header, to be called after each event

-    AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
+    if (fHeader) delete fHeader;
+    fHeader = new AliGenPythiaEventHeader("Pythia");

 //
 // Event type  
 //
 // Event type  

-    ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
+    ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));

 //
 // Number of trials
 //
 // Number of trials

-    ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials);
+    ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);

 //
 // Event Vertex 
 //
 // Event Vertex 

-    header->SetPrimaryVertex(fEventVertex);
+    fHeader->SetPrimaryVertex(fVertex);

 //
 // Jets that have triggered
     if (fProcess == kPyJets)
 //
 // Jets that have triggered
     if (fProcess == kPyJets)


@@ -786,11 +637,11 @@ void AliGenPythia::MakeHeader()
        GetJets(njet, ntrig, jets);
        
        for (Int_t i = 0; i < ntrig; i++) {
        GetJets(njet, ntrig, jets);
        
        for (Int_t i = 0; i < ntrig; i++) {

-           ((AliGenPythiaEventHeader*) header)->AddJet(jets[0][i], jets[1][i], jets[2][i], 
+           ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i], 

                                                        jets[3][i]);
        }
     }
                                                        jets[3][i]);
        }
     }

-    gAlice->SetGenEventHeader(header);
+    gAlice->SetGenEventHeader(fHeader);

 }
        
 
 }
        
 


@@ -888,19 +739,6 @@ void AliGenPythia::RecJetsUA1(Float_t eCellMin, Float_t eCellSeed, Float_t eMin,
 //  Calls the Pythia jet finding algorithm to find jets in the current event
 //
 //
 //  Calls the Pythia jet finding algorithm to find jets in the current event
 //
 //

-//  Configure detector (EMCAL like)
-//
-    fPythia->SetPARU(51,2.);
-    fPythia->SetMSTU(51,Int_t(96 * 2./0.7));
-    fPythia->SetMSTU(52,3 * 144);
-//
-//  Configure Jet Finder
-//  
-    fPythia->SetPARU(58, eCellMin);
-    fPythia->SetPARU(52, eCellSeed);
-    fPythia->SetPARU(53, eMin);
-    fPythia->SetPARU(54, rMax);
-    fPythia->SetMSTU(54, 2);

 //
 //  Save jets
     Int_t n     = fPythia->GetN();
 //
 //  Save jets
     Int_t n     = fPythia->GetN();


@@ -946,20 +784,7 @@ void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 //  Loading jets from common block
 //
     } else {
 //  Loading jets from common block
 //
     } else {

-//
-//  Configure detector (EMCAL like)
-//
-       fPythia->SetPARU(51,2.);
-       fPythia->SetMSTU(51,Int_t(96 * 2./0.7));
-       fPythia->SetMSTU(52,3 * 144);
-//
-//  Configure Jet Finder
-//  
-       fPythia->SetPARU(58,  0.0);
-       fPythia->SetPARU(52,  4.0);
-       fPythia->SetPARU(53, 10.0);
-       fPythia->SetPARU(54,  1.0);
-       fPythia->SetMSTU(54,  2);
+

 //
 //  Run Jet Finder
        fPythia->Pycell(nJets);
 //
 //  Run Jet Finder
        fPythia->Pycell(nJets);


@@ -979,7 +804,7 @@ void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 
        if (
            eta > fEtaMinJet && eta < fEtaMaxJet && 
 
        if (
            eta > fEtaMinJet && eta < fEtaMaxJet && 

-           phi > fPhiMinJet && eta < fPhiMaxJet &&
+           phi > fPhiMinJet && phi < fPhiMaxJet &&

            et  > fEtMinJet  && et  < fEtMaxJet     
            ) 
        {
            et  > fEtMinJet  && et  < fEtMaxJet     
            ) 
        {