#!/bin/bash
# Run this script to compute diphoton production cross sections 
# at the LHC. The procedure for the computation is as follows. 
# The configuration files (resbos.in) 
# for each subchannel are copied from the subdirectory ./ins, 
# and the key variables are modified according to the values specified 
# in this script. Computation in ResBos and matching of resummed cross sections
# is performed for each subchannel. The final histograms 
# or ntuples for all subchannels are summed up at the end  of 
# all runs to produce the final result. 
#
# When editing the script, separate new commands by semicolons and 
# pay attention to spaces. They can make a huge difference!
#
# Author:  Pavel Nadolsky
##############################################################################

# ===========================================================================
# A. User-specified parameters
# ===========================================================================

# 1. Type of the cross section: resummed (res) or NLO (nlo)
TypeXSec=res;

# 2. Compute the qq + qg -> gamma gamma X cross section (yes/no)?
ComputeQQ=yes;

# 3. If ComputeQQ=no, go to section 4
# Computation of the qq resummed cross section requires two input grids, 
# 'w321_aa_'$NameBaseQQ'.out' and 'y_aa_'$NameBaseQQ'.out', where NameBaseQQ 
# is the base of the grid's name (e.g., 'lhc_ct6m' for
# 'w321_aa_lhc_ct6m.out'). Computation of the NLO cross section requires
# one input grid, 'dsi_aa_'$NameBaseQQ'.out'. Specify the name base here.
NameBaseQQ="lhc_ct6m";

# 4. Compute the gg -> gamma gamma X cross section (yes/no)?
ComputeGG=no;

# 5. Computation of the resummed cross section requires two input grids, 
# 'w321_ag_'$NameBaseGG'.out' and 'y_ag_'$NameBaseGG'.out', where NameBaseGG 
# is the base of the grid's name (e.g., 'lhc_ct6m' for 
# 'w321_ag_lhc_ct6m.out'). Computation of the NLO cross section requires 
# one input grid, 'dsi_ag_'$NameBaseGG'.out'. Specify the name base here.
NameBaseGG="lhc_ct6m";

#6. Combine qq and gg channels if both are computed (yes/no)? 
#Not active if only one channel is computed.
CombineQQGG=yes;

 #Output type:
# "Gbook"  -- GBOOK (ASCII tables) 
# "HBook"  -- PAW HBOOK histograms 
# "PAWNT1" -- PAW weighted ntuple, with kinematics of each event 
#             given in terms of Q, QT, y, cos(theta_*), phi_*
# "PAWNT2" -- PAW weighted ntuple, with kinematics of each event 
#              given in terms of 4-momenta of two decay particles
# "ROOTNT1"-- ROOT weighted ntuple, with kinematics of each event 
#             given in terms of Q, QT, y, cos(theta_*), phi_*
# "ROOTNT2"  -- ROOT weighted ntuple, with kinematics of each event 
#              given in terms of 4-momenta of two decay particles
OutputType="ROOTNT2";             


# 7. Cuts on the final-state particles
#--------------------------------------------------
yGammaMax=2.5;               #Maximal absolute value for the rapidity 
                             #of the photon 
pTMin1=40.0;                 #Minimal p_T of the photon with a larger p_T
pTMin2=25.0;                 #Minimal p_T of the photon with a smaller P_T


DelR=0.3;                    #Minimal Del R separation between the photons 

iQQT=0;                      #Set iQQT=1 (2) to impose the cut Q > QT (Q < QT)
DelRiso=0.4;                 #The size (Del R) of the isolation cone 
ETMaxIso=15.0;               #Maximal hadronic transverse energy allowed 
                             #inside the isolation cone

QMin=35;     #Minimal invariant mass of the diphoton (>= 40 GeV)
QMax=250;    #Maximal invariant mass of the diphoton (<= 200 GeV)
qTMin=0.0;   #Minimal transverse momentum of the diphoton 
qTMax=300.;  #Maximal transverse momentum of the diphoton (> 20 GeV)
yMin=-9;     #Minimal rapidity of the diphoton
yMax=9;      #Maximal rapidity of the diphoton


# 8. Other parameters
#--------------------------------------------------
GridDir="./grids/";     #Directory with the input grids
InsDir="./ins/lhc/"; #Directory with the configuration (.in) files
LOGFILE="/dev/stdout";   #Name of the log file or /dev/stdout to print
                         #to the screen

# End of user-specified parameters =====================================

#=======================================================================
# B. Internal parameters (do not modify if not sure)
#=======================================================================

# The value of QT below which the resummed cross section is approximated
# entirely by the W piece, or by Delta Sigma piece (in GeV)
qTSepQQ=1.0;  #For the qq channel
qTSepGG=2.0;  #For the gg channel  

