Pythia.java

Go to the documentation of this file.

package cedar.jetweb.generator;

import java.util.HashMap;
import java.util.Map;
import java.util.Vector;
import java.util.Map.Entry;
import java.util.Iterator;

import java.io.BufferedWriter;
import java.io.IOException;

import cedar.jetweb.JetWebException;

import cedar.jetweb.model.MCProcessType;
import cedar.jetweb.model.Model;
import cedar.jetweb.db.DBGeneratorManager;
import cedar.jetweb.job.Cut;
import cedar.jetweb.job.CutCollection;

public class Pythia extends Generator {

    private static boolean particleNameMapInit_ = true;
    private static HashMap<String, String> particleNameMap_ =
    new HashMap<String, String>();

    private boolean processNamesInit_ = true;

    private static HashMap<Integer, String> processNameMap_ = null;
    private static HashMap<String, Vector<Integer>> 
    processSelectorMap_ = null;
    private static HashMap<String, Integer> MSTPMap_ = null;

    private static HashMap<String, Vector<Object>> 
    allowedCuts_ = new HashMap<String, Vector<Object>>();
    private static boolean allowedCutsInit_ = true;

    protected Pythia(String version) {
    generatorName = "pythia";
    generatorVersion = version;

    // Here we should read the default parameters from the DB.

    }

    public String getProcessSelector(MCProcessType proc) {

    StringBuffer sel = new StringBuffer("");

    if(processNamesInit_) initProcessNames();

    Vector<Integer> PIds =  proc.getCompatibleProcesses();

    System.out.println("getProcessSelector: PIds = "+
               PIds);

    Vector<Integer> MSTP14 = new Vector<Integer>(PIds.size());

    boolean found = false;

    for(Integer Id: PIds){
        
        MCProcessType subProc = new MCProcessType(Id);
        Vector<Integer> msubs = 
        processSelectorMap_.get(subProc.getMCProcessName());

        Integer mstp = MSTPMap_.get(subProc.getMCProcessName());
        if(mstp!=null &&
           !MSTP14.contains(mstp)){
        MSTP14.add(mstp);
        }
        if(msubs!=null){
        found = true;
        for(Integer msub: msubs){
            sel.append("MSUB " + msub + " = 1\n");
        }
        sel.append("C \n");
        }   
    }
    if(proc.isPhotoproduction()){
        if(MSTP14.size()!=0){
        sel.append("HZMSTP14 ");
        
        for(Integer mstp: MSTP14){
            sel.append(mstp+" ");
        }
        sel.append("\n");
        sel.append("C\n");
        }
        return sel.toString();
    }
    boolean charm = proc.getMCProcessName().equals(MCProcessType.CHARM);
    boolean beauty = proc.getMCProcessName().equals(MCProcessType.BEAUTY);

    if(charm ||beauty){
        
        int kfCode = charm? 4: 5;
        int antiKFCode = -kfCode;
        int jj=0;
        for(int ii=-40; ii!=antiKFCode; ++ii){
        jj = ii+41;
        sel.append("KFIN2 " + jj + " = 0\n");
        }
        jj = antiKFCode + 41;
        sel.append("KFIN2 "+ jj + " = 1\n");
        for(int ii=antiKFCode + 1; ii!=kfCode; ++ii){
        jj = ii + 41;
        sel.append("KFIN2 " + jj +" =0\n");
        }
        jj = kfCode + 41;
        sel.append("KFIN2 " + jj + " = 1\n");
        for(int ii=kfCode +1; ii!=41; ++ii){
        jj = ii+41;
        sel.append("KFIN2 " + jj +" =0\n");
        }
    }

    if (proc.getMCProcessName().equals(MCProcessType.MINBIAS)){
        if (proc.isPP()){       
        //StringBuffer sel = new StringBuffer("MSEL  1\n");
        return sel.toString();
        } else if (proc.isGammaGamma()){
        //StringBuffer sel = new StringBuffer("MSEL  1\n");
        sel.append("HZMSTP14 0  4  1\n");
        return sel.toString();
        } else {
        return "C Unknown process type";
        }

    } else if (proc.getMCProcessName().equals(MCProcessType.HIGHET)||
           proc.getMCProcessName().equals(MCProcessType.CHARM)||
           proc.getMCProcessName().equals(MCProcessType.BEAUTY)){
        
        if (proc.isPhotoproduction()){      
        //StringBuffer sel = new StringBuffer("HZMSTP14 0  1\n");

        //don't think we need hzmstp now because
        //direct and resolved in 2 separate runs

        //sel.append("HZMSTP14 0  1\n");
        //sel.append("MSEL  1\n");
        return sel.toString();
        } else if (proc.isPP()){        
        //StringBuffer sel = new StringBuffer("MSEL  1\n");

        return sel.toString();
        } else if (proc.isGammaGamma()){
        //StringBuffer sel = new StringBuffer("MSEL  1\n");
        sel.append("HZMSTP14 0  4  1\n");
        return sel.toString();
        } else if(!found){
        return "C Unknown process type";
        }
    } else if(!found){
        return "C Unknown process type";
    }else{
        return sel.toString();
    }

    return(sel.toString());
    }

