Fork me on GitHub

Changeset 971db5b in git for readers


Ignore:
Timestamp:
May 25, 2016, 6:00:55 PM (8 years ago)
Author:
Alexandre Mertens <alexandre.mertens@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
2acc2e9
Parents:
fb98f40
Message:

first version working with particles (not jet) and fixed value of Energy

File:
1 edited

Legend:

Unmodified
Added
Removed
  • readers/DelphesPythia8.cpp

    rfb98f40 r971db5b  
    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, 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  // Select particle mass; where relevant according to Breit-Wigner.
     161  double mm = pdt.mSel(id);
     162  double pp = Pythia8::sqrtpos(ee*ee - mm*mm);
     163
     164  // Special case when particle is supposed to be at rest.
     165  if (atRest) {
     166    ee = mm;
     167    pp = 0.;
     168  }
     169
     170  // Angles as input or uniform in solid angle.
     171  double cThe, sThe, phi;
     172  if (thetaIn >= 0.) {
     173    cThe = cos(thetaIn);
     174    sThe = sin(thetaIn);
     175    phi  = phiIn;
     176  } else {
     177    cThe = 2. * rndm.flat() - 1.;
     178    sThe = Pythia8::sqrtpos(1. - cThe * cThe);
     179    phi = 2. * M_PI * rndm.flat();
     180  }
     181
     182  // Store the particle in the event record.
     183  event.append( id, 1, 0, 0, pp * sThe * cos(phi), pp * sThe * sin(phi), pp * cThe, ee, mm);
     184
     185}
     186
     187
    146188
    147189//---------------------------------------------------------------------------
     
    213255    // Initialize Pythia
    214256    pythia = new Pythia8::Pythia;
     257    //Pythia8::Event& event = pythia->event;
     258    //Pythia8::ParticleData& pdt = pythia->particleData;
    215259
    216260    if(pythia == NULL)
     
    230274    timesAllowErrors = pythia->mode("Main:timesAllowErrors");
    231275
    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");
     276    // Check if particle gun
     277    if (!pythia->flag("Main:spareFlag1"))
     278    {
     279      inputFile = fopen(pythia->word("Beams:LHEF").c_str(), "r");
     280      if(inputFile)
     281      {
     282        reader = new DelphesLHEFReader;
     283        reader->SetInputFile(inputFile);
     284
     285        branchEventLHEF = treeWriter->NewBranch("EventLHEF", LHEFEvent::Class());
     286        branchWeightLHEF = treeWriter->NewBranch("WeightLHEF", LHEFWeight::Class());
     287
     288        allParticleOutputArrayLHEF = modularDelphes->ExportArray("allParticlesLHEF");
     289        stableParticleOutputArrayLHEF = modularDelphes->ExportArray("stableParticlesLHEF");
     290        partonOutputArrayLHEF = modularDelphes->ExportArray("partonsLHEF");
     291      }
    244292    }
    245293
     
    260308      while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF,
    261309        stableParticleOutputArrayLHEF, partonOutputArrayLHEF) && !reader->EventReady());
     310
     311      if (pythia->flag("Main:spareFlag1"))
     312      {
     313        fillParticle( pythia->mode("Main:spareMode1"), 30, -1., 0.,pythia->event, pythia->particleData, pythia->rndm, 0);
     314      }
    262315
    263316      if(!pythia->next())
Note: See TracChangeset for help on using the changeset viewer.