Changes in readers/DelphesPythia8.cpp [c0f9b5b:2acc2e9] in git

File:

• readers/DelphesPythia8.cpp (modified) (4 diffs)

readers/DelphesPythia8.cpp

@@ -144 +144 @@
   interrupted = true;
 }
+
+
+// Single-particle gun. The particle must be a colour singlet.
+// Input: flavour, energy, direction (theta, phi).
+// If theta < 0 then random choice over solid angle.
+// Optional final argument to put particle at rest => E = m.
+// from pythia8 example 21
+
+void fillParticle(int id, double ee_max, double thetaIn, double phiIn,
+  Pythia8::Event& event, Pythia8::ParticleData& pdt, Pythia8::Rndm& rndm, bool atRest = false) {
+
+  // Reset event record to allow for new event.
+  event.reset();
+
+  // Angles uniform in solid angle.
+  double cThe, sThe, phi, ee;
+  cThe = 2. * rndm.flat() - 1.;
+  sThe = Pythia8::sqrtpos(1. - cThe * cThe);
+  phi = 2. * M_PI * rndm.flat();
+  ee = pow(10,1+(log10(ee_max)-1)*rndm.flat());
+  double mm = pdt.mSel(id);
+  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
+
+  // Store the particle in the event record.
+  event.append( id, 1, 0, 0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
+
+}
+
+void fillPartons(int type, double ee_max, Pythia8::Event& event, Pythia8::ParticleData& pdt,
+  Pythia8::Rndm& rndm) {
+
+  // Reset event record to allow for new event.
+  event.reset();
+
+  // Angles uniform in solid angle.
+  double cThe, sThe, phi, ee;
+
+  // Information on a q qbar system, to be hadronized.
+
+  cThe = 2. * rndm.flat() - 1.;
+  sThe = Pythia8::sqrtpos(1. - cThe * cThe);
+  phi = 2. * M_PI * rndm.flat();
+  ee = pow(10,1+(log10(ee_max)-1)*rndm.flat());
+  double mm = pdt.m0(type);
+  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
+  if (type == 21)
+  {
+    event.append( 21, 23, 101, 102, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee);
+    event.append( 21, 23, 102, 101, -pp * sThe * cos(phi), -pp * sThe * sin(phi), -pp * cThe, ee);
+  }
+  else
+  {
+    event.append(  type, 23, 101,   0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
+    event.append( -type, 23,   0, 101, -pp * sThe * cos(phi), -pp * sThe * sin(phi), -pp * cThe, ee, mm);
+  }
+}
+
 
 //---------------------------------------------------------------------------
@@ -213 +270 @@
     // Initialize Pythia
     pythia = new Pythia8::Pythia;
+    //Pythia8::Event& event = pythia->event;
+    //Pythia8::ParticleData& pdt = pythia->particleData;
 
     if(pythia == NULL)
@@ -230 +289 @@
     timesAllowErrors = pythia->mode("Main:timesAllowErrors");
 
-    inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
-    if(inputFile)
-    {
-      reader = new DelphesLHEFReader;
-      reader->SetInputFile(inputFile);
-
-      branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
-      branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
-
-      allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
-      stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
-      partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
+    // Check if particle gun
+    if (!pythia->flag("Main:spareFlag1"))
+    {
+      inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
+      if(inputFile)
+      {
+        reader = new DelphesLHEFReader;
+        reader->SetInputFile(inputFile);
+
+        branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
+        branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
+
+        allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
+        stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
+        partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
+      }
     }
 
@@ -260 +323 @@
       while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF,
         stableParticleOutputArrayLHEF, partonOutputArrayLHEF) && !reader->EventReady());
+
+      if (pythia->flag("Main:spareFlag1"))
+      {
+        if (pythia->mode("Main:spareMode1") == 11 || pythia->mode("Main:spareMode1") == 13 || pythia->mode("Main:spareMode1") == 22) 
+        { 
+          fillParticle( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), -1., 0.,pythia->event, pythia->particleData, pythia->rndm, 0);
+        }
+        else fillPartons( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->event, pythia->particleData, pythia->rndm);
+      }
 
       if(!pythia->next())