    //    public String getMCProcessName(HashMap<String,HashMap<Integer,Integer>> parameters) 
    public String getMCProcessName(HashMap<String, Object> parameters)
    throws JetWebException{

    //System.out.println("Calling Pythia.getMCProcessName");


    if (processNamesInit_) initProcessNames();

    //  System.out.println("parameters = "+parameters);


    //Integer msel = intParameters.get("MSEL");
    HashMap<Integer, Integer> msub = 
        (HashMap<Integer, Integer>)parameters.get("MSUB");
    HashMap<Integer, HashMap<Integer, Integer>> kfin = 
        (HashMap<Integer, HashMap<Integer, Integer>>)
        parameters.get("KFIN");
    HashMap<Integer, Integer> mstp = 
        (HashMap<Integer, Integer>) parameters.get("MSTP");

    //  System.out.println("kfin = "+kfin);

    if(msub==null||
       //kfin==null||
       mstp==null){
        throw new JetWebException
       ("MSUB and/or KFIN not present in realParameterArrays",
    "cedar.jetweb.generator.pythia: getMCProcessType()");
    }

    Vector<Integer> isubs = new Vector<Integer>();

    Iterator it  = msub.entrySet().iterator();

    while(it.hasNext()){
        Entry ent = (Entry)it.next();
        Integer val = (Integer)ent.getValue();
        if(val!=0){
        Integer isub = (Integer)ent.getKey();
        isubs.add(isub);
        }

    }

    boolean cOnly = false;
    boolean bOnly = false;
    boolean anythingElse = false;

    it = kfin.entrySet().iterator();

    while(it.hasNext()){
        Entry ent = (Entry)it.next();
        
        Integer key = (Integer) ent.getKey();

        if(key.equals(2)){
        Iterator row = 
            ((HashMap<Integer, Integer>) ent.getValue()).entrySet().iterator();
        while(row.hasNext()){
            Entry rowEnt = (Entry)row.next();
            Integer col = (Integer)rowEnt.getKey();
            Integer val = (Integer)rowEnt.getValue();
            
            if(col.equals(4) || col.equals(-4)){
            if(!val.equals(0)) cOnly = true;
            }else if(col.equals(5) || col.equals(-5)){
            if(!val.equals(0)) bOnly = true;
            }else{
            if(!val.equals(0)) anythingElse = true;
            }
            
        }
        }
    }
    
    anythingElse = anythingElse || (cOnly && bOnly);
    cOnly = cOnly && !anythingElse;
    bOnly = bOnly && !anythingElse;

    String processName = processNameMap_.get(isubs.get(0));
    
    System.out.println("Pythia: found process type "+
               processName);

    for(Integer ii=1; ii!= isubs.size(); ++ii){
        if(!processNameMap_.get(isubs.get(ii)).equals(processName)){
        throw new JetWebException
            ("More than one process type present",
             "cedar.jetweb.generator.Pythia: getMCProcessType()");
        }
    }

    if(cOnly){ 
        processName = "Charm "+processName;
    }else if(bOnly){
        processName = "Beauty "+processName;
    }
    
    if(processName.equals(MCProcessType.HIGHET)){
        if(cOnly){
        processName = MCProcessType.CHARM;
        }else if(bOnly){
        processName = MCProcessType.BEAUTY;
        }
    }else if(processName.equals(MCProcessType.DIS)){
        if(cOnly){
        processName = MCProcessType.CHARMDIS;
        }else if(bOnly){
        processName = MCProcessType.BEAUTYDIS;
        }
    }
    return processName;

    }

