#include <FastSimulation/PileUpProducer/plugins/PileUpProducer.h>
Public Member Functions | |
virtual void | beginJob (const edm::EventSetup &c) |
virtual void | endJob () |
PileUpProducer (edm::ParameterSet const &p) | |
virtual void | produce (edm::Event &e, const edm::EventSetup &c) |
virtual | ~PileUpProducer () |
Private Member Functions | |
bool | read (std::string inputFile) |
Read former minbias configuration (from previous run). | |
void | save () |
Save current minbias configuration (for later use). | |
Private Attributes | |
double | averageNumber_ |
std::string | inputFile |
unsigned | myOutputBuffer |
std::ofstream | myOutputFile |
const RandomEngine * | random |
std::vector< TBranch * > | theBranches |
std::vector< unsigned > | theCurrentEntry |
std::vector< unsigned > | theCurrentMinBiasEvt |
std::vector< std::string > | theFileNames |
std::vector< TFile * > | theFiles |
std::vector< unsigned > | theNumberOfEntries |
unsigned | theNumberOfFiles |
std::vector< unsigned > | theNumberOfMinBiasEvts |
std::vector< PUEvent * > | thePUEvents |
std::vector< TTree * > | theTrees |
PrimaryVertexGenerator * | theVertexGenerator |
Definition at line 22 of file PileUpProducer.h.
PileUpProducer::PileUpProducer | ( | edm::ParameterSet const & | p | ) | [explicit] |
Definition at line 27 of file PileUpProducer.cc.
References averageNumber_, Exception, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), inputFile, random, theBranches, theCurrentEntry, theCurrentMinBiasEvt, theFileNames, theFiles, theNumberOfEntries, theNumberOfFiles, theNumberOfMinBiasEvts, thePUEvents, theTrees, and theVertexGenerator.
00028 { 00029 00030 // This producer produces a HepMCProduct, with all pileup vertices/particles 00031 produces<edm::HepMCProduct>("PileUpEvents"); 00032 00033 // Initialize the random number generator service 00034 edm::Service<edm::RandomNumberGenerator> rng; 00035 if ( ! rng.isAvailable() ) { 00036 throw cms::Exception("Configuration") 00037 << "PileUpProducer requires the RandomGeneratorService\n" 00038 "which is not present in the configuration file.\n" 00039 "You must add the service in the configuration file\n" 00040 "or remove the module that requires it"; 00041 } 00042 random = new RandomEngine(&(*rng)); 00043 00044 // The pile-up event generation condition 00045 const edm::ParameterSet& pu = p.getParameter<edm::ParameterSet>("PileUpSimulator"); 00046 averageNumber_ = pu.getParameter<double>("averageNumber"); 00047 theFileNames = pu.getUntrackedParameter<std::vector<std::string> >("fileNames"); 00048 inputFile = pu.getUntrackedParameter<std::string>("inputFile"); 00049 theNumberOfFiles = theFileNames.size(); 00050 theFiles.resize(theNumberOfFiles); 00051 theTrees.resize(theNumberOfFiles); 00052 theBranches.resize(theNumberOfFiles); 00053 thePUEvents.resize(theNumberOfFiles); 00054 theCurrentEntry.resize(theNumberOfFiles); 00055 theCurrentMinBiasEvt.resize(theNumberOfFiles); 00056 theNumberOfEntries.resize(theNumberOfFiles); 00057 theNumberOfMinBiasEvts.resize(theNumberOfFiles); 00058 00059 // Initialize the primary vertex generator 00060 const edm::ParameterSet& vtx = p.getParameter<edm::ParameterSet>("VertexGenerator"); 00061 std::string vtxType = vtx.getParameter<std::string>("type"); 00062 if ( vtxType == "Gaussian" ) 00063 theVertexGenerator = new GaussianPrimaryVertexGenerator(vtx,random); 00064 else if ( vtxType == "Flat" ) 00065 theVertexGenerator = new FlatPrimaryVertexGenerator(vtx,random); 00066 else if ( vtxType == "BetaFunc" ) 00067 theVertexGenerator = new BetaFuncPrimaryVertexGenerator(vtx,random); 00068 else 00069 theVertexGenerator = new NoPrimaryVertexGenerator(); 00070 00071 }
PileUpProducer::~PileUpProducer | ( | ) | [virtual] |
Definition at line 73 of file PileUpProducer.cc.
References theVertexGenerator.
00073 { 00074 00075 delete theVertexGenerator; 00076 00077 }
void PileUpProducer::beginJob | ( | const edm::EventSetup & | c | ) | [virtual] |
Reimplemented from edm::EDProducer.
Definition at line 79 of file PileUpProducer.cc.
References GenMuonPlsPt100GeV_cfg::cout, lat::endl(), Exception, file, RandomEngine::flatShoot(), edm::FileInPath::fullPath(), iggi_31X_cfg::input, inputFile, myOutputBuffer, myOutputFile, random, read(), theBranches, theCurrentEntry, theCurrentMinBiasEvt, theFileNames, theFiles, theNumberOfEntries, theNumberOfFiles, theNumberOfMinBiasEvts, thePUEvents, and theTrees.
00080 { 00081 00082 std::cout << " PileUpProducer initializing " << std::endl; 00083 gROOT->cd(); 00084 00085 std::string fullPath; 00086 00087 // Read the information from a previous run (to keep reproducibility) 00088 bool input = this->read(inputFile); 00089 if ( input ) 00090 std::cout << "***WARNING*** You are reading pile-up information from the file " 00091 << inputFile << " created in an earlier run." 00092 << std::endl; 00093 00094 // Open the file for saving the information of the current run 00095 myOutputFile.open ("PileUpOutputFile.txt"); 00096 myOutputBuffer = 0; 00097 00098 // Open the root files 00099 std::cout << "Opening minimum-bias event files ... " << std::endl; 00100 for ( unsigned file=0; file<theNumberOfFiles; ++file ) { 00101 00102 edm::FileInPath myDataFile("FastSimulation/PileUpProducer/data/"+theFileNames[file]); 00103 fullPath = myDataFile.fullPath(); 00104 // theFiles[file] = TFile::Open(theFileNames[file].c_str()); 00105 theFiles[file] = TFile::Open(fullPath.c_str()); 00106 if ( !theFiles[file] ) throw cms::Exception("FastSimulation/PileUpProducer") 00107 << "File " << theFileNames[file] << " " << fullPath << " not found "; 00108 // 00109 theTrees[file] = (TTree*) theFiles[file]->Get("MinBiasEvents"); 00110 if ( !theTrees[file] ) throw cms::Exception("FastSimulation/PileUpProducer") 00111 << "Tree with name MinBiasEvents not found in " << theFileNames[file]; 00112 // 00113 theBranches[file] = theTrees[file]->GetBranch("puEvent"); 00114 if ( !theBranches[file] ) throw cms::Exception("FastSimulation/PileUpProducer") 00115 << "Branch with name puEvent not found in " << theFileNames[file]; 00116 // 00117 thePUEvents[file] = new PUEvent(); 00118 theBranches[file]->SetAddress(&thePUEvents[file]); 00119 // 00120 theNumberOfEntries[file] = theTrees[file]->GetEntries(); 00121 // Add some randomness (if there was no input file) 00122 if ( !input ) 00123 theCurrentEntry[file] 00124 = (unsigned) (theNumberOfEntries[file] * random->flatShoot()); 00125 00126 theTrees[file]->GetEntry(theCurrentEntry[file]); 00127 unsigned NMinBias = thePUEvents[file]->nMinBias(); 00128 theNumberOfMinBiasEvts[file] = NMinBias; 00129 // Add some randomness (if there was no input file) 00130 if ( !input ) 00131 theCurrentMinBiasEvt[file] = 00132 (unsigned) (theNumberOfMinBiasEvts[file] * random->flatShoot()); 00133 00134 /* 00135 std::cout << "File " << theFileNames[file] 00136 << " is opened with " << theNumberOfEntries[file] 00137 << " entries and will be read from Entry/Event " 00138 << theCurrentEntry[file] << "/" << theCurrentMinBiasEvt[file] 00139 << std::endl; 00140 */ 00141 } 00142 00143 // Return Loot in the same state as it was when entering. 00144 gROOT->cd(); 00145 00146 }
Reimplemented from edm::EDProducer.
Definition at line 148 of file PileUpProducer.cc.
References GenMuonPlsPt100GeV_cfg::cout, lat::endl(), file, myOutputFile, and theFiles.
00149 { 00150 std::cout << " PileUpProducer terminating " << std::endl; 00151 // Close all local files 00152 // Among other things, this allows the TROOT destructor to end up 00153 // without crashing, while trying to close these files from outside 00154 std::cout << "Closing minimum-bias event files... " << std::endl; 00155 for ( unsigned file=0; file<theFiles.size(); ++file ) { 00156 00157 // std::cout << "Closing " << theFileNames[file] << std::endl; 00158 theFiles[file]->Close(); 00159 00160 } 00161 00162 // Close the output file 00163 myOutputFile.close(); 00164 00165 // And return Loot in the same state as it was when entering. 00166 gROOT->cd(); 00167 00168 }
void PileUpProducer::produce | ( | edm::Event & | e, | |
const edm::EventSetup & | c | |||
) | [virtual] |
Implements edm::EDProducer.
Definition at line 170 of file PileUpProducer.cc.
References averageNumber_, PrimaryVertexGenerator::boost(), funct::cos(), relval_parameters_module::energy, file, PUEvent::PUMinBiasEvt::first, RandomEngine::flatShoot(), PrimaryVertexGenerator::generate(), PUEvent::PUParticle::id, int, PUEvent::PUParticle::mass, RandomEngine::poissonShoot(), edm::Event::put(), PUEvent::PUParticle::px, PUEvent::PUParticle::py, PUEvent::PUParticle::pz, random, save(), funct::sin(), PUEvent::PUMinBiasEvt::size, funct::sqrt(), theCurrentEntry, theCurrentMinBiasEvt, theNumberOfEntries, theNumberOfFiles, theNumberOfMinBiasEvts, thePUEvents, theTrees, and theVertexGenerator.
00171 { 00172 00173 // Create the GenEvent and the HepMCProduct 00174 std::auto_ptr<edm::HepMCProduct> pu_product(new edm::HepMCProduct()); 00175 HepMC::GenEvent* evt = new HepMC::GenEvent(); 00176 00177 // How many pile-up events? 00178 int PUevts = (int) random->poissonShoot(averageNumber_); 00179 00180 // Get N events from random files 00181 for ( int ievt=0; ievt<PUevts; ++ievt ) { 00182 00183 // Draw a file in a ramdom manner 00184 unsigned file = (unsigned) (theNumberOfFiles * random->flatShoot()); 00185 /* 00186 if ( debug ) 00187 std::cout << "The file chosen for event " << ievt 00188 << " is the file number " << file << std::endl; 00189 */ 00190 00191 // Smear the primary vertex and express it in mm (stupid GenEvent convention...) 00192 theVertexGenerator->generate(); 00193 HepMC::FourVector smearedVertex = 00194 HepMC::FourVector(theVertexGenerator->X()*10., 00195 theVertexGenerator->Y()*10., 00196 theVertexGenerator->Z()*10., 00197 0.); 00198 HepMC::GenVertex* aVertex = new HepMC::GenVertex(smearedVertex); 00199 evt->add_vertex(aVertex); 00200 00201 // Some rotation around the z axis, for more randomness 00202 double theAngle = random->flatShoot() * 2. * 3.14159265358979323; 00203 double cAngle = std::cos(theAngle); 00204 double sAngle = std::sin(theAngle); 00205 00206 /* 00207 if ( debug ) 00208 std::cout << "File chosen : " << file 00209 << " Current entry in this file " << theCurrentEntry[file] 00210 << " Current minbias in this chunk= " << theCurrentMinBiasEvt[file] 00211 << " Total number of minbias in this chunk = " << theNumberOfMinBiasEvts[file] << std::endl; 00212 */ 00213 00214 // theFiles[file]->cd(); 00215 // gDirectory->ls(); 00216 // Check we are not either at the end of a minbias bunch 00217 // or at the end of a file 00218 if ( theCurrentMinBiasEvt[file] == theNumberOfMinBiasEvts[file] ) { 00219 // if ( debug ) std::cout << "End of MinBias bunch ! "; 00220 ++theCurrentEntry[file]; 00221 // if ( debug) std::cout << "Read the next entry " << theCurrentEntry[file] << std::endl; 00222 theCurrentMinBiasEvt[file] = 0; 00223 if ( theCurrentEntry[file] == theNumberOfEntries[file] ) { 00224 theCurrentEntry[file] = 0; 00225 // if ( debug ) std::cout << "End of file - Rewind! " << std::endl; 00226 } 00227 // if ( debug ) std::cout << "The PUEvent is reset ... "; 00228 thePUEvents[file]->reset(); 00229 unsigned myEntry = theCurrentEntry[file]; 00230 /* 00231 if ( debug ) std::cout << "The new entry " << myEntry 00232 << " is read ... in TTree " << theTrees[file] << " "; 00233 */ 00234 theTrees[file]->GetEntry(myEntry); 00235 /* 00236 if ( debug ) 00237 std::cout << "The number of interactions in the new entry is ... "; 00238 */ 00239 theNumberOfMinBiasEvts[file] = thePUEvents[file]->nMinBias(); 00240 // if ( debug ) std::cout << theNumberOfMinBiasEvts[file] << std::endl; 00241 } 00242 00243 // Read a minbias event chunk 00244 const PUEvent::PUMinBiasEvt& aMinBiasEvt 00245 = thePUEvents[file]->thePUMinBiasEvts()[theCurrentMinBiasEvt[file]]; 00246 00247 // Find corresponding particles 00248 unsigned firstTrack = aMinBiasEvt.first; 00249 unsigned trackSize = firstTrack + aMinBiasEvt.size; 00250 /* 00251 if ( debug ) std::cout << "First and last+1 tracks are " 00252 << firstTrack << " " << trackSize << std::endl; 00253 */ 00254 00255 // Loop on particles 00256 for ( unsigned iTrack=firstTrack; iTrack<trackSize; ++iTrack ) { 00257 00258 const PUEvent::PUParticle& aParticle 00259 = thePUEvents[file]->thePUParticles()[iTrack]; 00260 /* 00261 if ( debug) 00262 std::cout << "Track " << iTrack 00263 << " id/px/py/pz/mass " 00264 << aParticle.id << " " 00265 << aParticle.px << " " 00266 << aParticle.py << " " 00267 << aParticle.pz << " " 00268 << aParticle.mass << " " << std::endl; 00269 */ 00270 00271 // Create a FourVector, with rotation 00272 double energy = std::sqrt( aParticle.px*aParticle.px 00273 + aParticle.py*aParticle.py 00274 + aParticle.pz*aParticle.pz 00275 + aParticle.mass*aParticle.mass ); 00276 00277 HepMC::FourVector myPart(cAngle * aParticle.px + sAngle * aParticle.py, 00278 -sAngle * aParticle.px + cAngle * aParticle.py, 00279 aParticle.pz, energy); 00280 00281 // Add a GenParticle 00282 HepMC::GenParticle* aGenParticle = new HepMC::GenParticle(myPart,aParticle.id); 00283 aVertex->add_particle_out(aGenParticle); 00284 00285 } 00286 // End of particle loop 00287 00288 // Increment for next time 00289 ++theCurrentMinBiasEvt[file]; 00290 00291 } 00292 // End of pile-up event loop 00293 00294 // evt->print(); 00295 00296 // Fill the HepMCProduct from the GenEvent 00297 if ( evt ) { 00298 pu_product->addHepMCData( evt ); 00299 // Boost in case of beam crossing angle 00300 TMatrixD* boost = theVertexGenerator->boost(); 00301 if ( boost ) pu_product->boostToLab(boost,"momentum"); 00302 } 00303 00304 // Put the HepMCProduct onto the event 00305 iEvent.put(pu_product,"PileUpEvents"); 00306 // delete evt; 00307 00308 // Save the current location in each pile-up event files 00309 this->save(); 00310 00311 }
bool PileUpProducer::read | ( | std::string | inputFile | ) | [private] |
Read former minbias configuration (from previous run).
Definition at line 336 of file PileUpProducer.cc.
References size, theCurrentEntry, and theCurrentMinBiasEvt.
Referenced by beginJob().
00336 { 00337 00338 ifstream myInputFile; 00339 struct stat results; 00340 unsigned size1 = theCurrentEntry.size()*sizeof(unsigned); 00341 unsigned size2 = theCurrentMinBiasEvt.size()*sizeof(unsigned); 00342 unsigned size = 0; 00343 00344 00345 // Open the file (if any) 00346 myInputFile.open (inputFile.c_str()); 00347 if ( myInputFile.is_open() ) { 00348 00349 // Get the size of the file 00350 if ( stat(inputFile.c_str(), &results) == 0 ) size = results.st_size; 00351 else return false; // Something is wrong with that file ! 00352 00353 // Position the pointer just before the last record 00354 myInputFile.seekg(size-size1-size2); 00355 00356 // Read the information 00357 myInputFile.read((char*)(&theCurrentEntry.front()),size1); 00358 myInputFile.read((char*)&theCurrentMinBiasEvt.front(),size2); 00359 myInputFile.close(); 00360 00361 return true; 00362 00363 } 00364 00365 return false; 00366 00367 }
void PileUpProducer::save | ( | ) | [private] |
Save current minbias configuration (for later use).
Definition at line 314 of file PileUpProducer.cc.
References myOutputBuffer, myOutputFile, theCurrentEntry, and theCurrentMinBiasEvt.
Referenced by produce().
00314 { 00315 00316 // Size of buffer 00317 ++myOutputBuffer; 00318 00319 // Periodically close the current file and open a new one 00320 if ( myOutputBuffer/1000*1000 == myOutputBuffer ) { 00321 myOutputFile.close(); 00322 myOutputFile.open ("PileUpOutputFile.txt"); 00323 // myOutputFile.seekp(0); // No need to rewind in that case 00324 } 00325 00326 // Save the current position to file 00327 myOutputFile.write((const char*)(&theCurrentEntry.front()), 00328 theCurrentEntry.size()*sizeof(unsigned)); 00329 myOutputFile.write((const char*)&theCurrentMinBiasEvt.front(), 00330 theCurrentMinBiasEvt.size()*sizeof(unsigned)); 00331 myOutputFile.flush(); 00332 00333 }
double PileUpProducer::averageNumber_ [private] |
std::string PileUpProducer::inputFile [private] |
unsigned PileUpProducer::myOutputBuffer [private] |
std::ofstream PileUpProducer::myOutputFile [private] |
const RandomEngine* PileUpProducer::random [private] |
Definition at line 46 of file PileUpProducer.h.
Referenced by beginJob(), PileUpProducer(), and produce().
std::vector<TBranch*> PileUpProducer::theBranches [private] |
std::vector<unsigned> PileUpProducer::theCurrentEntry [private] |
Definition at line 55 of file PileUpProducer.h.
Referenced by beginJob(), PileUpProducer(), produce(), read(), and save().
std::vector<unsigned> PileUpProducer::theCurrentMinBiasEvt [private] |
Definition at line 56 of file PileUpProducer.h.
Referenced by beginJob(), PileUpProducer(), produce(), read(), and save().
std::vector<std::string> PileUpProducer::theFileNames [private] |
std::vector<TFile*> PileUpProducer::theFiles [private] |
Definition at line 51 of file PileUpProducer.h.
Referenced by beginJob(), endJob(), and PileUpProducer().
std::vector<unsigned> PileUpProducer::theNumberOfEntries [private] |
Definition at line 57 of file PileUpProducer.h.
Referenced by beginJob(), PileUpProducer(), and produce().
unsigned PileUpProducer::theNumberOfFiles [private] |
Definition at line 49 of file PileUpProducer.h.
Referenced by beginJob(), PileUpProducer(), and produce().
std::vector<unsigned> PileUpProducer::theNumberOfMinBiasEvts [private] |
Definition at line 58 of file PileUpProducer.h.
Referenced by beginJob(), PileUpProducer(), and produce().
std::vector<PUEvent*> PileUpProducer::thePUEvents [private] |
Definition at line 54 of file PileUpProducer.h.
Referenced by beginJob(), PileUpProducer(), and produce().
std::vector<TTree*> PileUpProducer::theTrees [private] |
Definition at line 52 of file PileUpProducer.h.
Referenced by beginJob(), PileUpProducer(), and produce().
Definition at line 43 of file PileUpProducer.h.
Referenced by PileUpProducer(), produce(), and ~PileUpProducer().