#To optimize Monte-Carlo integration, the W-A grid is integrated in 2 steps,
# first over the range qT < qTInt ~ 10-20 GeV, then over the range qT > qTInt.
#Specify the values of qTInt for the qq and gg channels here
qTIntQQ=20;
qTIntGG=20;

#Type of regularization at small Q_T not affected by the isolation 
#(Q_T < ETMaxIso)
# isub=12 : resummation-like subtraction inside the isolation cone 
#           (preferred)
#     =-2 : smooth-cone isolation (reasonable alternative form) 
isub=12

#The line containing qT_sep (to be inserted later and omitted here),  PDF's, 
#type of the process, etc.
LineqTSep=', 61, 0, 5, 4, '$isub', 0, 5.0                  > qTSep,PDF,Prc,Sc,ACt,iSb,iHQA';

#=======================================================================
# C. Subroutines
#=======================================================================
repline1()
# Replaces a line LINE1 by a line LINE2 in the list of files FILELIST
{

[ -z "$1" ] &&
{
echo "EXIT in the subroutine repline1: no parameters provided";
exit;
}
FILELIST=$1;

# If LINE1 and LINE2 are empty, exit
 [ -z "$2" ] &&
 { echo "EXIT in the subroutine repline1: no lines to process"
   exit;
 }

 LINE1=$2
 LINE2=$3


# Names of temporary files
 TEM_NAME="df381uz412aj"
 TEMFILE1=$TEM_NAME"1"
 TEMFILE2=$TEM_NAME"2"
 TEMFILE3=$TEM_NAME"3"
[ -f $TEMFILE1 ] && rm $TEMFILE1;
[ -f $TEMFILE2 ] && rm $TEMFILE2;
[ -f $TEMFILE3 ] && rm $TEMFILE3;

# Loop over the file names
 I1=0
 J1=0
for i in $FILELIST ; do

# FILELEN1 and FILELEN2 are lengths of TEMFILE1 and TEMFILE2
  FILELEN1=0
  FILELEN2=0

# Find the positions of the LINE1 in the file $i
  LINENUM=$( nl -ba -w6 $i | grep -i "$LINE1" | cut -c1-6)

 cp $i $TEMFILE2
 cp $i $i"~"
 # Replace all the occurencies of LINE1 in the file $i by LINE2
 for j in $LINENUM; do
    [ -f $TEMFILE1 ] && 
     FILELEN1=($(wc -l $TEMFILE1))
     FILELEN2=($(wc -l $TEMFILE2))
     HEADLEN=$((j-FILELEN1-1));
     TAILLEN=$((FILELEN2-HEADLEN-1));
     
    head -n$HEADLEN $TEMFILE2 >> $TEMFILE1
    echo "$LINE2" >> $TEMFILE1
    tail -n$TAILLEN $TEMFILE2 > $TEMFILE3
    mv $TEMFILE3 $TEMFILE2  

     J1=$((J1+1))     
 done # for j
 
 cat $TEMFILE2 >> $TEMFILE1

 [ -f $TEMFILE1 ] && mv $TEMFILE1 $i;
 [ -f $TEMFILE2 ] && rm $TEMFILE2;
 [ -f $TEMFILE3 ] && rm $TEMFILE3;

  I1=$((I1+1))
done  # for i

} #repline1 ->

determine_book_type()
# Set parameters according to the type of the output
# (gbook, hbook, or ntuple)
{
  if [ $OutputType == "Gbook" ]; then
    BOOKEXT=".gbook";
    RESBOSBOOK="fort.36";
    BOOKPLUS="gbook+";      #The tool to add gbooks
  elif [[ $OutputType == "HBook" || $OutputType == "PAWNT1" ||  $OutputType == "PAWNT2" ]]; then
    BOOKEXT=".hbook";
    RESBOSBOOK="resbos.hbook";
    BOOKPLUS="hbook+";      #The tool to add PAW hbooks and ntuple
  elif [[ $OutputType == "ROOTNT1" || $OutputType == "ROOTNT2" ]]; then 
    BOOKEXT=".root"
    RESBOSBOOK="resbos.root";
    BOOKPLUS=$TOOLDIR/"treemerge"      #The tool to add ROOT ntuples
  else 
    echo determine_book_type: unknown output format = $OutputType;
    exit
  fi

  [ -e $BOOKPLUS ] || make $BOOKPLUS;


} #determine_book_type ->

