- Timestamp:
- May 25, 2016, 6:00:55 PM (8 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- 2acc2e9
- Parents:
- fb98f40
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
readers/DelphesPythia8.cpp
rfb98f40 r971db5b 144 144 interrupted = true; 145 145 } 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 154 void 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 146 188 147 189 //--------------------------------------------------------------------------- … … 213 255 // Initialize Pythia 214 256 pythia = new Pythia8::Pythia; 257 //Pythia8::Event& event = pythia->event; 258 //Pythia8::ParticleData& pdt = pythia->particleData; 215 259 216 260 if(pythia == NULL) … … 230 274 timesAllowErrors = pythia->mode("Main:timesAllowErrors"); 231 275 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 } 244 292 } 245 293 … … 260 308 while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF, 261 309 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 } 262 315 263 316 if(!pythia->next())
Note:
See TracChangeset
for help on using the changeset viewer.