Improvements to RICH examples

Author: preghenella

examples/scripts/rich.sh

@@ -1,27 +1,31 @@
 #! /usr/bin/env bash
 
-NJOBS=5         # number of max parallel runs
-NRUNS=5         # number of runs
-NEVENTS=10000   # number of events in a run
+NJOBS=5                # number of max parallel runs
+NRUNS=5                # number of runs
+NEVENTS=10000          # number of events in a run
 
-BFIELD=5.       # magnetic field   [kG]
-RICHRAD=100.    # RICH radius      [cm]
-RICHLEN=200.    # RICH half length [cm]
-RICHETA=1.443   # RICH max pseudorapidity
+LUTTAG="werner"        # LUT tag name
+PY8CFG="pythia8_ccbar" # pythia8 configuration
+BFIELD=5.              # magnetic field   [kG]
+RICHRAD=100.           # RICH radius      [cm]
+RICHLEN=200.           # RICH half length [cm]
+RICHETA=1.443          # RICH max pseudorapidity
 
 ### calculate max eta from geometry
 RICHETA=`awk -v a=$RICHRAD -v b=$RICHLEN 'BEGIN {th=atan2(a,b)*0.5; sth=sin(th); cth=cos(th); print -log(sth/cth)}'`
 echo "maxEta = $RICHETA"
 
 ### copy relevant files in the working directory
 cp $DELPHESO2_ROOT/examples/cards/propagate.2kG.tcl propagate.tcl
+cp $DELPHESO2_ROOT/examples/pythia8/$PY8CFG.cfg pythia8.cfg
 cp $DELPHESO2_ROOT/examples/smearing/rich.C .
-cp $DELPHESO2_ROOT/examples/pythia8/pythia8_inel.cfg pythia8.cfg
 
 ### adjust pythia8 configuration
+echo "" >> pythia8.cfg
+echo "### run time configuration" >> pythia8.cfg
 echo "Main:numberOfEvents $NEVENTS" >> pythia8.cfg
-#echo "Beams:allowVertexSpread on " >> pythia8.cfg
-#echo "Beams:sigmaTime 60." >> pythia8.cfg
+echo "Beams:allowVertexSpread off " >> pythia8.cfg
+echo "Beams:sigmaTime 60." >> pythia8.cfg
 
 ### set magnetic field
 sed -i -e "s/set barrel_Bz .*$/set barrel_Bz ${BFIELD}e\-1/" propagate.tcl
@@ -36,7 +40,7 @@ sed -i -e "s/set barrel_Acceptance .*$/set barrel_Acceptance \{ 0.0 + 1.0 * fabs
 
 ### create LUTs
 BFIELDT=`awk -v a=$BFIELD 'BEGIN {print a*0.1}'`
-$DELPHESO2_ROOT/examples/scripts/create_luts.sh werner $BFIELDT $TOFRAD
+$DELPHESO2_ROOT/examples/scripts/create_luts.sh $LUTTAG $BFIELDT $TOFRAD
 
 ### loop over runs
 rm -f .running.* delphes.*.root
@@ -54,8 +58,9 @@ for I in $(seq 1 $NRUNS); do
 
     ### copy pythia8 configuration and adjust it
     cp pythia8.cfg pythia8.$I.cfg
+    echo "Random:setSeed on" >> pythia8.$I.cfg
     echo "Random:seed = $I" >> pythia8.$I.cfg
-
+    
     ### run Delphes and analysis    
     DelphesPythia8 propagate.tcl pythia8.$I.cfg delphes.$I.root  &> delphes.$I.log &&
 	root -b -q -l "rich.C(\"delphes.$I.root\", \"rich.$I.root\")" &> rich.$I.log &&

examples/smearing/rich.C

@@ -39,7 +39,7 @@ rich(const char *inputFile = "delphes.root",
   smearer.loadTable(2212, "lutCovm.pr.dat");
 
   // logx binning
-  const Int_t nbins = 80;
+  const Int_t nbins = 200;
   double xmin = 1.e-2;
   double xmax = 1.e2;
   double logxmin = std::log10(xmin);
@@ -72,11 +72,11 @@ rich(const char *inputFile = "delphes.root",
     hRecP[i] = new TH1F(Form("hRecP_%s", pname[i]), ";#it{p} (GeV/#it{c})", nbins, xbins);
     hRecPt[i] = new TH1F(Form("hRecPt_%s", pname[i]), ";#it{p} (GeV/#it{c})", nbins, xbins);
     hAngleP_true[i] = new TH2F(Form("hAngleP_true_%s", pname[i]), ";#it{p} (GeV/#it{c});#theta (rad)", nbins, xbins, 250, 0., 0.25);
-    hNsigmaP[i] = new TH2F(Form("hNsigmaP_%s", pname[i]), Form(";#it{p} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 200, -10., 10.);
-    hNsigmaPt[i] = new TH2F(Form("hNsigmaPt_%s", pname[i]), Form(";#it{p_{T}} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 200, -10., 10.);
+    hNsigmaP[i] = new TH2F(Form("hNsigmaP_%s", pname[i]), Form(";#it{p} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 500, -25., 25.);
+    hNsigmaPt[i] = new TH2F(Form("hNsigmaPt_%s", pname[i]), Form(";#it{p_{T}} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 500, -25., 25.);
     for (int j = 0; j < 5; ++j) {
-      hNsigmaP_true[i][j] = new TH2F(Form("hNsigmaP_%s_true_%s", pname[i], pname[j]), Form(";#it{p} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 200, -10., 10.);
-      hNsigmaPt_true[i][j] = new TH2F(Form("hNsigmaPt_%s_true_%s", pname[i], pname[j]), Form(";#it{p_{T}} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 200, -10., 10.);
+      hNsigmaP_true[i][j] = new TH2F(Form("hNsigmaP_%s_true_%s", pname[i], pname[j]), Form(";#it{p} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 500, -25., 25.);
+      hNsigmaPt_true[i][j] = new TH2F(Form("hNsigmaPt_%s_true_%s", pname[i], pname[j]), Form(";#it{p_{T}} (GeV/#it{c});n#sigma_{%s}", plabel[i]), nbins, xbins, 500, -25., 25.);
     }
   }