This commit was generated by cvs2svn to compensate for changes in r7641,

which included commits to RCS files with non-trunk default branches.

diff --git a/DPMJET/ConfigDPM.C b/DPMJET/ConfigDPM.C
new file mode 100644 (file)
index 0000000..71acf3e
--- /dev/null
+++ b/DPMJET/ConfigDPM.C
@@ -0,0 +1,380 @@
+static Int_t    eventsPerRun = 100;
+enum PprGeo_t 
+{
+    kHoles, kNoHoles
+};
+static PprGeo_t geo = kHoles;
+
+void Config()
+{
+    // 7-DEC-2000 09:00
+    // Switch on Transition Radiation simulation. 6/12/00 18:00
+    // iZDC=1  7/12/00 09:00
+    // ThetaRange is (0., 180.). It was (0.28,179.72) 7/12/00 09:00
+    // Theta range given through pseudorapidity limits 22/6/2001
+
+    // Set Random Number seed
+    gRandom->SetSeed(123456);
+
+
+    // libraries required by geant321
+    gSystem->Load("libgeant321");
+
+    new TGeant3("C++ Interface to Geant3");
+
+    if (!gSystem->Getenv("CONFIG_FILE"))
+    {
+        TFile  *rootfile = new TFile("galice.root", "recreate");
+
+        rootfile->SetCompressionLevel(2);
+    }
+
+    TGeant3 *geant3 = (TGeant3 *) gMC;
+    AliIonPDGCodes *PDGcodes = new AliIonPDGCodes();
+    PDGcodes->AddParticlesToPdgDataBase();
+    geant3->DefineIons();
+
+    //
+    // Set External decayer
+    TVirtualMCDecayer *decayer = new AliDecayerPythia();
+
+    decayer->SetForceDecay(kAll);
+    decayer->Init();
+    gMC->SetExternalDecayer(decayer);
+    //
+    //
+    //=======================================================================
+    // ******* GEANT STEERING parameters FOR ALICE SIMULATION *******
+    geant3->SetTRIG(1);         //Number of events to be processed 
+    geant3->SetSWIT(4, 10);
+    geant3->SetDEBU(0, 0, 1);
+    //geant3->SetSWIT(2,2);
+    geant3->SetDCAY(1);
+    geant3->SetPAIR(1);
+    geant3->SetCOMP(1);
+    geant3->SetPHOT(1);
+    geant3->SetPFIS(0);
+    geant3->SetDRAY(0);
+    geant3->SetANNI(1);
+    geant3->SetBREM(1);
+    geant3->SetMUNU(1);
+    geant3->SetCKOV(1);
+    geant3->SetHADR(0);         //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3)
+    geant3->SetLOSS(2);
+    geant3->SetMULS(1);
+    geant3->SetRAYL(1);
+    geant3->SetAUTO(1);         //Select automatic STMIN etc... calc. (AUTO 1) or manual (AUTO 0)
+    geant3->SetABAN(0);         //Restore 3.16 behaviour for abandoned tracks
+    geant3->SetOPTI(2);         //Select optimisation level for GEANT geometry searches (0,1,2)
+    geant3->SetERAN(5.e-7);
+
+    Float_t cut = 1.e-3;        // 1MeV cut by default
+    Float_t tofmax = 1.e10;
+
+    //             GAM ELEC NHAD CHAD MUON EBREM MUHAB EDEL MUDEL MUPA TOFMAX
+    geant3->SetCUTS(cut, cut, cut, cut, cut, cut, cut, cut, cut, cut,
+                    tofmax);
+    //
+    //=======================================================================
+    // ************* STEERING parameters FOR ALICE SIMULATION **************
+    // --- Specify event type to be tracked through the ALICE setup
+    // --- All positions are in cm, angles in degrees, and P and E in GeV
+    // ####  AliGenDPMjet generation   ######################################
+    AliGenDPMjet *gener = new AliGenDPMjet(1);
+    gener->SetBeamEnergy(2700.);
+    gener->SetEnergyCMS(5400.);
+    gener->SetProjectile(208,82);
+    gener->SetTarget(208,82);
+    gener->SetImpactParameterRange(12.,16.); 
+    gener->SetTrackingFlag(1);
+    gener->Init();
+    // 
+    // Activate this line if you want the vertex smearing to happen
+    // track by track
+    //
+    //gener->SetVertexSmear(perTrack); 
+    // Field (L3 0.4 T)
+    AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., 1);
+    rootfile->cd();
+    gAlice->SetField(field);    
+
+
+    Int_t   iABSO  =  1;
+    Int_t   iDIPO  =  1;
+    Int_t   iFMD   =  0;
+    Int_t   iFRAME =  0;
+    Int_t   iHALL  =  0;
+    Int_t   iITS   =  0;
+    Int_t   iMAG   =  0;
+    Int_t   iMUON  =  0;
+    Int_t   iPHOS  =  0;
+    Int_t   iPIPE  =  1;
+    Int_t   iPMD   =  0;
+    Int_t   iRICH  =  0;
+    Int_t   iSHIL  =  1;
+    Int_t   iSTART =  0;
+    Int_t   iTOF   =  0;
+    Int_t   iTPC   =  0;
+    Int_t   iTRD   =  0;
+    Int_t   iZDC   =  1;
+    Int_t   iEMCAL =  0;
+    Int_t   iCRT   =  0;
+    Int_t   iVZERO =  0;
+
+    //=================== Alice BODY parameters =============================
+    AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
+
+    if (iMAG)
+    {
+        //=================== MAG parameters ============================
+        // --- Start with Magnet since detector layouts may be depending ---
+        // --- on the selected Magnet dimensions ---
+        AliMAG *MAG = new AliMAG("MAG", "Magnet");
+    }
+
+
+    if (iABSO)
+    {
+        //=================== ABSO parameters ============================
+        AliABSO *ABSO = new AliABSOv0("ABSO", "Muon Absorber");
+    }
+
+    if (iDIPO)
+    {
+        //=================== DIPO parameters ============================
+
+        AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2");
+    }
+
+    if (iHALL)
+    {
+        //=================== HALL parameters ============================
+
+        AliHALL *HALL = new AliHALL("HALL", "Alice Hall");
+    }
+
+
+    if (iFRAME)
+    {
+        //=================== FRAME parameters ============================
+
+        AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
+       if (geo == kHoles) {
+           FRAME->SetHoles(1);
+       } else {
+           FRAME->SetHoles(0);
+       }
+    }
+
+    if (iSHIL)
+    {
+        //=================== SHIL parameters ============================
+
+        AliSHIL *SHIL = new AliSHILv2("SHIL", "Shielding Version 2");
+    }
+
+
+    if (iPIPE)
+    {
+        //=================== PIPE parameters ============================
+
+        AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe");
+    }
+ 
+    if(iITS) {
+
+    //=================== ITS parameters ============================
+    //
+    // As the innermost detector in ALICE, the Inner Tracking System "impacts" on
+    // almost all other detectors. This involves the fact that the ITS geometry
+    // still has several options to be followed in parallel in order to determine
+    // the best set-up which minimizes the induced background. All the geometries
+    // available to date are described in the following. Read carefully the comments
+    // and use the default version (the only one uncommented) unless you are making
+    // comparisons and you know what you are doing. In this case just uncomment the
+    // ITS geometry you want to use and run Aliroot.
+    //
+    // Detailed geometries:         
+    //
+    //
+    //AliITS *ITS  = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services");
+    //
+    //AliITS *ITS  = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services");
+    //
+       AliITSvPPRasymm *ITS  = new AliITSvPPRasymm("ITS","New ITS PPR detailed version with asymmetric services");
+       ITS->SetMinorVersion(2);                                         // don't touch this parameter if you're not an ITS developer
+       ITS->SetReadDet(kFALSE);                                         // don't touch this parameter if you're not an ITS developer
+    //    ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymm2.det");  // don't touch this parameter if you're not an ITS developer
+       ITS->SetThicknessDet1(200.);   // detector thickness on layer 1 must be in the range [100,300]
+       ITS->SetThicknessDet2(200.);   // detector thickness on layer 2 must be in the range [100,300]
+       ITS->SetThicknessChip1(200.);  // chip thickness on layer 1 must be in the range [150,300]
+       ITS->SetThicknessChip2(200.);  // chip thickness on layer 2 must be in the range [150,300]
+       ITS->SetRails(0);            // 1 --> rails in ; 0 --> rails out
+       ITS->SetCoolingFluid(1);   // 1 --> water ; 0 --> freon
+       //
+    //AliITSvPPRsymm *ITS  = new AliITSvPPRsymm("ITS","New ITS PPR detailed version with symmetric services");
+    //ITS->SetMinorVersion(2);                                       // don't touch this parameter if you're not an ITS developer
+    //ITS->SetReadDet(kFALSE);                                       // don't touch this parameter if you're not an ITS developer
+    //ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRsymm2.det"); // don't touch this parameter if you're not an ITS developer
+    //ITS->SetThicknessDet1(200.);   // detector thickness on layer 1 must be in the range [100,300]
+    //ITS->SetThicknessDet2(200.);   // detector thickness on layer 2 must be in the range [100,300]
+    //ITS->SetThicknessChip1(200.);  // chip thickness on layer 1 must be in the range [150,300]
+    //ITS->SetThicknessChip2(200.);  // chip thickness on layer 2 must be in the range [150,300]
+    //ITS->SetRails(0);              // 1 --> rails in ; 0 --> rails out
+    //ITS->SetCoolingFluid(1);       // 1 --> water ; 0 --> freon
+    //
+    //
+    // Coarse geometries (warning: no hits are produced with these coarse geometries and they unuseful 
+    // for reconstruction !):
+    //                                                     
+    //
+    //AliITSvPPRcoarseasymm *ITS  = new AliITSvPPRcoarseasymm("ITS","New ITS PPR coarse version with asymmetric services");
+    //ITS->SetRails(0);                // 1 --> rails in ; 0 --> rails out
+    //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
+    //
+    //AliITS *ITS  = new AliITSvPPRcoarsesymm("ITS","New ITS PPR coarse version with symmetric services");
+    //ITS->SetRails(0);                // 1 --> rails in ; 0 --> rails out
+    //ITS->SetSupportMaterial(0);      // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
+    //                      
+    //
+    //
+    // Geant3 <-> EUCLID conversion
+    // ============================
+    //
+    // SetEUCLID is a flag to output (=1) or not to output (=0) both geometry and
+    // media to two ASCII files (called by default ITSgeometry.euc and
+    // ITSgeometry.tme) in a format understandable to the CAD system EUCLID.
+    // The default (=0) means that you dont want to use this facility.
+    //
+       ITS->SetEUCLID(0);  
+    }
+
+    if (iTPC)
+    {
+        //============================ TPC parameters ================================
+        // --- This allows the user to specify sectors for the SLOW (TPC geometry 2)
+        // --- Simulator. SecAL (SecAU) <0 means that ALL lower (upper)
+        // --- sectors are specified, any value other than that requires at least one 
+        // --- sector (lower or upper)to be specified!
+        // --- Reminder: sectors 1-24 are lower sectors (1-12 -> z>0, 13-24 -> z<0)
+        // ---           sectors 25-72 are the upper ones (25-48 -> z>0, 49-72 -> z<0)
+        // --- SecLows - number of lower sectors specified (up to 6)
+        // --- SecUps - number of upper sectors specified (up to 12)
+        // --- Sens - sensitive strips for the Slow Simulator !!!
+        // --- This does NOT work if all S or L-sectors are specified, i.e.
+        // --- if SecAL or SecAU < 0
+        //
+        //
+        //-----------------------------------------------------------------------------
+
+        //  gROOT->LoadMacro("SetTPCParam.C");
+        //  AliTPCParam *param = SetTPCParam();
+        AliTPC *TPC = new AliTPCv2("TPC", "Default");
+
+        // All sectors included 
+        TPC->SetSecAL(-1);
+        TPC->SetSecAU(-1);
+
+    }
+
+
+    if (iTOF) {
+       if (geo == kHoles) {
+        //=================== TOF parameters ============================
+           AliTOF *TOF = new AliTOFv2FHoles("TOF", "TOF with Holes");
+       } else {
+           AliTOF *TOF = new AliTOFv4T0("TOF", "normal TOF");
+       }
+    }
+
+
+    if (iRICH)
+    {
+        //=================== RICH parameters ===========================
+        AliRICH *RICH = new AliRICHv3("RICH", "normal RICH");
+
+    }
+
+
+    if (iZDC)
+    {
+        //=================== ZDC parameters ============================
+
+        AliZDC *ZDC = new AliZDCv2("ZDC", "normal ZDC");
+       ZDC->NoShower();
+    }
+
+    if (iTRD)
+    {
+        //=================== TRD parameters ============================
+
+        AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
+
+        // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 90% Xe + 10% CO2)
+        TRD->SetGasMix(1);
+       if (geo == kHoles) {
+           // With hole in front of PHOS
+           TRD->SetPHOShole();
+           // With hole in front of RICH
+           TRD->SetRICHhole();
+       }
+           // Switch on TR
+           AliTRDsim *TRDsim = TRD->CreateTR();
+    }
+
+    if (iFMD)
+    {
+        //=================== FMD parameters ============================
+       AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
+        FMD->SetRingsSi1(256);
+        FMD->SetRingsSi2(128);
+        FMD->SetSectorsSi1(20);
+        FMD->SetSectorsSi2(40);      
+   }
+
+    if (iMUON)
+    {
+        //=================== MUON parameters ===========================
+
+        AliMUON *MUON = new AliMUONv1("MUON", "default");
+    }
+    //=================== PHOS parameters ===========================
+
+    if (iPHOS)
+    {
+        AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");
+    }
+
+
+    if (iPMD)
+    {
+        //=================== PMD parameters ============================
+        AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
+    }
+
+    if (iSTART)
+    {
+        //=================== START parameters ============================
+        AliSTART *START = new AliSTARTv1("START", "START Detector");
+    }
+
+    if (iEMCAL)
+    {
+        //=================== EMCAL parameters ============================
+        AliEMCAL *EMCAL = new AliEMCALv1("EMCAL", "G56_2_55_19_104_14");
+    }
+
+     if (iCRT)
+    {
+        //=================== CRT parameters ============================
+        AliCRT *CRT = new AliCRTv0("CRT", "normal ACORDE");
+    }
+
+     if (iVZERO)
+    {
+        //=================== CRT parameters ============================
+        AliVZERO *VZERO = new AliVZEROv2("VZERO", "normal VZERO");
+    }
+
+}
+

