#!/usr/bin/perl -w

# Perl script (sorry) for generating rattests for all radioactive decay isotopes.
# A directory is made for each isotope and the relevant .mac, .C and .config files are written to it. 
# Some tweaking of ranges is required for the timing plots (namely BiPo and Kr85)

@decay0_iso = ("Ra228","Ac228","Pb212","Bi212","Tl208","Po212","BiPo212","Th234","Pa234m","Bi214","Po214","Pb210","Bi210",
               "BiPo214","C14","K40","Kr85","Ar39","Y88","Sb126","Ar42","Na22","Rh106","K42","Y90","Co60","Bi212Tl208","C11","C10","Be11","He6");
@decaychain_iso_alpha = ("238U","234U","230Th","226Ra","222Rn","218Po","210Po","232Th","228Th","224Ra","220Rn","216Po");
@decaychain_iso_gamma = ("88Zr","102Rhm","44Ti","121TeGS","123TeMS_IT","125TeMS_IT");
@decaychain_iso_beta = ("214Pb","210Tl","124Sb","110Agm","110Ag","126SntoMS","MSto126Te","102RhChain","26AlChain","56CoChain","58CoChain",
"68GaChain","82RbChain","84RbChain","214BiTlChain","106Ru","44ScChain","90Sr","16N","57CoChain","24Na","46Sc","47Sc","48Sc","121TeMSChain","127TeMSChain","127TeGS","129TeMSChain","129TeGS");

# no geant4 isotopes at the moment, if you add something here, don't forget to uncomment the geant4 section at the end of this file
#@geant4_iso = ("ISO");

foreach $d0 (@decay0_iso){
    if(!-e $d0){
        $command = "mkdir $d0";
        system($command);
    }
    $config = "$d0/rattest.config";
    $macro = "$d0/$d0".".mac";
    $rootscr = "$d0/$d0".".C";
    open (CONFIG, ">$config");
        print CONFIG "# -*- python-*-\n";
        print CONFIG "description = '''A simple decay of $d0 using decay0.'''\n";
        print CONFIG "rat_macro = '$d0".".mac'\n";
        print CONFIG "root_macro = '$d0".".C'\n";
    close (CONFIG);
    open (MAC, ">$macro");
        print MAC "#/**\n#*	Rattest Macro for $d0 Decay using decay0,\n#* automatically generated by j.r.wilson\@qmul.ac.uk\n#*/\n";
        print MAC "/rat/physics_list/OmitMuonicProcesses true\n";
        print MAC "/rat/physics_list/OmitHadronicProcesses true\n\n";
        print MAC "#Disable photons and Cerenkov for speed\n";
        print MAC "/rat/physics_list/OmitCerenkov true\n\n";
        print MAC "/run/initialize\n\n";
        print MAC "/rat/physics/scintillation/off\n";

        print MAC "/rat/proc prune\n";
        print MAC "/rat/procset prune \"evs,mc.hits,mcevs,mc.pmts\"\n";

        print MAC "/rat/proc outroot\n";
        print MAC "/rat/procset file \"$d0.root\"\n\n";

        print MAC "/generator/add combo decay0:fill:uniform\n";
        print MAC "/generator/vtx/set backg $d0 \n";
        print MAC "/generator/pos/set inner_av\n\n";

        print MAC "/rat/run/start 10000\n";
        print MAC "exit\n";
    close (MAC);
    open (ROOT, ">$rootscr");
        print ROOT "void make_plots(TFile *event_file, TTree *T, TFile *out_file)\n";
        print ROOT "{\n";
        print ROOT "  TH1F *e_ke = new TH1F(\"e_ke\", \"e- (or e+) Kinetic Energy\",1000,0,3.0);\n";
        print ROOT "  e_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  e_ke->SetXTitle(\"Energy (MeV)\");\n\n";

        print ROOT "  TH1F *gamma_ke = new TH1F(\"gamma_ke\", \"Gamma Kinetic Energy\",1000,0,3.0);\n";
        print ROOT "  gamma_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  gamma_ke->SetXTitle(\"Gamma Energy (MeV)\");\n\n";
  
        print ROOT "  TH1F *alpha_ke = new TH1F(\"alpha_ke\", \"Alpha Kinetic Energy\",1000,0,10.0);\n";
        print ROOT "  alpha_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  alpha_ke->SetXTitle(\"Alpha Energy (MeV)\");\n\n";
  
        print ROOT "  TH1F *time = new TH1F(\"time\", \"Initial Time of Particles\",1000,0,100.0);\n";
        print ROOT "  time->SetYTitle(\"Events per bin\");\n";
        print ROOT "  time->SetXTitle(\"time (ns)\");\n\n";

        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>e_ke\", \"ds.mc.particles.pdgCode==11||ds.mc.particles.pdgCode==-11\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>gamma_ke\", \"ds.mc.particles.pdgCode==22\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>alpha_ke\", \"ds.mc.particles.pdgCode==1000020040\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.mcTime>>time\",\"\",\"goff\");\n\n";

        print ROOT "  e_ke->Write();\n";
        print ROOT "  gamma_ke->Write();\n";
        print ROOT "  alpha_ke->Write();\n";
        print ROOT "  time->Write();\n";  
        print ROOT "}\n";

    close (ROOT);
    
}

