]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliGenPythia.cxx
- Updates necessary for newio.
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
index 8743ba83b7c97a5c0d93d25eb9d17df2f35c6a70..549c6ca3066b9f1208020415b74b8f0d0c22e4fc 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.5  2003/04/04 08:12:56  morsch
-Boost for assymmetric systems added.
-
-Revision 1.4  2003/04/01 20:39:00  morsch
-Pycell configurable.
-
-Revision 1.3  2003/03/31 15:47:16  morsch
-Unhandled cases added.
-
-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)
@@ -247,9 +36,12 @@ Introduction of the Copyright and cvs Log
 #include "AliGenPythia.h"
 #include "AliGenPythiaEventHeader.h"
 #include "AliPythia.h"
+#include "AliPythiaRndm.h"
 #include "AliRun.h"
+#include "AliStack.h"
+#include "AliRunLoader.h"
 
- ClassImp(AliGenPythia)
+ClassImp(AliGenPythia)
 
 AliGenPythia::AliGenPythia()
                  :AliGenMC()
@@ -257,6 +49,7 @@ AliGenPythia::AliGenPythia()
 // Default Constructor
   fParticles = 0;
   fPythia    = 0;
+  fHeader = 0;
   fDecayer   = new AliDecayerPythia();
   SetEventListRange();
   SetJetPhiRange();
@@ -265,7 +58,11 @@ AliGenPythia::AliGenPythia()
   SetGammaPhiRange();
   SetGammaEtaRange();
   SetPtKick();
+  SetQuench();
+  
   fSetNuclei = kFALSE;
+  if (!AliPythiaRndm::GetPythiaRandom()) 
+    AliPythiaRndm::SetPythiaRandom(GetRandom());
 }
 
 AliGenPythia::AliGenPythia(Int_t npart)
@@ -286,11 +83,12 @@ AliGenPythia::AliGenPythia(Int_t npart)
     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);
-    fEventVertex.Set(3);
+    fHeader = 0;
     SetEventListRange();
     SetJetPhiRange();
     SetJetEtaRange();
@@ -298,6 +96,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetGammaPhiRange();
     SetGammaEtaRange();
     SetJetReconstructionMode();
+    SetQuench();
     SetPtKick();
     // Options determining what to keep in the stack (Heavy flavour generation)
     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
@@ -312,6 +111,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
 }
 
 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
+    :AliGenMC(Pythia)
 {
 // copy constructor
     Pythia.Copy(*this);
@@ -351,7 +151,7 @@ void AliGenPythia::Init()
 // Initialisation
     
     SetMC(AliPythia::Instance());
-    fPythia=(AliPythia*) fgMCEvGen;
+    fPythia=(AliPythia*) fMCEvGen;
 //
     fParentWeight=1./Float_t(fNpart);
 //
@@ -483,6 +283,7 @@ void AliGenPythia::Init()
 void AliGenPythia::Generate()
 {
 // Generate one event
+    
     fDecayer->ForceDecay();
 
     Float_t polar[3]   =   {0,0,0};
@@ -508,7 +309,19 @@ void AliGenPythia::Generate()
 //  event loop    
     while(1)
     {
+       if (fQuench) {
+           fPythia->SetMSTJ(1, 0);
+       }
+
        fPythia->Pyevnt();
+
+       if (fQuench) {
+           fPythia->Quench();
+           fPythia->SetMSTJ(1, 1);
+           fPythia->Pyexec();
+       }
+       
+       
        if (gAlice->GetEvNumber()>=fDebugEventFirst &&
            gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
        fTrials++;
@@ -528,9 +341,9 @@ void AliGenPythia::Generate()
        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];
@@ -647,7 +460,7 @@ void AliGenPythia::Generate()
 // 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;
@@ -674,11 +487,11 @@ void AliGenPythia::Generate()
                    Float_t tof   = kconv*iparticle->T();
                    Int_t ipa     = iparticle->GetFirstMother()-1;
                    Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
-                   SetTrack(fTrackIt*trackIt[i] ,
+                   PushTrack(fTrackIt*trackIt[i] ,
                                     iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
                    pParent[i] = nt;
                    KeepTrack(nt); 
-               } //  SetTrack loop
+               } //  PushTrack loop
            }
        } else {
            nc = GenerateMB();
@@ -719,7 +532,7 @@ void AliGenPythia::Generate()
     } // event loop
     SetHighWaterMark(nt);
 //  adjust weight due to kinematic selection
-    AdjustWeights();
+//    AdjustWeights();
 //  get cross-section
     fXsection=fPythia->GetPARI(1);
 }
@@ -770,7 +583,7 @@ Int_t  AliGenPythia::GenerateMB()
            origin[1] = fOrigin[1]+iparticle->Vy()/10.;
            origin[2] = fOrigin[2]+iparticle->Vz()/10.;
            Float_t tof=kconv*iparticle->T();
-           SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
+           PushTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
                     tof, kPPrimary, nt, 1., ks);
            KeepTrack(nt);
            pParent[i] = nt;
@@ -822,16 +635,17 @@ void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 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  
-    ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
+    ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
 //
 // Number of trials
-    ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials);
+    ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
 //
 // Event Vertex 
-    header->SetPrimaryVertex(fEventVertex);
+    fHeader->SetPrimaryVertex(fVertex);
 //
 // Jets that have triggered
     if (fProcess == kPyJets)
@@ -841,11 +655,11 @@ void AliGenPythia::MakeHeader()
        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]);
        }
     }
-    gAlice->SetGenEventHeader(header);
+    gAlice->SetGenEventHeader(fHeader);
 }
        
 
@@ -900,6 +714,7 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 {
 // Assignment operator
+    rhs.Copy(*this);
     return *this;
 }
 
@@ -909,12 +724,12 @@ void  AliGenPythia::LoadEvent()
 // Load event into Pythia Common Block
 //
  
-
-    Int_t npart = (Int_t) (gAlice->TreeK())->GetEntries(); 
+    AliRunLoader* rl = AliRunLoader::GetRunLoader();
+    Int_t npart = (rl->Stack())-> GetNprimary();
    (fPythia->GetPyjets())->N = npart;
 
     for (Int_t part = 0; part < npart; part++) {
-       TParticle *MPart = gAlice->Particle(part);
+       TParticle *MPart = (rl->Stack())->Particle(part);
        Int_t kf     = MPart->GetPdgCode();
        Int_t ks     = MPart->GetStatusCode();
        Float_t px = MPart->Px();
@@ -936,8 +751,7 @@ void  AliGenPythia::LoadEvent()
     }
 }
 
-void AliGenPythia::RecJetsUA1(Float_t eCellMin, Float_t eCellSeed, Float_t eMin, Float_t rMax, 
-                             Int_t& njets, Float_t jets [4][50])
+void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
 {
 //
 //  Calls the Pythia jet finding algorithm to find jets in the current event
@@ -1008,7 +822,7 @@ void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 
        if (
            eta > fEtaMinJet && eta < fEtaMaxJet && 
-           phi > fPhiMinJet && eta < fPhiMaxJet &&
+           phi > fPhiMinJet && phi < fPhiMaxJet &&
            et  > fEtMinJet  && et  < fEtMaxJet     
            ) 
        {