    private void initAllowedCuts(){

    // Initialises the names of the allowed cuts and maps them onto the
    // equivalent Pythia parameters
    // Only PTMin for the moment

    if(!allowedCutsInit_) return;
    allowedCutsInit_ = false;

    Vector<Object> tmp = new Vector<Object>(2);
    tmp.add(new String("CKIN"));
    tmp.add(new Integer(3));

    allowedCuts_.put("PTMin", tmp);

    return;
    }

    private void initProcessNames(){

    if(!processNamesInit_) return;
    processNamesInit_ = false;

    processSelectorMap_ = new HashMap<String, Vector<Integer>>();
    MSTPMap_ = new HashMap<String, Integer>();

    Vector<Integer> PIds = new Vector<Integer>();

    PIds.add(11);
    PIds.add(12);
    PIds.add(13);
    PIds.add(28);
    PIds.add(53);
    PIds.add(68);

    //processSelectorMap_.put(MCProcessType.HIGHET, PIds);
    processSelectorMap_.put(MCProcessType.PARTONPARTON, PIds);
    processSelectorMap_.put(MCProcessType.CHARMPARTONPARTON, PIds);
    processSelectorMap_.put(MCProcessType.BEAUTYPARTONPARTON, PIds);
    MSTPMap_.put(MCProcessType.PARTONPARTON, 1);
    MSTPMap_.put(MCProcessType.CHARMPARTONPARTON, 1);

    PIds = new Vector<Integer>();
    
    PIds.add(33);
    PIds.add(34);
    PIds.add(54);
    //PIds.add(58);
    PIds.add(80);
    PIds.add(84);
    //PIds.add(85);

    processSelectorMap_.put(MCProcessType.PHOTONPARTON, PIds);
    processSelectorMap_.put(MCProcessType.CHARMPHOTONPARTON, PIds);
    processSelectorMap_.put(MCProcessType.BEAUTYPHOTONPARTON, PIds);
    MSTPMap_.put(MCProcessType.PHOTONPARTON, 0);
    MSTPMap_.put(MCProcessType.CHARMPHOTONPARTON, 0);
    MSTPMap_.put(MCProcessType.BEAUTYPHOTONPARTON, 0);

    PIds = new Vector<Integer>();
    PIds.add(58);
    PIds.add(85);
    processSelectorMap_.put(MCProcessType.PHOTONPHOTON, PIds);
    processSelectorMap_.put(MCProcessType.CHARMPHOTONPHOTON, PIds);
    processSelectorMap_.put(MCProcessType.BEAUTYPHOTONPHOTON, PIds);
    MSTPMap_.put(MCProcessType.PHOTONPHOTON, 4);
    MSTPMap_.put(MCProcessType.CHARMPHOTONPHOTON, 4);
    MSTPMap_.put(MCProcessType.BEAUTYPHOTONPHOTON, 4);

    PIds = new Vector<Integer>();

    PIds.add(10);
    PIds.add(83);
    PIds.add(99);

    for(Integer ii=131; ii!= 137; ++ii){
        PIds.add(ii);
    }

    processSelectorMap_.put(MCProcessType.DIS, PIds);

    PIds = new Vector<Integer>();
    for(Integer ii=91; ii!=96; ++ii){
        PIds.add(ii);
    }

    processSelectorMap_.put(MCProcessType.MINBIAS, PIds);
    MSTPMap_.put(MCProcessType.MINBIAS, 1);
    //  PIds = new Vector<Integer>();

    processNameMap_ = new HashMap<Integer, String>();

    for(Map.Entry entry: processSelectorMap_.entrySet()){
        String procName = (String)entry.getKey();
        if(!procName.startsWith("Charm") &&
           !procName.startsWith("Beauty")){
        PIds = (Vector<Integer>)entry.getValue();
        for(Integer Id: PIds){
            processNameMap_.put(Id, procName);
        }
        }
        
    }
    /*
    processNameMap_.put(11, MCProcessType.HIGHET);
    processNameMap_.put(12, MCProcessType.HIGHET);
    processNameMap_.put(13, MCProcessType.HIGHET);
    processNameMap_.put(28, MCProcessType.HIGHET);
    processNameMap_.put(53, MCProcessType.HIGHET);
    processNameMap_.put(68, MCProcessType.HIGHET);

    processNameMap_.put(33, MCProcessType.DIRECT);
    processNameMap_.put(34, MCProcessType.DIRECT);
    processNameMap_.put(54, MCProcessType.DIRECT);
    processNameMap_.put(58, MCProcessType.DIRECT);
    processNameMap_.put(80, MCProcessType.DIRECT);
    processNameMap_.put(84, MCProcessType.DIRECT);
    processNameMap_.put(85, MCProcessType.DIRECT);

    processNameMap_.put(10, MCProcessType.DIS);
    processNameMap_.put(83, MCProcessType.DIS);
    processNameMap_.put(99, MCProcessType.DIS);

    for(Integer ii=131; ii !=137; ++ii){
        processNameMap_.put(ii, MCProcessType.DIS);
    }

    processNameMap_.put(91, MCProcessType.MINBIAS);
    processNameMap_.put(95, MCProcessType.MINBIAS);
    */
    return;
    }