foreach $d0 (@decaychain_iso_alpha){
    if(!-e $d0){
        $command = "mkdir $d0";
        system($command);
    }
    $config = "$d0/rattest.config";
    $macro = "$d0/$d0".".mac";
    $rootscr = "$d0/$d0".".C";
    open (CONFIG, ">$config");
        print CONFIG "# -*- python-*-\n";
        print CONFIG "description = '''A simple decay of $d0 using decaychain.'''\n";
        print CONFIG "rat_macro = '$d0".".mac'\n";
        print CONFIG "root_macro = '$d0".".C'\n";
    close (CONFIG);
    open (MAC, ">$macro");
        print MAC "#/**\n#*	Rattest Macro for $d0 Decay using decaychain,\n#* automatically generated by j.r.wilson\@qmul.ac.uk\n#*/\n";
        print MAC "/rat/physics_list/OmitMuonicProcesses true\n";
        print MAC "/rat/physics_list/OmitHadronicProcesses true\n\n";
        print MAC "#Disable photons and Cerenkov for speed\n";
        print MAC "/rat/physics_list/OmitCerenkov true\n\n";
        print MAC "/run/initialize\n\n";
        print MAC "/rat/physics/scintillation/off\n";

        print MAC "/rat/proc prune\n";
        print MAC "/rat/procset prune \"evs,mc.hits,mcevs,mc.pmts\"\n";

        print MAC "/rat/proc outroot\n";
        print MAC "/rat/procset file \"$d0.root\"\n\n";

        print MAC "/generator/add decaychain $d0:fill:uniform:alpha\n";
        print MAC "/generator/pos/set inner_av\n\n";

        print MAC "/rat/run/start 10000\n";
        print MAC "exit\n";
    close (MAC);
    open (ROOT, ">$rootscr");
        print ROOT "void make_plots(TFile *event_file, TTree *T, TFile *out_file)\n";
        print ROOT "{\n";
        print ROOT "  TH1F *e_ke = new TH1F(\"e_ke\", \"e- (or e+) Kinetic Energy\",1000,0,3.0);\n";
        print ROOT "  e_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  e_ke->SetXTitle(\"Energy (MeV)\");\n\n";

        print ROOT "  TH1F *gamma_ke = new TH1F(\"gamma_ke\", \"Gamma Kinetic Energy\",1000,0,3.0);\n";
        print ROOT "  gamma_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  gamma_ke->SetXTitle(\"Gamma Energy (MeV)\");\n\n";
  
        print ROOT "  TH1F *alpha_ke = new TH1F(\"alpha_ke\", \"Alpha Kinetic Energy\",1000,0,10.0);\n";
        print ROOT "  alpha_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  alpha_ke->SetXTitle(\"Alpha Energy (MeV)\");\n\n";
  
        print ROOT "  TH1F *time = new TH1F(\"time\", \"Initial Time of Particles\",1000,0,100.0);\n";
        print ROOT "  time->SetYTitle(\"Events per bin\");\n";
        print ROOT "  time->SetXTitle(\"time (ns)\");\n\n";

        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>e_ke\", \"ds.mc.particles.pdgCode==11||ds.mc.particles.pdgCode==-11\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>gamma_ke\", \"ds.mc.particles.pdgCode==22\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>alpha_ke\", \"ds.mc.particles.pdgCode==1000020040\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.mcTime>>time\",\"\",\"goff\");\n\n";

        print ROOT "  e_ke->Write();\n";
        print ROOT "  gamma_ke->Write();\n";
        print ROOT "  alpha_ke->Write();\n";
        print ROOT "  time->Write();\n";  
        print ROOT "}\n";

    close (ROOT);    
}
foreach $d0 (@decaychain_iso_gamma){
    if(!-e $d0){
        $command = "mkdir $d0";
        system($command);
    }
    $config = "$d0/rattest.config";
    $macro = "$d0/$d0".".mac";
    $rootscr = "$d0/$d0".".C";
    open (CONFIG, ">$config");
        print CONFIG "# -*- python-*-\n";
        print CONFIG "description = '''A simple decay of $d0 using decaychain.'''\n";
        print CONFIG "rat_macro = '$d0".".mac'\n";
        print CONFIG "root_macro = '$d0".".C'\n";
    close (CONFIG);
    open (MAC, ">$macro");
        print MAC "#/**\n#*	Rattest Macro for $d0 Decay using decaychain,\n#* automatically generated by j.r.wilson\@qmul.ac.uk\n#*/\n";
        print MAC "/rat/physics_list/OmitMuonicProcesses true\n";
        print MAC "/rat/physics_list/OmitHadronicProcesses true\n\n";
        print MAC "#Disable photons and Cerenkov for speed\n";
        print MAC "/rat/physics_list/OmitCerenkov true\n\n";
        print MAC "/run/initialize\n\n";
        print MAC "/rat/physics/scintillation/off\n";

        print MAC "/rat/proc prune\n";
        print MAC "/rat/procset prune \"evs,mc.hits,mcevs,mc.pmts\"\n";

        print MAC "/rat/proc outroot\n";
        print MAC "/rat/procset file \"$d0.root\"\n\n";

        print MAC "/generator/add decaychain $d0:fill:uniform:gamma\n";
        print MAC "/generator/pos/set inner_av\n\n";

        print MAC "/rat/run/start 10000\n";
        print MAC "exit\n";
    close (MAC);
    open (ROOT, ">$rootscr");
        print ROOT "void make_plots(TFile *event_file, TTree *T, TFile *out_file)\n";
        print ROOT "{\n";
        print ROOT "  TH1F *e_ke = new TH1F(\"e_ke\", \"e- (or e+) Kinetic Energy\",1000,0,3.0);\n";
        print ROOT "  e_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  e_ke->SetXTitle(\"Energy (MeV)\");\n\n";

        print ROOT "  TH1F *gamma_ke = new TH1F(\"gamma_ke\", \"Gamma Kinetic Energy\",1000,0,3.0);\n";
        print ROOT "  gamma_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  gamma_ke->SetXTitle(\"Gamma Energy (MeV)\");\n\n";
  
        print ROOT "  TH1F *alpha_ke = new TH1F(\"alpha_ke\", \"Alpha Kinetic Energy\",1000,0,10.0);\n";
        print ROOT "  alpha_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  alpha_ke->SetXTitle(\"Alpha Energy (MeV)\");\n\n";
  
        print ROOT "  TH1F *time = new TH1F(\"time\", \"Initial Time of Particles\",1000,0,100.0);\n";
        print ROOT "  time->SetYTitle(\"Events per bin\");\n";
        print ROOT "  time->SetXTitle(\"time (ns)\");\n\n";

        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>e_ke\", \"ds.mc.particles.pdgCode==11||ds.mc.particles.pdgCode==-11\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>gamma_ke\", \"ds.mc.particles.pdgCode==22\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>alpha_ke\", \"ds.mc.particles.pdgCode==1000020040\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.mcTime>>time\",\"\",\"goff\");\n\n";

        print ROOT "  e_ke->Write();\n";
        print ROOT "  gamma_ke->Write();\n";
        print ROOT "  alpha_ke->Write();\n";
        print ROOT "  time->Write();\n";  
        print ROOT "}\n";

    close (ROOT);    
}
foreach $d0 (@decaychain_iso_beta){
    if(!-e $d0){
        $command = "mkdir $d0";
        system($command);
    }
    $config = "$d0/rattest.config";
    $macro = "$d0/$d0".".mac";
    $rootscr = "$d0/$d0".".C";
    open (CONFIG, ">$config");
        print CONFIG "# -*- python-*-\n";
        print CONFIG "description = '''A simple decay of $d0 using decaychain.'''\n";
        print CONFIG "rat_macro = '$d0".".mac'\n";
        print CONFIG "root_macro = '$d0".".C'\n";
    close (CONFIG);
    open (MAC, ">$macro");
        print MAC "#/**\n#*	Rattest Macro for $d0 Decay using decaychain,\n#* automatically generated by j.r.wilson\@qmul.ac.uk\n#*/\n";
        print MAC "/rat/physics_list/OmitMuonicProcesses true\n";
        print MAC "/rat/physics_list/OmitHadronicProcesses true\n\n";
        print MAC "#Disable photons and Cerenkov for speed\n";
        print MAC "/rat/physics_list/OmitCerenkov true\n\n";
        print MAC "/run/initialize\n\n";
        print MAC "/rat/physics/scintillation/off\n";

        print MAC "/rat/proc prune\n";
        print MAC "/rat/procset prune \"evs,mc.hits,mcevs,mc.pmts\"\n";

        print MAC "/rat/proc outroot\n";
        print MAC "/rat/procset file \"$d0.root\"\n\n";

        print MAC "/generator/add decaychain $d0:fill:uniform\n";
        print MAC "/generator/pos/set inner_av\n\n";

        print MAC "/rat/run/start 10000\n";
        print MAC "exit\n";
    close (MAC);
    open (ROOT, ">$rootscr");
        print ROOT "void make_plots(TFile *event_file, TTree *T, TFile *out_file)\n";
        print ROOT "{\n";
        print ROOT "  TH1F *e_ke = new TH1F(\"e_ke\", \"e- (or e+) Kinetic Energy\",1000,0,3.0);\n";
        print ROOT "  e_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  e_ke->SetXTitle(\"Energy (MeV)\");\n\n";

        print ROOT "  TH1F *gamma_ke = new TH1F(\"gamma_ke\", \"Gamma Kinetic Energy\",1000,0,3.0);\n";
        print ROOT "  gamma_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  gamma_ke->SetXTitle(\"Gamma Energy (MeV)\");\n\n";
  
        print ROOT "  TH1F *alpha_ke = new TH1F(\"alpha_ke\", \"Alpha Kinetic Energy\",1000,0,10.0);\n";
        print ROOT "  alpha_ke->SetYTitle(\"Events per bin\");\n";
        print ROOT "  alpha_ke->SetXTitle(\"Alpha Energy (MeV)\");\n\n";
  
        print ROOT "  TH1F *time = new TH1F(\"time\", \"Initial Time of Particles\",1000,0,100.0);\n";
        print ROOT "  time->SetYTitle(\"Events per bin\");\n";
        print ROOT "  time->SetXTitle(\"time (ns)\");\n\n";

        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>e_ke\", \"ds.mc.particles.pdgCode==11||ds.mc.particles.pdgCode==-11\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>gamma_ke\", \"ds.mc.particles.pdgCode==22\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.kineticEnergy>>alpha_ke\", \"ds.mc.particles.pdgCode==1000020040\",\"goff\");\n";
        print ROOT "  T->Draw(\"ds.mc.particles.mcTime>>time\",\"\",\"goff\");\n\n";

        print ROOT "  e_ke->Write();\n";
        print ROOT "  gamma_ke->Write();\n";
        print ROOT "  alpha_ke->Write();\n";
        print ROOT "  time->Write();\n";  
        print ROOT "}\n";

    close (ROOT);    
}
#########################################################################
# since there are no geant4 isotopes at the moment this is commented, just uncomment and add isotop at the beginning of the file
#foreach $d0 (@geant4_iso){
#    if(!-e $d0){
#        $command = "mkdir $d0";
#        system($command);
#    }
#    $config = "$d0/rattest.config";
#    $macro = "$d0/$d0".".mac";
#    $rootscr = "$d0/$d0".".C";
#    open (CONFIG, ">$config");
#        print CONFIG "# -*- python-*-\n";
#        print CONFIG "description = '''A simple decay of $d0 using geant4.'''\n";
#        print CONFIG "rat_macro = '$d0".".mac'\n";
#        print CONFIG "root_macro = '$d0".".C'\n";
#    close (CONFIG);
#    open (MAC, ">$macro");
#        print MAC "#/**\n#*	Rattest Macro for $d0 Decay using geant4,\n#* automatically generated by j.r.wilson\@qmul.ac.uk\n#*/\n";
#        print MAC "/rat/physics_list/OmitMuonicProcesses true\n";
#        print MAC "/rat/physics_list/OmitHadronicProcesses true\n\n";
#        print MAC "#Disable photons and Cerenkov for speed\n";
#        print MAC "/rat/physics_list/OmitCerenkov true\n\n";
#        print MAC "/run/initialize\n\n";
#        print MAC "/rat/physics/scintillation/off\n";
#
#        print MAC "/rat/proc prune\n";
#        print MAC "/rat/procset prune \"evs,mc.hits,mcevs,mc.pmts\"\n";
#
#        print MAC "/rat/proc outroot\n";
#        print MAC "/rat/procset file \"$d0.root\"\n\n";
#
#        print MAC "/generator/add combo gun:fill:uniform\n";
#        print MAC "/generator/vtx/set $d0 0. 0. 0.\n";
#        print MAC "/generator/pos/set inner_av\n\n";
#
#        print MAC "/rat/run/start 10000\n";
#        print MAC "exit\n";
#    close (MAC);
#    open (ROOT, ">$rootscr");
#        print ROOT "void make_plots(TFile *event_file, TTree *T, TFile *out_file)\n";
#        print ROOT "{\n";
#        print ROOT "  TH1F *e_gen = new TH1F(\"e_gen\", \"Total energy deposited in Scintillator\",1000,0,3.0);\n";
#        print ROOT "  e_gen->SetYTitle(\"Events per bin\");\n";
#        print ROOT "  e_gen->SetXTitle(\"Energy (MeV)\");\n\n";
#
#        print ROOT "  TH1F *e_quenched = new TH1F(\"e_quenched\", \"Quenched energy deposited in Scintillator\",1000,0,3.0);\n";
#        print ROOT "  e_quenched->SetYTitle(\"Events per bin\");\n";
#        print ROOT "  e_quenched->SetXTitle(\"Energy (MeV)\");\n\n";
#   
#        print ROOT "  TH1F *time = new TH1F(\"time\", \"Time of energy deposit\",1000,0,100.0);\n";
#        print ROOT "  time->SetYTitle(\"Events per bin\");\n";
#        print ROOT "  time->SetXTitle(\"time (ns)\");\n\n";
#
#        print ROOT "  T->Draw(\"ds.mc.scintEnergyDeposit>>e_gen\", \"\",\"goff\");\n";
#        print ROOT "  T->Draw(\"ds.mc.scintQuenchedEnergyDeposit>>e_quenched\", \"\",\"goff\");\n";
#        print ROOT "  T->Draw(\"ds.mc.intialScintTime>>time\",\"\",\"goff\");\n\n";
#
#        print ROOT "  e_gen->Write();\n";
#        print ROOT "  e_quenched->Write();\n";
#        print ROOT "  time->Write();\n";  
#        print ROOT "}\n";
#
#    close (ROOT);    
#}