diff --git a/DPMJET/doc/dpmjet3.manual b/DPMJET/doc/dpmjet3.manual
new file mode 100644 (file)
index 0000000..805651c
--- /dev/null
+++ b/DPMJET/doc/dpmjet3.manual
@@ -0,0 +1,943 @@
+ 
+     +-------------------------------------------------------------+
+     |                                                             |
+     |                                                             |
+     |                        DPMJET 3.0                           |
+     |                                                             |
+     |                                                             |
+     |         S. Roesler+), R. Engel#), J. Ranft*)                |
+     |                                                             |
+     |         +) CERN, TIS-RP                                     |
+     |            CH-1211 Geneva 23, Switzerland                   |
+     |            Email: Stefan.Roesler@cern.ch                    |
+     |                                                             |
+     |         #) University of Delaware, BRI                      |
+     |            Newark, DE 19716, USA                            |
+     |                                                             |
+     |         *) University of Siegen, Dept. of Physics           |
+     |            D-57068 Siegen, Germany                          |
+     |                                                             |
+     |                                                             |
+     |         http://home.cern.ch/sroesler/dpmjet3.html           |
+     |                                                             |
+     |                                                             |
+     |       Monte Carlo models used for event generation:         |
+     |          PHOJET 1.12, PYTHIA 6.115 and LEPTO 6.5.1          |
+     |                                                             |
+     +-------------------------------------------------------------+
+ 
+ 
+            
+                            INPUT OPTIONS
+                           _______________
+
+
+ The input of DPMJET consists of option cards. Option cards have all the
+ same structure and have to be given in fixed format except for the
+ section enclosed by PHOINPUT and ENDINPUT which can be given in free format.
+
+      CODEWD, (WHAT(I),I=1,6), SDUM   (default fixed format is
+                                       A10,6E10.0,A8)
+
+   where:
+
+   - CODEWD is the option keyword
+   - The WHAT-parameters are numerical data
+   - SDUM may contain literal data
+
+ The order of the input cards is free, with the exception of the START card
+ which initiates event generation and the BEAM card (see below).
+
+ Most definitions have some default values. If these are acceptable,
+ it is not compulsory that the corresponding option card appear
+ explicitly in the input sequence. Similarly, for most WHAT and/or SDUM
+ parameters a default value is applied if the corresponding field is left 
+ blank (or set = 0.0) in the input card.
+
+ Several option cards may appear more than once in the input sequence.
+ In most cases, each of such additional cards obviously adds more
+ definitions to those already given, provided they are different and not
+ contradictory. In case of conflict, the last given generally overrides
+ the previous one(s).
+
+ Any line starting with "*" is treated as a comment line.
+
+
+ 1) List of input options
+ ------------------------
+
+ *** general options ***
+
+ TITLE       title of run    
+ START       start of event generation
+ STOP        stop of the event generation                      
+ RNDMINIT    initialization of random number generator         
+
+ *** interacting particles ***
+
+ PROJPAR     projectile parameters
+ TARPAR      target parameters
+ EMULSION    definition of nuclear target emulsions or composite targets 
+
+ *** collision energy ***
+
+ ENERGY      interaction energy (per nucleon, lab)
+ MOMENTUM    interaction momentum (per nucleon, lab)
+ CMENERGY    interaction energy (nucleon-nucleon cms)
+ BEAM        definition of beam parameters
+
+ *** model for hadron / lepton / photon - nucleon interactions ***
+
+ MODEL       model to be used to treat nucleon-nucleon interactions        
+ PHOINPUT    start of PHOJET-specific input
+ ENDINPUT    end of PHOJET-specific input
+ HADRIN      HADRIN module                             
+ LEPTO-CUT   parameter CUT in LEPTO-common /LEPTOU/
+ LEPTO-LST   parameter LST in LEPTO-common /LEPTOU/
+ LEPTO-PARL  parameter PARL in LEPTO-common /LEPTOU/
+
+ *** Glauber formalism - cross sections ***
+
+ GLAUB-PAR   parameters in Glauber-formalism                    
+ GLAUB-INI   pre-initialization of profile function
+ FLUCTUAT    treatment of cross section fluctuations            
+ VDM-PAR1    parameters in gamma-nucleus cross section calculation     
+ VDM-PAR2    parameters in gamma-nucleus cross section calculation       
+ XS-TABLE    output of cross section table for requested interaction      
+ CENTRAL     biasing in impact parameter                       
+
+ *** parameters in DPM two-chain approximation ***
+
+ RECOMBIN    chain recombination                        
+ COMBIJET    chain fusion
+ CRONINPT    Cronin multiple scattering of partons
+ DIQUARKS    sea-diquark/antidiquark-pairs
+
+ *** hadronization and JETSET-parameters ***
+
+ LUND-MSTU   set parameter MSTU in JETSET-common /LUDAT1/          
+ LUND-PARJ   set parameter PARJ in JETSET-common /LUDAT1/           
+ LUND-PARU   set parameter PARJ in JETSET-common /LUDAT1/         
+ POPCORN     Popcorn-effect in fragmentation                   
+ PARDECAY    decay of Sigma0, Asigma0, pion^0
+
+ *** nuclear fragmentation ***
+
+ FERMI       Fermi momentum of nucleons
+ TAUFOR      formation time suppressed intranuclear cascade            
+ PAULI       treatment of Pauli's principle
+ COULOMB     treatment of Coulomb force
+ EVAP        evaporation / fragmentation module                             
+
+ *** output and checks ***
+
+ FRAME       Lorentz-frame in which final state is given in DTEVT1         
+ HISTOGRAM   activate different classes of histograms         
+ EMCCHECK    extended energy-momentum / quantum-number conservation check   
+
+ *** lepton tagger ***
+
+ L-TAG       lepton tagger (lepton-nucleus interactions only)            
+ L-ETAG      lepton tagger (lepton-nucleus interactions only)
+ ECMS-CUT    lepton tagger (lepton-nucleus interactions only)
+
+ *** for code development only ***
+
+ INTPT       intrinsic transverse momenta of partons
+ OUTLEVEL    output control switches                     
+ RESONANC    treatment of low mass chains                      
+ SEASU3      treatment of strange-quarks at chain ends             
+ XCUTS       thresholds for x-sampling                        
+
+
+ 2) Description of input options
+ -------------------------------
+            
+
+ *** general options ***
+
+______________________________________________________________________
+
+                control card:  codewd = TITLE
+
+   what (1..6), sdum   no meaning
+
+   Note:  The control-card following this card must consist of
+          a string of characters usually giving the title of
+          the run.
+______________________________________________________________________
+
+                control card:  codewd = START
+
+   what (1) =   number of events                default: 100.
+______________________________________________________________________
+
+                control card:  codewd = STOP
+
+                stop of the event generation
+
+   what (1..6)  no meaning
+______________________________________________________________________
+
+                control card:  RNDMINIT
+
+            initialization of random number generator
+
+   what (1..4)    values for initialization (= 1..168)
+   what (5..6), sdum    no meaning
+______________________________________________________________________
+
+
+ *** interacting particles ***
+
+______________________________________________________________________
+
+                control card:  codewd = PROJPAR
+
+              definition of projectile properties
+
+   what (1) =  (nucleus)   mass number 
+               (photon)    virtuality Q^2
+               (lepton, PHOJET)  
+                           maximum virtuality Q^2 of emitted photon
+               (otherwise) no meaning
+   what (2) =  (nucleus)   charge number
+               (otherwise) no meaning
+   what (3..6) no meaning
+   sdum        (hadrons,photons,leptons) particle code word
+
+   Note: In general, projectile nuclei are defined by what (1) and 
+         what (2). All other projectiles are defined by sdum.
+______________________________________________________________________
+
+                control card:  codewd = TARPAR
+
+                definition of target properties
+
+   what (1) =  (nucleus)   mass number
+               (otherwise) no meaning
+   what (2) =  (nucleus)   charge number
+               (otherwise) no meaning
+   what (3..6)   no meaning
+   sdum        (hadrons) particle code word
+
+   Note: In general, target nuclei are defined by what (1) and
+         what (2). Target hadrons are defined by sdum.
+______________________________________________________________________
+
+
+                control card:  codewd = EMULSION
+
+      definition of nuclear target emulsions or composite targets
+
+   what(1)      mass number of emulsion component
+   what(2)      charge of emulsion component
+   what(3)      fraction of events with this target
+   what(4,5,6)  as what(1,2,3) but for a further component
+                                           default: no emulsion
+   sdum         no meaning
+
+   Note: If this input-card is once used with valid parameters
+         TARPAR is obsolete. Not the absolute values of the fractions 
+         are important but only relative values.
+         This control card can be repeatedly used to define
+         emulsions / composite targets consisting of up to 10
+         elements.
+______________________________________________________________________
+
+
+ *** collision energy ***
+
+______________________________________________________________________
+
+                control card:  codewd = ENERGY
+
+                definition of laboratory energy
+
+   what (1) >  0:  what (1)  = total energy per nucleon (GeV)                
+            <  0:  |what(1)| = kinetic energy per nucleon (GeV)
+                                            default: 200 GeV
+               if |what(2)| > 0: min. total/kinetic energy per nucleon 
+                                 for variable energy runs
+   what (2)    max. energy per nucleon for variable energy runs
+            >  0:  what (2) = total energy per nucleon (GeV)
+            <  0:  |what(1)| = kinetic energy per nucleon (GeV)
+______________________________________________________________________
+
+                control card:  codewd = MOMENTUM
+
+                definition of laboratory momentum
+
+   what (1) =  momentum per nucleon (GeV/c) of projectile in Lab.
+                                            default: 200 GeV/c
+   what (2..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = CMENERGY
+
+   what (1) =  nucleon-nucleon c.m. energy           default: none
+   what (2..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = BEAM
+
+                definition of beam parameters
+              
+   what (1/2)  > 0 : energy per nucleon of beam 1/2 (GeV)           
+               < 0 : |what(1/2)| = energy per charge of 
+                     beam 1/2 (GeV)                      
+               (beam 1 is directed into positive z-direction)
+   what (3)    beam crossing angle, defined as 2x angle between
+               one beam and the z-axis (micro rad)            
+   what (4)    angle with x-axis defining the collision plane 
+   what (5..6), sdum   no meaning                            
+                                                             
+   Note: This card requires previously defined projectile and
+         target identities (PROJPAR, TARPAR) !              
+______________________________________________________________________
+
+
+ *** model for hadron / lepton / photon - nucleon interactions ***
+
+______________________________________________________________________
+
+                control card:  codewd = MODEL
+
+   Model used to describe nucleon(hadron,photon,lepton)-nucleon 
+                                                   interactions
+
+   what (1)    (only if sdum = LEPTO) 
+               variable INTER (see LEPTO-manual)
+                    = 1  gamma exchange                         
+                    = 2  W+-   exchange                         
+                    = 3  Z0    exchange                         
+                    = 4  gamma/Z0 exchange                      
+
+   sdum = DTUNUC    two-chain model as for versions 1.xx
+                    (nucleon/hadron-nucleon interactions only)
+        = PHOJET    multiple chains including minijets
+        = LEPTO     DIS
+                                                   default: PHOJET
+   what (2..6)      no meaning
+______________________________________________________________________
+
+                control card:  codewd = PHOINPUT
+
+   Start of PHOJET-specific input.
+   For details and a list of PHOJET input cards see the PHOJET-manual 
+   available at
+       http://lepton.bartol.udel.edu/~eng/phojet.html
+   Note:  This part of the input has to be closed by the ENDINPUT-card
+
+   what (1..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = ENDINPUT
+
+   End of PHOJET-specific input.
+   what (1..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = HADRIN
+
+                        HADRIN module
+
+   what (1) = 0. elastic/inelastic interactions with probab.
+                 as defined by cross-sections
+            = 1. inelastic interactions forced
+            = 2. elastic interactions forced               default: 1
+   what (2)      upper threshold in total energy (GeV) below
+                 which interactions are sampled by HADRIN
+                    default: steady transition btw. HADRIN and
+                             DPM in the range 4-7 GeV
+   what (3..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = LEPTO-CUT                   
+                                                                    
+           set parameter CUT in LEPTO-common /LEPTOU/               
+                                                                    
+   what (1) =  index in CUT-array                            
+   what (2) =  new value of CUT( int(what(1)) )            
+   what (3), what(4) and what (5), what(6) further        
+               parameter in the same way as what (1) and what (2) 
+                         default: default-LEPTO parameters  
+
+   Note: see LEPTO-manual.
+______________________________________________________________________
+                                                            
+                control card:  codewd = LEPTO-LST           
+                                                            
+           set parameter LST in LEPTO-common /LEPTOU/       
+                                                            
+   what (1) =  index in LST-array                           
+   what (2) =  new value of LST( int(what(1)) )             
+   what (3), what(4) and what (5), what(6) further          
+               parameter in the same way as what (1) and  what (2)
+                         default: default-LEPTO parameters       
+                                                                
+   Note: see LEPTO-manual.
+______________________________________________________________________
+                                                       
+                control card:  codewd = LEPTO-PARL     
+                                                       
+           set parameter PARL in LEPTO-common /LEPTOU/ 
+                                                       
+   what (1) =  index in PARL-array                
+   what (2) =  new value of PARL( int(what(1)) )  
+   what (3), what(4) and what (5), what(6) further     
+               parameter in the same way as what (1) and what (2) 
+                         default: default-LEPTO parameters    
+
+   Note: see LEPTO-manual.
+______________________________________________________________________
+
+
+ *** Glauber formalism - cross sections ***
+
+______________________________________________________________________
+
+                control card:  codewd = GLAUB-PAR
+
+                 parameters in Glauber-formalism
+
+   what (1)  number of nucleon configurations sampled in integration
+             over nuclear density                      default: 1000
+   what (2)  number of bins for integration over impact-parameter and
+             for profile-function calculation          default: 49
+   what (3)  = 1 calculation of tot., el. and qel. cross sections
+             otherwise calculation of production cross sections only
+                                                       default: 0
+   what (4)  = 1   read pre-calculated impact-parameter distribution
+                   from "sdum".glb for fixed projectile/target/energy 
+                   runs
+             = -1  dump calculated impact-parameter distribution
+                   into "sdum".glb for fixed or variable projectile/
+                   target/energy runs
+             = 100 read pre-calculated impact-parameter distribution
+                   from "sdum".glb for variable projectile/target/
+                   energy runs
+                                                       default: 0
+   what (5..6)   no meaning
+   sdum      if |what (4)| = 1 name of in/output-file (sdum.glb)
+______________________________________________________________________
+
+                control card:  codewd = GLAUB-INI
+
+             pre-initialization of profile function
+
+   what (1)      lower energy limit for initialization
+            > 0  Lab. frame
+            < 0  nucleon-nucleon cms
+   what (2)      upper energy limit for initialization
+            > 0  Lab. frame
+            < 0  nucleon-nucleon cms
+   what (3) > 0  # of equidistant lin. bins in E
+            < 0  # of equidistant log. bins in E
+   what (4)      maximum projectile mass number for which the Glauber
+                 data are initialized for each projectile mass number
+                 (if <= mass given with the PROJPAR-card)
+                                          default: 18
+   what (5)      steps in mass number starting from what (4)
+                 up to mass number defined with PROJPAR-card
+                 for which Glauber data are initialized
+                                          default: 5 
+   what (6)      no meaning
+   sdum          no meaning
+______________________________________________________________________
+
+                control card:  codewd = FLUCTUAT
+
+            Treatment of cross section fluctuations
+
+   what (1) = 1  cross section fluctuations treated    default: 0.
+   what (2..6), sdum   no meaning
+______________________________________________________________________
+
+
+                control card:  codewd = VDM-PAR1
+
+       parameters in gamma-nucleus cross section calculation
+
+   what (1) =  Lambda^2                       default: 2.
+   what (2)    lower limit in M^2 integration
+            =  1  (3m_pi)^2
+            =  2  (m_rho0)^2
+            =  3  (m_phi)^2                   default: 1
+   what (3)    upper limit in M^2 integration
+            =  1   s/2
+            =  2   s/4
+            =  3   s                          default: 3
+   what (4)    CKMT F_2 structure function
+            =  2212  proton
+            =  100   deuteron                 default: 2212
+   what (5)    calculation of gamma-nucleon xsections
+            =  1  according to CKMT-parametrization of F_2
+            =  2  integrating SIGVP over M^2
+            =  3  using SIGGA
+            =  4  PHOJET cross sections       default:  4
+
+   what (6), sdum    no meaning
+______________________________________________________________________
+
+                control card:  codewd = VDM-PAR2
+
+       parameters in gamma-nucleus cross section calculation
+
+   what (1) = 0 no suppression of shadowing by direct photon
+                processes
+            = 1 suppression ..                   default: 1
+   what (2) = 0 no suppression of shadowing by anomalous
+                component if photon-F_2
+            = 1 suppression ..                   default: 1
+   what (3) = 0 no suppression of shadowing by coherence
+                length of the photon
+            = 1 suppression ..                   default: 1
+   what (4) = 1 longitudinal polarized photons are taken into
+                account
+                eps*R*Q^2/M^2 = what(4)*Q^2/M^2  default: 0
+   what (5..6), sdum    no meaning
+______________________________________________________________________
+
+                control card:  codewd = XS-TABLE
+
+     output of cross section table for requested interaction
+               - particle production deactivated ! -
+
+   what (1)      lower energy limit for tabulation
+            > 0  Lab. frame
+            < 0  nucleon-nucleon cms
+   what (2)      upper energy limit for tabulation
+            > 0  Lab. frame
+            < 0  nucleon-nucleon cms
+   what (3) > 0  # of equidistant lin. bins in E
+            < 0  # of equidistant log. bins in E
+   what (4)      lower limit of particle virtuality (photons)
+   what (5)      upper limit of particle virtuality (photons)
+   what (6) > 0  # of equidistant lin. bins in Q^2
+            < 0  # of equidistant log. bins in Q^2
+______________________________________________________________________
+
+                control card:  codewd = CENTRAL
+
+                  Biasing in impact parameter
+
+   what (1) = 1.  central production
+                  (not recommended, has to be updated)
+   what (1) < 0  and > -100
+         what (2) = min. impact parameter                     
+         what (3) = max. impact parameter                       
+   what (1) < -99
+         what (2) = fraction of cross section
+   what (4..6), sdum   no meaning
+
+   Note: if what (1) = -1 : evaporation is suppressed
+         if what (1) < -1 : evaporation is allowed
+______________________________________________________________________
+
+
+ *** parameters in DPM two-chain approximation ***
+
+______________________________________________________________________
+
+                control card:  codewd = RECOMBIN
+
+                      Chain recombination
+         (recombine S-S and V-V chains to V-S chains)
+
+   what (1) = -1. recombination switched off    default: 1
+   what (2..6), sdum   no meaning
+
+   Note: Limited applicability for MODEL = PHOJET.
+______________________________________________________________________
+
+                control card:  codewd = COMBIJET
+
+                chain fusion (2 q-aq --> qq-aqaq)
+
+   what (1) = 1   fusion treated                     default: 0.
+   what (2)       minimum number of uncombined chains from
+                  single projectile or target nucleons
+                                                     default: 0.
+   what (3..6), sdum   no meaning
+
+   Note: Limited applicability for MODEL = PHOJET.
+______________________________________________________________________
+
+                control card:  codewd = CRONINPT
+
+     Cronin effect (multiple scattering of partons at chain ends)
+
+   what (1) = -1  Cronin effect not treated
+                  default (h+A): 1   (A+A): 0
+   what (2)       scattering parameter          default: 0.64
+   what (3..6), sdum   no meaning
+
+   Note: The Cronin-treatment should not be invoked for A+A int.
+______________________________________________________________________
+
+                control card:  codewd = DIQUARKS
+
+   what (1) = -1.  sea-diquark/antidiquark-pairs not treated
+                                                     default: -1.
+   what (2..6), sdum   no meaning
+
+   Note: Limited applicability for MODEL = PHOJET.
+______________________________________________________________________
+
+
+ *** hadronization and JETSET-parameters ***
+
+______________________________________________________________________
+
+                control card:  codewd = LUND-MSTU
+
+            parameter MSTU in JETSET-common /LUDAT1/
+
+   what (1) =  index according to LUND-common block
+   what (2) =  new value of MSTU( int(what(1)) )
+   what (3), what(4) and what (5), what(6) further
+   parameter used in the same way as what (1) and  what (2)
+
+   Note: The use of this card is not recommended. Some parameters
+         can presently not be changed with this card anyway. 
+         Limited applicability for MODEL = PHOJET.
+______________________________________________________________________
+
+                control card:  codewd = LUND-PARJ
+
+            parameter PARJ in JETSET-common /LUDAT1/
+
+   what (1) =  index according to LUND-common block
+   what (2) =  new value of PARJ( int(what(1)) )
+   what (3), what(4) and what (5), what(6) further
+   parameter used in the same way as what (1) and what (2)
+
+   Note: The use of this card is not recommended. Some parameters
+         can presently not be changed with this card anyway. 
+         Limited applicability for MODEL = PHOJET.
+______________________________________________________________________
+
+                control card:  codewd = LUND-PARU
+
+            parameter PARJ in JETSET-common /LUDAT1/
+
+   what (1) =  index according to LUND-common block
+   what (2) =  new value of PARU( int(what(1)) )
+   what (3), what(4) and what (5), what(6) further
+   parameter used in the same way as what (1) and what (2)
+
+   Note: The use of this card is not recommended. Some parameters
+         can presently not be changed with this card anyway. 
+         Limited applicability for MODEL = PHOJET.
+______________________________________________________________________
+
+                control card:  codewd = POPCORN
+
+               "Popcorn-effect" in fragmentation
+
+   what (1) < 0  Popcorn-effect switched off (MSTJ(12) = 1)
+            >=0  Popcorn-effect treated (PARJ(5) = what (1))
+                                                    default: 0.15
+   what (2..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = PARDECAY
+
+   what (1) = 1.  Sigma0/Asigma0 decay treated by JETSET
+            = 2.  pion^0 decay after intranucl. cascade
+                                        default: 0 (no such decays)
+   what (2..6), sdum   no meaning
+______________________________________________________________________
+
+
+ *** nuclear fragmentation ***
+
+______________________________________________________________________
+
+                control card:  codewd = FERMI
+
+   what (1) = -1 Fermi-motion of nucleons not treated    default: 1
+   what (2) =    scale factor for Fermi-momentum         default: 0.68
+   what (3..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = TAUFOR
+
+           formation time suppressed intranuclear cascade
+
+   what (1)      formation time (in fm/c)             default: 3.1 fm/c
+   what (2)      number of generations followed       default: 25
+   what (3) = 1. p_t-dependent formation zone
+            = 2. constant formation zone              default: 1
+   what (4)      modus of selection of nucleus where the
+                 cascade if followed first
+            = 1.  proj./target-nucleus with probab. 1/2
+            = 2.  nucleus with highest mass
+            = 3.  proj. nucleus if particle is moving in pos. z
+                  targ. nucleus if particle is moving in neg. z
+                                                      default: 1
+   what (5..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = PAULI
+
+   what (1) =  -1  Pauli's principle for secondary
+                   interactions not treated           default: 1
+   what (2..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = COULOMB
+
+   what (1) = -1. Coulomb-energy treatment switched off   default: 1
+   what (2..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = EVAP
+
+                 evaporation module of FLUKA
+
+   The following options and defaults apply only if the code is linked
+   to the FLUKA-library (see README file).
+
+   what (1) =< -1 ==> evaporation is switched off
+            >=  1 ==> evaporation is performed
+
+          what (1) = i1 + i2*10 + i3*100 + i4*10000
+                     (i1, i2, i3, i4 >= 0 )
+
+    i1 is the flag for selecting the T=0 level density option used
+       =  1: standard EVAP level densities with Cook pairing
+             energies
+       =  2: Z,N-dependent Gilbert & Cameron level densities
+                                                         (default)
+       =  3: Julich A-dependent level densities
+       =  4: Z,N-dependent Brancazio & Cameron level densities
+
+    i2 >= 1: high energy fission activated
+             (default high energy fission is activated)
+
+    i3 =  0: No energy dependence for level densities
+       =  1: Standard Ignyatuk (1975, 1st) energy dependence
+             for level densities (default)
+       =  2: Standard Ignyatuk (1975, 1st) energy dependence
+             for level densities with NOT used set of parameters
+       =  3: Standard Ignyatuk (1975, 1st) energy dependence
+             for level densities with NOT used set of parameters
+       =  4: Second   Ignyatuk (1975, 2nd) energy dependence
+             for level densities
+       =  5: Second   Ignyatuk (1975, 2nd) energy dependence
+             for level densities with fit 1 Iljinov & Mebel set of
+             parameters
+       =  6: Second   Ignyatuk (1975, 2nd) energy dependence
+             for level densities with fit 2 Iljinov & Mebel set of
+             parameters
+       =  7: Second   Ignyatuk (1975, 2nd) energy dependence
+             for level densities with fit 3 Iljinov & Mebel set of
+             parameters
+       =  8: Second   Ignyatuk (1975, 2nd) energy dependence
+             for level densities with fit 4 Iljinov & Mebel set of
+             parameters
+
+    i4 >= 1: Original Gilbert and Cameron pairing energies used
+             (default Cook's modified pairing energies)
+
+   what (2) = ig + 10 * if   (ig and if must have the same sign)
+
+    ig =< -1 ==> deexcitation gammas are not produced
+                 (if the evaporation step is not performed
+                  they are never produced)
+    if =< -1 ==> Fermi Break Up is not invoked
+                 (if the evaporation step is not performed
+                  it is never invoked)
+    The default is: deexcitation gamma are produced and
+                    Fermi break up is activated
+   what (3..6), sdum   no meaning
+______________________________________________________________________
+
+
+ *** output and checks ***
+
+______________________________________________________________________
+
+                control card:  codewd = FRAME
+
+           frame in which final state is given in DTEVT1
+
+   what (1) = 1  target rest frame (laboratory)
+            = 2  nucleon-nucleon cms                     default: 1
+______________________________________________________________________
+
+
+                control card:  codewd = HISTOGRAM
+
+            activate different classes of histograms
+
+                                 default: no histograms
+______________________________________________________________________
+
+                control card:  codewd = EMCCHECK
+
+     extended energy-momentum / quantum-number conservation check
+
+   what (1) = -1   extended check not performed        default: -1.
+   what (2..6), sdum   no meaning
+______________________________________________________________________
+
+
+ *** lepton tagger ***
+
+______________________________________________________________________
+
+                control card:  codewd = L-TAG
+        (lepton-nucleus interactions with MODEL=PHOJET only)
+
+                         lepton tagger:
+    definition of kinematic cuts for radiated photon and
+    outgoing lepton detection in lepton-nucleus interactions
+
+   what (1) = y_min
+   what (2) = y_max
+   what (3) = Q^2_min
+   what (4) = Q^2_max
+   what (5) = theta_min  (Lab)
+   what (6) = theta_max  (Lab)                        default: no cuts
+   sdum    no meaning
+______________________________________________________________________
+
+
+                control card:  codewd = L-ETAG
+        (lepton-nucleus interactions with MODEL=PHOJET only)
+
+                         lepton tagger:
+   what (1) = min. outgoing lepton energy  (in Lab)
+   what (2) = min. photon energy           (in Lab)
+   what (3) = max. photon energy           (in Lab)   default: no cuts
+   what (2..6), sdum    no meaning 
+______________________________________________________________________
+
+                control card:  codewd = ECMS-CUT
+        (lepton-nucleus interactions with MODEL=PHOJET only)
+
+   what (1) = min. c.m. energy to be sampled
+   what (2) = max. c.m. energy to be sampled
+   what (3) = min x_Bj         to be sampled          default: no cuts
+   what (3..6), sdum    no meaning
+______________________________________________________________________
+
+
+ *** for code development only ***
+
+______________________________________________________________________
+
+                control card:  codewd = INTPT
+
+   what (1) = -1   intrinsic transverse momenta of partons
+                   not treated                default: 1
+   what (2..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = OUTLEVEL
+
+                     output control switches
+
+   what (1) =  internal rejection informations             default: 0
+   what (2) =  energy-momentum conservation check output   default: 0
+   what (3..6) internal warning messages                   default: 0
+______________________________________________________________________
+
+                control card:  codewd = RESONANC
+
+                  treatment of low mass chains
+
+   what (1) = -1 low chain masses are not corrected for resonance
+                 masses                                    default: 1.
+   what (2) = -1 massless partons                 default: 1. (massive)
+   what (3) = -1 chain-system containing chain of too small
+                 mass is rejected (note: this does not fully
+                 apply to S-S chains)                      default: 0.
+   what (4..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = SEASU3
+
+           Treatment of strange-quarks at chain ends
+
+   what (1)   (SEASQ)  strange-quark suppression factor
+              iflav = 1.+rndm*(2.+SEASQ)                   default: 1.
+   what (2..6), sdum   no meaning
+______________________________________________________________________
+
+                control card:  codewd = XCUTS
+
+                  thresholds for x-sampling
+
+   what (1)    defines lower threshold for val.-q x-value (CVQ)
+                                                       default: 1.
+   what (2)    defines lower threshold for val.-qq x-value (CDQ)
+                                                       default: 2.
+   what (3)    defines lower threshold for sea-q x-value (CSEA)
+                                                       default: 0.2
+   what (4)    sea-q x-values in S-S chains (SSMIMA)   default: 0.14
+   what (5)    not used                                default: 2.
+   what (6), sdum   no meaning
+
+   Note: Lower thresholds (what(1..3)) are def. as x_thr=CXXX/ECM
+______________________________________________________________________
+
+
+
+                     EVENT HISTORY - COMMON /DTEVT1/
+                    _________________________________
+
+
+ 1) Common block DTEVT1
+ ----------------------
+
+ The complete history of each event can be found in COMMON /DTEVT1/ 
+ and COMMON /EXTEVT/.
+
+ COMMON /DTEVT1/ :
+         NHKK         number of entries in common block
+         NEVHKK       number of the event
+         ISTHKK(i)    status code for entry i
+         IDHKK(i)     identifier for the entry
+                      (for particles: identifier according to the PDG scheme)
+         JMOHKK(1,i)  pointer to the entry of the first mother of entry i
+         JMOHKK(2,i)  pointer to the entry of the second mother of entry i
+         JDAHKK(1,i)  pointer to the entry of the first daughter of entry i
+         JDAHKK(2,i)  pointer to the entry of the second daughter of entry i
+         PHKK(1..3,i) 3-momentum
+         PHKK(4,i)    energy
+         PHKK(5,i)    mass
+         VHKK / WHKK  spatial position of particle in target / projectile
+                      rest frame
+
+ 2) Final state particles
+ ------------------------
+
+ The final state particles from the actual event (number NEVHKK)
+ can be found in DTEVT1 and identified by their status:
+
+    ISTHKK(i) = 1    final state particle produced in
+                     photon-/hadron-/nucleon-nucleon collisions or
+                     in intranuclear cascade processes
+               -1    nucleons, deuterons, H-3, He-3, He-4 evaporated
+                     from excited nucleus, fragmentation and fission
+                     products (A > 4) and photons produced in nuclear 
+                     deexcitation processes
+               1001  residual nucleus (ground state)
+
+ The types of these particles/nuclei are given in IDHKK as follows
+
+    all final state particles except nuclei :
+      IDHKK(i)=particle identifier according to PDG numbering scheme
+    nuclei (A > 1: evaporation / fragmentation / fission products, and 
+            residual nuclei) :
+      IDHKK(i)=80000, IDRES(i)=mass number, IDXRES(i)=charge number
+
+ The 4-momenta and masses can be found in PHKK (target nucleus rest frame
+ unless defined by the FRAME-card):
+                  PHKK(1..3,i) 3-momentum (p_x,p_y,p_z)
+                  PHKK(4,i)    energy
+                  PHKK(5,i)    mass

diff --git a/DPMJET/dpmjet3.0-4.f b/DPMJET/dpmjet3.0-4.f
new file mode 100644 (file)
index 0000000..382eb94
--- /dev/null
+++ b/DPMJET/dpmjet3.0-4.f
@@ -0,0 +1,37139 @@
+*
+*    +-------------------------------------------------------------+
+*    |                                                             |
+*    |                                                             |
+*    |                        DPMJET 3.0                           |
+*    |                                                             |
+*    |                                                             |
+*    |         S. Roesler+), R. Engel#), J. Ranft*)                |
+*    |                                                             |
+*    |         +) CERN, TIS-RP                                     |
+*    |            CH-1211 Geneva 23, Switzerland                   |
+*    |            Email: Stefan.Roesler@cern.ch                    |
+*    |                                                             |
+*    |         #) University of Delaware, BRI                      |
+*    |            Newark, DE 19716, USA                            |
+*    |                                                             |
+*    |         *) University of Siegen, Dept. of Physics           |
+*    |            D-57068 Siegen, Germany                          |
+*    |                                                             |
+*    |                                                             |
+*    |       http://home.cern.ch/sroesler/dpmjet3.html             |
+*    |                                                             |
+*    |                                                             |
+*    |       Monte Carlo models used for event generation:         |
+*    |          PHOJET 1.12, JETSET 7.4 and LEPTO 6.5.1            |
+*    |                                                             |
+*    +-------------------------------------------------------------+
+*
+*
+*===init===============================================================*
+*
+CDECK  ID>, DT_INIT
+      SUBROUTINE DT_INIT(NCASES,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
+     &                                             IDP,IGLAU)
+
+************************************************************************
+* Initialization of event generation                                   *
+* This version dated  7.4.98  is written by S. Roesler.                *
+************************************************************************
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      PARAMETER ( LINP = 5 ,
+     &            LOUT = 6 ,
+     &            LDAT = 9 )
+
+      PARAMETER (ZERO=0.0D0,ONE=1.0D0)
+
+* particle properties (BAMJET index convention)
+      CHARACTER*8  ANAME
+      COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
+     &                IICH(210),IIBAR(210),K1(210),K2(210)
+* names of hadrons used in input-cards
+      CHARACTER*8 BTYPE
+      COMMON /DTPAIN/ BTYPE(30)
+
+      INCLUDE './flukapro/(DIMPAR)'
+      INCLUDE './flukapro/(PAREVT)'
+      INCLUDE './flukapro/(EVAPAR)'
+      INCLUDE './flukapro/(FRBKCM)'
+
+      PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
+
+* emulsion treatment
+      COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
+     &                NCOMPO,IEMUL
+* Glauber formalism: parameters
+      COMMON /DTGLAM/ RASH(NCOMPX),RBSH(NCOMPX),
+     &                BMAX(NCOMPX),BSTEP(NCOMPX),
+     &                SIGSH,ROSH,GSH,BSITE(0:NEB,NQB,NCOMPX,KSITEB),
+     &                NSITEB,NSTATB
+* Glauber formalism: cross sections
+      COMMON /DTGLXS/ ECMNN(NEB),Q2G(NQB),ECMNOW,Q2,
+     &                XSTOT(NEB,NQB,NCOMPX),XSELA(NEB,NQB,NCOMPX),
+     &                XSQEP(NEB,NQB,NCOMPX),XSQET(NEB,NQB,NCOMPX),
+     &                XSQE2(NEB,NQB,NCOMPX),XSPRO(NEB,NQB,NCOMPX),
+     &                XSDEL(NEB,NQB,NCOMPX),XSDQE(NEB,NQB,NCOMPX),
+     &                XETOT(NEB,NQB,NCOMPX),XEELA(NEB,NQB,NCOMPX),
+     &                XEQEP(NEB,NQB,NCOMPX),XEQET(NEB,NQB,NCOMPX),
+     &                XEQE2(NEB,NQB,NCOMPX),XEPRO(NEB,NQB,NCOMPX),
+     &                XEDEL(NEB,NQB,NCOMPX),XEDQE(NEB,NQB,NCOMPX),
+     &                BSLOPE,NEBINI,NQBINI
+* interface HADRIN-DPM
+      COMMON /HNTHRE/ EHADTH,EHADLO,EHADHI,INTHAD,IDXTA
+* central particle production, impact parameter biasing
+      COMMON /DTIMPA/ BIMIN,BIMAX,XSFRAC,ICENTR
+* parameter for intranuclear cascade
+      LOGICAL LPAULI
+      COMMON /DTFOTI/ TAUFOR,KTAUGE,ITAUVE,INCMOD,LPAULI
+* various options for treatment of partons (DTUNUC 1.x)
+* (chain recombination, Cronin,..)
+      LOGICAL LCO2CR,LINTPT
+      COMMON /DTCHAI/ SEASQ,CRONCO,CUTOF,MKCRON,ISICHA,IRECOM,
+     &                LCO2CR,LINTPT
+* threshold values for x-sampling (DTUNUC 1.x)
+      COMMON /DTXCUT/ XSEACU,UNON,UNOM,UNOSEA,CVQ,CDQ,CSEA,SSMIMA,
+     &                SSMIMQ,VVMTHR
+* flags for input different options
+      LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
+      COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
+     &                LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
+* nuclear potential
+      LOGICAL LFERMI
+      COMMON /DTNPOT/ PFERMP(2),PFERMN(2),FERMOD,
+     &                EBINDP(2),EBINDN(2),EPOT(2,210),
+     &                ETACOU(2),ICOUL,LFERMI
+* n-n cross section fluctuations
+      PARAMETER (NBINS = 1000)
+      COMMON /DTXSFL/ FLUIXX(NBINS),IFLUCT
+* flags for particle decays
+      COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
+     &                IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
+     &                NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
+* diquark-breaking mechanism
+      COMMON /DTDIQB/ DBRKR(3,8),DBRKA(3,8),CHAM1,CHAM3,CHAB1,CHAB3
+* nucleon-nucleon event-generator
+      CHARACTER*8 CMODEL
+      LOGICAL LPHOIN
+      COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
+* properties of interacting particles
+      COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
+* properties of photon/lepton projectiles
+      COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
+* flags for diffractive interactions (DTUNUC 1.x)
+      COMMON /DTFLG3/ ISINGD,IDOUBD,IFLAGD,IDIFF
+* parameters for hA-diffraction
+      COMMON /DTDIHA/ DIBETA,DIALPH
+* Lorentz-parameters of the current interaction
+      COMMON /DTLTRA/ GACMS(2),BGCMS(2),GALAB,BGLAB,BLAB,
+     &                UMO,PPCM,EPROJ,PPROJ
+* kinematical cuts for lepton-nucleus interactions
+      COMMON /DTLCUT/ ECMIN,ECMAX,XBJMIN,ELMIN,EGMIN,EGMAX,YMIN,YMAX,
+     &                Q2MIN,Q2MAX,THMIN,THMAX,Q2LI,Q2HI,ECMLI,ECMHI
+* VDM parameter for photon-nucleus interactions
+      COMMON /DTVDMP/ RL2,EPSPOL,INTRGE(2),IDPDF,MODEGA,ISHAD(3)
+* Glauber formalism: flags and parameters for statistics
+      LOGICAL LPROD
+      CHARACTER*8 CGLB
+      COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
+* cuts for variable energy runs
+      COMMON /DTVARE/ VARELO,VAREHI,VARCLO,VARCHI
+* flags for activated histograms
+      COMMON /DTHIS3/ IHISPP(50),IHISXS(50),IXSTBL
+
+      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
+
+      COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
+
+* LEPTO
+**LUND single / double precision
+      REAL CUT,PARL,TMPX,TMPY,TMPW2,TMPQ2,TMPU
+      COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),
+     &                TMPX,TMPY,TMPW2,TMPQ2,TMPU
+* LEPTO
+      REAL RPPN
+      COMMON /LEPTOI/ RPPN,LEPIN,INTER
+* steering flags for qel neutrino scattering modules
+      COMMON /QNEUTO/ DSIGSU,DSIGMC,NDSIG,NEUTYP,NEUDEC
+* event flag
+      COMMON /DTEVNO/ NEVENT,ICASCA
+
+      INTEGER PYCOMP
+
+C     DIMENSION XPARA(5)
+      DIMENSION XDUMB(40),IPRANG(5)
+
+      PARAMETER (MXCARD=58)
+      CHARACTER*78 CLINE,CTITLE
+      CHARACTER*60 CWHAT
+      CHARACTER*8  BLANK,SDUM
+      CHARACTER*10 CODE,CODEWD
+      CHARACTER*72 HEADER
+      LOGICAL LSTART,LEINP,LXSTAB
+      DIMENSION WHAT(6),CODE(MXCARD)
+      DATA CODE/
+     &   'TITLE     ','PROJPAR   ','TARPAR    ','ENERGY    ',
+     &   'MOMENTUM  ','CMENERGY  ','EMULSION  ','FERMI     ',
+     &   'TAUFOR    ','PAULI     ','COULOMB   ','HADRIN    ',
+     &   'EVAP      ','EMCCHECK  ','MODEL     ','PHOINPUT  ',
+     &   'GLAUBERI  ','FLUCTUAT  ','CENTRAL   ','RECOMBIN  ',
+     &   'COMBIJET  ','XCUTS     ','INTPT     ','CRONINPT  ',
+     &   'SEADISTR  ','SEASU3    ','DIQUARKS  ','RESONANC  ',
+     &   'DIFFRACT  ','SINGLECH  ','NOFRAGME  ','HADRONIZE ',
+     &   'POPCORN   ','PARDECAY  ','BEAM      ','LUND-MSTU ',
+     &   'LUND-MSTJ ','LUND-MDCY ','LUND-PARJ ','LUND-PARU ',
+     &   'OUTLEVEL  ','FRAME     ','L-TAG     ','L-ETAG    ',
+     &   'ECMS-CUT  ','VDM-PAR1  ','HISTOGRAM ','XS-TABLE  ',
+     &   'GLAUB-PAR ','GLAUB-INI ','VDM-PAR2  ','XS-QELPRO ',
+     &   'RNDMINIT  ','LEPTO-CUT ','LEPTO-LST ','LEPTO-PARL',
+     &   'START     ','STOP      '/
+      DATA BLANK /'        '/
+
+      DATA LSTART,LXSTAB,IFIRST /.TRUE.,.FALSE.,1/
+      DATA CMEOLD /0.0D0/
+
+* --- Added by Chiara
+      
+      CHARACTER*100  ALIROOT  
+      CHARACTER*100 FILNAM
+      INTEGER*4 LNROOT
+      LOGICAL EXISTS
+      ALIROOT=' '
+
+*---------------------------------------------------------------------
+* at the first call of INIT: initialize event generation
+      EPNSAV = EPN
+      IF (LSTART) THEN
+         CALL DT_TITLE
+*   initialization and test of the random number generator
+         IF (ITRSPT.NE.1) THEN
+
+            CALL FL48UT (ISRM48,ISEED1,ISEED2)
+            CALL FL48IN (54217137,ISEED1,ISEED2)
+
+         ENDIF
+*   initialization of BAMJET, DECAY and HADRIN
+         CALL DT_DDATAR
+         CALL DT_DHADDE
+         CALL DT_DCHANT
+         CALL DT_DCHANH
+*   set default values for input variables
+         CALL DT_DEFAUL(EPN,PPN)
+         IGLAU  = 0
+         IXSQEL = 0
+*   flag for collision energy input
+         LEINP  = .FALSE.
+         LSTART = .FALSE.
+      ENDIF
+
+*---------------------------------------------------------------------
+   10 CONTINUE
+
+* bypass reading input cards (e.g. for use with Fluka)
+*  in this case Epn is expected to carry the beam momentum
+      IF (NCASES.EQ.-1) THEN
+         IP      = NPMASS
+         IPZ     = NPCHAR
+         PPN     = EPNSAV
+         EPN     = ZERO
+         CMENER  = ZERO
+         LEINP   = .TRUE.
+         MKCRON  = 0
+         WHAT(1) = 1
+         WHAT(2) = 0
+         CODEWD  = 'START     '
+         GOTO 900
+      ENDIF
+
+* read control card from input-unit LINP
+C      READ(LINP,'(A78)',END=9999) CLINE
+* ###   Read control card from specified file 
+* ### Changed by Chiara (original version LINP=5)
+*      OPEN(UNIT=7,
+*     + FILE='/home/oppedisa/AliRoot/new/DPMJET/inp/PbPbLHC.inp',
+*     + STATUS='OLD')
+
+      CALL GETENVF('ALICE_ROOT',ALIROOT)
+      LNROOT = LNBLNK(ALIROOT)
+
+      FILNAM=ALIROOT(1:LNROOT)//'/DPMJET/inp/PbPbLHC.inp'
+      OPEN(UNIT=7,FILE=FILNAM,STATUS='OLD')
+
+
+      READ(7,'(A78)',END=9999) CLINE
+
+      IF (CLINE(1:1).EQ.'*') THEN
+* comment-line
+C         WRITE(LOUT,'(A78)') CLINE
+         GOTO 10
+      ENDIF
+C     READ(CLINE,1000,END=9999) CODEWD,(WHAT(I),I=1,6),SDUM
+C1000 FORMAT(A10,6E10.0,A8)
+      DO 1008 I=1,6
+         WHAT(I) = ZERO
+ 1008 CONTINUE
+      READ(CLINE,1006,END=9999) CODEWD,CWHAT,SDUM
+ 1006 FORMAT(A10,A60,A8)
+      READ(CWHAT,*,END=1007) (WHAT(I),I=1,6)
+ 1007 CONTINUE
+      WRITE(LOUT,1001) CODEWD,(WHAT(I),I=1,6),SDUM
+ 1001 FORMAT(A10,6G10.3,A8)
+
+  900 CONTINUE
+
+* check for valid control card and get card index
+      ICW = 0
+      DO 11 I=1,MXCARD
+         IF (CODEWD.EQ.CODE(I)) ICW = I
+   11 CONTINUE
+      IF (ICW.EQ.0) THEN
+         WRITE(LOUT,1002) CODEWD
+ 1002    FORMAT(/,1X,'---> ',A10,': invalid control-card !',/)
+         GOTO 10
+      ENDIF
+
+      GOTO(
+*------------------------------------------------------------
+*       TITLE   ,  PROJPAR ,  TARPAR  ,  ENERGY  ,  MOMENTUM,
+     &  100     ,  110     ,  120     ,  130     ,  140     ,
+*
+*------------------------------------------------------------
+*       CMENERGY,  EMULSION,  FERMI   ,  TAUFOR  ,  PAULI   ,
+     &  150     ,  160     ,  170     ,  180     ,  190     ,
+*
+*------------------------------------------------------------
+*       COULOMB ,  HADRIN  ,  EVAP    ,  EMCCHECK,  MODEL   ,
+     &  200     ,  210     ,  220     ,  230     ,  240     ,
+*
+*------------------------------------------------------------
+*       PHOINPUT,  GLAUBERI,  FLUCTUAT,  CENTRAL ,  RECOMBIN,
+     &  250     ,  260     ,  270     ,  280     ,  290     ,
+*
+*------------------------------------------------------------
+*       COMBIJET,  XCUTS   ,  INTPT   ,  CRONINPT,  SEADISTR,
+     &  300     ,  310     ,  320     ,  330     ,  340     ,
+*
+*------------------------------------------------------------
+*       SEASU3  ,  DIQUARKS,  RESONANC,  DIFFRACT,  SINGLECH,
+     &  350     ,  360     ,  370     ,  380     ,  390     ,
+*
+*------------------------------------------------------------
+*       NOFRAGME, HADRONIZE,  POPCORN ,  PARDECAY,  BEAM    ,
+     &  400     ,  410     ,  420     ,  430     ,  440     ,
+*
+*------------------------------------------------------------
+*      LUND-MSTU, LUND-MSTJ, LUND-MDCY, LUND-PARJ, LUND-PARU,
+     &  450     ,  451     ,  452     ,  460     ,  470     ,
+*
+*------------------------------------------------------------
+*       OUTLEVEL,  FRAME   , L-TAG    ,  L-ETAG  ,  ECMS-CUT,
+     &  480     ,  490     ,  500     ,  510     ,  520     ,
+*
+*------------------------------------------------------------
+*       VDM-PAR1, HISTOGRAM, XS-TABLE , GLAUB-PAR, GLAUB-INI,
+     &  530     ,  540     ,  550     ,  560     ,  565     ,
+*
+*------------------------------------------------------------
+*               ,          ,  VDM-PAR2, XS-QELPRO, RNDMINIT ,
+     &                        570     ,  580     ,  590     ,
+*
+*------------------------------------------------------------
+*      LEPTO-CUT, LEPTO-LST,LEPTO-PARL,  START   ,  STOP    )
+     &  600     ,  610     ,  620     ,  630     ,  640     ) , ICW
+*
+*------------------------------------------------------------
+
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = TITLE                       *
+*                                                                   *
+*       what (1..6), sdum   no meaning                              *
+*                                                                   *
+*       Note:  The control-card following this must consist of      *
+*              a string of characters usually giving the title of   *
+*              the run.                                             *
+*                                                                   *
+*********************************************************************
+
+  100 CONTINUE
+C      READ(LINP,'(A78)') CTITLE
+* ###   Read control card from specified file 
+* ### Changed by Chiara (original version LINP=5)
+      READ(7,'(A78)') CTITLE
+
+      WRITE(LOUT,'(//,5X,A78,//)') CTITLE
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = PROJPAR                     *
+*                                                                   *
+*       what (1) =  mass number of projectile nucleus  default: 1   *
+*       what (2) =  charge of projectile nucleus       default: 1   *
+*       what (3..6)   no meaning                                    *
+*       sdum        projectile particle code word                   *
+*                                                                   *
+*       Note: If sdum is defined what (1..2) have no meaning.       *
+*                                                                   *
+*********************************************************************
+
+  110 CONTINUE
+      IF (SDUM.EQ.BLANK) THEN
+         IP     = INT(WHAT(1))
+         IPZ    = INT(WHAT(2))
+         IJPROJ = 1
+         IBPROJ = 1
+      ELSE
+         IJPROJ = 0
+         DO 111 II=1,30
+            IF (SDUM.EQ.BTYPE(II)) THEN
+               IP     = 1
+               IPZ    = 1
+               IF (II.EQ.26) THEN
+                  IJPROJ = 135
+               ELSEIF (II.EQ.27) THEN
+                  IJPROJ = 136
+               ELSEIF (II.EQ.28) THEN
+                  IJPROJ = 133
+               ELSEIF (II.EQ.29) THEN
+                  IJPROJ = 134
+               ELSE
+                  IJPROJ = II
+               ENDIF
+               IBPROJ = IIBAR(IJPROJ)
+* photon
+               IF ((IJPROJ.EQ.7).AND.(WHAT(1).GT.ZERO)) VIRT = WHAT(1)
+* lepton
+               IF (((IJPROJ.EQ. 3).OR.(IJPROJ.EQ. 4).OR.
+     &              (IJPROJ.EQ.10).OR.(IJPROJ.EQ.11)).AND.
+     &                              (WHAT(1).GT.ZERO)) Q2HI = WHAT(1)
+            ENDIF
+  111    CONTINUE
+         IF (IJPROJ.EQ.0) THEN
+            WRITE(LOUT,1110)
+ 1110       FORMAT(/,1X,'invalid PROJPAR card !',/)
+            GOTO 9999
+         ENDIF
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = TARPAR                      *
+*                                                                   *
+*       what (1) =  mass number of target nucleus      default: 1   *
+*       what (2) =  charge of target nucleus           default: 1   *
+*       what (3..6)   no meaning                                    *
+*       sdum        target particle code word                       *
+*                                                                   *
+*       Note: If sdum is defined what (1..2) have no meaning.       *
+*                                                                   *
+*********************************************************************
+
+  120 CONTINUE
+      IF (SDUM.EQ.BLANK) THEN
+         IT     = INT(WHAT(1))
+         ITZ    = INT(WHAT(2))
+         IJTARG = 1
+         IBTARG = 1
+      ELSE
+         IJTARG = 0
+         DO 121 II=1,30
+            IF (SDUM.EQ.BTYPE(II)) THEN
+               IT     = 1
+               ITZ    = 1
+               IJTARG = II
+               IBTARG = IIBAR(IJTARG)
+            ENDIF
+  121    CONTINUE
+         IF (IJTARG.EQ.0) THEN
+            WRITE(LOUT,1120)
+ 1120       FORMAT(/,1X,'invalid TARPAR card !',/)
+            GOTO 9999
+         ENDIF
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = ENERGY                      *
+*                                                                   *
+*       what (1) =  energy (GeV) of projectile in Lab.              *
+*                   if what(1) < 0:  |what(1)| = kinetic energy     *
+*                                                default: 200 GeV   *
+*                   if |what(2)| > 0: min. energy for variable      *
+*                                     energy runs                   *
+*       what (2) =  max. energy for variable energy runs            *
+*                   if what(2) < 0:  |what(2)| = kinetic energy     *
+*                                                                   *
+*********************************************************************
+
+  130 CONTINUE
+      EPN    = WHAT(1)
+      PPN    = ZERO
+      CMENER = ZERO
+      IF ((ABS(WHAT(2)).GT.ZERO).AND.
+     &    (ABS(WHAT(2)).GT.ABS(WHAT(1)))) THEN
+         VARELO = WHAT(1)
+         VAREHI = WHAT(2)
+         EPN    = VAREHI
+      ENDIF
+      LEINP  = .TRUE.
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = MOMENTUM                    *
+*                                                                   *
+*       what (1) =  momentum (GeV/c) of projectile in Lab.          *
+*                                                default: 200 GeV/c *
+*       what (2..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  140 CONTINUE
+      EPN    = ZERO
+      PPN    = WHAT(1)
+      CMENER = ZERO
+      LEINP  = .TRUE.
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = CMENERGY                    *
+*                                                                   *
+*       what (1) =  energy in nucleon-nucleon cms.                  *
+*                                                default: none      *
+*       what (2..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  150 CONTINUE
+      EPN    = ZERO
+      PPN    = ZERO
+      CMENER = WHAT(1)
+      LEINP  = .TRUE.
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = EMULSION                    *
+*                                                                   *
+*               definition of nuclear emulsions                     *
+*                                                                   *
+*     what(1)      mass number of emulsion component                *
+*     what(2)      charge of emulsion component                     *
+*     what(3)      fraction of events in which a scattering on a    *
+*                  nucleus of this properties is performed          *
+*     what(4,5,6)  as what(1,2,3) but for another component         *
+*                                             default: no emulsion  *
+*     sdum         no meaning                                       *
+*                                                                   *
+*     Note: If this input-card is once used with valid parameters   *
+*           TARPAR is obsolete.                                     *
+*           Not the absolute values of the fractions are important  *
+*           but only the ratios of fractions of different comp.     *
+*           This control card can be repeatedly used to define      *
+*           emulsions consisting of up to 10 elements.              *
+*                                                                   *
+*********************************************************************
+
+  160 CONTINUE
+      IF ((WHAT(1).GT.ZERO).AND.(WHAT(2).GT.ZERO)
+     &                     .AND.(ABS(WHAT(3)).GT.ZERO)) THEN
+         NCOMPO = NCOMPO+1
+         IF (NCOMPO.GT.NCOMPX) THEN
+            WRITE(LOUT,1600)
+            STOP
+         ENDIF
+         IEMUMA(NCOMPO) = INT(WHAT(1))
+         IEMUCH(NCOMPO) = INT(WHAT(2))
+         EMUFRA(NCOMPO) = WHAT(3)
+         IEMUL = 1
+C        CALL SHMAKF(IDUM,IDUM,IEMUMA(NCOMPO),IEMUCH(NCOMPO))
+      ENDIF
+      IF ((WHAT(4).GT.ZERO).AND.(WHAT(5).GT.ZERO)
+     &                     .AND.(ABS(WHAT(6)).GT.ZERO)) THEN
+         NCOMPO = NCOMPO+1
+         IF (NCOMPO.GT.NCOMPX) THEN
+            WRITE(LOUT,1001)
+            STOP
+         ENDIF
+         IEMUMA(NCOMPO) = INT(WHAT(4))
+         IEMUCH(NCOMPO) = INT(WHAT(5))
+         EMUFRA(NCOMPO) = WHAT(6)
+C        CALL SHMAKF(IDUM,IDUM,IEMUMA(NCOMPO),IEMUCH(NCOMPO))
+      ENDIF
+ 1600 FORMAT(1X,'too many emulsion components - program stopped')
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = FERMI                       *
+*                                                                   *
+*       what (1) = -1 Fermi-motion of nucleons not treated          *
+*                                                 default: 1        *
+*       what (2) =    scale factor for Fermi-momentum               *
+*                                                 default: 0.75     *
+*       what (3..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  170 CONTINUE
+      IF (WHAT(1).EQ.-1.0D0) THEN
+         LFERMI = .FALSE.
+      ELSE
+         LFERMI = .TRUE.
+      ENDIF
+      XMOD = WHAT(2)
+      IF (XMOD.GE.ZERO) FERMOD = XMOD
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = TAUFOR                      *
+*                                                                   *
+*          formation time supressed intranuclear cascade            *
+*                                                                   *
+*    what (1)      formation time (in fm/c)                         *
+*                  note: what(1)=10. corresponds roughly to an      *
+*                        average formation time of 1 fm/c           *
+*                                                 default: 5. fm/c  *
+*    what (2)      number of generations followed                   *
+*                                                 default: 25       *
+*    what (3) = 1. p_t-dependent formation zone                     *
+*             = 2. constant formation zone                          *
+*                                                 default: 1        *
+*    what (4)      modus of selection of nucleus where the          *
+*                  cascade if followed first                        *
+*             = 1.  proj./target-nucleus with probab. 1/2           *
+*             = 2.  nucleus with highest mass                       *
+*             = 3.  proj. nucleus if particle is moving in pos. z   *
+*                   targ. nucleus if particle is moving in neg. z   *
+*                                                 default: 1        *
+*    what (5..6), sdum   no meaning                                 *
+*                                                                   *
+*********************************************************************
+
+  180 CONTINUE
+      TAUFOR = WHAT(1)
+      KTAUGE = INT(WHAT(2))
+      INCMOD = 1
+      IF ((WHAT(3).GE.1.0D0).AND.(WHAT(3).LE.2.0D0))
+     &                                    ITAUVE = INT(WHAT(3))
+      IF ((WHAT(4).GE.1.0D0).AND.(WHAT(4).LE.3.0D0))
+     &                                    INCMOD = INT(WHAT(4))
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = PAULI                       *
+*                                                                   *
+*       what (1) =  -1  Pauli's principle for secondary             *
+*                       interactions not treated                    *
+*                                                    default: 1     *
+*       what (2..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  190 CONTINUE
+      IF (WHAT(1).EQ.-1.0D0) THEN
+         LPAULI = .FALSE.
+      ELSE
+         LPAULI = .TRUE.
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = COULOMB                     *
+*                                                                   *
+*       what (1) = -1. Coulomb-energy treatment switched off        *
+*                                                    default: 1     *
+*       what (2..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  200 CONTINUE
+      ICOUL = 1
+      IF (WHAT(1).EQ.-1.0D0) THEN
+         ICOUL = 0
+      ELSE
+         ICOUL = 1
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = HADRIN                      *
+*                                                                   *
+*                       HADRIN module                               *
+*                                                                   *
+*    what (1) = 0. elastic/inelastic interactions with probab.      *
+*                  as defined by cross-sections                     *
+*             = 1. inelastic interactions forced                    *
+*             = 2. elastic interactions forced                      *
+*                                                 default: 1        *
+*    what (2)      upper threshold in total energy (GeV) below      *
+*                  which interactions are sampled by HADRIN         *
+*                                                 default: 5. GeV   *
+*    what (3..6), sdum   no meaning                                 *
+*                                                                   *
+*********************************************************************
+
+  210 CONTINUE
+      IWHAT = INT(WHAT(1))
+      IF ((IWHAT.GE.0).AND.(IWHAT.LE.2)) INTHAD = IWHAT
+      IF ((WHAT(2).GT.ZERO).AND.(WHAT(2).LT.15.0D0)) EHADTH = WHAT(2)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = EVAP                        *
+*                                                                   *
+*                    evaporation module                             *
+*                                                                   *
+*  what (1) =< -1 ==> evaporation is switched off                   *
+*           >=  1 ==> evaporation is performed                      *
+*                                                                   *
+*         what (1) = i1 + i2*10 + i3*100 + i4*10000                 *
+*                    (i1, i2, i3, i4 >= 0 )                         *
+*                                                                   *
+*   i1 is the flag for selecting the T=0 level density option used  *
+*      =  1: standard EVAP level densities with Cook pairing        *
+*            energies                                               *
+*      =  2: Z,N-dependent Gilbert & Cameron level densities        *
+*                                                        (default)  *
+*      =  3: Julich A-dependent level densities                     *
+*      =  4: Z,N-dependent Brancazio & Cameron level densities      *
+*                                                                   *
+*   i2 >= 1: high energy fission activated                          *
+*            (default high energy fission activated)                *
+*                                                                   *
+*   i3 =  0: No energy dependence for level densities               *
+*      =  1: Standard Ignyatuk (1975, 1st) energy dependence        *
+*            for level densities (default)                          *
+*      =  2: Standard Ignyatuk (1975, 1st) energy dependence        *
+*            for level densities with NOT used set of parameters    *
+*      =  3: Standard Ignyatuk (1975, 1st) energy dependence        *
+*            for level densities with NOT used set of parameters    *
+*      =  4: Second   Ignyatuk (1975, 2nd) energy dependence        *
+*            for level densities                                    *
+*      =  5: Second   Ignyatuk (1975, 2nd) energy dependence        *
+*            for level densities with fit 1 Iljinov & Mebel set of  *
+*            parameters                                             *
+*      =  6: Second   Ignyatuk (1975, 2nd) energy dependence        *
+*            for level densities with fit 2 Iljinov & Mebel set of  *
+*            parameters                                             *
+*      =  7: Second   Ignyatuk (1975, 2nd) energy dependence        *
+*            for level densities with fit 3 Iljinov & Mebel set of  *
+*            parameters                                             *
+*      =  8: Second   Ignyatuk (1975, 2nd) energy dependence        *
+*            for level densities with fit 4 Iljinov & Mebel set of  *
+*            parameters                                             *
+*                                                                   *
+*   i4 >= 1: Original Gilbert and Cameron pairing energies used     *
+*            (default Cook's modified pairing energies)             *
+*                                                                   *
+*  what (2) = ig + 10 * if   (ig and if must have the same sign)    *
+*                                                                   *
+*   ig =< -1 ==> deexcitation gammas are not produced               *
+*                (if the evaporation step is not performed          *
+*                 they are never produced)                          *
+*   if =< -1 ==> Fermi Break Up is not invoked                      *
+*                (if the evaporation step is not performed          *
+*                 it is never invoked)                              *
+*   The default is: deexcitation gamma produced and Fermi break up  *
+*                   activated for the new  preequilibrium, not      *
+*                   activated otherwise.                            *
+*  what (3..6), sdum   no meaning                                   *
+*                                                                   *
+*********************************************************************
+
+ 220  CONTINUE
+
+      IF (WHAT(1).LE.-1.0D0) THEN
+         LEVPRT = .FALSE.
+         LDEEXG = .FALSE.
+         LHEAVY = .FALSE.
+         GOTO 10
+      ENDIF
+      WHTSAV = WHAT (1)
+      IF ( NINT (WHAT (1)) .GE. 10000 ) THEN
+         LLVMOD   = .FALSE.
+         JLVHLP   = NINT (WHAT (1)) / 10000
+         WHAT (1) = WHAT (1) - 10000.D+00 * JLVHLP
+      END IF
+      IF ( NINT (WHAT (1)) .GE. 100 ) THEN
+         JLVMOD   = NINT (WHAT (1)) / 100
+         WHAT (1) = WHAT (1) - 100.D+00 * JLVMOD
+      END IF
+      IF ( NINT (WHAT (1)) .GE. 10  ) THEN
+         IFISS    = 1
+         JLVHLP   = NINT (WHAT (1)) / 10
+         WHAT (1) = WHAT (1) - 10.D+00 * JLVHLP
+      ELSE IF ( NINT (WHTSAV) .NE. 0 ) THEN
+         IFISS    = 0
+      END IF
+      IF ( NINT (WHAT (1)) .GE. 0 ) THEN
+         LEVPRT = .TRUE.
+         ILVMOD = NINT (WHAT(1))
+         IF ( ABS (NINT (WHAT (2))) .GE. 10  ) THEN
+            LFRMBK   = .TRUE.
+            JLVHLP   = NINT (WHAT (2)) / 10
+            WHAT (2) = WHAT (2) - 10.D+00 * JLVHLP
+         ELSE IF ( NINT (WHAT (2)) .NE. 0 ) THEN
+            LFRMBK   = .FALSE.
+         END IF
+         IF ( NINT (WHAT (2)) .GE. 0 ) THEN
+            LDEEXG = .TRUE.
+         ELSE
+            LDEEXG = .FALSE.
+         END IF
+**sr heavies are always put to /FKFHVY/
+C        IF ( NINT (WHAT(3)) .GE. 1 ) THEN
+C           LHEAVY = .TRUE.
+C        ELSE
+C           LHEAVY = .FALSE.
+C        END IF
+         LHEAVY = .TRUE.
+      ELSE
+         LEVPRT = .FALSE.
+         LDEEXG = .FALSE.
+         LHEAVY = .FALSE.
+      END IF
+
+      LOLDEV = .FALSE.
+
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = EMCCHECK                    *
+*                                                                   *
+*    extended energy-momentum / quantum-number conservation check   *
+*                                                                   *
+*       what (1) = -1   extended check not performed                *
+*                                                    default: 1.    *
+*       what (2..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  230 CONTINUE
+      IF (WHAT(1).EQ.-1) THEN
+         LEMCCK = .FALSE.
+      ELSE
+         LEMCCK = .TRUE.
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = MODEL                       *
+*                                                                   *
+*     Model to be used to treat nucleon-nucleon interactions        *
+*                                                                   *
+*       sdum = DTUNUC    two-chain model                            *
+*            = PHOJET    multiple chains including minijets         *
+*            = LEPTO     DIS                                        *
+*            = QNEUTRIN  quasi-elastic neutrino scattering          *
+*                                                  default: PHOJET  *
+*                                                                   *
+*       if sdum = LEPTO:                                            *
+*       what (1)         (variable INTER)                           *
+*                        = 1  gamma exchange                        *
+*                        = 2  W+-   exchange                        *
+*                        = 3  Z0    exchange                        *
+*                        = 4  gamma/Z0 exchange                     *
+*                                                                   *
+*       if sdum = QNEUTRIN:                                         *
+*       what (1)         = 0  elastic scattering on nucleon and     *
+*                             tau does not decay (default)          *
+*                        = 1  decay of tau into mu..                *
+*                        = 2  decay of tau into e..                 *
+*                        = 10 CC events on p and n                  *
+*                        = 11 NC events on p and n                  *
+*                                                                   *
+*       what (2..6)      no meaning                                 *
+*                                                                   *
+*********************************************************************
+
+  240 CONTINUE
+      IF (SDUM.EQ.CMODEL(1)) THEN
+         MCGENE = 1
+      ELSEIF (SDUM.EQ.CMODEL(2)) THEN
+         MCGENE = 2
+      ELSEIF (SDUM.EQ.CMODEL(3)) THEN
+         MCGENE = 3
+         IF ((WHAT(1).GE.1.0D0).AND.(WHAT(1).LE.4.0D0))
+     &      INTER = INT(WHAT(1))
+      ELSEIF (SDUM.EQ.CMODEL(4)) THEN
+         MCGENE = 4
+         IWHAT  = INT(WHAT(1))
+         IF ((IWHAT.EQ.1 ).OR.(IWHAT.EQ.2 ).OR.
+     &       (IWHAT.EQ.10).OR.(IWHAT.EQ.11))
+     &      NEUDEC = IWHAT
+      ELSE
+         STOP ' Unknown model !'
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = PHOINPUT                    *
+*                                                                   *
+*       Start of input-section for PHOJET-specific input-cards      *
+*       Note:  This section will not be finished before giving      *
+*              ENDINPUT-card                                        *
+*       what (1..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  250 CONTINUE
+      IF (LPHOIN) THEN
+
+C         CALL PHO_INIT(LINP,IREJ1)
+* ###   Read control card from specified file 
+* ### Changed by Chiara (original version LINP=5)
+         CALL PHO_INIT(7,IREJ1)
+
+         IF (IREJ1.NE.0) THEN
+            WRITE(LOUT,'(1X,A)')'INIT:   reading PHOJET-input failed'
+            STOP
+         ENDIF
+         LPHOIN = .FALSE.
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = GLAUBERI                    *
+*                                                                   *
+*        Pre-initialization of impact parameter selection           *
+*                                                                   *
+*        what (1..6), sdum   no meaning                             *
+*                                                                   *
+*********************************************************************
+
+  260 CONTINUE
+      IF (IFIRST.NE.99) THEN
+         CALL DT_RNDMST(12,34,56,78)
+         CALL DT_RNDMTE(1)
+         OPEN(40,FILE='shm.out',STATUS='UNKNOWN')
+C        OPEN(11,FILE='shm.dbg',STATUS='UNKNOWN')
+         IFIRST = 99
+      ENDIF
+
+      IPPN = 8
+      PLOW = 10.0D0
+C     IPPN = 1
+C     PLOW = 100.0D0
+      PHI  = 1.0D5
+      APLOW = LOG10(PLOW)
+      APHI  = LOG10(PHI)
+      ADP   = (APHI-APLOW)/DBLE(IPPN)
+
+      IPLOW = 1
+      IDIP  = 1
+      IIP   = 5
+C     IPLOW = 1
+C     IDIP  = 1
+C     IIP   = 1
+      IPRANG(1) = 1
+      IPRANG(2) = 2
+      IPRANG(3) = 5
+      IPRANG(4) = 10
+      IPRANG(5) = 20
+
+      ITLOW = 30
+      IDIT  = 3
+      IIT   = 60
+C     IDIT  = 10
+C     IIT   = 21
+
+      DO 473 NCIT=1,IIT
+         IT   = ITLOW+(NCIT-1)*IDIT
+C        IPHI = IT
+C        IDIP = 10
+C        IIP  = (IPHI-IPLOW)/IDIP
+C        IF (IIP.EQ.0) IIP = 1
+C        IF (IT.EQ.IPLOW) IIP = 0
+
+         DO 472 NCIP=1,IIP
+            IP = IPRANG(NCIP)
+CC           IF (NCIP.LE.IIP) THEN
+C               IP = IPLOW+(NCIP-1)*IDIP
+CC           ELSE
+CC              IP = IT
+CC           ENDIF
+            IF (IP.GT.IT) GOTO 472
+
+            DO 471 NCP=1,IPPN+1
+               APPN = APLOW+DBLE(NCP-1)*ADP
+               PPN  = 10**APPN
+
+               OPEN(12,FILE='shm.sta',STATUS='UNKNOWN')
+               WRITE(12,'(1X,2I5,E15.3)') IP,IT,PPN
+               CLOSE(12)
+
+               XLIM1 = 0.0D0
+               XLIM2 = 50.0D0
+               XLIM3 = ZERO
+               IBIN  = 50
+               CALL DT_NEWHGR(XDUM,XDUM,XDUM,XDUMB,-1,IHDUM)
+               CALL DT_NEWHGR(XLIM1,XLIM2,XLIM3,XDUMB,IBIN,IHSHMA)
+
+               NEVFIT = 5
+C              IF ((IP.GT.10).OR.(IT.GT.10)) THEN
+C                 NEVFIT = 5
+C              ELSE
+C                 NEVFIT = 10
+C              ENDIF
+               SIGAV  = 0.0D0
+
+               DO 478 I=1,NEVFIT
+                  CALL DT_SHMAKI(IP,IDUM1,IT,IDUM1,IJPROJ,PPN,99)
+                  SIGAV = SIGAV+XSPRO(1,1,1)
+                  DO 479 J=1,50
+                     XC = DBLE(J)
+                     CALL DT_FILHGR(XC,BSITE(1,1,1,J),IHSHMA,I)
+  479             CONTINUE
+  478          CONTINUE
+
+               CALL DT_EVTHIS(IDUM)
+               HEADER = ' BSITE'
+C              CALL OUTGEN(IHSHMA,0,0,0,0,0,HEADER,0,NEVFIT,ONE,0,1,-1)
+
+C              CALL GENFIT(XPARA)
+C              WRITE(40,'(2I4,E11.3,F6.0,5E11.3)')
+C    &              IP,IT,PPN,SIGAV/DBLE(NEVFIT),XPARA
+
+  471       CONTINUE
+
+  472    CONTINUE
+
+  473 CONTINUE
+
+      STOP
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = FLUCTUAT                    *
+*                                                                   *
+*           Treatment of cross section fluctuations                 *
+*                                                                   *
+*       what (1) = 1  treat cross section fluctuations              *
+*                                                    default: 0.    *
+*       what (1..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+ 270  CONTINUE
+      IFLUCT = 0
+      IF (WHAT(1).EQ.ONE) THEN
+         IFLUCT = 1
+         CALL DT_FLUINI
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = CENTRAL                     *
+*                                                                   *
+*       what (1) = 1.  central production forced     default: 0     *
+*  if what (1) < 0 and > -100                                       *
+*       what (2) = min. impact parameter             default: 0     *
+*       what (3) = max. impact parameter             default: b_max *
+*  if what (1) < -99                                                *
+*       what (2) = fraction of cross section         default: 1     *
+*  if what (1) = -1 : evaporation/fzc suppressed                    *
+*  if what (1) < -1 : evaporation/fzc allowed                       *
+*                                                                   *
+*       what (4..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  280 CONTINUE
+      ICENTR = INT(WHAT(1))
+      IF (ICENTR.LT.0) THEN
+         IF (ICENTR.GT.-100) THEN
+            BIMIN = WHAT(2)
+            BIMAX = WHAT(3)
+         ELSE
+            XSFRAC = WHAT(2)
+         ENDIF
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = RECOMBIN                    *
+*                                                                   *
+*                     Chain recombination                           *
+*        (recombine S-S and V-V chains to V-S chains)               *
+*                                                                   *
+*       what (1) = -1. recombination switched off    default: 1     *
+*       what (2..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  290 CONTINUE
+      IRECOM = 1
+      IF (WHAT(1).EQ.-1.0D0) IRECOM = 0
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = COMBIJET                    *
+*                                                                   *
+*               chain fusion (2 q-aq --> qq-aqaq)                   *
+*                                                                   *
+*       what (1) = 1   fusion treated                               *
+*                                                    default: 0.    *
+*       what (2)       minimum number of uncombined chains from     *
+*                      single projectile or target nucleons         *
+*                                                    default: 0.    *
+*       what (3..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  300 CONTINUE
+      LCO2CR = .FALSE.
+      IF (INT(WHAT(1)).EQ.1) LCO2CR = .TRUE.
+      IF (WHAT(2).GE.ZERO) CUTOF = WHAT(2)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = XCUTS                       *
+*                                                                   *
+*                 thresholds for x-sampling                         *
+*                                                                   *
+*    what (1)    defines lower threshold for val.-q x-value (CVQ)   *
+*                                                 default: 1.       *
+*    what (2)    defines lower threshold for val.-qq x-value (CDQ)  *
+*                                                 default: 2.       *
+*    what (3)    defines lower threshold for sea-q x-value (CSEA)   *
+*                                                 default: 0.2      *
+*    what (4)    sea-q x-values in S-S chains (SSMIMA)              *
+*                                                 default: 0.14     *
+*    what (5)    not used                                           *
+*                                                 default: 2.       *
+*    what (6), sdum   no meaning                                    *
+*                                                                   *
+*    Note: Lower thresholds (what(1..3)) are def. as x_thr=CXXX/ECM *
+*                                                                   *
+*********************************************************************
+
+  310 CONTINUE
+      IF (WHAT(1).GE.0.5D0) CVQ    = WHAT(1)
+      IF (WHAT(2).GE.ONE)   CDQ    = WHAT(2)
+      IF (WHAT(3).GE.0.1D0) CSEA   = WHAT(3)
+      IF (WHAT(4).GE.ZERO) THEN
+         SSMIMA = WHAT(4)
+         SSMIMQ = SSMIMA**2
+      ENDIF
+      IF (WHAT(5).GT.2.0D0) VVMTHR = WHAT(5)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = INTPT                       *
+*                                                                   *
+*     what (1) = -1   intrinsic transverse momenta of partons       *
+*                     not treated                default: 1         *
+*     what (2..6), sdum   no meaning                                *
+*                                                                   *
+*********************************************************************
+
+  320 CONTINUE
+      IF (WHAT(1).EQ.-1.0D0) THEN
+         LINTPT = .FALSE.
+      ELSE
+         LINTPT = .TRUE.
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = CRONINPT                    *
+*                                                                   *
+*    Cronin effect (multiple scattering of partons at chain ends)   *
+*                                                                   *
+*       what (1) = -1  Cronin effect not treated     default: 1     *
+*       what (2) = 0   scattering parameter          default: 0.64  *
+*       what (3..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  330 CONTINUE
+      IF (WHAT(1).EQ.-1.0D0) THEN
+         MKCRON = 0
+      ELSE
+         MKCRON = 1
+      ENDIF
+      CRONCO = WHAT(2)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = SEADISTR                    *
+*                                                                   *
+*     what (1)  (XSEACO)  sea(x) prop. 1/x**what (1)   default: 1.  *
+*     what (2)  (UNON)                                 default: 2.  *
+*     what (3)  (UNOM)                                 default: 1.5 *
+*     what (4)  (UNOSEA)                               default: 5.  *
+*                        qdis(x) prop. (1-x)**what (1)  etc.        *
+*     what (5..6), sdum   no meaning                                *
+*                                                                   *
+*********************************************************************
+
+  340 CONTINUE
+      XSEACO = WHAT(1)
+      XSEACU = 1.05D0-XSEACO
+      UNON   = WHAT(2)
+      IF (UNON.LT.0.1D0) UNON = 2.0D0
+      UNOM   = WHAT(3)
+      IF (UNOM.LT.0.1D0) UNOM = 1.5D0
+      UNOSEA = WHAT(4)
+      IF (UNOSEA.LT.0.1D0) UNOSEA = 5.0D0
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = SEASU3                      *
+*                                                                   *
+*          Treatment of strange-quarks at chain ends                *
+*                                                                   *
+*       what (1)   (SEASQ)  strange-quark supression factor         *
+*                  iflav = 1.+rndm*(2.+SEASQ)                       *
+*                                                    default: 1.    *
+*       what (2..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  350 CONTINUE
+      SEASQ = WHAT(1)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = DIQUARKS                    *
+*                                                                   *
+*     what (1) = -1.  sea-diquark/antidiquark-pairs not treated     *
+*                                                    default: 1.    *
+*     what (2..6), sdum   no meaning                                *
+*                                                                   *
+*********************************************************************
+
+ 360  CONTINUE
+      IF (WHAT(1).EQ.-1.0D0) THEN
+         LSEADI = .FALSE.
+      ELSE
+         LSEADI = .TRUE.
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = RESONANC                    *
+*                                                                   *
+*                 treatment of low mass chains                      *
+*                                                                   *
+*    what (1) = -1 low chain masses are not corrected for resonance *
+*                  masses (obsolete for BAMJET-fragmentation)       *
+*                                       default: 1.                 *
+*    what (2) = -1 massless partons     default: 1. (massive)       *
+*                                       default: 1. (massive)       *
+*    what (3) = -1 chain-system containing chain of too small       *
+*                  mass is rejected (note: this does not fully      *
+*                  apply to S-S chains) default: 0.                 *
+*    what (4..6), sdum   no meaning                                 *
+*                                                                   *
+*********************************************************************
+
+  370 CONTINUE
+      IRESCO = 1
+      IMSHL  = 1
+      IRESRJ = 0
+      IF (WHAT(1).EQ.-ONE) IRESCO = 0
+      IF (WHAT(2).EQ.-ONE) IMSHL  = 0
+      IF (WHAT(3).EQ.-ONE) IRESRJ = 1
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = DIFFRACT                    *
+*                                                                   *
+*                Treatment of diffractive events                    *
+*                                                                   *
+*     what (1) = (ISINGD) 0  no single diffraction                  *
+*                         1  single diffraction included            *
+*                       +-2  single diffractive events only         *
+*                       +-3  projectile single diffraction only     *
+*                       +-4  target single diffraction only         *
+*                        -5  double pomeron exchange only           *
+*                      (neg. sign applies to PHOJET events)         *
+*                                                     default: 0.   *
+*                                                                   *
+*     what (2) = (IDOUBD) 0  no double diffraction                  *
+*                         1  double diffraction included            *
+*                         2  double diffractive events only         *
+*                                                     default: 0.   *
+*     what (3) = 1 projectile diffraction treated (2-channel form.) *
+*                                                     default: 0.   *
+*     what (4) = alpha-parameter in projectile diffraction          *
+*                                                     default: 0.   *
+*     what (5..6), sdum   no meaning                                *
+*                                                                   *
+*********************************************************************
+
+  380 CONTINUE
+      IF (ABS(WHAT(1)).GT.ZERO) ISINGD = INT(WHAT(1))
+      IF (ABS(WHAT(2)).GT.ZERO) IDOUBD = INT(WHAT(2))
+      IF ((ISINGD.GT.1).AND.(IDOUBD.GT.1)) THEN
+         WRITE(LOUT,1380)
+ 1380    FORMAT(1X,'INIT:   inconsistent DIFFRACT - input !',/,
+     &          11X,'IDOUBD is reset to zero')
+         IDOUBD = 0
+      ENDIF
+      IF (WHAT(3).GT.ZERO) DIBETA = WHAT(3)
+      IF (WHAT(4).GT.ZERO) DIALPH = WHAT(4)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = SINGLECH                    *
+*                                                                   *
+*       what (1) = 1.  Regge contribution (one chain) included      *
+*                                                   default: 0.     *
+*       what (2..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+ 390  CONTINUE
+      ISICHA = 0
+      IF (WHAT(1).EQ.ONE) ISICHA = 1
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = NOFRAGME                    *
+*                                                                   *
+*                 biased chain hadronization                        *
+*                                                                   *
+*       what (1..6) = -1  no of hadronizsation of S-S chains        *
+*                   = -2  no of hadronizsation of D-S chains        *
+*                   = -3  no of hadronizsation of S-D chains        *
+*                   = -4  no of hadronizsation of S-V chains        *
+*                   = -5  no of hadronizsation of D-V chains        *
+*                   = -6  no of hadronizsation of V-S chains        *
+*                   = -7  no of hadronizsation of V-D chains        *
+*                   = -8  no of hadronizsation of V-V chains        *
+*                   = -9  no of hadronizsation of comb. chains      *
+*                                  default:  complete hadronization *
+*       sdum   no meaning                                           *
+*                                                                   *
+*********************************************************************
+
+  400 CONTINUE
+      DO 401 I=1,6
+         ICHAIN = INT(WHAT(I))
+         IF ((ICHAIN.LE.-1).AND.(ICHAIN.GE.-9))
+     &      LHADRO(ABS(ICHAIN)) = .FALSE.
+  401 CONTINUE
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = HADRONIZE                   *
+*                                                                   *
+*           hadronization model and parameter switch                *
+*                                                                   *
+*       what (1) = 1    hadronization via BAMJET                    *
+*                = 2    hadronization via JETSET                    *
+*                                                    default: 2     *
+*       what (2) = 1..3 parameter set to be used                    *
+*                       JETSET: 3 sets available                    *
+*                               ( = 3 default JETSET-parameters)    *
+*                       BAMJET: 1 set available                     *
+*                                                    default: 1     *
+*       what (3..6), sdum   no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  410 CONTINUE
+      IWHAT1 = INT(WHAT(1))
+      IWHAT2 = INT(WHAT(2))
+      IF ((IWHAT1.EQ.1).OR.(IWHAT1.EQ.2)) IFRAG(1) = IWHAT1
+      IF ((IWHAT1.EQ.2).AND.(IWHAT2.GE.1).AND.(IWHAT2.LE.3))
+     &                                    IFRAG(2) = IWHAT2
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = POPCORN                     *
+*                                                                   *
+*  "Popcorn-effect" in fragmentation and diquark breaking diagrams  *
+*                                                                   *
+*   what (1) = (PDB) frac. of diquark fragmenting directly into     *
+*                    baryons (PYTHIA/JETSET fragmentation)          *
+*                    (JETSET: = 0. Popcorn mechanism switched off)  *
+*                                                    default: 0.5   *
+*   what (2) = probability for accepting a diquark breaking         *
+*              diagram involving the generation of a u/d quark-     *
+*              antiquark pair                        default: 0.0   *
+*   what (3) = same a what (2), here for s quark-antiquark pair     *
+*                                                    default: 0.0   *
+*   what (4..6), sdum   no meaning                                  *
+*                                                                   *
+*********************************************************************
+
+  420 CONTINUE
+      IF (WHAT(1).GE.0.0D0) PDB = WHAT(1)
+      IF (WHAT(2).GE.0.0D0) THEN
+         PDBSEA(1) = WHAT(2)
+         PDBSEA(2) = WHAT(2)
+      ENDIF
+      IF (WHAT(3).GE.0.0D0) PDBSEA(3) = WHAT(3)
+      DO 421 I=1,8
+         DBRKA(1,I) = DBRKR(1,I)*PDBSEA(1)/(1.D0-PDBSEA(1))
+         DBRKA(2,I) = DBRKR(2,I)*PDBSEA(2)/(1.D0-PDBSEA(2))
+         DBRKA(3,I) = DBRKR(3,I)*PDBSEA(3)/(1.D0-PDBSEA(3))
+  421 CONTINUE
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = PARDECAY                    *
+*                                                                   *
+*      what (1) = 1.  Sigma0/Asigma0 are decaying within JETSET     *
+*               = 2.  pion^0 decay after intranucl. cascade         *
+*                                                default: no decay  *
+*      what (2..6), sdum   no meaning                               *
+*                                                                   *
+*********************************************************************
+
+ 430  CONTINUE
+      IF (WHAT(1).EQ.ONE)  ISIG0 = 1
+      IF (WHAT(1).EQ.2.0D0) IPI0 = 1
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = BEAM                        *
+*                                                                   *
+*              definition of beam parameters                        *
+*                                                                   *
+*      what (1/2)  > 0 : energy of beam 1/2 (GeV)                   *
+*                  < 0 : abs(what(1/2)) energy per charge of        *
+*                        beam 1/2 (GeV)                             *
+*                  (beam 1 is directed into positive z-direction)   *
+*      what (3)    beam crossing angle, defined as 2x angle between *
+*                  one beam and the z-axis (micro rad)              *
+*      what (4)    angle with x-axis defining the collision plane   *
+*      what (5..6), sdum   no meaning                               *
+*                                                                   *
+*      Note: this card requires previously defined projectile and   *
+*            target identities (PROJPAR, TARPAR)                    *
+*                                                                   *
+*********************************************************************
+
+  440 CONTINUE
+      CALL DT_BEAMPR(WHAT,PPN,1)
+      EPN    = ZERO
+      CMENER = ZERO
+      LEINP  = .TRUE.
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = LUND-MSTU                   *
+*                                                                   *
+*          set parameter MSTU in JETSET-common /LUDAT1/             *
+*                                                                   *
+*       what (1) =  index according to LUND-common block            *
+*       what (2) =  new value of MSTU( int(what(1)) )               *
+*       what (3), what(4) and what (5), what(6) further             *
+*                   parameter in the same way as what (1) and       *
+*                   what (2)                                        *
+*                        default: default-Lund or corresponding to  *
+*                                 the set given in HADRONIZE        *
+*                                                                   *
+*********************************************************************
+
+  450 CONTINUE
+      IF (WHAT(1).GT.ZERO) THEN
+         NMSTU = NMSTU+1
+         IMSTU(NMSTU) = INT(WHAT(1))
+         MSTUX(NMSTU) = INT(WHAT(2))
+      ENDIF
+      IF (WHAT(3).GT.ZERO) THEN
+         NMSTU = NMSTU+1
+         IMSTU(NMSTU) = INT(WHAT(3))
+         MSTUX(NMSTU) = INT(WHAT(4))
+      ENDIF
+      IF (WHAT(5).GT.ZERO) THEN
+         NMSTU = NMSTU+1
+         IMSTU(NMSTU) = INT(WHAT(5))
+         MSTUX(NMSTU) = INT(WHAT(6))
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = LUND-MSTJ                   *
+*                                                                   *
+*          set parameter MSTJ in JETSET-common /LUDAT1/             *
+*                                                                   *
+*       what (1) =  index according to LUND-common block            *
+*       what (2) =  new value of MSTJ( int(what(1)) )               *
+*       what (3), what(4) and what (5), what(6) further             *
+*                   parameter in the same way as what (1) and       *
+*                   what (2)                                        *
+*                        default: default-Lund or corresponding to  *
+*                                 the set given in HADRONIZE        *
+*                                                                   *
+*********************************************************************
+
+  451 CONTINUE
+      IF (WHAT(1).GT.ZERO) THEN
+         NMSTJ = NMSTJ+1
+         IMSTJ(NMSTJ) = INT(WHAT(1))
+         MSTJX(NMSTJ) = INT(WHAT(2))
+      ENDIF
+      IF (WHAT(3).GT.ZERO) THEN
+         NMSTJ = NMSTJ+1
+         IMSTJ(NMSTJ) = INT(WHAT(3))
+         MSTJX(NMSTJ) = INT(WHAT(4))
+      ENDIF
+      IF (WHAT(5).GT.ZERO) THEN
+         NMSTJ = NMSTJ+1
+         IMSTJ(NMSTJ) = INT(WHAT(5))
+         MSTJX(NMSTJ) = INT(WHAT(6))
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = LUND-MDCY                   *
+*                                                                   *
+*  set parameter MDCY(I,1) for particle decays in JETSET-common     *
+*                                                      /LUDAT3/     *
+*                                                                   *
+*       what (1-6) = PDG particle index of particle which should    *
+*                    not decay                                      *
+*                        default: default-Lund or forced in         *
+*                                 DT_INITJS                         *
+*                                                                   *
+*********************************************************************
+
+  452 CONTINUE
+      DO 4521 I=1,6
+         IF (WHAT(I).NE.ZERO) THEN
+
+            KC = PYCOMP(INT(WHAT(I)))
+
+            MDCY(KC,1) = 0
+         ENDIF
+ 4521 CONTINUE
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = LUND-PARJ                   *
+*                                                                   *
+*          set parameter PARJ in JETSET-common /LUDAT1/             *
+*                                                                   *
+*       what (1) =  index according to LUND-common block            *
+*       what (2) =  new value of PARJ( int(what(1)) )               *
+*       what (3), what(4) and what (5), what(6) further             *
+*                   parameter in the same way as what (1) and       *
+*                   what (2)                                        *
+*                        default: default-Lund or corresponding to  *
+*                                 the set given in HADRONIZE        *
+*                                                                   *
+*********************************************************************
+
+  460 CONTINUE
+      IF (WHAT(1).NE.ZERO) THEN
+         NPARJ = NPARJ+1
+         IPARJ(NPARJ) = INT(WHAT(1))
+         PARJX(NPARJ) = WHAT(2)
+      ENDIF
+      IF (WHAT(3).NE.ZERO) THEN
+         NPARJ = NPARJ+1
+         IPARJ(NPARJ) = INT(WHAT(3))
+         PARJX(NPARJ) = WHAT(4)
+      ENDIF
+      IF (WHAT(5).NE.ZERO) THEN
+         NPARJ = NPARJ+1
+         IPARJ(NPARJ) = INT(WHAT(5))
+         PARJX(NPARJ) = WHAT(6)
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = LUND-PARU                   *
+*                                                                   *
+*          set parameter PARJ in JETSET-common /LUDAT1/             *
+*                                                                   *
+*       what (1) =  index according to LUND-common block            *
+*       what (2) =  new value of PARU( int(what(1)) )               *
+*       what (3), what(4) and what (5), what(6) further             *
+*                   parameter in the same way as what (1) and       *
+*                   what (2)                                        *
+*                        default: default-Lund or corresponding to  *
+*                                 the set given in HADRONIZE        *
+*                                                                   *
+*********************************************************************
+
+  470 CONTINUE
+      IF (WHAT(1).GT.ZERO) THEN
+         NPARU = NPARU+1
+         IPARU(NPARU) = INT(WHAT(1))
+         PARUX(NPARU) = WHAT(2)
+      ENDIF
+      IF (WHAT(3).GT.ZERO) THEN
+         NPARU = NPARU+1
+         IPARU(NPARU) = INT(WHAT(3))
+         PARUX(NPARU) = WHAT(4)
+      ENDIF
+      IF (WHAT(5).GT.ZERO) THEN
+         NPARU = NPARU+1
+         IPARU(NPARU) = INT(WHAT(5))
+         PARUX(NPARU) = WHAT(6)
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = OUTLEVEL                    *
+*                                                                   *
+*                    output control switches                        *
+*                                                                   *
+*       what (1) =  internal rejection informations  default: 0     *
+*       what (2) =  energy-momentum conservation check output       *
+*                                                    default: 0     *
+*       what (3) =  internal warning messages        default: 0     *
+*       what (4..6), sdum    not yet used                           *
+*                                                                   *
+*********************************************************************
+
+  480 CONTINUE
+      DO 481 K=1,6
+         IOULEV(K) = INT(WHAT(K))
+  481 CONTINUE
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = FRAME                       *
+*                                                                   *
+*          frame in which final state is given in DTEVT1            *
+*                                                                   *
+*       what (1) = 1  target rest frame (laboratory)                *
+*                = 2  nucleon-nucleon cms                           *
+*                                                    default: 1     *
+*                                                                   *
+*********************************************************************
+
+  490 CONTINUE
+      KFRAME = INT(WHAT(1))
+      IF ((KFRAME.GE.1).AND.(KFRAME.LE.2)) IFRAME = KFRAME
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = L-TAG                       *
+*                                                                   *
+*                        lepton tagger:                             *
+*   definition of kinematical cuts for radiated photon and          *
+*   outgoing lepton detection in lepton-nucleus interactions        *
+*                                                                   *
+*       what (1) = y_min                                            *
+*       what (2) = y_max                                            *
+*       what (3) = Q^2_min                                          *
+*       what (4) = Q^2_max                                          *
+*       what (5) = theta_min  (Lab)                                 *
+*       what (6) = theta_max  (Lab)                                 *
+*                                       default: no cuts            *
+*       sdum    no meaning                                          *
+*                                                                   *
+*********************************************************************
+
+  500 CONTINUE
+      YMIN  = WHAT(1)
+      YMAX  = WHAT(2)
+      Q2MIN = WHAT(3)
+      Q2MAX = WHAT(4)
+      THMIN = WHAT(5)
+      THMAX = WHAT(6)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = L-ETAG                      *
+*                                                                   *
+*                        lepton tagger:                             *
+*       what (1) = min. outgoing lepton energy  (in Lab)            *
+*       what (2) = min. photon energy           (in Lab)            *
+*       what (3) = max. photon energy           (in Lab)            *
+*                                       default: no cuts            *
+*       what (2..6), sdum    no meaning                             *
+*                                                                   *
+*********************************************************************
+
+  510 CONTINUE
+      ELMIN = MAX(WHAT(1),ZERO)
+      EGMIN = MAX(WHAT(2),ZERO)
+      EGMAX = MAX(WHAT(3),ZERO)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = ECMS-CUT                    *
+*                                                                   *
+*     what (1) = min. c.m. energy to be sampled                     *
+*     what (2) = max. c.m. energy to be sampled                     *
+*     what (3) = min x_Bj         to be sampled                     *
+*                                       default: no cuts            *
+*     what (3..6), sdum    no meaning                               *
+*                                                                   *
+*********************************************************************
+
+  520 CONTINUE
+      ECMIN  = WHAT(1)
+      ECMAX  = WHAT(2)
+      IF (ECMIN.GT.ECMAX) ECMIN = ECMAX
+      XBJMIN = MAX(WHAT(3),ZERO)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = VDM-PAR1                    *
+*                                                                   *
+*      parameters in gamma-nucleus cross section calculation        *
+*                                                                   *
+*       what (1) =  Lambda^2                       default: 2.      *
+*       what (2)    lower limit in M^2 integration                  *
+*                =  1  (3m_pi)^2                                    *
+*                =  2  (m_rho0)^2                                   *
+*                =  3  (m_phi)^2                   default: 1       *
+*       what (3)    upper limit in M^2 integration                  *
+*                =  1   s/2                                         *
+*                =  2   s/4                                         *
+*                =  3   s                          default: 3       *
+*       what (4)    CKMT F_2 structure function                     *
+*                =  2212  proton                                    *
+*                =  100   deuteron                 default: 2212    *
+*       what (5)    calculation of gamma-nucleon xsections          *
+*                =  1  according to CKMT-parametrization of F_2     *
+*                =  2  integrating SIGVP over M^2                   *
+*                =  3  using SIGGA                                  *
+*                =  4  PHOJET cross sections       default:  4      *
+*                                                                   *
+*       what (6), sdum    no meaning                                *
+*                                                                   *
+*********************************************************************
+
+  530 CONTINUE
+      IF (WHAT(1).GE.ZERO) RL2 = WHAT(1)
+      IF ((WHAT(2).GE.1).AND.(WHAT(2).LE.3)) INTRGE(1) = INT(WHAT(2))
+      IF ((WHAT(3).GE.1).AND.(WHAT(3).LE.3)) INTRGE(2) = INT(WHAT(3))
+      IF ((WHAT(4).EQ.2212).OR.(WHAT(4).EQ.100)) IDPDF = INT(WHAT(4))
+      IF ((WHAT(5).GE.1).AND.(WHAT(5).LE.4)) MODEGA = INT(WHAT(5))
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = HISTOGRAM                   *
+*                                                                   *
+*           activate different classes of histograms                *
+*                                                                   *
+*                                default: no histograms             *
+*                                                                   *
+*********************************************************************
+
+  540 CONTINUE
+      DO 541 J=1,6
+         IF ((WHAT(J).GE.100).AND.(WHAT(J).LE.150)) THEN
+            IHISPP(INT(WHAT(J))-100) = 1
+         ELSEIF ((ABS(WHAT(J)).GE.200).AND.(ABS(WHAT(J)).LE.250)) THEN
+            IHISXS(INT(ABS(WHAT(J)))-200) = 1
+            IF (WHAT(J).LT.ZERO) IXSTBL = 1
+         ENDIF
+  541 CONTINUE
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = XS-TABLE                    *
+*                                                                   *
+*    output of cross section table for requested interaction        *
+*              - particle production deactivated ! -                *
+*                                                                   *
+*       what (1)      lower energy limit for tabulation             *
+*                > 0  Lab. frame                                    *
+*                < 0  nucleon-nucleon cms                           *
+*       what (2)      upper energy limit for tabulation             *
+*                > 0  Lab. frame                                    *
+*                < 0  nucleon-nucleon cms                           *
+*       what (3) > 0  # of equidistant lin. bins in E               *
+*                < 0  # of equidistant log. bins in E               *
+*       what (4)      lower limit of particle virtuality (photons)  *
+*       what (5)      upper limit of particle virtuality (photons)  *
+*       what (6) > 0  # of equidistant lin. bins in Q^2             *
+*                < 0  # of equidistant log. bins in Q^2             *
+*                                                                   *
+*********************************************************************
+
+  550 CONTINUE
+      IF (WHAT(1).EQ.99999.0D0) THEN
+         IRATIO = INT(WHAT(2))
+         GOTO 10
+      ENDIF
+      CMENER = ABS(WHAT(2))
+      IF (.NOT.LXSTAB) THEN
+
+         CALL BERTTP
+         CALL INCINI
+
+      ENDIF
+      IF ((.NOT.LXSTAB).OR.(CMENER.NE.CMEOLD)) THEN
+         CMEOLD = CMENER
+         IF (WHAT(2).GT.ZERO)
+     &      CMENER = SQRT(2.0D0*AAM(1)**2+2.0D0*WHAT(2)*AAM(1))
+         EPN = ZERO
+         PPN = ZERO
+C        WRITE(LOUT,*) 'CMENER = ',CMENER
+         CALL DT_LTINI(IJPROJ,IJTARG,EPN,PPN,CMENER,1)
+         CALL DT_PHOINI
+      ENDIF
+      CALL DT_XSTABL(WHAT,IXSQEL,IRATIO)
+      IXSQEL = 0
+      LXSTAB = .TRUE.
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = GLAUB-PAR                   *
+*                                                                   *
+*                parameters in Glauber-formalism                    *
+*                                                                   *
+*    what (1)  # of nucleon configurations sampled in integration   *
+*              over nuclear desity                default: 1000     *
+*    what (2)  # of bins for integration over impact-parameter and  *
+*              for profile-function calculation   default: 49       *
+*    what (3)  = 1 calculation of tot., el. and qel. cross sections *
+*                                                 default: 0        *
+*    what (4)  = 1   read pre-calculated impact-parameter distrib.  *
+*                    from "sdum".glb                                *
+*              =-1   dump pre-calculated impact-parameter distrib.  *
+*                    into "sdum".glb                                *
+*              = 100 read pre-calculated impact-parameter distrib.  *
+*                    for variable projectile/target/energy runs     *
+*                    from "sdum".glb                                *
+*                                                 default: 0        *
+*    what (5..6)   no meaning                                       *
+*    sdum      if |what (4)| = 1 name of in/output-file (sdum.glb)  *
+*                                                                   *
+*********************************************************************
+
+  560 CONTINUE
+      IF (WHAT(1).GT.ZERO) JSTATB = INT(WHAT(1))
+      IF (WHAT(2).GT.ZERO) JBINSB = INT(WHAT(2))
+      IF (WHAT(3).EQ.ONE) LPROD = .FALSE.
+      IF ((ABS(WHAT(4)).EQ.ONE).OR.(WHAT(4).EQ.100)) THEN
+         IOGLB = INT(WHAT(4))
+         CGLB  = SDUM
+      ENDIF
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = GLAUB-INI                   *
+*                                                                   *
+*             pre-initialization of profile function                *
+*                                                                   *
+*       what (1)      lower energy limit for initialization         *
+*                > 0  Lab. frame                                    *
+*                < 0  nucleon-nucleon cms                           *
+*       what (2)      upper energy limit for initialization         *
+*                > 0  Lab. frame                                    *
+*                < 0  nucleon-nucleon cms                           *
+*       what (3) > 0  # of equidistant lin. bins in E               *
+*                < 0  # of equidistant log. bins in E               *
+*       what (4)      maximum projectile mass number for which the  *
+*                     Glauber data are initialized for each         *
+*                     projectile mass number                        *
+*                     (if <= mass given with the PROJPAR-card)      *
+*                                              default: 18          *
+*       what (5)      steps in mass number starting from what (4)   *
+*                     up to mass number defined with PROJPAR-card   *
+*                     for which Glauber data are initialized        *
+*                                              default: 5           *
+*       what (6)      no meaning                                    *
+*       sdum          no meaning                                    *
+*                                                                   *
+*********************************************************************
+
+  565 CONTINUE
+      IOGLB = -100
+      CALL DT_GLBINI(WHAT)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = VDM-PAR2                    *
+*                                                                   *
+*      parameters in gamma-nucleus cross section calculation        *
+*                                                                   *
+*      what (1) = 0 no suppression of shadowing by direct photon    *
+*                   processes                                       *
+*               = 1 suppression ..                   default: 1     *
+*      what (2) = 0 no suppression of shadowing by anomalous        *
+*                   component if photon-F_2                         *
+*               = 1 suppression ..                   default: 1     *
+*      what (3) = 0 no suppression of shadowing by coherence        *
+*                   length of the photon                            *
+*               = 1 suppression ..                   default: 1     *
+*      what (4) = 1 longitudinal polarized photons are taken into   *
+*                   account                                         *
+*                   eps*R*Q^2/M^2 = what(4)*Q^2/M^2  default: 0     *
+*      what (5..6), sdum    no meaning                              *
+*                                                                   *
+*********************************************************************
+
+  570 CONTINUE
+      IF ((WHAT(1).EQ.ZERO).OR.(WHAT(1).EQ.ONE)) ISHAD(1) = INT(WHAT(1))
+      IF ((WHAT(2).EQ.ZERO).OR.(WHAT(2).EQ.ONE)) ISHAD(2) = INT(WHAT(2))
+      IF ((WHAT(3).EQ.ZERO).OR.(WHAT(3).EQ.ONE)) ISHAD(3) = INT(WHAT(3))
+      EPSPOL  = WHAT(4)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  XS-QELPRO                            *
+*                                                                   *
+*     what (1..6), sdum    no meaning                               *
+*                                                                   *
+*********************************************************************
+
+  580 CONTINUE
+      IXSQEL = ABS(WHAT(1))
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  RNDMINIT                             *
+*                                                                   *
+*           initialization of random number generator               *
+*                                                                   *
+*     what (1..4)    values for initialization (= 1..168)           *
+*     what (5..6), sdum    no meaning                               *
+*                                                                   *
+*********************************************************************
+
+  590 CONTINUE
+      IF ((WHAT(1).LT.1.0D0).OR.(WHAT(1).GT.168.0D0)) THEN
+         NA1 = 22
+      ELSE
+         NA1 = WHAT(1)
+      ENDIF
+      IF ((WHAT(2).LT.1.0D0).OR.(WHAT(2).GT.168.0D0)) THEN
+         NA2 = 54
+      ELSE
+         NA2 = WHAT(2)
+      ENDIF
+      IF ((WHAT(3).LT.1.0D0).OR.(WHAT(3).GT.168.0D0)) THEN
+         NA3 = 76
+      ELSE
+         NA3 = WHAT(3)
+      ENDIF
+      IF ((WHAT(4).LT.1.0D0).OR.(WHAT(4).GT.168.0D0)) THEN
+         NA4 = 92
+      ELSE
+         NA4 = WHAT(4)
+      ENDIF
+      CALL DT_RNDMST(NA1,NA2,NA3,NA4)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = LEPTO-CUT                   *
+*                                                                   *
+*          set parameter CUT in LEPTO-common /LEPTOU/               *
+*                                                                   *
+*       what (1) =  index in CUT-array                              *
+*       what (2) =  new value of CUT( int(what(1)) )                *
+*       what (3), what(4) and what (5), what(6) further             *
+*                   parameter in the same way as what (1) and       *
+*                   what (2)                                        *
+*                        default: default-LEPTO parameters          *
+*                                                                   *
+*********************************************************************
+
+  600 CONTINUE
+      IF (WHAT(1).GT.ZERO) CUT(INT(WHAT(1))) = WHAT(2)
+      IF (WHAT(3).GT.ZERO) CUT(INT(WHAT(3))) = WHAT(4)
+      IF (WHAT(5).GT.ZERO) CUT(INT(WHAT(5))) = WHAT(6)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = LEPTO-LST                   *
+*                                                                   *
+*          set parameter LST in LEPTO-common /LEPTOU/               *
+*                                                                   *
+*       what (1) =  index in LST-array                              *
+*       what (2) =  new value of LST( int(what(1)) )                *
+*       what (3), what(4) and what (5), what(6) further             *
+*                   parameter in the same way as what (1) and       *
+*                   what (2)                                        *
+*                        default: default-LEPTO parameters          *
+*                                                                   *
+*********************************************************************
+
+  610 CONTINUE
+      IF (WHAT(1).GT.ZERO) LST(INT(WHAT(1))) = INT(WHAT(2))
+      IF (WHAT(3).GT.ZERO) LST(INT(WHAT(3))) = INT(WHAT(4))
+      IF (WHAT(5).GT.ZERO) LST(INT(WHAT(5))) = INT(WHAT(6))
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = LEPTO-PARL                  *
+*                                                                   *
+*          set parameter PARL in LEPTO-common /LEPTOU/              *
+*                                                                   *
+*       what (1) =  index in PARL-array                             *
+*       what (2) =  new value of PARL( int(what(1)) )               *
+*       what (3), what(4) and what (5), what(6) further             *
+*                   parameter in the same way as what (1) and       *
+*                   what (2)                                        *
+*                        default: default-LEPTO parameters          *
+*                                                                   *
+*********************************************************************
+
+  620 CONTINUE
+      IF (WHAT(1).GT.ZERO) PARL(INT(WHAT(1))) = WHAT(2)
+      IF (WHAT(3).GT.ZERO) PARL(INT(WHAT(3))) = WHAT(4)
+      IF (WHAT(5).GT.ZERO) PARL(INT(WHAT(5))) = WHAT(6)
+      GOTO 10
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = START                       *
+*                                                                   *
+*       what (1) =   number of events                default: 100.  *
+*       what (2) = 0 Glauber initialization follows                 *
+*                = 1 Glauber initialization supressed, fitted       *
+*                    results are used instead                       *
+*                    (this does not apply if emulsion-treatment     *
+*                     is requested)                                 *
+*                = 2 Glauber initialization is written to           *
+*                    output-file shmakov.out                        *
+*                = 3 Glauber initialization is read from input-file *
+*                    shmakov.out                     default: 0     *
+*       what (3..6)  no meaning                                     *
+*       what (3..6)  no meaning                                     *
+*                                                                   *
+*********************************************************************
+
+  630 CONTINUE
+
+* check for cross-section table output only
+      IF (LXSTAB) STOP
+
+      NCASES = INT(WHAT(1))
+      IF (NCASES.LE.0) NCASES = 100
+      IGLAU = INT(WHAT(2))
+      IF ((IGLAU.NE.1).AND.(IGLAU.NE.2).AND.(IGLAU.NE.3))
+     &                                            IGLAU = 0
+
+      NPMASS = IP
+      NPCHAR = IPZ
+      NTMASS = IT
+      NTCHAR = ITZ
+      IDP    = IJPROJ
+      IDT    = IJTARG
+      IF (IDP.LE.0) IDP = 1
+* muon neutrinos: temporary (missing index)
+* (new patch in projpar: therefore the following this is probably not
+*  necessary anymore..)
+C     IF (IDP.EQ.26) IDP = 5
+C     IF (IDP.EQ.27) IDP = 6
+
+* redefine collision energy
+      IF (LEINP) THEN
+         IF (ABS(VAREHI).GT.ZERO) THEN
+            PDUM = ZERO
+            IF (VARELO.LT.EHADLO) VARELO = EHADLO
+            CALL DT_LTINI(IDP,IDT,VARELO,PDUM,VARCLO,1)
+            PDUM = ZERO
+            CALL DT_LTINI(IDP,IDT,VAREHI,PDUM,VARCHI,1)
+         ENDIF
+         CALL DT_LTINI(IDP,IDT,EPN,PPN,CMENER,1)
+      ELSE
+         WRITE(LOUT,1003)
+ 1003    FORMAT(1X,'INIT:   collision energy not defined!',/,
+     &          1X,'              -program stopped-      ')
+         STOP
+      ENDIF
+
+* switch off evaporation (even if requested) if central coll. requ.
+      IF ((ICENTR.EQ.-1).OR.(ICENTR.GT.0).OR.(XSFRAC.LT.0.5D0)) THEN
+         IF (LEVPRT) THEN
+            WRITE(LOUT,1004)
+ 1004       FORMAT(1X,/,'Warning!  Evaporation request rejected since',
+     &             ' central collisions forced.')
+            LEVPRT = .FALSE.
+            LDEEXG = .FALSE.
+            LHEAVY = .FALSE.
+         ENDIF
+      ENDIF
+
+* initialization of evaporation-module
+
+*  initialize evaporation if the code is not used as Fluka event generator
+      IF (ITRSPT.NE.1) THEN
+         CALL BERTTP
+         CALL INCINI
+      ENDIF
+      IF (LEVPRT) LHEAVY = .TRUE.
+
+
+* save the default JETSET-parameter
+      CALL DT_JSPARA(0)
+
+* force use of phojet for g-A
+      IF ((IDP.EQ.7).AND.(MCGENE.NE.3)) MCGENE = 2
+* initialization of nucleon-nucleon event generator
+      IF (MCGENE.EQ.2) CALL DT_PHOINI
+* initialization of LEPTO event generator
+      IF (MCGENE.EQ.3) THEN
+
+         STOP ' This version does not contain LEPTO !'
+
+      ENDIF
+
+* initialization of quasi-elastic neutrino scattering
+      IF (MCGENE.EQ.4) THEN
+         IF (IJPROJ.EQ.5) THEN
+            NEUTYP = 1
+         ELSEIF (IJPROJ.EQ.6) THEN
+            NEUTYP = 2
+         ELSEIF (IJPROJ.EQ.135) THEN
+            NEUTYP = 3
+         ELSEIF (IJPROJ.EQ.136) THEN
+            NEUTYP = 4
+         ELSEIF (IJPROJ.EQ.133) THEN
+            NEUTYP = 5
+         ELSEIF (IJPROJ.EQ.134) THEN
+            NEUTYP = 6
+         ENDIF
+      ENDIF
+
+* normalize fractions of emulsion components
+      IF (NCOMPO.GT.0) THEN
+         SUMFRA = ZERO
+         DO 491 I=1,NCOMPO
+            SUMFRA = SUMFRA+EMUFRA(I)
+  491    CONTINUE
+         IF (SUMFRA.GT.ZERO) THEN
+            DO 492 I=1,NCOMPO
+               EMUFRA(I) = EMUFRA(I)/SUMFRA
+  492       CONTINUE
+         ENDIF
+      ENDIF
+
+* disallow Cronin's multiple scattering for nucleus-nucleus interactions
+      IF ((IP.GT.1).AND.(MKCRON.GT.0)) THEN
+         WRITE(LOUT,1005)
+ 1005    FORMAT(/,1X,'INIT:  multiple scattering disallowed',/)
+         MKCRON = 0
+      ENDIF
+
+* initialization of Glauber-formalism (moved to xAEVT, sr 26.3.96)
+C     IF (NCOMPO.LE.0) THEN
+C        CALL DT_SHMAKI(IP,IPZ,IT,ITZ,IDP,PPN,IGLAU)
+C     ELSE
+C        DO 493 I=1,NCOMPO
+C           CALL DT_SHMAKI(IP,IPZ,IEMUMA(I),IEMUCH(I),IDP,PPN,0)
+C 493    CONTINUE
+C     ENDIF
+
+* pre-tabulation of elastic cross-sections
+      CALL DT_SIGTBL(JDUM,JDUM,DUM,DUM,-1)
+
+      CALL DT_XTIME
+
+      RETURN
+
+*********************************************************************
+*                                                                   *
+*               control card:  codewd = STOP                        *
+*                                                                   *
+*               stop of the event generation                        *
+*                                                                   *
+*       what (1..6)  no meaning                                     *
+*                                                                   *
+*********************************************************************
+
+ 9999 CONTINUE
+      WRITE(LOUT,9000)
+ 9000 FORMAT(1X,'---> unexpected end of input !')
+
+  640 CONTINUE
+      STOP
+
+      END
+*
+*===kkinc==============================================================*
+*
+CDECK  ID>, DT_KKINC
+      SUBROUTINE DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,KKMAT,
+     &                                                         IREJ)
+
+************************************************************************
+* Treatment of complete nucleus-nucleus or hadron-nucleus scattering   *
+* This subroutine is an update of the previous version written         *
+* by J. Ranft/ H.-J. Moehring.                                         *
+* This version dated 19.11.95 is written by S. Roesler                 *
+************************************************************************
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      PARAMETER ( LINP = 5 ,
+     &            LOUT = 6 ,
+     &            LDAT = 9 )
+
+      PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY5=1.0D-5,
+     &           TINY2=1.0D-2,TINY3=1.0D-3)
+
+      LOGICAL LFZC
+
+* event history
+
+      PARAMETER (NMXHKK=200000)
+
+      COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
+     &                JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
+     &                PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
+* extended event history
+      COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
+     &                IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
+     &                IHIST(2,NMXHKK)
+* particle properties (BAMJET index convention)
+      CHARACTER*8  ANAME
+      COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
+     &                IICH(210),IIBAR(210),K1(210),K2(210)
+* properties of interacting particles
+      COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
+* Lorentz-parameters of the current interaction
+      COMMON /DTLTRA/ GACMS(2),BGCMS(2),GALAB,BGLAB,BLAB,
+     &                UMO,PPCM,EPROJ,PPROJ
+* flags for input different options
+      LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
+      COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
+     &                LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
+* flags for particle decays
+      COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
+     &                IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
+     &                NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
+* cuts for variable energy runs
+      COMMON /DTVARE/ VARELO,VAREHI,VARCLO,VARCHI
+* Glauber formalism: flags and parameters for statistics
+      LOGICAL LPROD
+      CHARACTER*8 CGLB
+      COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
+
+      DIMENSION WHAT(6)
+
+      IREJ  = 0
+      ILOOP = 0
+  100 CONTINUE
+      IF (ILOOP.EQ.4) THEN
+         WRITE(LOUT,1000) NEVHKK
+ 1000    FORMAT(1X,'KKINC: event ',I8,' rejected!')
+         GOTO 9999
+      ENDIF
+      ILOOP = ILOOP+1
+
+* variable energy-runs, recalculate parameters for LT's
+      IF ((ABS(VAREHI).GT.ZERO).OR.(IOGLB.EQ.100)) THEN
+         PDUM = ZERO
+         CDUM = ZERO
+         CALL DT_LTINI(IDP,1,EPN,PDUM,CDUM,1)
+      ENDIF
+      IF (EPN.GT.EPROJ) THEN
+         WRITE(LOUT,'(A,E9.3,2A,E9.3,A)')
+     &      ' Requested energy (',EPN,'GeV) exceeds',
+     &      ' initialization energy (',EPROJ,'GeV) !'
+         STOP
+      ENDIF
+
+* re-initialize /DTPRTA/
+      IP  = NPMASS
+      IPZ = NPCHAR
+      IT  = NTMASS
+      ITZ = NTCHAR
+      IJPROJ = IDP
+      IBPROJ = IIBAR(IJPROJ)
+
+* calculate nuclear potentials (common /DTNPOT/)
+      CALL DT_NCLPOT(IPZ,IP,ITZ,IT,ZERO,ZERO,0)
+
+* initialize treatment for residual nuclei
+      CALL DT_RESNCL(EPN,NLOOP,1)
+
+* sample hadron/nucleus-nucleus interaction
+      CALL DT_KKEVNT(KKMAT,IREJ1)
+      IF (IREJ1.GT.0) THEN
+         IF (IOULEV(1).GT.0) WRITE(LOUT,*) 'rejected 1 in KKINC'
+         GOTO 9999
+      ENDIF
+
+      IF ((NPMASS.GT.1).OR.(NTMASS.GT.1)) THEN
+
+* intranuclear cascade of final state particles for KTAUGE generations
+* of secondaries
+         CALL DT_FOZOCA(LFZC,IREJ1)
+         IF (IREJ1.GT.0) THEN
+            IF (IOULEV(1).GT.0) WRITE(LOUT,*) 'rejected 2 in KKINC'
+            GOTO 9999
+         ENDIF
+
+* baryons unable to escape the nuclear potential are treated as
+* excited nucleons (ISTHKK=15,16)
+         CALL DT_SCN4BA
+
+* decay of resonances produced in intranuclear cascade processes
+**sr 15-11-95 should be obsolete
+C        IF (LFZC) CALL DT_DECAY1
+
+  101    CONTINUE
+* treatment of residual nuclei
+         CALL DT_RESNCL(EPN,NLOOP,2)
+
+* evaporation / fission / fragmentation
+* (if intranuclear cascade was sampled only)
+         IF (LFZC) THEN
+            CALL DT_FICONF(IJPROJ,IP,IPZ,IT,ITZ,NLOOP,IREJ1)
+            IF (IREJ1.GT.1) GOTO 101
+            IF (IREJ1.EQ.1) GOTO 100
+         ENDIF
+
+      ENDIF
+
+* transform finale state into Lab.
+      IFLAG = 2
+      CALL DT_BEAMPR(WHAT,DUM,IFLAG)
+      IF ((IFRAME.EQ.1).AND.(IFLAG.EQ.-1)) CALL DT_LT2LAB
+
+      IF (IPI0.EQ.1) CALL DT_DECPI0
+
+C     IF (NEVHKK.EQ.5) CALL DT_EVTOUT(4)
+
+      RETURN
+ 9999 CONTINUE
+      IREJ = 1
+      RETURN
+      END
+*
+*===defaul=============================================================*
+*
+CDECK  ID>, DT_DEFAUL
+      SUBROUTINE DT_DEFAUL(EPN,PPN)
+
+************************************************************************
+* Variables are set to default values.                                 *
+* This version dated 8.5.95 is written by S. Roesler.                  *
+************************************************************************
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+      PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY10=1.0D-10)
+      PARAMETER (TWOPI  = 6.283185307179586454D+00)
+
+* particle properties (BAMJET index convention)
+      CHARACTER*8  ANAME
+      COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
+     &                IICH(210),IIBAR(210),K1(210),K2(210)
+* nuclear potential
+      LOGICAL LFERMI
+      COMMON /DTNPOT/ PFERMP(2),PFERMN(2),FERMOD,
+     &                EBINDP(2),EBINDN(2),EPOT(2,210),
+     &                ETACOU(2),ICOUL,LFERMI
+* interface HADRIN-DPM
+      COMMON /HNTHRE/ EHADTH,EHADLO,EHADHI,INTHAD,IDXTA
+* central particle production, impact parameter biasing
+      COMMON /DTIMPA/ BIMIN,BIMAX,XSFRAC,ICENTR
+* properties of interacting particles
+      COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
+* properties of photon/lepton projectiles
+      COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
+
+      PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
+
+* emulsion treatment
+      COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
+     &                NCOMPO,IEMUL
+* parameter for intranuclear cascade
+      LOGICAL LPAULI
+      COMMON /DTFOTI/ TAUFOR,KTAUGE,ITAUVE,INCMOD,LPAULI
+* various options for treatment of partons (DTUNUC 1.x)
+* (chain recombination, Cronin,..)
+      LOGICAL LCO2CR,LINTPT
+      COMMON /DTCHAI/ SEASQ,CRONCO,CUTOF,MKCRON,ISICHA,IRECOM,
+     &                LCO2CR,LINTPT
+* threshold values for x-sampling (DTUNUC 1.x)
+      COMMON /DTXCUT/ XSEACU,UNON,UNOM,UNOSEA,CVQ,CDQ,CSEA,SSMIMA,
+     &                SSMIMQ,VVMTHR
+* flags for input different options
+      LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
+      COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
+     &                LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
+* n-n cross section fluctuations
+      PARAMETER (NBINS = 1000)
+      COMMON /DTXSFL/ FLUIXX(NBINS),IFLUCT
+* flags for particle decays
+      COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
+     &                IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
+     &                NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
+* diquark-breaking mechanism
+      COMMON /DTDIQB/ DBRKR(3,8),DBRKA(3,8),CHAM1,CHAM3,CHAB1,CHAB3
+* nucleon-nucleon event-generator
+      CHARACTER*8 CMODEL
+      LOGICAL LPHOIN
+      COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
+* flags for diffractive interactions (DTUNUC 1.x)
+      COMMON /DTFLG3/ ISINGD,IDOUBD,IFLAGD,IDIFF
+* VDM parameter for photon-nucleus interactions
+      COMMON /DTVDMP/ RL2,EPSPOL,INTRGE(2),IDPDF,MODEGA,ISHAD(3)
+* Glauber formalism: flags and parameters for statistics
+      LOGICAL LPROD
+      CHARACTER*8 CGLB
+      COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
+* kinematical cuts for lepton-nucleus interactions
+      COMMON /DTLCUT/ ECMIN,ECMAX,XBJMIN,ELMIN,EGMIN,EGMAX,YMIN,YMAX,
+     &                Q2MIN,Q2MAX,THMIN,THMAX,Q2LI,Q2HI,ECMLI,ECMHI
+* flags for activated histograms
+      COMMON /DTHIS3/ IHISPP(50),IHISXS(50),IXSTBL
+* cuts for variable energy runs
+      COMMON /DTVARE/ VARELO,VAREHI,VARCLO,VARCHI
+* parameters for hA-diffraction
+      COMMON /DTDIHA/ DIBETA,DIALPH
+* LEPTO
+      REAL RPPN
+      COMMON /LEPTOI/ RPPN,LEPIN,INTER
+* steering flags for qel neutrino scattering modules
+      COMMON /QNEUTO/ DSIGSU,DSIGMC,NDSIG,NEUTYP,NEUDEC
+* event flag
+      COMMON /DTEVNO/ NEVENT,ICASCA
+
+      DATA POTMES /0.002D0/
+
+* common /DTNPOT/
+      DO 10 I=1,2
+         PFERMP(I) = ZERO
+         PFERMN(I) = ZERO
+         EBINDP(I) = ZERO
+         EBINDN(I) = ZERO
+         DO 11 J=1,210
+            EPOT(I,J) = ZERO
+   11    CONTINUE
+* nucleus independent meson potential
+         EPOT(I,13) = POTMES
+         EPOT(I,14) = POTMES
+         EPOT(I,15) = POTMES
+         EPOT(I,16) = POTMES
+         EPOT(I,23) = POTMES
+         EPOT(I,24) = POTMES
+         EPOT(I,25) = POTMES
+   10 CONTINUE
+**sr 7.4.98: changed after corrected B-sampling
+C     FERMOD    = 0.55D0
+      FERMOD    = 0.68D0
+      ETACOU(1) = ZERO
+      ETACOU(2) = ZERO
+      ICOUL     = 1
+      LFERMI    = .TRUE.
+
+* common /HNTHRE/
+      EHADTH = -99.0D0
+      EHADLO = 4.06D0
+      EHADHI = 6.0D0
+      INTHAD = 1
+      IDXTA  = 2
+
+* common /DTIMPA/
+      ICENTR = 0
+      BIMIN  = ZERO
+      BIMAX  = 1.0D10
+      XSFRAC = 1.0D0
+
+* common /DTPRTA/
+      IP  = 1
+      IPZ = 1
+      IT  = 1
+      ITZ = 1
+      IJPROJ = 1
+      IBPROJ = 1
+      IJTARG = 1
+      IBTARG = 1
+* common /DTGPRO/
+      VIRT = ZERO
+      DO 14 I=1,4
+         PGAMM(I)  = ZERO
+         PLEPT0(I) = ZERO
+         PLEPT1(I) = ZERO
+         PNUCL(I)  = ZERO
+   14 CONTINUE
+      IDIREC   = 0
+
+* common /DTFOTI/
+**sr 7.4.98: changed after corrected B-sampling
+C     TAUFOR = 4.4D0
+      TAUFOR = 3.1D0
+      KTAUGE = 25
+      ITAUVE = 1
+      INCMOD = 1
+      LPAULI = .TRUE.
+
+* common /DTCHAI/
+      SEASQ  = ONE
+      MKCRON = 1
+      CRONCO = 0.64D0
+      ISICHA = 0
+      CUTOF  = 100.0D0
+      LCO2CR = .FALSE.
+      IRECOM = 1
+      LINTPT = .TRUE.
+
+* common /DTXCUT/
+*  definition of soft quark distributions
+      XSEACU = 0.05D0
+      UNON   = 2.0D0
+      UNOM   = 1.5D0
+      UNOSEA = 5.0D0
+*  cutoff parameters for x-sampling
+      CVQ    = 1.0D0
+      CDQ    = 2.0D0
+C     CSEA   = 0.3D0
+      CSEA   = 0.1D0
+      SSMIMA = 1.2D0
+      SSMIMQ = SSMIMA**2
+      VVMTHR = 2.0D0
+
+* common /DTXSFL/
+      IFLUCT = 0
+
+* common /DTFRPA/
+      PDB = 0.15D0
+      PDBSEA(1) = 0.0D0
+      PDBSEA(2) = 0.0D0
+      PDBSEA(3) = 0.0D0
+      ISIG0 = 0
+      IPI0  = 0
+      NMSTU = 0
+      NPARU = 0
+      NMSTJ = 0
+      NPARJ = 0
+
+* common /DTDIQB/
+      DO 15 I=1,8
+         DBRKR(1,I) = 5.0D0
+         DBRKR(2,I) = 5.0D0
+         DBRKR(3,I) = 10.0D0
+         DBRKA(1,I) = ZERO
+         DBRKA(2,I) = ZERO
+         DBRKA(3,I) = ZERO
+   15 CONTINUE
+      CHAM1 = 0.2D0
+      CHAM3 = 0.5D0
+      CHAB1 = 0.7D0
+      CHAB3 = 1.0D0
+
+* common /DTFLG3/
+      ISINGD = 0
+      IDOUBD = 0
+      IFLAGD = 0
+      IDIFF  = 0
+
+* common /DTMODL/
+      MCGENE    = 2
+      CMODEL(1) = 'DTUNUC  '
+      CMODEL(2) = 'PHOJET  '
+      CMODEL(3) = 'LEPTO   '
+      CMODEL(4) = 'QNEUTRIN'
+      LPHOIN    = .TRUE.
+      ELOJET    = 5.0D0
+
+* common /DTLCUT/
+      ECMIN  = 3.5D0
+      ECMAX  = 1.0D10
+      XBJMIN = ZERO
+      ELMIN = ZERO
+      EGMIN = ZERO
+      EGMAX = 1.0D10
+      YMIN  = TINY10
+      YMAX  = 0.999D0
+      Q2MIN = TINY10
+      Q2MAX = 10.0D0
+      THMIN = ZERO
+      THMAX = TWOPI
+      Q2LI  = ZERO
+      Q2HI  = 1.0D10
+      ECMLI = ZERO
+      ECMHI = 1.0D10
+
+* common /DTVDMP/
+      RL2       = 2.0D0
+      INTRGE(1) = 1
+      INTRGE(2) = 3
+      IDPDF     = 2212
+      MODEGA    = 4
+      ISHAD(1)  = 1
+      ISHAD(2)  = 1
+      ISHAD(3)  = 1
+      EPSPOL    = ZERO
+
+* common /DTGLGP/
+      JSTATB = 1000
+      JBINSB = 49
+      CGLB   = '        '
+      IF (ITRSPT.EQ.1) THEN
+         IOGLB  = 100
+      ELSE
+         IOGLB  = 0
+      ENDIF
+      LPROD  = .TRUE.
+
+* common /DTHIS3/
+      DO 16 I=1,50
+         IHISPP(I) = 0
+         IHISXS(I) = 0
+   16 CONTINUE
+      IXSTBL = 0
+
+* common /DTVARE/
+      VARELO = ZERO
+      VAREHI = ZERO
+      VARCLO = ZERO
+      VARCHI = ZERO
+
+* common /DTDIHA/
+      DIBETA = -1.0D0
+      DIALPH = ZERO
+
+* common /LEPTOI/
+      RPPN  = 0.0
+      LEPIN = 0
+      INTER = 0
+
+* common /QNEUTO/
+      NEUTYP = 1
+      NEUDEC = 0
+
+* common /DTEVNO/
+      NEVENT = 1
+      IF (ITRSPT.EQ.1) THEN
+         ICASCA = 1
+      ELSE
+         ICASCA = 0
+      ENDIF
+
+* default Lab.-energy
+      EPN = 200.0D0
+      PPN = SQRT((EPN-AAM(IJPROJ))*(EPN+AAM(IJPROJ)))
+
+      RETURN
+      END
+*
+*===aaevt==============================================================*
+*
+CDECK  ID>, DT_AAEVT
+      SUBROUTINE DT_AAEVT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
+     &                                             IDP,IGLAU)
+
+************************************************************************
+* This version dated 22.03.96 is written by S. Roesler.                *
+************************************************************************
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      PARAMETER ( LINP = 5 ,
+     &            LOUT = 6 ,
+     &            LDAT = 9 )
+
+      PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
+
+* emulsion treatment
+      COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
+     &                NCOMPO,IEMUL
+* event flag
+      COMMON /DTEVNO/ NEVENT,ICASCA
+
+      CHARACTER*8 DATE,HHMMSS
+      DIMENSION IDMNYR(3)
+
+      KKMAT  = 1
+      NMSG   = MAX(NEVTS/100,1)
+
+* initialization of run-statistics and histograms
+      CALL DT_STATIS(1)
+
+      CALL PHO_PHIST(1000,DUM)
+
+* initialization of Glauber-formalism
+      IF (NCOMPO.LE.0) THEN
+         CALL DT_SHMAKI(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,IGLAU)
+      ELSE
+         DO 1 I=1,NCOMPO
+            CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
+    1    CONTINUE
+      ENDIF
+      CALL DT_SIGEMU
+
+      CALL IDATE(IDMNYR)
+      WRITE(DATE,'(I2,''/'',I2,''/'',I2)')
+     &   IDMNYR(1),IDMNYR(2),MOD(IDMNYR(3),100)
+      CALL ITIME(IDMNYR)
+      WRITE(HHMMSS,'(I2,'':'',I2,'':'',I2)')
+     &   IDMNYR(1),IDMNYR(2),IDMNYR(3)
+      WRITE(LOUT,1001) DATE,HHMMSS
+ 1001 FORMAT(/,' DT_AAEVT: Initialisation finished. ( Date: ',A8,
+     &       '   Time: ',A8,' )')
+
+* generate NEVTS events
+      DO 2 IEVT=1,NEVTS
+
+*  print run-status message
+         IF (MOD(IEVT,NMSG).EQ.0) THEN
+            CALL IDATE(IDMNYR)
+            WRITE(DATE,'(I2,''/'',I2,''/'',I2)')
+     &         IDMNYR(1),IDMNYR(2),MOD(IDMNYR(3),100)
+            CALL ITIME(IDMNYR)
+            WRITE(HHMMSS,'(I2,'':'',I2,'':'',I2)')
+     &         IDMNYR(1),IDMNYR(2),IDMNYR(3)
+            WRITE(LOUT,1000) IEVT-1,NEVTS,DATE,HHMMSS
+ 1000       FORMAT(/,1X,I8,' out of ',I8,' events sampled ( Date: ',A,
+     &             '   Time: ',A,' )',/)
+C           WRITE(LOUT,1000) IEVT-1
+C1000       FORMAT(1X,I8,' events sampled')
+         ENDIF
+         NEVENT = IEVT
+*  treat nuclear emulsions
+         IF (IEMUL.GT.0) CALL DT_GETEMU(NTMASS,NTCHAR,KKMAT,0)
+*  composite targets only
+         KKMAT = -KKMAT
+*  sample this event
+         CALL DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,KKMAT,IREJ)
+
+         CALL PHO_PHIST(2000,DUM)
+
+    2 CONTINUE
+
+* print run-statistics and histograms to output-unit 6
+
+      CALL PHO_PHIST(3000,DUM)
+
+      CALL DT_STATIS(2)
+
+      RETURN
+      END
+*
+*===laevt==============================================================*
+*
+CDECK  ID>, DT_LAEVT
+      SUBROUTINE DT_LAEVT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
+     &                                             IDP,IGLAU)
+
+************************************************************************
+* Interface to run DPMJET for lepton-nucleus interactions.             *
+* Kinematics is sampled using the equivalent photon approximation      *
+* Based on GPHERA-routine by R. Engel.                                 *
+* This version dated 23.03.96 is written by S. Roesler.                *
+************************************************************************
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      PARAMETER ( LINP = 5 ,
+     &            LOUT = 6 ,
+     &            LDAT = 9 )
+
+      PARAMETER (TINY10=1.0D-10,TINY4=1.0D-4,
+     &           ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,THREE=3.0D0)
+      PARAMETER (TWOPI  = 6.283185307179586454D+00,
+     &           PI     = TWOPI/TWO,
+     &           ALPHEM = ONE/137.0D0)
+
+C     CHARACTER*72 HEADER
+
+* particle properties (BAMJET index convention)
+      CHARACTER*8  ANAME
+      COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
+     &                IICH(210),IIBAR(210),K1(210),K2(210)
+* event history
+
+      PARAMETER (NMXHKK=200000)
+
+      COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
+     &                JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
+     &                PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
+* extended event history
+      COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
+     &                IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
+     &                IHIST(2,NMXHKK)
+* kinematical cuts for lepton-nucleus interactions
+      COMMON /DTLCUT/ ECMIN,ECMAX,XBJMIN,ELMIN,EGMIN,EGMAX,YMIN,YMAX,
+     &                Q2MIN,Q2MAX,THMIN,THMAX,Q2LI,Q2HI,ECMLI,ECMHI
+* properties of interacting particles
+      COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
+* properties of photon/lepton projectiles
+      COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
+* kinematics at lepton-gamma vertex
+      COMMON /DTLGVX/ PPL0(4),PPL1(4),PPG(4),PPA(4)
+* flags for activated histograms
+      COMMON /DTHIS3/ IHISPP(50),IHISXS(50),IXSTBL
+
+      PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
+
+* emulsion treatment
+      COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
+     &                NCOMPO,IEMUL
+* Glauber formalism: cross sections
+      COMMON /DTGLXS/ ECMNN(NEB),Q2G(NQB),ECMNOW,Q2,
+     &                XSTOT(NEB,NQB,NCOMPX),XSELA(NEB,NQB,NCOMPX),
+     &                XSQEP(NEB,NQB,NCOMPX),XSQET(NEB,NQB,NCOMPX),
+     &                XSQE2(NEB,NQB,NCOMPX),XSPRO(NEB,NQB,NCOMPX),
+     &                XSDEL(NEB,NQB,NCOMPX),XSDQE(NEB,NQB,NCOMPX),
+     &                XETOT(NEB,NQB,NCOMPX),XEELA(NEB,NQB,NCOMPX),
+     &                XEQEP(NEB,NQB,NCOMPX),XEQET(NEB,NQB,NCOMPX),
+     &                XEQE2(NEB,NQB,NCOMPX),XEPRO(NEB,NQB,NCOMPX),
+     &                XEDEL(NEB,NQB,NCOMPX),XEDQE(NEB,NQB,NCOMPX),
+     &                BSLOPE,NEBINI,NQBINI
+* nucleon-nucleon event-generator
+      CHARACTER*8 CMODEL
+      LOGICAL LPHOIN
+      COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
+* flags for input different options
+      LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
+      COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
+     &                LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
+* event flag
+      COMMON /DTEVNO/ NEVENT,ICASCA
+
+      DIMENSION XDUMB(40),BGTA(4)
+
+* LEPTO
+      IF (MCGENE.EQ.3) THEN
+
+         STOP ' This version does not contain LEPTO !'
+
+      ENDIF
+
+      KKMAT  = 1
+      NMSG   = MAX(NEVTS/10,1)
+
+* mass of incident lepton
+      AMLPT  = AAM(IDP)
+      AMLPT2 = AMLPT**2
+      IDPPDG = IDT_IPDGHA(IDP)
+
+* consistency of kinematical limits
+      Q2MIN  = MAX(Q2MIN,TINY10)
+      Q2MAX  = MAX(Q2MAX,TINY10)
+      YMIN   = MIN(MAX(YMIN,TINY10),0.999D0)
+      YMAX   = MIN(MAX(YMAX,TINY10),0.999D0)
+
+* total energy of the lepton-nucleon system
+      PTOTLN = SQRT( (PLEPT0(1)+PNUCL(1))**2+(PLEPT0(2)+PNUCL(2))**2
+     &                                      +(PLEPT0(3)+PNUCL(3))**2 )
+      ETOTLN = PLEPT0(4)+PNUCL(4)
+      ECMLN  = SQRT((ETOTLN-PTOTLN)*(ETOTLN+PTOTLN))
+      ECMAX  = MIN(ECMAX,ECMLN)
+      WRITE(LOUT,1003) ECMIN,ECMAX,YMIN,YMAX,Q2MIN,Q2MAX,EGMIN,
+     &                 THMIN,THMAX,ELMIN
+ 1003 FORMAT(1X,'LAEVT:',16X,'kinematical cuts',/,22X,
+     &       '------------------',/,9X,'W (min)   =',
+     &       F7.1,' GeV    (max) =',F7.1,' GeV',/,9X,'y (min)   =',
+     &       F7.3,8X,'(max) =',F7.3,/,9X,'Q^2 (min) =',F7.1,
+     &       ' GeV^2  (max) =',F7.1,' GeV^2',/,' (Lab)   E_g (min) ='
+     &       ,F7.1,' GeV',/,' (Lab) theta (min) =',F7.4,8X,'(max) =',
+     &       F7.4,'   for E_lpt >',F7.1,' GeV',/)
+
+* Lorentz-parameter for transf. into Lab
+      BGTA(1) = PNUCL(1)/AAM(1)
+      BGTA(2) = PNUCL(2)/AAM(1)
+      BGTA(3) = PNUCL(3)/AAM(1)
+      BGTA(4) = PNUCL(4)/AAM(1)
+* LT of incident lepton into Lab and dump it in DTEVT1
+      CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
+     &            PLEPT0(1),PLEPT0(2),PLEPT0(3),PLEPT0(4),
+     &            PLTOT,PPL0(1),PPL0(2),PPL0(3),PPL0(4))
+      CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
+     &            PNUCL(1),PNUCL(2),PNUCL(3),PNUCL(4),
+     &            PLTOT,PPA(1),PPA(2),PPA(3),PPA(4))
+* maximum energy of photon nucleon system
+      PTOTGN = SQRT((YMAX*PPL0(1)+PPA(1))**2+(YMAX*PPL0(2)+PPA(2))**2
+     &                                      +(YMAX*PPL0(3)+PPA(3))**2)
+      ETOTGN = YMAX*PPL0(4)+PPA(4)
+      EGNMAX = SQRT((ETOTGN-PTOTGN)*(ETOTGN+PTOTGN))
+      EGNMAX = MIN(EGNMAX,ECMAX)
+* minimum energy of photon nucleon system
+      PTOTGN = SQRT((YMIN*PPL0(1)+PPA(1))**2+(YMIN*PPL0(2)+PPA(2))**2
+     &                                      +(YMIN*PPL0(3)+PPA(3))**2)
+      ETOTGN = YMIN*PPL0(4)+PPA(4)
+      EGNMIN = SQRT((ETOTGN-PTOTGN)*(ETOTGN+PTOTGN))
+      EGNMIN = MAX(EGNMIN,ECMIN)
+
+* limits for Glauber-initialization
+      Q2LI  = Q2MIN
+      Q2HI  = MAX(Q2LI,MIN(Q2HI,Q2MAX))
+      ECMLI = MAX(EGNMIN,THREE)
+      ECMHI = EGNMAX
+      WRITE(LOUT,1004) EGNMIN,EGNMAX,ECMLI,ECMHI,Q2LI,Q2HI
+ 1004 FORMAT(1X,'resulting limits:',/,9X,'W (min)   =',F7.1,
+     &       ' GeV    (max) =',F7.1,' GeV',/,/,' limits for ',
+     &       'Glauber-initialization:',/,9X,'W (min)   =',F7.1,
+     &       ' GeV    (max) =',F7.1,' GeV',/,9X,'Q^2 (min) =',F7.1,
+     &       ' GeV^2  (max) =',F7.1,' GeV^2',/)
+* initialization of Glauber-formalism
+      IF (NCOMPO.LE.0) THEN
+         CALL DT_SHMAKI(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,IGLAU)
+      ELSE
+         DO 9 I=1,NCOMPO
+            CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
+    9    CONTINUE
+      ENDIF
+      CALL DT_SIGEMU
+
+* initialization of run-statistics and histograms
+      CALL DT_STATIS(1)
+
+      CALL PHO_PHIST(1000,DUM)
+
+* maximum photon-nucleus cross section
+      I1  = 1
+      I2  = 1
+      RAT = ONE
+      IF (EGNMAX.GE.ECMNN(NEBINI)) THEN
+         I1  = NEBINI
+         I2  = NEBINI
+         RAT = ONE
+      ELSEIF (EGNMAX.GT.ECMNN(1)) THEN
+         DO 5 I=2,NEBINI
+            IF (EGNMAX.LT.ECMNN(I)) THEN
+               I1  = I-1
+               I2  = I
+               RAT = (EGNMAX-ECMNN(I1))/(ECMNN(I2)-ECMNN(I1))
+               GOTO 6
+            ENDIF
+    5    CONTINUE
+    6    CONTINUE
+      ENDIF
+      SIGMAX = XSTOT(I1,1,1)+RAT*(XSTOT(I2,1,1)-XSTOT(I1,1,1))
+      EGNXX  = EGNMAX
+      I1  = 1
+      I2  = 1
+      RAT = ONE
+      IF (EGNMIN.GE.ECMNN(NEBINI)) THEN
+         I1  = NEBINI
+         I2  = NEBINI
+         RAT = ONE
+      ELSEIF (EGNMIN.GT.ECMNN(1)) THEN
+         DO 7 I=2,NEBINI
+            IF (EGNMIN.LT.ECMNN(I)) THEN
+               I1  = I-1
+               I2  = I
+               RAT = (EGNMIN-ECMNN(I1))/(ECMNN(I2)-ECMNN(I1))
+               GOTO 8
+            ENDIF
+    7    CONTINUE
+    8    CONTINUE
+      ENDIF
+      SIGXX = XSTOT(I1,1,1)+RAT*(XSTOT(I2,1,1)-XSTOT(I1,1,1))
+      IF (SIGXX.GT.SIGMAX) EGNXX = EGNMIN
+      SIGMAX = MAX(SIGMAX,SIGXX)
+      WRITE(LOUT,'(9X,A,F8.3,A)') 'Sigma_tot (max) =',SIGMAX,' mb'
+
+* plot photon flux table
+      AYMIN = LOG(YMIN)
+      AYMAX = LOG(YMAX)
+      AYRGE = AYMAX-AYMIN
+      MAXTAB = 50
+      ADY    = LOG(YMAX/YMIN)/DBLE(MAXTAB-1)
+C     WRITE(LOUT,'(/,1X,A)') 'LAEVT:   photon flux '
+      DO 1 I=1,MAXTAB
+         Y     = EXP(AYMIN+ADY*DBLE(I-1))
+         Q2LOW = MAX(Q2MIN,AMLPT2*Y**2/(ONE-Y))
+         FF1   = ALPHEM/TWOPI * ((ONE+(ONE-Y)**2)/Y*LOG(Q2MAX/Q2LOW)
+     &                           -TWO*AMLPT2*Y*(ONE/Q2LOW-ONE/Q2MAX))
+         FF2   = ALPHEM/TWOPI * ((ONE+(ONE-Y)**2)/Y*LOG(Q2MAX/Q2LOW)
+     &                           -TWO*(ONE-Y)/Y*(ONE-Q2LOW/Q2MAX))
+C        WRITE(LOUT,'(5X,3E15.4)') Y,FF1,FF2
+    1 CONTINUE
+
+* maximum residual weight for flux sampling (dy/y)
+      YY     = YMIN
+      Q2LOW  = MAX(Q2MIN,AMLPT2*YY**2/(ONE-YY))
+      WGHMAX = (ONE+(ONE-YY)**2)*LOG(Q2MAX/Q2LOW)
+     &         -TWO*AMLPT2*YY*(ONE/Q2LOW-ONE/Q2MAX)*YY
+
+      CALL DT_NEWHGR(YMIN,YMAX,ZERO,XDUMB,49,IHFLY0)
+      CALL DT_NEWHGR(YMIN,YMAX,ZERO,XDUMB,49,IHFLY1)
+      CALL DT_NEWHGR(YMIN,YMAX,ZERO,XDUMB,49,IHFLY2)
+      CALL DT_NEWHGR(Q2LOW,Q2MAX,ZERO,XDUMB,20,IHFLQ0)
+      CALL DT_NEWHGR(Q2LOW,Q2MAX,ZERO,XDUMB,20,IHFLQ1)
+      CALL DT_NEWHGR(Q2LOW,Q2MAX,ZERO,XDUMB,20,IHFLQ2)
+      CALL DT_NEWHGR(EGNMIN,EGNMAX,ZERO,XDUMB,20,IHFLE0)
+      CALL DT_NEWHGR(EGNMIN,EGNMAX,ZERO,XDUMB,20,IHFLE1)
+      CALL DT_NEWHGR(EGNMIN,EGNMAX,ZERO,XDUMB,20,IHFLE2)
+      CALL DT_NEWHGR(ZERO,EGMAX,ZERO,XDUMB,20,IHFLU0)
+      CALL DT_NEWHGR(ZERO,EGMAX,ZERO,XDUMB,20,IHFLU1)
+      CALL DT_NEWHGR(ZERO,EGMAX,ZERO,XDUMB,20,IHFLU2)
+      XBLOW = 0.001D0
+      CALL DT_NEWHGR(XBLOW,ONE,ZERO,XDUMB,-40,IHFLX0)
+      CALL DT_NEWHGR(XBLOW,ONE,ZERO,XDUMB,-40,IHFLX1)
+      CALL DT_NEWHGR(XBLOW,ONE,ZERO,XDUMB,-40,IHFLX2)
+
+      ITRY = 0
+      ITRW = 0
+      NC0  = 0
+      NC1  = 0
+
+* generate events
+      DO 2 IEVT=1,NEVTS
+         IF (MOD(IEVT,NMSG).EQ.0) THEN
+C           OPEN(LDAT,FILE='/scrtch3/hr/sroesler/statusd5.out',
+C    &                                         STATUS='UNKNOWN')
+            WRITE(LOUT,'(1X,I8,A)') IEVT-1,' events sampled'
+C           CLOSE(LDAT)
+         ENDIF
+         NEVENT = IEVT
+
+  100    CONTINUE
+         ITRY = ITRY+1
+
+*  sample y
+  101    CONTINUE
+         ITRW  = ITRW+1
+         YY    = EXP(AYRGE*DT_RNDM(RAT)+AYMIN)
+         Q2LOW = MAX(Q2MIN,AMLPT2*YY**2/(ONE-YY))
+         Q2LOG = LOG(Q2MAX/Q2LOW)
+         WGH   = (ONE+(ONE-YY)**2)*Q2LOG
+     &           -TWO*AMLPT2*YY*(ONE/Q2LOW-ONE/Q2MAX)*YY
+         IF (WGHMAX.LT.WGH) WRITE(LOUT,1000) YY,WGHMAX,WGH
+ 1000    FORMAT(1X,'LAEVT:   weight error!',3E12.5)
+         IF (DT_RNDM(YY)*WGHMAX.GT.WGH) GOTO 101
+
+*  sample Q2
+         YEFF = ONE+(ONE-YY)**2
+  102    CONTINUE
+         Q2  = Q2LOW*EXP(Q2LOG*DT_RNDM(YY))
+         WGH = (YEFF-TWO*(ONE-YY)*Q2LOW/Q2)/YEFF
+         IF (WGH.LT.DT_RNDM(Q2)) GOTO 102
+
+c        NC0 = NC0+1
+c        CALL DT_FILHGR(YY,ONE,IHFLY0,NC0)
+c        CALL DT_FILHGR(Q2,ONE,IHFLQ0,NC0)
+
+*  kinematics at lepton-photon vertex
+*   scattered electron
+         YQ2 = SQRT((ONE-YY)*Q2)
+         Q2E = Q2/(4.0D0*PLEPT0(4))
+         E1Y = (ONE-YY)*PLEPT0(4)
+         CALL DT_DSFECF(SIF,COF)
+         PLEPT1(1) = YQ2*COF
+         PLEPT1(2) = YQ2*SIF
+         PLEPT1(3) = E1Y-Q2E
+         PLEPT1(4) = E1Y+Q2E
+C        THETA = ACOS( (E1Y-Q2E)/(E1Y+Q2E) )
+*   radiated photon
+         PGAMM(1) = -PLEPT1(1)
+         PGAMM(2) = -PLEPT1(2)
+         PGAMM(3) = PLEPT0(3)-PLEPT1(3)
+         PGAMM(4) = PLEPT0(4)-PLEPT1(4)
+*   E_cm cut
+         PTOTGN = SQRT( (PGAMM(1)+PNUCL(1))**2+(PGAMM(2)+PNUCL(2))**2
+     &                                        +(PGAMM(3)+PNUCL(3))**2 )
+         ETOTGN = PGAMM(4)+PNUCL(4)
+         ECMGN  = (ETOTGN-PTOTGN)*(ETOTGN+PTOTGN)
+         IF (ECMGN.LT.0.1D0) GOTO 101
+         ECMGN  = SQRT(ECMGN)
+         IF ((ECMGN.LT.ECMIN).OR.(ECMGN.GT.ECMAX)) GOTO 101
+
+*  Lorentz-transformation into nucleon-rest system
+         CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
+     &               PGAMM(1),PGAMM(2),PGAMM(3),PGAMM(4),
+     &               PGTOT,PPG(1),PPG(2),PPG(3),PPG(4))
+         CALL DT_DALTRA(BGTA(4),-BGTA(1),-BGTA(2),-BGTA(3),
+     &               PLEPT1(1),PLEPT1(2),PLEPT1(3),PLEPT1(4),
+     &               PLTOT,PPL1(1),PPL1(2),PPL1(3),PPL1(4))
+*  temporary checks..
+         Q2TMP = ABS(PPG(4)**2-PGTOT**2)
+         IF (ABS(Q2-Q2TMP).GT.0.01D0) WRITE(LOUT,1001) Q2,Q2TMP
+ 1001    FORMAT(1X,'LAEVT:    inconsistent kinematics (Q2,Q2TMP) ',
+     &          2F10.4)
+         ECMTMP = SQRT((PPG(4)+AAM(1)-PGTOT)*(PPG(4)+AAM(1)+PGTOT))
+         IF (ABS(ECMGN-ECMTMP).GT.TINY10) WRITE(LOUT,1002) ECMGN,ECMTMP
+ 1002    FORMAT(1X,'LAEVT:    inconsistent kinematics (ECMGN,ECMTMP) ',
+     &          2F10.2)
+         YYTMP = PPG(4)/PPL0(4)
+         IF (ABS(YY-YYTMP).GT.0.01D0) WRITE(LOUT,1005) YY,YYTMP
+ 1005    FORMAT(1X,'LAEVT:    inconsistent kinematics (YY,YYTMP) ',
+     &          2F10.4)
+
+*  lepton tagger (Lab)
+         THETA = ACOS( PPL1(3)/PLTOT )
+         IF (PPL1(4).GT.ELMIN) THEN
+            IF ((THETA.LT.THMIN).OR.(THETA.GT.THMAX)) GOTO 101
+         ENDIF
+*  photon energy-cut (Lab)
+         IF (PPG(4).LT.EGMIN) GOTO 101
+         IF (PPG(4).GT.EGMAX) GOTO 101
+*   x_Bj cut
+         XBJ = ABS(Q2/(1.876D0*PPG(4)))
+         IF (XBJ.LT.XBJMIN) GOTO 101
+
+         NC0 = NC0+1
+         CALL DT_FILHGR(    Q2,ONE,IHFLQ0,NC0)
+         CALL DT_FILHGR(    YY,ONE,IHFLY0,NC0)
+         CALL DT_FILHGR(   XBJ,ONE,IHFLX0,NC0)
+         CALL DT_FILHGR(PPG(4),ONE,IHFLU0,NC0)
+         CALL DT_FILHGR( ECMGN,ONE,IHFLE0,NC0)
+
+*  rotation angles against z-axis
+         COD = PPG(3)/PGTOT
+C        SID = SQRT((ONE-COD)*(ONE+COD))
+         PPT = SQRT(PPG(1)**2+PPG(2)**2)
+         SID = PPT/PGTOT
+         COF = ONE
+         SIF = ZERO
+         IF (PGTOT*SID.GT.TINY10) THEN
+            COF   = PPG(1)/(SID*PGTOT)
+            SIF   = PPG(2)/(SID*PGTOT)
+            ANORF = SQRT(COF*COF+SIF*SIF)
+            COF   = COF/ANORF
+            SIF   = SIF/ANORF
+         ENDIF
+
+         IF (IXSTBL.EQ.0) THEN
+*  change to photon projectile
+            IJPROJ = 7
+*  set virtuality
+            VIRT = Q2
+*  re-initialize LTs with new kinematics
+*  !!PGAMM ist set in cms (ECMGN) along z
+            EPN = ZERO
+            PPN = ZERO
+            CALL DT_LTINI(IJPROJ,IJTARG,EPN,PPN,ECMGN,0)
+*  Introduced by Chiara -> force CMS-system
+*            IFRAME = 2
+*  to force Lab-system
+            IFRAME = 1
+*  get emulsion component if requested
+            IF (IEMUL.GT.0) CALL DT_GETEMU(NTMASS,NTCHAR,KKMAT,0)
+*  convolute with cross section
+            CALL DT_SIGGAT(Q2LOW,EGNXX,STOTX,KKMAT)
+            CALL DT_SIGGAT(Q2,ECMGN,STOT,KKMAT)
+            IF (STOTX.LT.STOT) WRITE(LOUT,'(1X,A,/,6E12.3)')
+     &         'LAEVT: warning STOTX<STOT ! ',Q2LOW,EGNMAX,STOTX,
+     &                                        Q2,ECMGN,STOT
+            IF (DT_RNDM(Q2)*STOTX.GT.STOT) GOTO 100
+            NC1 = NC1+1
+            CALL DT_FILHGR(    Q2,ONE,IHFLQ1,NC1)
+            CALL DT_FILHGR(    YY,ONE,IHFLY1,NC1)
+            CALL DT_FILHGR(   XBJ,ONE,IHFLX1,NC1)
+            CALL DT_FILHGR(PPG(4),ONE,IHFLU1,NC1)
+            CALL DT_FILHGR( ECMGN,ONE,IHFLE1,NC1)
+*  composite targets only
+            KKMAT = -KKMAT
+*  sample this event
+            CALL DT_KKINC(NPMASS,NPCHAR,NTMASS,NTCHAR,IJPROJ,EPN,KKMAT,
+     &                                                            IREJ)
+*  rotate momenta of final state particles back in photon-nucleon syst.
+            DO 4 I=NPOINT(4),NHKK
+               IF ((ABS(ISTHKK(I)).EQ.1).OR.(ISTHKK(I).EQ.1000).OR.
+     &                                      (ISTHKK(I).EQ.1001)) THEN
+                  PX = PHKK(1,I)
+                  PY = PHKK(2,I)
+                  PZ = PHKK(3,I)
+                  CALL DT_MYTRAN(1,PX,PY,PZ,COD,SID,COF,SIF,
+     &                        PHKK(1,I),PHKK(2,I),PHKK(3,I))
+               ENDIF
+    4       CONTINUE
+         ENDIF
+
+         CALL DT_FILHGR(    Q2,ONE,IHFLQ2,NC1)
+         CALL DT_FILHGR(    YY,ONE,IHFLY2,NC1)
+         CALL DT_FILHGR(   XBJ,ONE,IHFLX2,NC1)
+         CALL DT_FILHGR(PPG(4),ONE,IHFLU2,NC1)
+         CALL DT_FILHGR( ECMGN,ONE,IHFLE2,NC1)
+
+*  dump this event to histograms
+
+         CALL PHO_PHIST(2000,DUM)
+
+    2 CONTINUE
+
+      WGY    = ALPHEM/TWOPI*WGHMAX*DBLE(ITRY)/DBLE(ITRW)
+      WGY    = WGY*LOG(YMAX/YMIN)
+      WEIGHT = WGY*SIGMAX*DBLE(NEVTS)/DBLE(ITRY)
+
+C     HEADER = ' LAEVT:  Q^2 distribution 0'
+C     CALL DT_OUTHGR(IHFLQ0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  Q^2 distribution 1'
+C     CALL DT_OUTHGR(IHFLQ1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  Q^2 distribution 2'
+C     CALL DT_OUTHGR(IHFLQ2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  y   distribution 0'
+C     CALL DT_OUTHGR(IHFLY0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  y   distribution 1'
+C     CALL DT_OUTHGR(IHFLY1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  y   distribution 2'
+C     CALL DT_OUTHGR(IHFLY2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  x   distribution 0'
+C     CALL DT_OUTHGR(IHFLX0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  x   distribution 1'
+C     CALL DT_OUTHGR(IHFLX1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  x   distribution 2'
+C     CALL DT_OUTHGR(IHFLX2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  E_g distribution 0'
+C     CALL DT_OUTHGR(IHFLU0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  E_g distribution 1'
+C     CALL DT_OUTHGR(IHFLU1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  E_g distribution 2'
+C     CALL DT_OUTHGR(IHFLU2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  E_c distribution 0'
+C     CALL DT_OUTHGR(IHFLE0,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  E_c distribution 1'
+C     CALL DT_OUTHGR(IHFLE1,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+C     HEADER = ' LAEVT:  E_c distribution 2'
+C     CALL DT_OUTHGR(IHFLE2,0,0,0,0,0,HEADER,0,NEVTS,ONE,1,1,-1)
+
+* print run-statistics and histograms to output-unit 6
+
+      CALL PHO_PHIST(3000,DUM)
+
+      IF (IXSTBL.EQ.0) CALL DT_STATIS(2)
+
+      RETURN
+      END
+*
+*===dtuini=============================================================*
+*
+CDECK  ID>, DT_DTUINI
+      SUBROUTINE DT_DTUINI(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,
+     &                                               IDP,IEMU)
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
+
+* emulsion treatment
+      COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
+     &                NCOMPO,IEMUL
+* Glauber formalism: flags and parameters for statistics
+      LOGICAL LPROD
+      CHARACTER*8 CGLB
+      COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
+
+      CALL DT_INIT(NEVTS,EPN,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,IGLAU)
+      CALL DT_STATIS(1)
+
+      CALL PHO_PHIST(1000,DUM)
+
+      IF (NCOMPO.LE.0) THEN
+         CALL DT_SHMAKI(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPN,IGLAU)
+      ELSE
+         DO 1 I=1,NCOMPO
+            CALL DT_SHMAKI(NPMASS,NPCHAR,IEMUMA(I),IEMUCH(I),IDP,EPN,0)
+    1    CONTINUE
+      ENDIF
+      IF (IOGLB.NE.100) CALL DT_SIGEMU
+      IEMU = IEMUL
+
+      RETURN
+      END
+*
+*===dtuout=============================================================*
+*
+CDECK  ID>, DT_DTUOUT
+      SUBROUTINE DT_DTUOUT
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      CALL PHO_PHIST(3000,DUM)
+
+      CALL DT_STATIS(2)
+
+      RETURN
+      END
+*
+*===beam===============================================================*
+*
+CDECK  ID>, DT_BEAMPR
+      SUBROUTINE DT_BEAMPR(WHAT,PLAB,MODE)
+
+************************************************************************
+* Initialization of event generation                                   *
+* This version dated  7.4.98  is written by S. Roesler.                *
+************************************************************************
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      PARAMETER ( LINP = 5 ,
+     &            LOUT = 6 ,
+     &            LDAT = 9 )
+
+      PARAMETER (ZERO=0.0D0,ONE=1.0D0,TINY10=1.0D-10)
+      PARAMETER (TWOPI=6.283185307D0,BOG=TWOPI/360.0D0)
+
+      LOGICAL LBEAM
+
+* event history
+
+      PARAMETER (NMXHKK=200000)
+
+      COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
+     &                JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
+     &                PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
+* extended event history
+      COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
+     &                IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
+     &                IHIST(2,NMXHKK)
+* properties of interacting particles
+      COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
+* particle properties (BAMJET index convention)
+      CHARACTER*8  ANAME
+      COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
+     &                IICH(210),IIBAR(210),K1(210),K2(210)
+* beam momenta
+      COMMON /DTBEAM/ P1(4),P2(4)
+
+C     DIMENSION WHAT(6),P1(4),P2(4),P1CMS(4),P2CMS(4)
+      DIMENSION WHAT(6),P1CMS(4),P2CMS(4)
+
+      DATA LBEAM /.FALSE./
+
+      GOTO (1,2) MODE
+
+    1 CONTINUE
+
+      E1  = WHAT(1)
+      IF (E1.LT.ZERO) E1 = DBLE(IPZ)/DBLE(IP)*ABS(WHAT(1))
+      E2  = WHAT(2)
+      IF (E2.LT.ZERO) E2 = DBLE(ITZ)/DBLE(IT)*ABS(WHAT(2))
+      PP1 = SQRT( (E1+AAM(IJPROJ))*(E1-AAM(IJPROJ)) )
+      PP2 = SQRT( (E2+AAM(IJTARG))*(E2-AAM(IJTARG)) )
+      TH  = 1.D-6*WHAT(3)/2.D0
+      PH  = WHAT(4)*BOG
+      P1(1) = PP1*SIN(TH)*COS(PH)
+      P1(2) = PP1*SIN(TH)*SIN(PH)
+      P1(3) = PP1*COS(TH)
+      P1(4) = E1
+      P2(1) = PP2*SIN(TH)*COS(PH)
+      P2(2) = PP2*SIN(TH)*SIN(PH)
+      P2(3) = -PP2*COS(TH)
+      P2(4) = E2
+      ECM  = SQRT( (P1(4)+P2(4))**2-(P1(1)+P2(1))**2-(P1(2)+P2(2))**2
+     &                                              -(P1(3)+P2(3))**2 )
+      ELAB = (ECM**2-AAM(IJPROJ)**2-AAM(IJTARG)**2)/(2.0D0*AAM(IJTARG))
+      PLAB = SQRT( (ELAB+AAM(IJPROJ))*(ELAB-AAM(IJPROJ)) )
+      BGX  = (P1(1)+P2(1))/ECM
+      BGY  = (P1(2)+P2(2))/ECM
+      BGZ  = (P1(3)+P2(3))/ECM
+      BGE  = (P1(4)+P2(4))/ECM
+      CALL DT_DALTRA(BGE,-BGX,-BGY,-BGZ,P1(1),P1(2),P1(3),P1(4),
+     &            P1TOT,P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4))
+      CALL DT_DALTRA(BGE,-BGX,-BGY,-BGZ,P2(1),P2(2),P2(3),P2(4),
+     &            P2TOT,P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4))
+      COD = P1CMS(3)/P1TOT
+C     SID = SQRT((ONE-COD)*(ONE+COD))
+      PPT = SQRT(P1CMS(1)**2+P1CMS(2)**2)
+      SID = PPT/P1TOT
+      COF = ONE
+      SIF = ZERO
+      IF (P1TOT*SID.GT.TINY10) THEN
+         COF   = P1CMS(1)/(SID*P1TOT)
+         SIF   = P1CMS(2)/(SID*P1TOT)
+         ANORF = SQRT(COF*COF+SIF*SIF)
+         COF   = COF/ANORF
+         SIF   = SIF/ANORF
+      ENDIF
+**check
+C     WRITE(LOUT,'(4E15.4)') P1(1),P1(2),P1(3),P1(4)
+C     WRITE(LOUT,'(4E15.4)') P2(1),P2(2),P2(3),P2(4)
+C     WRITE(LOUT,'(5E15.4)') P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4),P1TOT
+C     WRITE(LOUT,'(5E15.4)') P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4),P2TOT
+C     PAX = ZERO
+C     PAY = ZERO
+C     PAZ = P1TOT
+C     PAE = SQRT(AAM(IJPROJ)**2+PAZ**2)
+C     PBX = ZERO
+C     PBY = ZERO
+C     PBZ = -P2TOT
+C     PBE = SQRT(AAM(IJTARG)**2+PBZ**2)
+C     WRITE(LOUT,'(4E15.4)') PAX,PAY,PAZ,PAE
+C     WRITE(LOUT,'(4E15.4)') PBX,PBY,PBZ,PBE
+C     CALL DT_MYTRAN(1,PAX,PAY,PAZ,COD,SID,COF,SIF,
+C    &            P1CMS(1),P1CMS(2),P1CMS(3))
+C     CALL DT_MYTRAN(1,PBX,PBY,PBZ,COD,SID,COF,SIF,
+C    &            P2CMS(1),P2CMS(2),P2CMS(3))
+C     WRITE(LOUT,'(4E15.4)') P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4)
+C     WRITE(LOUT,'(4E15.4)') P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4)
+C     CALL DT_DALTRA(BGE,BGX,BGY,BGZ,P1CMS(1),P1CMS(2),P1CMS(3),P1CMS(4),
+C    &            P1TOT,P1(1),P1(2),P1(3),P1(4))
+C     CALL DT_DALTRA(BGE,BGX,BGY,BGZ,P2CMS(1),P2CMS(2),P2CMS(3),P2CMS(4),
+C    &            P2TOT,P2(1),P2(2),P2(3),P2(4))
+C     WRITE(LOUT,'(4E15.4)') P1(1),P1(2),P1(3),P1(4)
+C     WRITE(LOUT,'(4E15.4)') P2(1),P2(2),P2(3),P2(4)
+C     STOP
+**
+
+      LBEAM = .TRUE.
+
+      RETURN
+
+    2 CONTINUE
+
+      IF (LBEAM) THEN
+         IF ( (NPOINT(4).EQ.0).OR.(NHKK.LT.NPOINT(4)) ) RETURN
+         DO 20 I=NPOINT(4),NHKK
+            IF ((ABS(ISTHKK(I)).EQ.1).OR.(ISTHKK(I).EQ.1000).OR.
+     &                                   (ISTHKK(I).EQ.1001)) THEN
+               CALL DT_MYTRAN(1,PHKK(1,I),PHKK(2,I),PHKK(3,I),
+     &                     COD,SID,COF,SIF,PXCMS,PYCMS,PZCMS)
+               PECMS = PHKK(4,I)
+               CALL DT_DALTRA(BGE,BGX,BGY,BGZ,PXCMS,PYCMS,PZCMS,PECMS,
+     &                     PTOT,PHKK(1,I),PHKK(2,I),PHKK(3,I),PHKK(4,I))
+            ENDIF
+   20    CONTINUE
+      ELSE
+         MODE = -1
+      ENDIF
+
+      RETURN
+      END
+*
+*===eventb=============================================================*
+*
+CDECK  ID>, DT_EVENTB
+      SUBROUTINE DT_EVENTB(NCSY,IREJ)
+
+************************************************************************
+* Treatment of nucleon-nucleon interactions with full two-component    *
+* Dual Parton Model.                                                   *
+*          NCSY     number of nucleon-nucleon interactions             *
+*          IREJ     rejection flag                                     *
+* This version dated 14.01.2000 is written by S. Roesler               *
+************************************************************************
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      PARAMETER ( LINP = 5 ,
+     &            LOUT = 6 ,
+     &            LDAT = 9 )
+
+      PARAMETER (TINY10=1.0D-10,ZERO=0.0D0,ONE=1.0D0)
+
+* event history
+
+      PARAMETER (NMXHKK=200000)
+
+      COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
+     &                JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
+     &                PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
+* extended event history
+      COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
+     &                IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
+     &                IHIST(2,NMXHKK)
+*! uncomment this line for internal phojet-fragmentation
+C #include "dtu_dtevtp.inc"
+* particle properties (BAMJET index convention)
+      CHARACTER*8  ANAME
+      COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
+     &                IICH(210),IIBAR(210),K1(210),K2(210)
+* flags for input different options
+      LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
+      COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
+     &                LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
+* rejection counter
+      COMMON /DTREJC/ IRPT,IRHHA,IRRES(2),LOMRES,LOBRES,
+     &                IRCHKI(2),IRFRAG,IRCRON(3),IREVT,
+     &                IREXCI(3),IRDIFF(2),IRINC
+* properties of interacting particles
+      COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
+* properties of photon/lepton projectiles
+      COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
+* various options for treatment of partons (DTUNUC 1.x)
+* (chain recombination, Cronin,..)
+      LOGICAL LCO2CR,LINTPT
+      COMMON /DTCHAI/ SEASQ,CRONCO,CUTOF,MKCRON,ISICHA,IRECOM,
+     &                LCO2CR,LINTPT
+* statistics
+      COMMON /DTSTA1/ ICREQU,ICSAMP,ICCPRO,ICDPR,ICDTA,
+     &                ICRJSS,ICVV2S,ICCHAI(2,9),ICRES(9),ICDIFF(5),
+     &                ICEVTG(8,0:30)
+* DTUNUC-PHOJET interface, Lorentz-param. of n-n subsystem
+      COMMON /DTLTSU/ BGX,BGY,BGZ,GAM
+* Glauber formalism: collision properties
+      COMMON /DTGLCP/ RPROJ,RTARG,BIMPAC,
+     &                NWTSAM,NWASAM,NWBSAM,NWTACC,NWAACC,NWBACC
+* flags for diffractive interactions (DTUNUC 1.x)
+      COMMON /DTFLG3/ ISINGD,IDOUBD,IFLAGD,IDIFF
+* statistics: double-Pomeron exchange
+      COMMON /DTFLG2/ INTFLG,IPOPO
+* flags for particle decays
+      COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
+     &                IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
+     &                NMSTU,NPARU,NMSTJ,NPARJ,PDB,PDBSEA(3),ISIG0,IPI0
+* nucleon-nucleon event-generator
+      CHARACTER*8 CMODEL
+      LOGICAL LPHOIN
+      COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
+C  nucleon-nucleus / nucleus-nucleus interface to DPMJET
+      INTEGER IDEQP,IDEQB,IHFLD,IHFLS
+      DOUBLE PRECISION ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB
+      COMMON /POHDFL/ ECMN,PCMN,SECM,SPCM,XPSUB,XTSUB,
+     &                IDEQP(2),IDEQB(2),IHFLD(2,2),IHFLS(2)
+C  model switches and parameters
+      CHARACTER*8 MDLNA
+      INTEGER ISWMDL,IPAMDL
+      DOUBLE PRECISION PARMDL
+      COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
+C  initial state parton radiation (internal part)
+      INTEGER MXISR3,MXISR4
+      PARAMETER ( MXISR3 = 50, MXISR4 = 100 )
+      INTEGER IFL1,IFL2,IBRA,IFANO,ISH,NACC
+      DOUBLE PRECISION Q2SH,PT2SH,XPSH,ZPSH,THSH,SHAT
+      COMMON /POINT6/ Q2SH(2,MXISR3),PT2SH(2,MXISR3),XPSH(2,MXISR3),
+     &                ZPSH(2,MXISR3),THSH(2,MXISR3),SHAT(MXISR3),
+     &                IFL1(2,MXISR3),IFL2(2,MXISR3),
+     &                IBRA(2,MXISR4),IFANO(2),ISH(2),NACC
+C  event debugging information
+      INTEGER NMAXD
+      PARAMETER (NMAXD=100)
+      INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
+     &        KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
+      COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
+     &                KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
+C  general process information
+      INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
+      COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
+
+      DIMENSION PP(4),PT(4),PTOT(4),PP1(4),PP2(4),PT1(4),PT2(4),
+     &          PPNN(4),PTNN(4),PTOTNN(4),PPSUB(4),PTSUB(4),
+     &          PPTCMS(4),PTTCMS(4),PPTMP(4),PTTMP(4),
+     &          KPRON(15),ISINGL(2000)
+
+* initial values for max. number of phojet scatterings and dtunuc chains
+* to be fragmented with one pyexec call
+      DATA MXPHFR,MXDTFR /10,100/
+
+      IREJ      = 0
+* pointer to first parton of the first chain in dtevt common
+      NPOINT(3) = NHKK+1
+* special flag for double-Pomeron statistics
+      IPOPO = 1
+* counter for low-mass (DTUNUC) interactions
+      NDTUSC = 0
+* counter for interactions treated by PHOJET
+      NPHOSC = 0
+
+* scan interactions for single nucleon-nucleon interactions
+* (this has to be checked here because Cronin modifies parton momenta)
+      NC = NPOINT(2)
+      IF (NCSY.GT.2000) STOP ' DT_EVENTB: NCSY > 2000 ! '
+      DO 8 I=1,NCSY
+         ISINGL(I) = 0
+         MOP = JMOHKK(1,NC)
+         MOT = JMOHKK(1,NC+1)
+         DIFF1 = ABS(PHKK(4,MOP)-PHKK(4,  NC)-PHKK(4,NC+2))
+         DIFF2 = ABS(PHKK(4,MOT)-PHKK(4,NC+1)-PHKK(4,NC+3))
+         IF ((DIFF1.LT.TINY10).AND.(DIFF2.LT.TINY10)) ISINGL(I) = 1
+         NC = NC+4
+    8 CONTINUE
+
+* multiple scattering of chain ends
+      IF ((IP.GT.1).AND.(MKCRON.NE.0)) CALL DT_CRONIN(1)
+      IF ((IT.GT.1).AND.(MKCRON.NE.0)) CALL DT_CRONIN(2)
+
+* switch to PHOJET-settings for JETSET parameter
+      CALL DT_INITJS(1)
+
+* loop over nucleon-nucleon interaction
+      NC = NPOINT(2)
+      DO 2 I=1,NCSY
+*
+*   pick up one nucleon-nucleon interaction from DTEVT1
+*     ppnn  / ptnn   - momenta of the interacting nucleons (cms)
+*     ptotnn         - total momentum of the interacting nucleons (cms)
+*     pp1,2 / pt1,2  - momenta of the four partons
+*     pp    / pt     - total momenta of the proj / targ partons
+*     ptot           - total momentum of the four partons
+         MOP = JMOHKK(1,NC)
+         MOT = JMOHKK(1,NC+1)
+         DO 3 K=1,4
+            PPNN(K)   = PHKK(K,MOP)
+            PTNN(K)   = PHKK(K,MOT)
+            PTOTNN(K) = PPNN(K)+PTNN(K)
+            PP1(K)    = PHKK(K,NC)
+            PT1(K)    = PHKK(K,NC+1)
+            PP2(K)    = PHKK(K,NC+2)
+            PT2(K)    = PHKK(K,NC+3)
+            PP(K)     = PP1(K)+PP2(K)
+            PT(K)     = PT1(K)+PT2(K)
+            PTOT(K)   = PP(K)+PT(K)
+    3    CONTINUE
+*
+*-----------------------------------------------------------------------
+*   this is a complete nucleon-nucleon interaction
+*
+         IF (ISINGL(I).EQ.1) THEN
+*
+*     initialize PHOJET-variables for remnant/valence-partons
+            IHFLD(1,1) = 0
+            IHFLD(1,2) = 0
+            IHFLD(2,1) = 0
+            IHFLD(2,2) = 0
+            IHFLS(1) = 1
+            IHFLS(2) = 1
+*     save current settings of PHOJET process and min. bias flags
+            DO 9 K=1,11
+               KPRON(K) = IPRON(K,1)
+    9       CONTINUE
+            ISWSAV   = ISWMDL(2)
+*
+*     check if forced sampling of diffractive interaction requested
+            IF (ISINGD.LT.-1) THEN
+               DO 90 K=1,11
+                  IPRON(K,1) = 0
+   90          CONTINUE
+               IF ((ISINGD.EQ.-2).OR.(ISINGD.EQ.-3)) IPRON(5,1) = 1
+               IF ((ISINGD.EQ.-2).OR.(ISINGD.EQ.-4)) IPRON(6,1) = 1
+               IF (ISINGD.EQ.-5) IPRON(4,1) = 1
+            ENDIF
+*
+*     for photons: a direct/anomalous interaction is not sampled
+*     in PHOJET but already in Glauber-formalism. Here we check if such
+*     an interaction is requested
+            IF (IJPROJ.EQ.7) THEN
+*       first switch off direct interactions
+               IPRON(8,1) = 0
+*       this is a direct interactions
+               IF (IDIREC.EQ.1) THEN
+                  DO 12 K=1,11
+                     IPRON(K,1) = 0
+   12             CONTINUE
+                  IPRON(8,1) = 1
+*       this is an anomalous interactions
+*         (iswmdl(2) = 0 only hard int. generated ( = 1 min. bias) )
+               ELSEIF (IDIREC.EQ.2) THEN
+                  ISWMDL(2) = 0
+               ENDIF
+            ELSE
+               IF (IDIREC.NE.0) STOP ' DT_EVENTB: IDIREC > 0 ! '
+            ENDIF
+*
+*     make sure that total momenta of partons, pp and pt, are on mass
+*     shell (Cronin may have srewed this up..)
+            CALL DT_MASHEL(PP,PT,PHKK(5,MOP),PHKK(5,MOT),PPNN,PTNN,IR1)
+            IF (IR1.NE.0) THEN
+               IF (IOULEV(1).GT.0) WRITE(LOUT,'(1X,A)')
+     &              'EVENTB:  mass shell correction rejected'
+               GOTO 9999
+            ENDIF
+*
+*     initialize the incoming particles in PHOJET
+            IF ((IP.EQ.1).AND.(IJPROJ.EQ.7)) THEN
+
+               CALL PHO_SETPAR(1,22,0,VIRT)
+
+            ELSE
+
+               CALL PHO_SETPAR(1,IDHKK(MOP),0,ZERO)
+
+            ENDIF
+
+            CALL PHO_SETPAR(2,IDHKK(MOT),0,ZERO)
+
+*
+*     initialize rejection loop counter for anomalous processes
+            IRJANO = 0
+  800       CONTINUE
+            IRJANO = IRJANO+1
+*
+*     temporary fix for ifano problem
+            IFANO(1) = 0
+            IFANO(2) = 0
+*
+*     generate complete hadron/nucleon/photon-nucleon event with PHOJET
+
+            CALL PHO_EVENT(2,PPNN,PTNN,DUM,IREJ1)
+
+*
+*     for photons: special consistency check for anomalous interactions
+            IF (IJPROJ.EQ.7) THEN
+               IF (IRJANO.LT.30) THEN
+                  IF (IFANO(1).NE.0) THEN
+*       here, an anomalous interaction was generated. Check if it
+*       was also requested. Otherwise reject this event.
+                     IF (IDIREC.EQ.0) GOTO 800
+                  ELSE
+*       here, an anomalous interaction was not generated. Check if it
+*       was requested in which case we need to reject this event.
+                     IF (IDIREC.EQ.2) GOTO 800
+                  ENDIF
+               ELSE
+                  WRITE(LOUT,*) ' DT_EVENTB: Warning! IRJANO > 30 ',
+     &                          IRJANO,IDIREC,NEVHKK
+               ENDIF
+            ENDIF
+*
+*     copy back original settings of PHOJET process and min. bias flags
+            DO 10 K=1,11
+               IPRON(K,1) = KPRON(K)
+   10       CONTINUE
+            ISWMDL(2) = ISWSAV
+*
+*     check if PHOJET has rejected this event
+            IF (IREJ1.NE.0) THEN
+C              IF (IOULEV(1).GT.0) WRITE(LOUT,'(1X,A,I4)')
+               WRITE(LOUT,'(1X,A,I4)')
+     &            'EVENTB:  chain system rejected',IDIREC
+
+               CALL PHO_PREVNT(0)
+
+               GOTO 9999
+            ENDIF
+*
+*     copy partons and strings from PHOJET common back into DTEVT for
+*     external fragmentation
+            MO1 = NC
+            MO2 = NC+3
+*!      uncomment this line for internal phojet-fragmentation
+C           CALL DT_GETFSP(MO1,MO2,PPNN,PTNN,-1)
+            NPHOSC = NPHOSC+1
+            CALL DT_GETPJE(MO1,MO2,PPNN,PTNN,-1,NPHOSC,IREJ1)
+            IF (IREJ1.NE.0) THEN
+               IF (IOULEV(1).GT.0)
+     &         WRITE(LOUT,'(1X,A,I4)') 'EVENTB: chain system rejected 1'
+               GOTO 9999
+            ENDIF
+*
+*     update statistics counter
+            ICEVTG(IDCH(NC),29) = ICEVTG(IDCH(NC),29)+1
+*
+*-----------------------------------------------------------------------
+*   this interaction involves "remnants"
+*
+         ELSE
+*
+*     total mass of this system
+            PPTOT  = SQRT(PTOT(1)**2+PTOT(2)**2+PTOT(3)**2)
+            AMTOT2 = (PTOT(4)-PPTOT)*(PTOT(4)+PPTOT)
+            IF (AMTOT2.LT.ZERO) THEN
+               AMTOT = ZERO
+            ELSE
+               AMTOT = SQRT(AMTOT2)
+            ENDIF
+*
+*     systems with masses larger than elojet are treated with PHOJET
+            IF (AMTOT.GT.ELOJET) THEN
+*
+*     initialize PHOJET-variables for remnant/valence-partons
+*       projectile parton flavors and valence flag
+               IHFLD(1,1) = IDHKK(NC)
+               IHFLD(1,2) = IDHKK(NC+2)
+               IHFLS(1)   = 0
+               IF ((IDCH(NC).EQ.6).OR.(IDCH(NC).EQ.7)
+     &                            .OR.(IDCH(NC).EQ.8)) IHFLS(1) = 1
+*       target parton flavors and valence flag
+               IHFLD(2,1) = IDHKK(NC+1)
+               IHFLD(2,2) = IDHKK(NC+3)
+               IHFLS(2)   = 0
+               IF ((IDCH(NC).EQ.4).OR.(IDCH(NC).EQ.5)
+     &                            .OR.(IDCH(NC).EQ.8)) IHFLS(2) = 1
+*       flag signalizing PHOJET how to treat the remnant:
+*         iremn = -1 sea-quark remnant: PHOJET takes flavors from ihfld
+*         iremn > -1 valence remnant: PHOJET assumes flavors according
+*                    to mother particle
+               IREMN1 = IHFLS(1)-1
+               IREMN2 = IHFLS(2)-1
+*
+*     initialize the incoming particles in PHOJET
+               IF ((IP.EQ.1).AND.(IJPROJ.EQ.7)) THEN
+
+                  CALL PHO_SETPAR(1,22,IREMN1,VIRT)
+
+               ELSE
+
+                  CALL PHO_SETPAR(1,IDHKK(MOP),IREMN1,ZERO)
+
+               ENDIF
+
+               CALL PHO_SETPAR(2,IDHKK(MOT),IREMN2,ZERO)
+
+*
+*     calculate Lorentz parameter of the nucleon-nucleon cm-system
+               PPTOTN = SQRT(PTOTNN(1)**2+PTOTNN(2)**2+PTOTNN(3)**2)
+               AMNN   = SQRT( (PTOTNN(4)-PPTOTN)*(PTOTNN(4)+PPTOTN) )
+               BGX    = PTOTNN(1)/AMNN
+               BGY    = PTOTNN(2)/AMNN
+               BGZ    = PTOTNN(3)/AMNN
+               GAM    = PTOTNN(4)/AMNN
+*     transform interacting nucleons into nucleon-nucleon cm-system
+               CALL DT_DALTRA(GAM,-BGX,-BGY,-BGZ,
+     &                     PPNN(1),PPNN(2),PPNN(3),PPNN(4),PPCMS,
+     &                     PPTCMS(1),PPTCMS(2),PPTCMS(3),PPTCMS(4))
+               CALL DT_DALTRA(GAM,-BGX,-BGY,-BGZ,
+     &                     PTNN(1),PTNN(2),PTNN(3),PTNN(4),PTCMS,
+     &                     PTTCMS(1),PTTCMS(2),PTTCMS(3),PTTCMS(4))
+*     transform (total) momenta of the proj and targ partons into
+*     nucleon-nucleon cm-system
+               CALL DT_DALTRA(GAM,-BGX,-BGY,-BGZ,
+     &                     PP(1),PP(2),PP(3),PP(4),
+     &                     PPTSUB,PPSUB(1),PPSUB(2),PPSUB(3),PPSUB(4))
+               CALL DT_DALTRA(GAM,-BGX,-BGY,-BGZ,
+     &                     PT(1),PT(2),PT(3),PT(4),
+     &                     PTTSUB,PTSUB(1),PTSUB(2),PTSUB(3),PTSUB(4))
+*     energy fractions of the proj and targ partons
+               XPSUB = MIN(PPSUB(4)/PPTCMS(4),ONE)
+               XTSUB = MIN(PTSUB(4)/PTTCMS(4),ONE)
+***
+* testprint
+c              PTOTCM = SQRT( (PPTCMS(1)+PTTCMS(1))**2 +
+c    &                        (PPTCMS(2)+PTTCMS(2))**2 +
+c    &                        (PPTCMS(3)+PTTCMS(3))**2 )
+c              EOLDCM = SQRT( (PPTCMS(4)+PTTCMS(4)-PTOTCM) *
+c    &                        (PPTCMS(4)+PTTCMS(4)+PTOTCM) )
+c              PTOTSU = SQRT( (PPSUB(1)+PTSUB(1))**2 +
+c    &                        (PPSUB(2)+PTSUB(2))**2 +
+c    &                        (PPSUB(3)+PTSUB(3))**2 )
+c              EOLDSU = SQRT( (PPSUB(4)+PTSUB(4)-PTOTSU) *
+c    &                        (PPSUB(4)+PTSUB(4)+PTOTSU) )
+***
+*
+*     save current settings of PHOJET process and min. bias flags
+               DO 7 K=1,11
+                  KPRON(K) = IPRON(K,1)
+    7          CONTINUE
+*     disallow direct photon int. (does not make sense here anyway)
+               IPRON(8,1) = 0
+*     disallow double pomeron processes (due to technical problems
+*     in PHOJET, needs to be solved sometime)
+               IPRON(4,1) = 0
+*     disallow diffraction for sea-diquarks
+               IF ((IABS(IHFLD(1,1)).GT.1100).AND.
+     &             (IABS(IHFLD(1,2)).GT.1100)) THEN
+                  IPRON(3,1) = 0
+                  IPRON(6,1) = 0
+               ENDIF
+               IF ((IABS(IHFLD(2,1)).GT.1100).AND.
+     &             (IABS(IHFLD(2,2)).GT.1100)) THEN
+                  IPRON(3,1) = 0
+                  IPRON(5,1) = 0
+               ENDIF
+*
+*     we need massless partons: transform them on mass shell
+               XMP = ZERO
+               XMT = ZERO
+               DO 6 K=1,4
+                  PPTMP(K) = PPSUB(K)
+                  PTTMP(K) = PTSUB(K)
+    6          CONTINUE
+               CALL DT_MASHEL(PPTMP,PTTMP,XMP,XMT,PPSUB,PTSUB,IREJ1)
+               PPSUTO  = SQRT(PPSUB(1)**2+PPSUB(2)**2+PPSUB(3)**2)
+               PTSUTO  = SQRT(PTSUB(1)**2+PTSUB(2)**2+PTSUB(3)**2)
+               PSUTOT = SQRT((PPSUB(1)+PTSUB(1))**2+
+     &                  (PPSUB(2)+PTSUB(2))**2+(PPSUB(3)+PTSUB(3))**2)
+*     total energy of the subsysten after mass transformation
+*      (should be the same as before..)
+               SECM = SQRT( (PPSUB(4)+PTSUB(4)-PSUTOT)*
+     &                      (PPSUB(4)+PTSUB(4)+PSUTOT) )
+*
+*     after mass shell transformation the x_sub - relation has to be
+*     corrected. We therefore create "pseudo-momenta" of mother-nucleons.
+*
+*     The old version was to scale based on the original x_sub and the
+*     4-momenta of the subsystem. At very high energy this could lead to
+*     "pseudo-cm energies" of the parent system considerably exceeding
+*     the true cm energy. Now we keep the true cm energy and calculate
+*     new x_sub instead.
+C old version  PPTCMS(4) = PPSUB(4)/XPSUB
+               PPTCMS(4) = MAX(PPTCMS(4),PPSUB(4))
+               XPSUB = PPSUB(4)/PPTCMS(4)
+               IF (IJPROJ.EQ.7) THEN
+                  AMP2  = PHKK(5,MOT)**2
+                  PTOT1 = SQRT(PPTCMS(4)**2-AMP2)
+               ELSE
+*???????
+                  PTOT1 = SQRT((PPTCMS(4)-PHKK(5,MOP))
+     &                        *(PPTCMS(4)+PHKK(5,MOP)))
+C                 PTOT1 = SQRT((PPTCMS(4)-PHKK(5,MOT))
+C    &                        *(PPTCMS(4)+PHKK(5,MOT)))
+               ENDIF
+C old version  PTTCMS(4) = PTSUB(4)/XTSUB
+               PTTCMS(4) = MAX(PTTCMS(4),PTSUB(4))
+               XTSUB = PTSUB(4)/PTTCMS(4)
+               PTOT2 = SQRT((PTTCMS(4)-PHKK(5,MOT))
+     &                     *(PTTCMS(4)+PHKK(5,MOT)))
+               DO 4 K=1,3
+                  PPTCMS(K) = PTOT1*PPSUB(K)/PPSUTO
+                  PTTCMS(K) = PTOT2*PTSUB(K)/PTSUTO
+    4          CONTINUE
+***
+* testprint
+*
+*     ppnn  / ptnn   - momenta of the int. nucleons (cms, negl. Fermi)
+*     ptotnn         - total momentum of the int. nucleons (cms, negl. Fermi)
+*     pptcms/ pttcms - momenta of the interacting nucleons (cms)
+*     pp1,2 / pt1,2  - momenta of the four partons
+*
+*     pp    / pt     - total momenta of the pr/ta partons (cms, negl. Fermi)
+*     ptot           - total momentum of the four partons (cms, negl. Fermi)
+*     ppsub / ptsub  - total momenta of the proj / targ partons (cms)
+*
+c              PTOTCM = SQRT( (PPTCMS(1)+PTTCMS(1))**2 +
+c    &                        (PPTCMS(2)+PTTCMS(2))**2 +
+c    &                        (PPTCMS(3)+PTTCMS(3))**2 )
+c              ENEWCM = SQRT( (PPTCMS(4)+PTTCMS(4)-PTOTCM) *
+c    &                        (PPTCMS(4)+PTTCMS(4)+PTOTCM) )
+c              PTOTSU = SQRT( (PPSUB(1)+PTSUB(1))**2 +
+c    &                        (PPSUB(2)+PTSUB(2))**2 +
+c    &                        (PPSUB(3)+PTSUB(3))**2 )
+c              ENEWSU = SQRT( (PPSUB(4)+PTSUB(4)-PTOTSU) *
+c    &                        (PPSUB(4)+PTSUB(4)+PTOTSU) )
+c              IF (ENEWCM/EOLDCM.GT.1.1D0) THEN
+c                 WRITE(*,*) ' EOLDCM, ENEWCM : ',EOLDCM,ENEWCM
+c                 WRITE(*,*) ' EOLDSU, ENEWSU : ',EOLDSU,ENEWSU
+c                 WRITE(*,*) ' XPSUB,  XTSUB  : ',XPSUB,XTSUB
+c              ENDIF
+c              BBGX = (PPTCMS(1)+PTTCMS(1))/ENEWCM
+c              BBGY = (PPTCMS(2)+PTTCMS(2))/ENEWCM
+c              BBGZ = (PPTCMS(3)+PTTCMS(3))/ENEWCM
+c              BGAM = (PPTCMS(4)+PTTCMS(4))/ENEWCM
+*     transform interacting nucleons into nucleon-nucleon cm-system
+c              CALL DT_DALTRA(BGAM,-BBGX,-BBGY,-BBGZ,
+c    &                    PPTCMS(1),PPTCMS(2),PPTCMS(3),PPTCMS(4),PPTOT,
+c    &                     PPNEW1,PPNEW2,PPNEW3,PPNEW4)
+c              CALL DT_DALTRA(BGAM,-BBGX,-BBGY,-BBGZ,
+c    &                    PTTCMS(1),PTTCMS(2),PTTCMS(3),PTTCMS(4),PTTOT,
+c    &                     PTNEW1,PTNEW2,PTNEW3,PTNEW4)
+c              CALL DT_DALTRA(BGAM,-BBGX,-BBGY,-BBGZ,
+c    &                     PPSUB(1),PPSUB(2),PPSUB(3),PPSUB(4),PPTOT,
+c    &                     PPSUB1,PPSUB2,PPSUB3,PPSUB4)
+c              CALL DT_DALTRA(BGAM,-BBGX,-BBGY,-BBGZ,
+c    &                     PTSUB(1),PTSUB(2),PTSUB(3),PTSUB(4),PTTOT,
+c    &                     PTSUB1,PTSUB2,PTSUB3,PTSUB4)
+c              PTSTCM = SQRT( (PPNEW1+PTNEW1)**2 +
+c    &                        (PPNEW2+PTNEW2)**2 +
+c    &                        (PPNEW3+PTNEW3)**2 )
+c              ETSTCM = SQRT( (PPNEW4+PTNEW4-PTSTCM) *
+c    &                        (PPNEW4+PTNEW4+PTSTCM) )
+c              PTSTSU = SQRT( (PPSUB1+PTSUB1)**2 +
+c    &                        (PPSUB2+PTSUB2)**2 +
+c    &                        (PPSUB3+PTSUB3)**2 )
+c              ETSTSU = SQRT( (PPSUB4+PTSUB4-PTSTSU) *
+c    &                        (PPSUB4+PTSUB4+PTSTSU) )
+C              WRITE(*,*) ' mother cmE :'
+C              WRITE(*,*) ETSTCM,ENEWCM
+C              WRITE(*,*) ' subsystem cmE :'
+C              WRITE(*,*) ETSTSU,ENEWSU
+C              WRITE(*,*) ' projectile mother :'
+C              WRITE(*,*) PPNEW1,PPNEW2,PPNEW3,PPNEW4
+C              WRITE(*,*) ' target mother :'
+C              WRITE(*,*) PTNEW1,PTNEW2,PTNEW3,PTNEW4
+C              WRITE(*,*) ' projectile subsystem:'
+C              WRITE(*,*) PPSUB1,PPSUB2,PPSUB3,PPSUB4
+C              WRITE(*,*) ' target subsystem:'
+C              WRITE(*,*) PTSUB1,PTSUB2,PTSUB3,PTSUB4
+C              WRITE(*,*) ' projectile subsystem should be:'
+C              WRITE(*,*) ZERO,ZERO,XPSUB*ETSTCM/2.0D0,
+C    &                    XPSUB*ETSTCM/2.0D0
+C              WRITE(*,*) ' target subsystem should be:'
+C              WRITE(*,*) ZERO,ZERO,-XTSUB*ETSTCM/2.0D0,
+C    &                    XTSUB*ETSTCM/2.0D0
+C              WRITE(*,*) ' subsystem cmE should be: '
+C              WRITE(*,*) SQRT(XPSUB*XTSUB)*ETSTCM,XPSUB,XTSUB
+***
+*
+*     generate complete remnant - nucleon/remnant event with PHOJET
+
+               CALL PHO_EVENT(3,PPTCMS,PTTCMS,DUM,IREJ1)
+
+*
+*     copy back original settings of PHOJET process flags
+               DO 11 K=1,11
+                  IPRON(K,1) = KPRON(K)
+   11          CONTINUE
+*
+*     check if PHOJET has rejected this event
+               IF (IREJ1.NE.0) THEN
+                  IF (IOULEV(1).GT.0)
+     &            WRITE(LOUT,'(1X,A)') 'EVENTB:  chain system rejected'
+                  WRITE(LOUT,*)
+     &                 'XPSUB,XTSUB,SECM ',XPSUB,XTSUB,SECM,AMTOT
+
+                  CALL PHO_PREVNT(0)
+
+                  GOTO 9999
+               ENDIF
+*
+*     copy partons and strings from PHOJET common back into DTEVT for
+*     external fragmentation
+               MO1 = NC
+               MO2 = NC+3
+*!      uncomment this line for internal phojet-fragmentation
+C              CALL DT_GETFSP(MO1,MO2,PP,PT,1)
+               NPHOSC = NPHOSC+1
+               CALL DT_GETPJE(MO1,MO2,PP,PT,1,NPHOSC,IREJ1)
+               IF (IREJ1.NE.0) THEN
+                  IF (IOULEV(1).GT.0) WRITE(LOUT,'(1X,A,I4)')
+     &               'EVENTB: chain system rejected 2'
+                  GOTO 9999
+               ENDIF
+*
+*     update statistics counter
+               ICEVTG(IDCH(NC),2) = ICEVTG(IDCH(NC),2)+1
+*
+*-----------------------------------------------------------------------
+* two-chain approx. for smaller systems
+*
+            ELSE
+*
+               NDTUSC = NDTUSC+1
+*   special flag for double-Pomeron statistics
+               IPOPO = 0
+*
+*   pick up flavors at the ends of the two chains
+               IFP1 = IDHKK(NC)
+               IFT1 = IDHKK(NC+1)
+               IFP2 = IDHKK(NC+2)
+               IFT2 = IDHKK(NC+3)
+*   ..and the indices of the mothers
+               MOP1 = NC
+               MOT1 = NC+1
+               MOP2 = NC+2
+               MOT2 = NC+3
+               CALL DT_GETCSY(IFP1,PP1,MOP1,IFP2,PP2,MOP2,
+     &                     IFT1,PT1,MOT1,IFT2,PT2,MOT2,IREJ1)
+*
+*   check if this chain system was rejected
+               IF (IREJ1.GT.0) THEN
+                  IF (IOULEV(1).GT.0) THEN
+                     WRITE(LOUT,*) 'rejected 1 in EVENTB'
+                     WRITE(LOUT,'(1X,4(I6,4E12.3,/),E12.3)')
+     &                  IFP1,PP1,IFT1,PT1,IFP2,PP2,IFT2,PT2,AMTOT
+                  ENDIF
+                  IRHHA = IRHHA+1
+                  GOTO 9999
+               ENDIF
+*   the following lines are for sea-sea chains rejected in GETCSY
+               IF (IREJ1.EQ.-1) NDTUSC = NDTUSC-1
+               ICEVTG(IDCH(NC),1) = ICEVTG(IDCH(NC),1)+1
+            ENDIF
+*
+         ENDIF
+*
+*     update statistics counter
+         ICEVTG(IDCH(NC),0) = ICEVTG(IDCH(NC),0)+1
+*
+         NC = NC+4
+*
+    2 CONTINUE
+*
+*-----------------------------------------------------------------------
+* treatment of low-mass chains (if there are any)
+*
+      IF (NDTUSC.GT.0) THEN
+*
+*   correct chains of very low masses for possible resonances
+         IF (IRESCO.EQ.1) THEN
+            CALL DT_EVTRES(IREJ1)
+            IF (IREJ1.GT.0) THEN
+               IF (IOULEV(1).GT.0) WRITE(LOUT,*) 'rejected 2a in EVENTB'
+               IRRES(1) = IRRES(1)+1
+               GOTO 9999
+            ENDIF
+         ENDIF
+*   fragmentation of low-mass chains
+*!  uncomment this line for internal phojet-fragmentation
+*   (of course it will still be fragmented by DPMJET-routines but it
+*    has to be done here instead of further below)
+C        CALL DT_EVTFRA(IREJ1)
+C        IF (IREJ1.GT.0) THEN
+C           IF (IOULEV(1).GT.0) WRITE(LOUT,*) 'rejected 2b in EVENTB'
+C           IRFRAG = IRFRAG+1
+C           GOTO 9999
+C        ENDIF
+      ELSE
+*! uncomment this line for internal phojet-fragmentation
+C        NPOINT(4) = NHKK+1
+         IF (NPOINT(4).LE.NPOINT(3)) NPOINT(4) = NHKK+1
+      ENDIF
+*
+*-----------------------------------------------------------------------
+* new di-quark breaking mechanisms
+*
+      MXLEFT = 2
+      CALL DT_CHASTA(0)
+      IF ((PDBSEA(1).GT.0.0D0).OR.(PDBSEA(2).GT.0.0D0)
+     &                        .OR.(PDBSEA(3).GT.0.0D0)) THEN
+         CALL DT_DIQBRK
+         MXLEFT = 4
+      ENDIF
+*
+*-----------------------------------------------------------------------
+* hadronize this event
+*
+*   hadronize PHOJET chain systems
+      NPYMAX = 0
+      NPJE   = NPHOSC/MXPHFR
+      IF (MXPHFR.LT.MXLEFT) MXLEFT = 2
+      IF (NPJE.GT.1) THEN
+         NLEFT = NPHOSC-NPJE*MXPHFR
+         DO 20 JFRG=1,NPJE
+            NFRG = JFRG*MXPHFR
+            IF ((JFRG.EQ.NPJE).AND.(NLEFT.LE.MXLEFT)) THEN
+               CALL DT_EVTFRG(1,NPHOSC,NPYMEM,IREJ1)
+               IF (IREJ1.GT.0) GOTO 22
+               NLEFT = 0
+            ELSE
+               CALL DT_EVTFRG(1,NFRG,NPYMEM,IREJ1)
+               IF (IREJ1.GT.0) GOTO 22
+            ENDIF
+            IF (NPYMEM.GT.NPYMAX) NPYMAX = NPYMEM
+   20    CONTINUE
+         IF (NLEFT.GT.0) THEN
+            CALL DT_EVTFRG(1,NPHOSC,NPYMEM,IREJ1)
+            IF (IREJ1.GT.0) GOTO 22
+            IF (NPYMEM.GT.NPYMAX) NPYMAX = NPYMEM
+         ENDIF
+      ELSE
+         CALL DT_EVTFRG(1,NPHOSC,NPYMEM,IREJ1)
+         IF (IREJ1.GT.0) GOTO 22
+         IF (NPYMEM.GT.NPYMAX) NPYMAX = NPYMEM
+      ENDIF
+*
+*   check max. filling level of jetset common and
+*   reduce mxphfr if necessary
+      IF (NPYMAX.GT.3000) THEN
+         IF (NPYMAX.GT.3500) THEN
+            MXPHFR = MAX(1,MXPHFR-2)
+         ELSE
+            MXPHFR = MAX(1,MXPHFR-1)
+         ENDIF
+C        WRITE(LOUT,*) ' EVENTB: Mxphfr reduced to ',MXPHFR
+      ENDIF
+*
+*   hadronize DTUNUC chain systems
+   23 CONTINUE
+      IBACK = MXDTFR
+      CALL DT_EVTFRG(2,IBACK,NPYMEM,IREJ2)
+      IF (IREJ2.GT.0) GOTO 22
+*
+*   check max. filling level of jetset common and
+*   reduce mxdtfr if necessary
+      IF (NPYMEM.GT.3000) THEN
+         IF (NPYMEM.GT.3500) THEN
+            MXDTFR = MAX(1,MXDTFR-20)
+         ELSE
+            MXDTFR = MAX(1,MXDTFR-10)
+         ENDIF
+C        WRITE(LOUT,*) ' EVENTB: Mxdtfr reduced to ',MXDTFR
+      ENDIF
+*
+      IF (IBACK.EQ.-1) GOTO 23
+*
+   22 CONTINUE
+C     CALL DT_EVTFRG(1,IREJ1)
+C     CALL DT_EVTFRG(2,IREJ2)
+      IF ((IREJ1.GT.0).OR.(IREJ2.GT.0)) THEN
+         IF (IOULEV(1).GT.0) WRITE(LOUT,*) 'rejected 1 in EVENTB'
+         IRFRAG = IRFRAG+1
+         GOTO 9999
+      ENDIF
+*
+* get final state particles from /DTEVTP/
+*! uncomment this line for internal phojet-fragmentation
+C     CALL DT_GETFSP(IDUM,IDUM,PP,PT,2)
+
+      IF (IJPROJ.NE.7)
+     &   CALL DT_EMC2(9,10,0,0,0,3,1,0,0,0,0,3,4,88,IREJ3)
+C     IF (IREJ3.NE.0) GOTO 9999
+
+      RETURN
+
+ 9999 CONTINUE
+      IREVT = IREVT+1
+      IREJ  = 1
+      RETURN
+      END
+*
+*===getpje=============================================================*
+*
+CDECK  ID>, DT_GETPJE
+      SUBROUTINE DT_GETPJE(MO1,MO2,PP,PT,MODE,IPJE,IREJ)
+
+************************************************************************
+* This subroutine copies PHOJET partons and strings from POEVT1 into   *
+* DTEVT1.                                                              *
+*      MO1,MO2   indices of first and last mother-parton in DTEVT1     *
+*      PP,PT     4-momenta of projectile/target being handled by       *
+*                PHOJET                                                *
+* This version dated 11.12.99 is written by S. Roesler                 *
+************************************************************************
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      PARAMETER ( LINP = 5 ,
+     &            LOUT = 6 ,
+     &            LDAT = 9 )
+
+      PARAMETER (TINY10=1.0D-10,TINY1=1.0D-1,
+     &           ZERO=0.0D0,ONE=1.0D0,OHALF=0.5D0)
+
+      LOGICAL LFLIP
+
+* event history
+
+      PARAMETER (NMXHKK=200000)
+
+      COMMON /DTEVT1/ NHKK,NEVHKK,ISTHKK(NMXHKK),IDHKK(NMXHKK),
+     &                JMOHKK(2,NMXHKK),JDAHKK(2,NMXHKK),
+     &                PHKK(5,NMXHKK),VHKK(4,NMXHKK),WHKK(4,NMXHKK)
+* extended event history
+      COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
+     &                IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
+     &                IHIST(2,NMXHKK)
+* Lorentz-parameters of the current interaction
+      COMMON /DTLTRA/ GACMS(2),BGCMS(2),GALAB,BGLAB,BLAB,
+     &                UMO,PPCM,EPROJ,PPROJ
+* DTUNUC-PHOJET interface, Lorentz-param. of n-n subsystem
+      COMMON /DTLTSU/ BGX,BGY,BGZ,GAM
+* flags for input different options
+      LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
+      COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
+     &                LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IFRAME,ITRSPT
+* statistics: double-Pomeron exchange
+      COMMON /DTFLG2/ INTFLG,IPOPO
+* statistics
+      COMMON /DTSTA1/ ICREQU,ICSAMP,ICCPRO,ICDPR,ICDTA,
+     &                ICRJSS,ICVV2S,ICCHAI(2,9),ICRES(9),ICDIFF(5),
+     &                ICEVTG(8,0:30)
+* rejection counter
+      COMMON /DTREJC/ IRPT,IRHHA,IRRES(2),LOMRES,LOBRES,
+     &                IRCHKI(2),IRFRAG,IRCRON(3),IREVT,
+     &                IREXCI(3),IRDIFF(2),IRINC
+
+C  standard particle data interface
+      INTEGER NMXHEP
+
+      PARAMETER (NMXHEP=4000)
+
+      INTEGER NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
+      DOUBLE PRECISION PHEP,VHEP
+      COMMON /POEVT1/ NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
+     &                JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),
+     &                VHEP(4,NMXHEP)
+C  extension to standard particle data interface (PHOJET specific)
+      INTEGER IMPART,IPHIST,ICOLOR
+      COMMON /POEVT2/ IMPART(NMXHEP),IPHIST(2,NMXHEP),ICOLOR(2,NMXHEP)
+
+C  color string configurations including collapsed strings and hadrons
+      INTEGER MSTR
+      PARAMETER (MSTR=500)
+      INTEGER NPOS,NCODE,IPAR1,IPAR2,IPAR3,IPAR4,NNCH,IBHAD
+      COMMON /POSTRG/ NPOS(4,MSTR),NCODE(MSTR),
+     &                IPAR1(MSTR),IPAR2(MSTR),IPAR3(MSTR),IPAR4(MSTR),
+     &                NNCH(MSTR),IBHAD(MSTR),ISTR
+C  general process information
+      INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
+      COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
+C  model switches and parameters
+      CHARACTER*8 MDLNA
+      INTEGER ISWMDL,IPAMDL
+      DOUBLE PRECISION PARMDL
+      COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
+C  event debugging information
+      INTEGER NMAXD
+      PARAMETER (NMAXD=100)
+      INTEGER IDEB,KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,KHTRG,KSLOO,
+     &        KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
+      COMMON /PODEBG/ IDEB(NMAXD),KSPOM,KHPOM,KSREG,KHDIR,KACCEP,KSTRG,
+     &                KHTRG,KSLOO,KHLOO,KSDPO,KHDPO,KEVENT,KSOFT,KHARD
+
+      DIMENSION PP(4),PT(4)
+      DATA MAXLOP /10000/
+
+      INHKK = NHKK
+      LFLIP = .TRUE.
+    1 CONTINUE
+      NPVAL = 0
+      NTVAL = 0
+      IREJ  = 0
+
+*   store initial momenta for energy-momentum conservation check
+      IF (LEMCCK) THEN
+         CALL DT_EVTEMC(PP(1),PP(2),PP(3),PP(4),1,IDUM1,IDUM2)
+         CALL DT_EVTEMC(PT(1),PT(2),PT(3),PT(4),2,IDUM1,IDUM2)
+      ENDIF
+* copy partons and strings from POEVT1 into DTEVT1
+      DO 11 I=1,ISTR
+C        IF ((NCODE(I).EQ.-99).AND.(IPAMDL(17).EQ.0)) THEN
+         IF (NCODE(I).EQ.-99) THEN
+            IDXSTG = NPOS(1,I)
+            IDSTG  = IDHEP(IDXSTG)
+            PX = PHEP(1,IDXSTG)
+            PY = PHEP(2,IDXSTG)
+            PZ = PHEP(3,IDXSTG)
+            PE = PHEP(4,IDXSTG)
+            IF (MODE.LT.0) THEN
+               ISTAT = 70000+IPJE
+               CALL DT_EVTPUT(2,ISTAT,MO1,MO2,PX,PY,PZ,PE,
+     &                        11,IDSTG,0)
+               IF (LEMCCK) THEN
+                  PX = -PX
+                  PY = -PY
+                  PZ = -PZ
+                  PE = -PE
+                  CALL DT_EVTEMC(PX,PY,PZ,PE,2,IDUM1,IDUM2)
+               ENDIF
+            ELSE
+               CALL DT_DALTRA(GAM,BGX,BGY,BGZ,PX,PY,PZ,PE,PTOTMP,
+     &                        PPX,PPY,PPZ,PPE)
+               ISTAT = 70000+IPJE
+               CALL DT_EVTPUT(2,ISTAT,MO1,MO2,PPX,PPY,PPZ,PPE,
+     &                        11,IDSTG,0)
+               IF (LEMCCK) THEN
+                  PX = -PPX
+                  PY = -PPY
+                  PZ = -PPZ
+                  PE = -PPE
+                  CALL DT_EVTEMC(PX,PY,PZ,PE,2,IDUM1,IDUM2)
+               ENDIF
+            ENDIF
+            NOBAM(NHKK)   = 0
+            IHIST(1,NHKK) = IPHIST(1,IDXSTG)
+            IHIST(2,NHKK) = 0
+         ELSEIF (NCODE(I).GE.0) THEN
+*   indices of partons and string in POEVT1
+            IDX1 = ABS(JMOHEP(1,NPOS(1,I)))
+            IDX2 = ABS(JMOHEP(2,NPOS(1,I)))
+            IF ((IDX1.GT.IDX2).OR.(JMOHEP(2,NPOS(1,I)).GT.0)) THEN
+               WRITE(LOUT,*) ' GETPJE: IDX1.GT.IDX2 ',IDX1,IDX2,
+     &         ' or JMOHEP(2,NPOS(1,I)).GT.0 ',JMOHEP(2,NPOS(1,I)),' ! '
+               STOP ' GETPJE 1'
+            ENDIF
+            IDXSTG = NPOS(1,I)
+*   find "mother" string of the string
+            IDXMS1 = ABS(JMOHEP(1,IDX1))
+            IDXMS2 = ABS(JMOHEP(1,IDX2))
+            IF (IDXMS1.NE.IDXMS2) THEN
+               IDXMS1 = IDXSTG
+               IDXMS2 = IDXSTG
+C              STOP ' GETPJE: IDXMS1.NE.IDXMS2 !'
+            ENDIF
+*   search POEVT1 for the original hadron of the parton
+            ILOOP = 0
+            IPOM1 = 0
+   14       CONTINUE
+            ILOOP = ILOOP+1
+
+            IF (IDHEP(IDXMS1).EQ.990) IPOM1 = 1
+
+            IDXMS1 = ABS(JMOHEP(1,IDXMS1))
+            IF ((IDXMS1.NE.1).AND.(IDXMS1.NE.2).AND.
+     &          (ILOOP.LT.MAXLOP)) GOTO 14
+            IF (ILOOP.EQ.MAXLOP) WRITE(LOUT,*) ' GETPJE: MAXLOP in 1 ! '
+            IPOM2 = 0
+            ILOOP = 0
+   15       CONTINUE
+            ILOOP = ILOOP+1
+
+            IF (IDHEP(IDXMS2).EQ.990) IPOM2 = 1
+
+            IF ((ILOOP.EQ.1).OR.(IDHEP(IDXMS2).GE.7777)) THEN
+               IDXMS2 = ABS(JMOHEP(2,IDXMS2))
+            ELSE
+               IDXMS2 = ABS(JMOHEP(1,IDXMS2))
+            ENDIF
+            IF ((IDXMS2.NE.1).AND.(IDXMS2.NE.2).AND.
+     &          (ILOOP.LT.MAXLOP)) GOTO 15
+            IF (ILOOP.EQ.MAXLOP) WRITE(LOUT,*) ' GETPJE: MAXLOP in 5 ! '
+*   parton 1
+            IF (IDXMS1.EQ.1) THEN
+               ISPTN1 = ISTHKK(MO1)
+               M1PTN1 = MO1
+               M2PTN1 = MO1+2
+            ELSE
+               ISPTN1 = ISTHKK(MO2)
+               M1PTN1 = MO2-2
+               M2PTN1 = MO2
+            ENDIF
+*   parton 2
+            IF (IDXMS2.EQ.1) THEN
+               ISPTN2 = ISTHKK(MO1)
+               M1PTN2 = MO1
+               M2PTN2 = MO1+2
+            ELSE
+               ISPTN2 = ISTHKK(MO2)
+               M1PTN2 = MO2-2
+               M2PTN2 = MO2
+            ENDIF
+*   check for mis-identified mothers and switch mother indices if necessary
+            IF ((IDXMS1.EQ.IDXMS2).AND.(IPROCE.NE.5).AND.(IPROCE.NE.6)
+     &          .AND.((IDHEP(IDX1).NE.21).OR.(IDHEP(IDX2).NE.21)).AND.
+     &          (LFLIP)) THEN
+               IF (PHEP(3,IDX1).GT.PHEP(3,IDX2)) THEN
+                  ISPTN1 = ISTHKK(MO1)
+                  M1PTN1 = MO1
+                  M2PTN1 = MO1+2
+                  ISPTN2 = ISTHKK(MO2)
+                  M1PTN2 = MO2-2
+                  M2PTN2 = MO2
+               ELSE
+                  ISPTN1 = ISTHKK(MO2)
+                  M1PTN1 = MO2-2
+                  M2PTN1 = MO2
+                  ISPTN2 = ISTHKK(MO1)
+                  M1PTN2 = MO1
+                  M2PTN2 = MO1+2
+               ENDIF
+            ENDIF
+*   register partons in temporary common
+*     parton at chain end
+            PX = PHEP(1,IDX1)
+            PY = PHEP(2,IDX1)
+            PZ = PHEP(3,IDX1)
+            PE = PHEP(4,IDX1)
+* flag only partons coming from Pomeron with 41/42
+C           IF ((IPOM1.NE.0).OR.(NPOS(4,I).GE.4)) THEN
+            IF (IPOM1.NE.0) THEN
+               ISTX = ABS(ISPTN1)/10
+               IMO  = ABS(ISPTN1)-10*ISTX
+               ISPTN1 = -(40+IMO)
+            ELSE
+               IF ((ICOLOR(2,IDX1).EQ.0).OR.(IDHEP(IDX1).EQ.21)) THEN
+                  ISTX = ABS(ISPTN1)/10
+                  IMO  = ABS(ISPTN1)-10*ISTX
+                  IF ((IDHEP(IDX1).EQ.21).OR.
+     &                (ABS(IPHIST(1,IDX1)).GE.100)) THEN
+                     ISPTN1 = -(60+IMO)
+                  ELSE
+                     ISPTN1 = -(50+IMO)
+                  ENDIF
+               ENDIF
+            ENDIF
+            IF (ISPTN1.EQ.-21) NPVAL = NPVAL+1
+            IF (ISPTN1.EQ.-22) NTVAL = NTVAL+1
+            IF (MODE.LT.0) THEN
+               CALL DT_EVTPUT(ISPTN1,IDHEP(IDX1),M1PTN1,M2PTN1,PX,PY,
+     &                        PZ,PE,0,0,0)
+            ELSE
+               CALL DT_DALTRA(GAM,BGX,BGY,BGZ,PX,PY,PZ,PE,PTOTMP,
+     &                        PPX,PPY,PPZ,PPE)
+               CALL DT_EVTPUT(ISPTN1,IDHEP(IDX1),M1PTN1,M2PTN1,PPX,PPY,
+     &                        PPZ,PPE,0,0,0)
+            ENDIF
+            IHIST(1,NHKK) = IPHIST(1,IDX1)
+            IHIST(2,NHKK) = 0
+            DO 19 KK=1,4
+               VHKK(KK,NHKK) = VHKK(KK,M2PTN1)
+               WHKK(KK,NHKK) = WHKK(KK,M1PTN1)
+   19       CONTINUE
+            VHKK(4,NHKK) = VHKK(3,M2PTN1)/BLAB-VHKK(3,M1PTN1)/BGLAB
+            WHKK(4,NHKK) = -WHKK(3,M1PTN1)/BLAB+WHKK(3,M2PTN1)/BGLAB
+            M1STRG = NHKK
+*     gluon kinks
+            NGLUON = IDX2-IDX1-1
+            IF (NGLUON.GT.0) THEN
+               DO 17 IGLUON=1,NGLUON
+                  IDX   = IDX1+IGLUON
+                  IDXMS = ABS(JMOHEP(1,IDX))
+                  IF ((IDXMS.NE.1).AND.(IDXMS.NE.2)) THEN
+                     ILOOP = 0
+   16                CONTINUE
+                     ILOOP = ILOOP+1
+                     IDXMS = ABS(JMOHEP(1,IDXMS))
+                     IF ((IDXMS.NE.1).AND.(IDXMS.NE.2).AND.
+     &                   (ILOOP.LT.MAXLOP)) GOTO 16
+                     IF (ILOOP.EQ.MAXLOP)
+     &                  WRITE(LOUT,*) ' GETPJE: MAXLOP in 3 ! '
+                  ENDIF
+                  IF (IDXMS.EQ.1) THEN
+                     ISPTN = ISTHKK(MO1)
+                     M1PTN = MO1
+                     M2PTN = MO1+2
+                  ELSE
+                     ISPTN = ISTHKK(MO2)
+                     M1PTN = MO2-2
+                     M2PTN = MO2
+                  ENDIF
+                  PX = PHEP(1,IDX)
+                  PY = PHEP(2,IDX)
+                  PZ = PHEP(3,IDX)
+                  PE = PHEP(4,IDX)
+                  IF ((ICOLOR(2,IDX).EQ.0).OR.(IDHEP(IDX).EQ.21)) THEN
+                     ISTX = ABS(ISPTN)/10
+                     IMO  = ABS(ISPTN)-10*ISTX
+                     IF ((IDHEP(IDX).EQ.21).OR.
+     &                   (ABS(IPHIST(1,IDX)).GE.100)) THEN
+                        ISPTN = -(60+IMO)
+                     ELSE
+                        ISPTN = -(50+IMO)
+                     ENDIF
+                  ENDIF
+                  IF (ISPTN.EQ.-21) NPVAL = NPVAL+1
+                  IF (ISPTN.EQ.-22) NTVAL = NTVAL+1
+                  IF (MODE.LT.0) THEN
+                     CALL DT_EVTPUT(ISPTN,IDHEP(IDX),M1PTN,M2PTN,
+     &                              PX,PY,PZ,PE,0,0,0)
+                  ELSE
+                     CALL DT_DALTRA(GAM,BGX,BGY,BGZ,PX,PY,PZ,PE,PTOTMP,
+     &                              PPX,PPY,PPZ,PPE)
+                     CALL DT_EVTPUT(ISPTN,IDHEP(IDX),M1PTN,M2PTN,
+     &                              PPX,PPY,PPZ,PPE,0,0,0)
+                  ENDIF
+                  IHIST(1,NHKK) = IPHIST(1,IDX)
+                  IHIST(2,NHKK) = 0
+                  DO 20 KK=1,4
+                     VHKK(KK,NHKK) = VHKK(KK,M2PTN)
+                     WHKK(KK,NHKK) = WHKK(KK,M1PTN)
+   20             CONTINUE
+                  VHKK(4,NHKK)= VHKK(3,M2PTN)/BLAB-VHKK(3,M1PTN)/BGLAB
+                  WHKK(4,NHKK)= -WHKK(3,M1PTN)/BLAB+WHKK(3,M2PTN)/BGLAB
+   17          CONTINUE
+            ENDIF
+*     parton at chain end
+            PX = PHEP(1,IDX2)
+            PY = PHEP(2,IDX2)
+            PZ = PHEP(3,IDX2)
+            PE = PHEP(4,IDX2)
+* flag only partons coming from Pomeron with 41/42
+C           IF ((IPOM2.NE.0).OR.(NPOS(4,I).GE.4)) THEN
+            IF (IPOM2.NE.0) THEN
+               ISTX = ABS(ISPTN2)/10
+               IMO  = ABS(ISPTN2)-10*ISTX
+               ISPTN2 = -(40+IMO)
+            ELSE
+               IF ((ICOLOR(2,IDX2).EQ.0).OR.(IDHEP(IDX2).EQ.21)) THEN
+                  ISTX = ABS(ISPTN2)/10
+                  IMO  = ABS(ISPTN2)-10*ISTX
+                  IF ((IDHEP(IDX2).EQ.21).OR.
+     &                (ABS(IPHIST(1,IDX2)).GE.100)) THEN
+                     ISPTN2 = -(60+IMO)
+                  ELSE
+                     ISPTN2 = -(50+IMO)
+                  ENDIF
+               ENDIF
+            ENDIF
+            IF (ISPTN2.EQ.-21) NPVAL = NPVAL+1
+            IF (ISPTN2.EQ.-22) NTVAL = NTVAL+1
+            IF (MODE.LT.0) THEN
+               CALL DT_EVTPUT(ISPTN2,IDHEP(IDX2),M1PTN2,M2PTN2,
+     &                        PX,PY,PZ,PE,0,0,0)
+            ELSE
+               CALL DT_DALTRA(GAM,BGX,BGY,BGZ,PX,PY,PZ,PE,PTOTMP,
+     &                        PPX,PPY,PPZ,PPE)
+               CALL DT_EVTPUT(ISPTN2,IDHEP(IDX2),M1PTN2,M2PTN2,
+     &                        PPX,PPY,PPZ,PPE,0,0,0)
+            ENDIF
+            IHIST(1,NHKK) = IPHIST(1,IDX2)
+            IHIST(2,NHKK) = 0
+            DO 21 KK=1,4
+               VHKK(KK,NHKK) = VHKK(KK,M2PTN2)
+               WHKK(KK,NHKK) = WHKK(KK,M1PTN2)
+   21       CONTINUE
+            VHKK(4,NHKK) = VHKK(3,M2PTN2)/BLAB-VHKK(3,M1PTN2)/BGLAB
+            WHKK(4,NHKK) = -WHKK(3,M1PTN2)/BLAB+WHKK(3,M2PTN2)/BGLAB
+            M2STRG = NHKK
+*   register string
+            JSTRG = 100*IPROCE+NCODE(I)
+            PX = PHEP(1,IDXSTG)
+            PY = PHEP(2,IDXSTG)
+            PZ = PHEP(3,IDXSTG)
+            PE = PHEP(4,IDXSTG)
+            IF (MODE.LT.0) THEN
+               ISTAT = 70000+IPJE
+               CALL DT_EVTPUT(JSTRG,ISTAT,M1STRG,M2STRG,
+     &                        PX,PY,PZ,PE,0,0,0)
+               IF (LEMCCK) THEN
+                  PX = -PX
+                  PY = -PY
+                  PZ = -PZ
+                  PE = -PE
+                  CALL DT_EVTEMC(PX,PY,PZ,PE,2,IDUM1,IDUM2)
+               ENDIF
+            ELSE
+               CALL DT_DALTRA(GAM,BGX,BGY,BGZ,PX,PY,PZ,PE,PTOTMP,
+     &                        PPX,PPY,PPZ,PPE)
+               ISTAT = 70000+IPJE
+               CALL DT_EVTPUT(JSTRG,ISTAT,M1STRG,M2STRG,
+     &                        PPX,PPY,PPZ,PPE,0,0,0)
+               IF (LEMCCK) THEN
+                  PX = -PPX
+                  PY = -PPY
+                  PZ = -PPZ
+                  PE = -PPE
+                  CALL DT_EVTEMC(PX,PY,PZ,PE,2,IDUM1,IDUM2)
+               ENDIF
+            ENDIF
+            NOBAM(NHKK)   = 0
+            IHIST(1,NHKK) = 0
+            IHIST(2,NHKK) = 0
+            DO 18 KK=1,4
+               VHKK(KK,NHKK) = VHKK(KK,MO2)
+               WHKK(KK,NHKK) = WHKK(KK,MO1)
+   18       CONTINUE
+            VHKK(4,NHKK) = VHKK(3,MO2)/BLAB-VHKK(3,MO1)/BGLAB
+            WHKK(4,NHKK) = -WHKK(3,MO1)/BLAB+WHKK(3,MO2)/BGLAB
+         ENDIF
+   11 CONTINUE
+
+      IF ( ((NPVAL.GT.2).OR.(NTVAL.GT.2)).AND.(LFLIP) ) THEN
+         NHKK  = INHKK
+         LFLIP = .FALSE.
+         GOTO 1
+      ENDIF
+
+      IF (LEMCCK) THEN
+         IF (UMO.GT.1.0D5) THEN
+            CHKLEV = 1.0D0
+         ELSE
+            CHKLEV = TINY1
+         ENDIF
+         CALL DT_EVTEMC(DUM1,DUM2,DUM3,CHKLEV,-1,1000,IREJ2)
+
+         IF (IREJ2.GT.ZERO) CALL PHO_PREVNT(0)
+
+      ENDIF
+
+* internal statistics
+*   dble-Po statistics.
+      IF (IPROCE.NE.4) IPOPO = 0
+
+      INTFLG = IPROCE
+      IDCHSY = IDCH(MO1)
+      IF ((IPROCE.GE.1).AND.(IPROCE.LE.8)) THEN
+         ICEVTG(IDCHSY,IPROCE+2) = ICEVTG(IDCHSY,IPROCE+2)+1
+      ELSE
+         WRITE(LOUT,1000) IPROCE,NEVHKK,MO1
+ 1000    FORMAT(1X,'GETFSP:   warning! incons. process id. (',I2,
+     &          ') at evt(chain) ',I6,'(',I2,')')
+      ENDIF
+      IF (IPROCE.EQ.5) THEN
+         IF ((IDIFR1.GE.1).AND.(IDIFR1.LE.3)) THEN
+            ICEVTG(IDCHSY,18+IDIFR1) = ICEVTG(IDCHSY,18+IDIFR1)+1
+         ELSE
+C           WRITE(LOUT,1001) IPROCE,IDIFR1,IDIFR2
+ 1001       FORMAT(1X,'GETFSP:   warning! incons. diffrac. id. ',
+     &             '(IPROCE,IDIFR1,IDIFR2=',3I3,')')
+         ENDIF
+      ELSEIF (IPROCE.EQ.6) THEN
+         IF ((IDIFR2.GE.1).AND.(IDIFR2.LE.3)) THEN
+            ICEVTG(IDCHSY,21+IDIFR2) = ICEVTG(IDCHSY,21+IDIFR2)+1
+         ELSE
+C           WRITE(LOUT,1001) IPROCE,IDIFR1,IDIFR2
+         ENDIF
+      ELSEIF (IPROCE.EQ.7) THEN
+         IF ((IDIFR1.GE.1).AND.(IDIFR1.LE.3).AND.
+     &       (IDIFR2.GE.1).AND.(IDIFR2.LE.3)) THEN
+            IF ((IDIFR1.EQ.1).AND.(IDIFR2.EQ.1))
+     &         ICEVTG(IDCHSY,25) = ICEVTG(IDCHSY,25)+1
+            IF ((IDIFR1.EQ.2).AND.(IDIFR2.EQ.2))
+     &         ICEVTG(IDCHSY,26) = ICEVTG(IDCHSY,26)+1
+            IF ((IDIFR1.EQ.1).AND.(IDIFR2.EQ.2))
+     &         ICEVTG(IDCHSY,27) = ICEVTG(IDCHSY,27)+1
+            IF ((IDIFR1.EQ.2).AND.(IDIFR2.EQ.1))
+     &         ICEVTG(IDCHSY,28) = ICEVTG(IDCHSY,28)+1
+         ELSE
+            WRITE(LOUT,1001) IPROCE,IDIFR1,IDIFR2
+         ENDIF
+      ENDIF
+      IF ((IDIFR1+IDIFR2.EQ.0).AND.(KHDIR.GE.1).AND.(KHDIR.LE.3))
+     &                                                       THEN
+         ICEVTG(IDCHSY,10+KHDIR) = ICEVTG(IDCHSY,10+KHDIR)+1
+         ICEVTG(IDCHSY,10+KHDIR) = ICEVTG(IDCHSY,10+KHDIR)+1
+         ICEVTG(IDCHSY,10+KHDIR) = ICEVTG(IDCHSY,10+KHDIR)+1
+      ENDIF
+      ICEVTG(IDCHSY,14) = ICEVTG(IDCHSY,14)+KSPOM
+      ICEVTG(IDCHSY,15) = ICEVTG(IDCHSY,15)+KHPOM
+      ICEVTG(IDCHSY,16) = ICEVTG(IDCHSY,16)+KSREG
+      ICEVTG(IDCHSY,17) = ICEVTG(IDCHSY,17)+(KSTRG+KHTRG)
+      ICEVTG(IDCHSY,18) = ICEVTG(IDCHSY,18)+(KSLOO+KHLOO)
+
+      RETURN
+
+ 9999 CONTINUE
+      IREJ = 1
+      RETURN
+      END
+*
+*===phoini=============================================================*
+*
+CDECK  ID>, DT_PHOINI
+      SUBROUTINE DT_PHOINI
+
+************************************************************************
+* Initialization PHOJET-event generator for nucleon-nucleon interact.  *
+* This version dated 16.11.95 is written by S. Roesler                 *
+* Last change: s.r. 21.01.01                                           *
+************************************************************************
+
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      SAVE
+
+      PARAMETER ( LINP = 5 ,
+     &            LOUT = 6 ,
+     &            LDAT = 9 )
+
+      PARAMETER (TINY10=1.0D-10,ZERO=0.0D0,ONE=1.0D0)
+
+* nucleon-nucleon event-generator
+      CHARACTER*8 CMODEL
+      LOGICAL LPHOIN
+      COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
+* particle properties (BAMJET index convention)
+      CHARACTER*8  ANAME
+      COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
+     &                IICH(210),IIBAR(210),K1(210),K2(210)
+* Lorentz-parameters of the current interaction
+      COMMON /DTLTRA/ GACMS(2),BGCMS(2),GALAB,BGLAB,BLAB,
+     &                UMO,PPCM,EPROJ,PPROJ
+* properties of interacting particles
+      COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
+* properties of photon/lepton projectiles
+      COMMON /DTGPRO/ VIRT,PGAMM(4),PLEPT0(4),PLEPT1(4),PNUCL(4),IDIREC
+
+      PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
+
+* emulsion treatment
+      COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
+     &                NCOMPO,IEMUL
+* VDM parameter for photon-nucleus interactions
+      COMMON /DTVDMP/ RL2,EPSPOL,INTRGE(2),IDPDF,MODEGA,ISHAD(3)
+* nuclear potential
+      LOGICAL LFERMI
+      COMMON /DTNPOT/ PFERMP(2),PFERMN(2),FERMOD,
+     &                EBINDP(2),EBINDN(2),EPOT(2,210),
+     &                ETACOU(2),ICOUL,LFERMI
+* Glauber formalism: flags and parameters for statistics
+      LOGICAL LPROD
+      CHARACTER*8 CGLB
+      COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
+*
+* parameters for cascade calculations:
+* maximum mumber of PDF's which can be defined in phojet (limited
+* by the dimension of ipdfs in pho_setpdf)
+      PARAMETER (MAXPDF = 20)
+* PDF parametrization and number of set for the first 30 hadrons in
+* the bamjet-code list
+*   negative numbers mean that the PDF is set in phojet,
+*   zero stands for "not a hadron"
+      DIMENSION IPARPD(30),ISETPD(30)
+* PDF parametrization
+      DATA IPARPD /
+     &  -5,-5, 0, 0, 0, 0,-5,-5,-5, 0, 0, 5,-5,-5, 5, 5, 5, 5, 5, 5,
+     &   5, 5,-5, 5, 5, 0, 0, 0, 0, 0/
+* number of set
+      DATA ISETPD /
+     &  -6,-6, 0, 0, 0, 0,-3,-6,-6, 0, 0, 2,-2,-2, 2, 2, 6, 6, 2, 6,
+     &   6, 6,-2, 2, 2, 0, 0, 0, 0, 0/
+
+**PHOJET105a
+C     COMMON /GLOCMS/ XECM,XPCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
+C     PARAMETER ( MAXPRO = 16 )
+C     PARAMETER ( MAXTAB = 20 )
+C     COMMON /HAXSEC/ XSECTA(4,-1:MAXPRO,4,MAXTAB),XSECT(6,-1:MAXPRO),
+C    &                MXSECT(0:4,-1:MAXPRO,4),ECMSH(4,MAXTAB),ISTTAB
+C     CHARACTER*8 MDLNA
+C     COMMON /MODELS/ MDLNA(50),ISWMDL(50),PARMDL(200),IPAMDL(100)
+C     COMMON /PROCES/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15)
+**PHOJET110
+C  global event kinematics and particle IDs
+      INTEGER IFPAP,IFPAB
+      DOUBLE PRECISION ECM,PCM,PMASS,PVIRT
+      COMMON /POGCMS/ ECM,PCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
+C  hard cross sections and MC selection weights
+      INTEGER Max_pro_2
+      PARAMETER ( Max_pro_2 = 16 )
+      INTEGER IHa_last,IHb_last,MH_pro_on,MH_tried,
+     &  MH_acc_1,MH_acc_2
+      DOUBLE PRECISION Hfac,HWgx,HSig,Hdpt,HEcm_last,HQ2a_last,HQ2b_last
+      COMMON /POHRCS/ Hfac(-1:Max_pro_2),HWgx(-1:Max_pro_2),
+     &  HSig(-1:Max_pro_2),Hdpt(-1:Max_pro_2),
+     &  HEcm_last,HQ2a_last,HQ2b_last,IHa_last,IHb_last,
+     &  MH_pro_on(-1:Max_pro_2,0:4),MH_tried(-1:Max_pro_2,0:4),
+     &  MH_acc_1(-1:Max_pro_2,0:4),MH_acc_2(-1:Max_pro_2,0:4)
+C  model switches and parameters
+      CHARACTER*8 MDLNA
+      INTEGER ISWMDL,IPAMDL
+      DOUBLE PRECISION PARMDL
+      COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
+C  general process information
+      INTEGER IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON
+      COMMON /POPRCS/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15,4)
+**
+      DIMENSION PP(4),PT(4)
+
+      LOGICAL LSTART
+      DATA LSTART /.TRUE./
+
+      IJP = IJPROJ
+      IJT = IJTARG
+      Q2  = VIRT
+* lepton-projectiles: initialize real photon instead
+      IF ((IJP.EQ.3).OR.(IJP.EQ.4).OR.(IJP.EQ.10).OR.(IJP.EQ.11)) THEN
+         IJP = 7
+         Q2  = ZERO
+      ENDIF
+
+      IF (LPHOIN) CALL PHO_INIT(-1,IDUM)
+
+* switch Reggeon off
+C     IPAMDL(3)= 0
+      IF (IP.EQ.1) THEN
+         IFPAP(1) = IDT_IPDGHA(IJP)
+         IFPAB(1) = IJP
+      ELSE
+         IFPAP(1) = 2212
+         IFPAB(1) = IDT_ICIHAD(IFPAP(1))
+      ENDIF
+      PMASS(1) = AAM(IFPAB(1))-SQRT(Q2)
+      PVIRT(1) = PMASS(1)**2
+      IF (IT.EQ.1) THEN
+         IFPAP(2) = IDT_IPDGHA(IJT)
+         IFPAB(2) = IJT
+      ELSE
+         IFPAP(2) = 2212
+         IFPAB(2) = IDT_ICIHAD(IFPAP(2))
+      ENDIF
+      PMASS(2) = AAM(IFPAB(2))
+      PVIRT(2) = ZERO
+      DO 1 K=1,4
+         PP(K) = ZERO
+         PT(K) = ZERO
+    1 CONTINUE
+* get max. possible momenta of incoming particles to be used for PHOJET ini.
+      PPF = ZERO
+      PTF = ZERO
+      SCPF= 1.5D0
+      IF (UMO.GE.1.E5) THEN
+         SCPF= 5.0D0
+      ENDIF
+      IF (NCOMPO.GT.0) THEN
+         DO 2 I=1,NCOMPO
+            IF (IT.GT.1) THEN
+               CALL DT_NCLPOT(IEMUCH(I),IEMUMA(I),ITZ,IT,ZERO,ZERO,0)
+            ELSE
+               CALL DT_NCLPOT(IPZ,IP,IEMUCH(I),IEMUMA(I),ZERO,ZERO,0)
+            ENDIF
+            PPFTMP = MAX(PFERMP(1),PFERMN(1))
+            PTFTMP = MAX(PFERMP(2),PFERMN(2))
+            IF (PPFTMP.GT.PPF) PPF = PPFTMP
+            IF (PTFTMP.GT.PTF) PTF = PTFTMP
+    2    CONTINUE
+      ELSE
+         CALL DT_NCLPOT(IPZ,IP,ITZ,IT,ZERO,ZERO,0)
+         PPF = MAX(PFERMP(1),PFERMN(1))
+         PTF = MAX(PFERMP(2),PFERMN(2))
+      ENDIF
+      PTF = -PTF
+      PPF = SCPF*PPF
+      PTF = SCPF*PTF
+      IF (IJP.EQ.7) THEN
+         AMP2  = SIGN(PMASS(1)**2,PMASS(1))
+         PP(3) = PPCM
+         PP(4) = SQRT(AMP2+PP(3)**2)
+      ELSE
+         EPF = SQRT(PPF**2+PMASS(1)**2)
+         CALL DT_LTNUC(PPF,EPF,PP(3),PP(4),2)
+      ENDIF
+      ETF = SQRT(PTF**2+PMASS(2)**2)
+      CALL DT_LTNUC(PTF,ETF,PT(3),PT(4),3)
+      ECMINI = SQRT((PP(4)+PT(4))**2-(PP(1)+PT(1))**2-
+     &              (PP(2)+PT(2))**2-(PP(3)+PT(3))**2)
+      IF (LSTART) THEN
+C *** Commented by Chiara
+C         WRITE(LOUT,1001) IP,IPZ,SCPF,PPF,PP
+ 1001    FORMAT(
+     &      ' DT_PHOINI:    PHOJET initialized for projectile A,Z = ',
+     &      I3,',',I2,/,F4.1,'xp_F(max) = ',E10.3,'  p(max) = ',4E10.3)
+C *** Commented by Chiara
+C         IF (NCOMPO.GT.0) THEN
+C            WRITE(LOUT,1002) SCPF,PTF,PT
+C         ELSE
+C            WRITE(LOUT,1003) IT,ITZ,SCPF,PTF,PT
+C         ENDIF
+ 1002    FORMAT(
+     &      ' DT_PHOINI:    PHOJET initialized for target emulsion  ',
+     &          /,F4.1,'xp_F(max) = ',E10.3,'  p(max) = ',4E10.3)
+ 1003    FORMAT(
+     &      ' DT_PHOINI:    PHOJET initialized for target     A,Z = ',
+     &      I3,',',I2,/,F4.1,'xp_F(max) = ',E10.3,'  p(max) = ',4E10.3)
+C *** Commented by Chiara
+C         WRITE(LOUT,1004) ECMINI
+ 1004    FORMAT(' E_cm = ',E10.3)
+         IF (IJP.EQ.8) WRITE(LOUT,1005)
+ 1005    FORMAT(
+     &      ' DT_PHOINI: warning! proton parameters used for neutron',
+     &          ' projectile')
+         LSTART = .FALSE.
+      ENDIF
+* switch off new diffractive cross sections at low energies for nuclei
+* (temporary solution)
+      IF ((ISWMDL(30).NE.0).AND.((IP.GT.1).OR.(IT.GT.1))) THEN
+         WRITE(LOUT,'(1X,A)')
+     &      ' DT_PHOINI: model-switch 30 for nuclei re-set !'
+         CALL PHO_SETMDL(30,0,1)
+      ENDIF
+*
+C     IF (IJP.EQ.7) THEN
+C        AMP2  = SIGN(PMASS(1)**2,PMASS(1))
+C        PP(3) = PPCM
+C        PP(4) = SQRT(AMP2+PP(3)**2)
+C     ELSE
+C        PFERMX = ZERO
+C        IF (IP.GT.1) PFERMX = 0.5D0
+C        EFERMX = SQRT(PFERMX**2+PMASS(1)**2)
+C        CALL DT_LTNUC(PFERMX,EFERMX,PP(3),PP(4),2)
+C     ENDIF
+C     PFERMX = ZERO
+C     IF ((IT.GT.1).OR.(NCOMPO.GT.0)) PFERMX = -0.5D0
+C     EFERMX = SQRT(PFERMX**2+PMASS(2)**2)
+C     CALL DT_LTNUC(PFERMX,EFERMX,PT(3),PT(4),3)
+**sr 26.10.96
+      ISAV = IPAMDL(13)
+      IF ((ISHAD(2).EQ.1).AND.
+     &   ((IJPROJ.EQ. 7).OR.(IJPROJ.EQ.3).OR.(IJPROJ.EQ.4).OR.
+     &    (IJPROJ.EQ.10).OR.(IJPROJ.EQ.11))) IPAMDL(13) = 1
+**
+
+      CALL PHO_EVENT(-1,PP,PT,SIGMAX,IREJ1)
+
+**sr 26.10.96
+      IPAMDL(13) = ISAV
+**
+*
+* patch for cascade calculations:
+* define parton distribution functions for other hadrons, i.e. other
+* then defined already in phojet
+      IF (IOGLB.EQ.100) THEN
+         WRITE(LOUT,1006)
+ 1006    FORMAT(/,1X,'PHOINI: additional parton distribution functions',
+     &          ' assiged (ID,IPAR,ISET)',/)
+         NPDF = 0
+         DO 3 I=1,30
+            IF (IPARPD(I).NE.0) THEN
+               NPDF = NPDF+1
+               IF (NPDF.GT.MAXPDF) STOP ' PHOINI: npdf > maxpdf !'
+               IF ((IPARPD(I).GT.0).AND.(ISETPD(I).GT.0)) THEN
+