Modify_ResBos_In1()
# Writes modified parameters into the ResBos.In file, part 1
{
repline1 resbos.in "> Main data grid" $PrimaryGrid"                                   > Main data grid";

repline1 resbos.in "> Y piece grid" $SecondaryGrid"                                             > Y piece grid (same PDF or -)";

} # Modify_ResBos_In1 ->

Modify_ResBos_In2()
# Writes modified parameters into the ResBos.In file, part 2
{
repline1 resbos.in "> Cuts(1)" \
$pTMin1", "$yGammaMax", "$DelR", "$pTMin2", "$yGammaMax"		> Cuts(1): pT1,y1,DelR,pT2,y2";

repline1 resbos.in "> Cuts(2)" \
$QMin", "$QMax", "$qTMin", "$qTMax", "$yMin", "$yMax"          > Cuts(2): QMn,QMx,QTMn,QTMx,yWmn,yWmx"

repline1 resbos.in "> Cuts(3)" $iQQT", "$DelRiso", "$ETMaxIso"                                    > Cuts(3):iQQT,iso DelR,iso max ET";

repline1 resbos.in "> Output" $OutputType"                                        > Output fromat";

repline1 resbos.in "> qT" $"$qTSep$LineqTSep";

} # Modify_ResBos_In2 ->


#=======================================================================
# C. Main calculation
#=======================================================================

[ -e resbos ] || make resbos;
determine_book_type;

[ -f $LOGFILE ] && [ $LOGFILE != /dev/stdout ] && rm $$LOGFILE
qTMinFull=$qTMin;
qTMaxFull=$qTMax;

                                                 #Compute the qq channel
[ $ComputeQQ == "yes" ] && 
{
    qTSep=$qTSepQQ;
    qTInt=$qTIntQQ;

    if [ $TypeXSec == "res" ]; then
	{
	    qTMin=$qTSep;
	    qTMax=$qTInt;     
	    echo "Computing the resummed qq cross section at qT < "$qTMax" GeV)...";
	    cp $InsDir"aa_w2211_ct6m.in" resbos.in;
	    
	    PrimaryGrid=$GridDir"w321_aa_"$NameBaseQQ".out";
#       SecondaryGrid="-";
	    SecondaryGrid=$GridDir"y_aa_"$NameBaseQQ".out"; 

	    QQBOOK="res321_aa_"$NameBaseQQ$BOOKEXT;

	    Modify_ResBos_In1; Modify_ResBos_In2;
	    nice resbos >> $LOGFILE;
	    mv $RESBOSBOOK "tmp01"$BOOKEXT;
	    
	    qTMin=$qTInt;
	    qTMax=$qTMaxFull;

	    echo "Computing the resummed qq cross section at "$qTMin" < qT < "$qTMax" GeV)...";

	    Modify_ResBos_In1; Modify_ResBos_In2;
	    nice resbos >> $LOGFILE;
	    $BOOKPLUS $RESBOSBOOK "tmp01"$BOOKEXT "tmp1"$BOOKEXT;
	    cp "tmp1"$BOOKEXT "w321-a_aa_"$NameBaseQQ$BOOKEXT;
	} 
    else # $TypeXSec == "nlo"
	{
	    echo "Computing the small-pT NLO qq cross section...";
	    cp $InsDir"aa_dsi_ct6m.in" resbos.in;
	    
	    PrimaryGrid=$GridDir"dsi_aa_"$NameBaseQQ".out";
	    SecondaryGrid="-";
	    QQBOOK="nlo_aa_"$NameBaseQQ$BOOKEXT;
	    
	    qTMin=0.0;
	    qTMax=5.0;
	    Modify_ResBos_In1; 
	    Modify_ResBos_In2;
	    nice resbos >>  $LOGFILE;
	    cp $RESBOSBOOK "tmp1"$BOOKEXT;
	    mv $RESBOSBOOK "dsi_aa_"$NameBaseQQ$BOOKEXT;
	    
	}
    fi  # [ TypeXSec ->

    echo "Computing the large-pT NLO qq cross section";

    qTMin=0.0;
    qTMax=$qTMaxFull;

    cp $InsDir"aa_ps2_ct6m.in" resbos.in;
    Modify_ResBos_In2; #part 1 is not needed
    nice resbos >> $LOGFILE;

                                                  #Add the qq hbooks
    [ -f $QQBOOK ] && rm $QQBOOK;
    $BOOKPLUS $RESBOSBOOK "tmp1"$BOOKEXT $QQBOOK;
    mv $RESBOSBOOK "pert_aa_"$NameBaseQQ$BOOKEXT;

    [ -f $QQBOOK ] && echo The output for the qq channel is written to $QQBOOK;
 
} # [ $ComputeQQ ...

                                                 #Compute the ag channel
