Changeset a9c67b3a in git for readers/DelphesPythia8.cpp
- Timestamp:
- Jun 4, 2016, 8:10:24 PM (8 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- 4dd5358
- Parents:
- 3c46e17 (diff), 7b45ff5 (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. - File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
readers/DelphesPythia8.cpp
r3c46e17 ra9c67b3a 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_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 174 void 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 146 203 147 204 //--------------------------------------------------------------------------- … … 213 270 // Initialize Pythia 214 271 pythia = new Pythia8::Pythia; 272 //Pythia8::Event& event = pythia->event; 273 //Pythia8::ParticleData& pdt = pythia->particleData; 215 274 216 275 if(pythia == NULL) … … 230 289 timesAllowErrors = pythia->mode("Main:timesAllowErrors"); 231 290 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 } 244 307 } 245 308 … … 260 323 while(reader && reader->ReadBlock(factory, allParticleOutputArrayLHEF, 261 324 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 } 262 334 263 335 if(!pythia->next())
Note:
See TracChangeset
for help on using the changeset viewer.