    public boolean checkParameters(){
    return true;
    }

    protected String getParticleName(String pName){
    if(particleNameMapInit_) initParticleNames();

    String name = particleNameMap_.get(pName);
    
    if(name==null) return pName;

    return name;
    }

    private static void initParticleNames(){
    
    particleNameMapInit_ = false;

    particleNameMap_.put("e+","gamma/e+");
    particleNameMap_.put("e","gamma/e");

    return;
    }

    public void writeSteering(BufferedWriter ofile, Model model, 
                  MCProcessType proc) 
    throws JetWebException {

    try {
        // Common, except PDFs
        writeCommonFFKey(ofile,proc);
        // Turn on the appropriate HZTool routines
        writePapersFFKey(ofile,proc);
        // PDFs
        writePDFFFKey(ofile,model,proc);
        // All the other parameters
        writeGenFFKey(ofile,proc);

        
    } catch (IOException io) {
        throw new JetWebException(io,"Pythia.writeSteering failed in write to "+ofile);
    }

    }


    public boolean inModel(String pName){
    if (pName.equalsIgnoreCase("MSUB")) return false;
    if (pName.equalsIgnoreCase("KFIN")) return false;
    if (pName.equalsIgnoreCase("MSEL")) return false;
    if (pName.equalsIgnoreCase("CKIN")) return false;
    return true;
    }

