CMS 3D CMS Logo

PileUpProducer Class Reference

#include <FastSimulation/PileUpProducer/plugins/PileUpProducer.h>

Inheritance diagram for PileUpProducer:

edm::EDProducer edm::ProducerBase edm::ProductRegistryHelper

List of all members.

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 RandomEnginerandom
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
PrimaryVertexGeneratortheVertexGenerator


Detailed Description

Definition at line 22 of file PileUpProducer.h.


Constructor & Destructor Documentation

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 }


Member Function Documentation

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 }

void PileUpProducer::endJob ( void   )  [virtual]

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 }


Member Data Documentation

double PileUpProducer::averageNumber_ [private]

Definition at line 45 of file PileUpProducer.h.

Referenced by PileUpProducer(), and produce().

std::string PileUpProducer::inputFile [private]

Definition at line 48 of file PileUpProducer.h.

Referenced by beginJob(), and PileUpProducer().

unsigned PileUpProducer::myOutputBuffer [private]

Definition at line 61 of file PileUpProducer.h.

Referenced by beginJob(), and save().

std::ofstream PileUpProducer::myOutputFile [private]

Definition at line 60 of file PileUpProducer.h.

Referenced by beginJob(), endJob(), and save().

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]

Definition at line 53 of file PileUpProducer.h.

Referenced by beginJob(), and PileUpProducer().

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]

Definition at line 47 of file PileUpProducer.h.

Referenced by beginJob(), and PileUpProducer().

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().

PrimaryVertexGenerator* PileUpProducer::theVertexGenerator [private]

Definition at line 43 of file PileUpProducer.h.

Referenced by PileUpProducer(), produce(), and ~PileUpProducer().


The documentation for this class was generated from the following files:
Generated on Tue Jun 9 18:29:58 2009 for CMSSW by  doxygen 1.5.4