Fork me on GitHub

Changeset 1808db5 in git for readers


Ignore:
Timestamp:
May 28, 2016, 12:16:32 PM (9 years ago)
Author:
Pavel Demin <pavel-demin@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
e628a1d
Parents:
fb98f40 (diff), a889a47 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge pull request #22 from AlexandreMertens/particleGun

Particle gun

File:
1 edited

Legend:

Unmodified
Added
Removed
  • readers/DelphesPythia8.cpp

    rfb98f40 r1808db5  
    144144  interrupted = true;
    145145}
     146
     147
     148// Single-particle gun. The particle must be a colour singlet.
     149// Input: flavour, energy, direction (theta, phi).
     150// If theta < 0 then random choice over solid angle.
     151// Optional final argument to put particle at rest => E = m.
     152// from pythia8 example 21
     153
     154void fillParticle(int id, double ee_max, double thetaIn, double phiIn,
     155  Pythia8::Event& event, Pythia8::ParticleData& pdt, Pythia8::Rndm& rndm, bool atRest = false) {
     156
     157  // Reset event record to allow for new event.
     158  event.reset();
     159
     160  // Angles uniform in solid angle.
     161  double cThe, sThe, phi, ee;
     162  cThe = 2. * rndm.flat() - 1.;
     163  sThe = Pythia8::sqrtpos(1. - cThe * cThe);
     164  phi = 2. * M_PI * rndm.flat();
     165  ee = pow(10,1+(log10(ee_max)-1)*rndm.flat());
     166  double mm = pdt.mSel(id);
     167  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
     168
     169  // Store the particle in the event record.
     170  event.append( id, 1, 0, 0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
     171
     172}
     173
     174void fillPartons(int type, double ee_max, Pythia8::Event& event, Pythia8::ParticleData& pdt,
     175  Pythia8::Rndm& rndm) {
     176
     177  // Reset event record to allow for new event.
     178  event.reset();
     179
     180  // Angles uniform in solid angle.
     181  double cThe, sThe, phi, ee;
     182
     183  // Information on a q qbar system, to be hadronized.
     184
     185  cThe = 2. * rndm.flat() - 1.;
     186  sThe = Pythia8::sqrtpos(1. - cThe * cThe);
     187  phi = 2. * M_PI * rndm.flat();
     188  ee = pow(10,1+(log10(ee_max)-1)*rndm.flat());
     189  double mm = pdt.m0(type);
     190  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
     191  if (type == 21)
     192  {
     193    event.append( 21, 23, 101, 102, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee);
     194    event.append( 21, 23, 102, 101, -pp * sThe * cos(phi), -pp * sThe * sin(phi), -pp * cThe, ee);
     195  }
     196  else
     197  {
     198    event.append(  type, 23, 101,   0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
     199    event.append( -type, 23,   0, 101, -pp * sThe * cos(phi), -pp * sThe * sin(phi), -pp * cThe, ee, mm);
     200  }
     201}
     202
    146203
    147204//---------------------------------------------------------------------------
     
    213270    // Initialize Pythia
    214271    pythia = new Pythia8::Pythia;
     272    //Pythia8::Event& event = pythia->event;
     273    //Pythia8::ParticleData& pdt = pythia->particleData;
    215274
    216275    if(pythia == NULL)
     
    230289    timesAllowErrors = pythia->mode("Main:timesAllowErrors");
    231290
    232     inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
    233     if(inputFile)
    234     {
    235       reader = new DelphesLHEFReader;
    236       reader->SetInputFile(inputFile);
    237 
    238       branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
    239       branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
    240 
    241       allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
    242       stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
    243       partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
     291    // Check if particle gun
     292    if (!pythia->flag("Main:spareFlag1"))
     293    {
     294      inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
     295      if(inputFile)
     296      {
     297        reader = new DelphesLHEFReader;
     298        reader->SetInputFile(inputFile);
     299
     300        branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
     301        branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
     302
     303        allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
     304        stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
     305        partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
     306      }
    244307    }
    245308
     
    260323      while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF,
    261324        stableParticleOutputArrayLHEF, partonOutputArrayLHEF) && !reader->EventReady());
     325
     326      if (pythia->flag("Main:spareFlag1"))
     327      {
     328        if (pythia->mode("Main:spareMode1") == 11 || pythia->mode("Main:spareMode1") == 13 || pythia->mode("Main:spareMode1") == 22)
     329        {
     330          fillParticle( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), -1., 0.,pythia->event, pythia->particleData, pythia->rndm, 0);
     331        }
     332        else fillPartons( pythia->mode("Main:spareMode1"), pythia->parm("Main:spareParm1"), pythia->event, pythia->particleData, pythia->rndm);
     333      }
    262334
    263335      if(!pythia->next())
Note: See TracChangeset for help on using the changeset viewer.