    public boolean inModel(String aName, Integer index){
    
    if(!inModel(aName)) return false;

    //MSTP(14) controls photon structure (resolved/direct).
    //Jetweb steers this by the process type selected e.g. 
    //in photoproduction does 2 runs with MSTP(14)=0, 1
    //Therefore need to exclude it from the model

    if(aName.equalsIgnoreCase("MSTP") && index==14)return false;

    //MSTP(51-56) control PDFs.  Jetweb steers those itself via
    //lhapdf.  We therefore exclude these from the model as well

    if(aName.equalsIgnoreCase("MSTP")){
        if(index>50 && index <57) return false;
    }

    return true;
    }

    public boolean inProcessType(String pName){
    if (pName.equalsIgnoreCase("MSUB")) return true;
    if (pName.equalsIgnoreCase("KFIN")) return true;
    if (pName.equalsIgnoreCase("MSTP")) return true;
    if (pName.equalsIgnoreCase("MSEL")) return true;
    return false;
    }

    public boolean inCuts(String pName){
    if(allowedCutsInit_) initAllowedCuts();
    for(Map.Entry entry: allowedCuts_.entrySet()){
        Vector<Object> name = (Vector<Object>)entry.getValue();
        if(((String)name.get(0)).equalsIgnoreCase(pName)){
        return true;
        }
    }
    
    return false;
    }

    public Pythia addCut(Cut cut){
    if(allowedCutsInit_) initAllowedCuts();

    if(cuts_==null) cuts_ = new CutCollection();

    Vector<Object> cutName = allowedCuts_.get(cut.getName());
    //if cut not known by Pythia then return

        if(cutName==null) return this;
    
    // if existing cut in the collection supersedes this one 
    // then don't do anything
    if(!cuts_.addCut(cut)) return this;

    String arName = (String) cutName.get(0);
    Integer arIndex = (Integer) cutName.get(1);
    HashMap<Integer, Double> array = realParameterArrays.get(arName);
    if(array==null){
        array = new HashMap<Integer, Double>();
        realParameterArrays.put(arName, array);
    }

    //RealParameterArrays().get(arName).put(arIndex, cut.getValue());

    array.put(arIndex, cut.getValue());

    return this;
    }

    public CutCollection getCutCollection(){

    if(!dirtyCuts_) return cuts_;
    dirtyCuts_ = false;
    
    cuts_ = new CutCollection();

    if(allowedCutsInit_) initAllowedCuts();

    for(Map.Entry entry: allowedCuts_.entrySet()){
        Vector<Object> param = (Vector<Object>)entry.getValue();
        String name = (String)entry.getKey();
        String arName = (String)param.get(0);
        Integer arIndex = (Integer)param.get(1);
        
        Double value = realParameterArrays.get(arName).get(arIndex);

        Cut tmp = new Cut(name, value);
        cuts_.addCut(tmp);
    }
    return cuts_;
    }

    public Double getPthat(){

    HashMap<Integer,Double> ckin = realParameterArrays.get("CKIN");

    return  ckin.get(new Integer(3)); 

    }

    /*
    public boolean isRelevant(MCProcessType proc, String pName, 
                  Integer aIndex) 
    throws JetWebException {

    if (proc.hasPhoton()){
        if (DBGeneratorManager.checkPhotonRelevance(this,pName,aIndex)){
        return true;
        }
    }
    if (proc.hasProton()){
        if (DBGeneratorManager.checkProtonRelevance(this,pName,aIndex)){
        return true;
        }
    }
    return false;
    }
    */


    /*  
    if (getGenerator().equalsIgnoreCase("pythia")){

        // If the energy dependence on the ptmin is turned of,
        // parp(89) is irrelevant, so set it to the default.
        Double edepend = (Double)(getRunSeries().getParp().get(new Integer(90)));    
        if (edepend.doubleValue() == 0.0) { 
        getRunSeries().getParp().remove(new Integer(89));
        getRunSeries().getParp().put(new Integer(89),new Double(1000.0));
        } 
        
    }
    */