[ $ComputeGG == "yes" ] && 
{
    qTSep=$qTSepGG;
    qTInt=$qTIntGG;

    if [ $TypeXSec == "res" ]; then
	{
	    qTMin=$qTMinFull;
	    qTMax=$qTSep;     
	    echo "Computing the resummed gg cross section at qT < "$qTMax" GeV)...";
	    cp $InsDir"ag_w2211_ct6m.in" resbos.in;
	    
	    PrimaryGrid=$GridDir"w321_ag_"$NameBaseGG".out";
	    SecondaryGrid="-";
	    GGBOOK="res321_ag_"$NameBaseGG$BOOKEXT;
	    
	    Modify_ResBos_In1; Modify_ResBos_In2;
	    nice resbos >> $LOGFILE;
	    mv $RESBOSBOOK "tmp01"$BOOKEXT;

	    qTMin=$qTSep;
	    qTMax=$qTInt;     
	    echo "Computing the resummed gg cross section at "$qTMin" < qT < "$qTMax" GeV)...";
	    SecondaryGrid=$GridDir"y_ag_"$NameBaseGG".out";

	    Modify_ResBos_In1; Modify_ResBos_In2;
	    nice resbos >> $LOGFILE;
	    mv $RESBOSBOOK "tmp02"$BOOKEXT;

	    qTMin=$qTInt;
	    qTMax=$qTMaxFull;     
	    echo "Computing the resummed gg cross section at "$qTMin" < qT < "$qTMax" GeV)...";

	    Modify_ResBos_In1; Modify_ResBos_In2; 
	    nice resbos >> $LOGFILE;
	    $BOOKPLUS $RESBOSBOOK "tmp01"$BOOKEXT "tmp1"$BOOKEXT;
	    $BOOKPLUS "tmp1"$BOOKEXT "tmp02"$BOOKEXT "tmp2"$BOOKEXT;
	    cp "tmp2"$BOOKEXT "w321-a_ag_"$NameBaseGG$BOOKEXT;
   } 
  else # $TypeXSec == "nlo"
   {
       echo "Computing the small-pT NLO gg cross section...";
       cp $InsDir"ag_dsi_ct6m.in" resbos.in;
       PrimaryGrid=$GridDir"dsi_ag_"$NameBaseGG".out";
       SecondaryGrid="-";
       GGBOOK="nlo_ag_"$NameBaseGG$BOOKEXT;

       qTMin=0.0;
       qTMax=5.0;
       Modify_ResBos_In1; Modify_ResBos_In2;
       nice resbos >> $LOGFILE;
       cp $RESBOSBOOK "tmp2"$BOOKEXT;
       mv $RESBOSBOOK "dsi_ag_"$NameBaseGG$BOOKEXT;
       exit;
   }
  fi  # [ TypeXSec ->

    qTMin=$qTMinFull;
    qTMax=$qTMaxFull;

    echo "Computing the large-pT NLO gg cross section";
    cp $InsDir"ag_ps2_ct6m.in" resbos.in;
    Modify_ResBos_In2; #part 1 is not needed
    nice resbos >> $LOGFILE;

                                                  #Add the GG hbooks
    [ -f $GGBOOK ] && rm $GGBOOK;
    $BOOKPLUS $RESBOSBOOK "tmp2"$BOOKEXT $GGBOOK;
    mv $RESBOSBOOK "pert_ag_"$NameBaseGG$BOOKEXT;

    [ -f $GGBOOK ] && echo The output for the gg channel is written to $GGBOOK; 
} # [ $ComputeGG ...

[[ $ComputeQQ == "yes" && $ComputeGG == "yes" && $CombineQQGG == "yes" ]] &&
{
    if [ $TypeXSec == "res" ]; then
	FINALBOOK="res321_all_"$NameBaseQQ$BOOKEXT;
    else
	FINALBOOK="nlo_all_"$NameBaseQQ$BOOKEXT;
    fi

    [ -f $FINALBOOK ] && rm $FINALBOOK;
    $BOOKPLUS $QQBOOK $GGBOOK $FINALBOOK;
    
    [ -f $FINALBOOK ] && 
    echo The combined qq and gg output is written to $FINALBOOK;

} #[[ $ComputeQQ... ]]

echo Cleaning up...
for i in tmp*$BOOKEXT pert*$BOOKEXT w???-a*$BOOKEXT dsi*$BOOKEXT; do
    [ -f $i ] && rm $i;
done

echo The computation is finished.

#Produce a beep
ibeep=5;
while [ $ibeep -gt 0 ]; do
  echo  $'\a'; sleep 0.1; echo  $'\a'; sleep 0.5;
  ibeep=$((ibeep-1));
done