    public static void init(){

    /*
    //
    // MSTP
    //
    eeIrrelevantMstp = new Vector();
    epIrrelevantMstp = new Vector();
    ppIrrelevantMstp = new Vector();

    // Irrelevant for e+e-
    eeIrrelevantMstp.add(new Integer(14));
    eeIrrelevantMstp.add(new Integer(51));
    eeIrrelevantMstp.add(new Integer(52));
    eeIrrelevantMstp.add(new Integer(142));
    eeIrrelevantMstp.add(new Integer(191));

    // Irrelevant for ep
    epIrrelevantMstp.add(new Integer(14));
    epIrrelevantMstp.add(new Integer(48));
    epIrrelevantMstp.add(new Integer(142));
    epIrrelevantMstp.add(new Integer(191));

    // Irrelevant for proton beams
    ppIrrelevantMstp.add(new Integer(11));
    ppIrrelevantMstp.add(new Integer(12));
    ppIrrelevantMstp.add(new Integer(13));
    ppIrrelevantMstp.add(new Integer(14));
    ppIrrelevantMstp.add(new Integer(15));
    ppIrrelevantMstp.add(new Integer(16));
    ppIrrelevantMstp.add(new Integer(17));
    ppIrrelevantMstp.add(new Integer(18));
    ppIrrelevantMstp.add(new Integer(19));
    ppIrrelevantMstp.add(new Integer(20));
    ppIrrelevantMstp.add(new Integer(22));
    ppIrrelevantMstp.add(new Integer(23));
    ppIrrelevantMstp.add(new Integer(48));
    ppIrrelevantMstp.add(new Integer(55));
    ppIrrelevantMstp.add(new Integer(56));
    ppIrrelevantMstp.add(new Integer(59));
    ppIrrelevantMstp.add(new Integer(60));
    ppIrrelevantMstp.add(new Integer(66));
    ppIrrelevantMstp.add(new Integer(142));
    ppIrrelevantMstp.add(new Integer(191));

    //
    // PARP
    //
    epIrrelevantParp = new Vector();
    eeIrrelevantParp = new Vector();
    ppIrrelevantParp = new Vector();

    // Irrelevant for all
    // none

    // Irrelevant for e+e-
    eeIrrelevantParp.add(new Integer(91));
    eeIrrelevantParp.add(new Integer(92));
    eeIrrelevantParp.add(new Integer(93));


    // Irrelevant for ep

    // Irrelevant for pp
    ppIrrelevantParp.add(new Integer(13));
    ppIrrelevantParp.add(new Integer(15));
    ppIrrelevantParp.add(new Integer(16));
    ppIrrelevantParp.add(new Integer(17));
    ppIrrelevantParp.add(new Integer(18));
    ppIrrelevantParp.add(new Integer(99));
    ppIrrelevantParp.add(new Integer(100));


    //
    // MSTJ ??? Are any parameters irrelevant?
    //
    eeIrrelevantMstj = new Vector();
    epIrrelevantMstj = new Vector();
    ppIrrelevantMstj = new Vector();

    //
    // PARJ ??? Are any parameters irrelevant?
    //
    eeIrrelevantParj = new Vector();
    epIrrelevantParj = new Vector();
    ppIrrelevantParj = new Vector();
    */
    }
    /*


    // Everything below is PYTHIA specific and belong in the eventual
    // PYTHIA class.


    private static Vector epIrrelevantMstp;
    private static Vector eeIrrelevantMstp;
    private static Vector ppIrrelevantMstp;
    private static Vector epIrrelevantParp;
    private static Vector eeIrrelevantParp;
    private static Vector ppIrrelevantParp;
    private static Vector epIrrelevantMstj;
    private static Vector eeIrrelevantMstj;
    private static Vector ppIrrelevantMstj;
    private static Vector epIrrelevantParj;
    private static Vector eeIrrelevantParj;
    private static Vector ppIrrelevantParj;

Generated Wed Jan 17 09:14:27 GMT 2007