CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/SimG4Core/GFlash/TB/TreeProducerCalibSimul.cc

Go to the documentation of this file.
00001 #include "SimG4Core/GFlash/TB/TreeProducerCalibSimul.h"
00002 
00003 using namespace std;
00004 
00005 
00006 // -------------------------------------------------
00007 // contructor
00008 TreeProducerCalibSimul::TreeProducerCalibSimul( const edm::ParameterSet& iConfig )
00009 {
00010   // now do what ever initialization is needed
00011   xtalInBeam                 = iConfig.getUntrackedParameter<int>("xtalInBeam",-1000);
00012   rootfile_                  = iConfig.getUntrackedParameter<std::string>("rootfile","mySimMatrixTree.root");
00013   txtfile_                   = iConfig.getUntrackedParameter<std::string>("txtfile", "mySimMatrixTree.txt");
00014   EBRecHitCollection_        = iConfig.getParameter<std::string>("EBRecHitCollection");
00015   RecHitProducer_            = iConfig.getParameter<std::string>("RecHitProducer");
00016   hodoRecInfoCollection_     = iConfig.getParameter<std::string>("hodoRecInfoCollection");
00017   hodoRecInfoProducer_       = iConfig.getParameter<std::string>("hodoRecInfoProducer");
00018   tdcRecInfoCollection_      = iConfig.getParameter<std::string>("tdcRecInfoCollection");
00019   tdcRecInfoProducer_        = iConfig.getParameter<std::string>("tdcRecInfoProducer");
00020   eventHeaderCollection_     = iConfig.getParameter<std::string>("eventHeaderCollection");
00021   eventHeaderProducer_       = iConfig.getParameter<std::string>("eventHeaderProducer");
00022 
00023   // summary
00024   cout << endl;
00025   cout <<"Constructor" << endl;
00026   cout << endl;
00027   cout << "TreeProducerCalibSimul"    << endl;
00028   cout << "xtal in beam = " << xtalInBeam << endl;
00029   cout <<"Fetching hitCollection: "   << EBRecHitCollection_.c_str()     << " prod by " << RecHitProducer_.c_str()      <<endl; 
00030   cout <<"Fetching hodoCollection: "  << hodoRecInfoCollection_.c_str()  << " prod by " << hodoRecInfoProducer_.c_str() <<endl;
00031   cout <<"Fetching tdcCollection: "   << tdcRecInfoCollection_.c_str()   << " prod by " << tdcRecInfoProducer_.c_str()  <<endl;
00032   cout <<"Fetching evHeaCollection: " << eventHeaderCollection_.c_str()  << " prod by " << eventHeaderProducer_.c_str() <<endl;
00033   cout << endl;
00034 }
00035 
00036 
00037 // -------------------------------------------------
00038 // destructor
00039 TreeProducerCalibSimul::~TreeProducerCalibSimul()
00040 {
00041   cout << endl;
00042   cout << "Deleting" << endl;
00043   cout << endl;
00044 
00045   delete myTree;
00046 }
00047 
00048 
00049 
00050 // ------------------------------------------------------
00051 // initializations
00052 void TreeProducerCalibSimul::beginJob()
00053 {
00054   cout << endl;
00055   cout << "BeginJob" << endl;
00056   cout << endl;
00057 
00058   // tree
00059   myTree = new TreeMatrixCalib(rootfile_.c_str());
00060 
00061   // counters
00062   tot_events      = 0;
00063   tot_events_ok   = 0;
00064   noHits          = 0;
00065   noHodo          = 0;
00066   noTdc           = 0;
00067   noHeader        = 0;
00068 } 
00069 
00070 
00071 
00072 // -------------------------------------------
00073 // finalizing
00074 void TreeProducerCalibSimul::endJob() 
00075 {
00076   cout << endl;
00077   cout << "EndJob" << endl;
00078   cout << endl;
00079 
00080   ofstream *MyOut = new ofstream(txtfile_.c_str());   
00081   *MyOut << "total events: "                                   << tot_events      << endl;
00082   *MyOut << "events skipped because of no hits: "              << noHits          << endl;
00083   *MyOut << "events skipped because of no hodos: "             << noHodo          << endl;
00084   *MyOut << "events skipped because of no tdc: "               << noTdc           << endl;
00085   *MyOut << "events skipped because of no header: "            << noHeader        << endl;
00086   *MyOut << "total OK events (passing the basic selection): "  << tot_events_ok   << endl;
00087   MyOut->close();
00088   delete MyOut;
00089 }
00090 
00091 
00092 
00093 // -----------------------------------------------
00094 // my analysis
00095 void TreeProducerCalibSimul::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00096 {
00097   using namespace edm;
00098   using namespace cms;
00099 
00100   // counting events
00101   tot_events++;
00102 
00103   if ( tot_events%5000 == 0){ cout << "event " << tot_events << endl;}
00104   
00105 
00106   // ---------------------------------------------------------------------
00107   // taking what I need: hits
00108   Handle< EBRecHitCollection > pEBRecHits ;
00109   const EBRecHitCollection*  EBRecHits = 0 ;
00110   try {
00111     iEvent.getByLabel (RecHitProducer_, EBRecHitCollection_, pEBRecHits) ;
00112     EBRecHits = pEBRecHits.product(); 
00113   } catch ( std::exception& ex ) {
00114     std::cout << "Error! can't get the product " << EBRecHitCollection_.c_str () << std::endl ;
00115     std::cerr << "Error! can't get the product " << EBRecHitCollection_.c_str () << std::endl ;
00116   }
00117 
00118   // taking what I need: hodoscopes
00119   Handle<EcalTBHodoscopeRecInfo> pHodo;
00120   const EcalTBHodoscopeRecInfo* recHodo=0;
00121   try {
00122     iEvent.getByLabel( hodoRecInfoProducer_, hodoRecInfoCollection_, pHodo);
00123     recHodo = pHodo.product();
00124   } catch ( std::exception& ex ) {
00125     std::cout << "Error! can't get the product " << hodoRecInfoCollection_.c_str() << std::endl;
00126     std::cerr << "Error! can't get the product " << hodoRecInfoCollection_.c_str() << std::endl;
00127   }
00128   
00129   // taking what I need: tdc
00130   Handle<EcalTBTDCRecInfo> pTDC;
00131   const EcalTBTDCRecInfo* recTDC=0;
00132   try {
00133     iEvent.getByLabel( tdcRecInfoProducer_, tdcRecInfoCollection_, pTDC);
00134     recTDC = pTDC.product(); 
00135   } catch ( std::exception& ex ) {
00136     std::cout << "Error! can't get the product " << tdcRecInfoCollection_.c_str() << std::endl;
00137     std::cerr << "Error! can't get the product " << tdcRecInfoCollection_.c_str() << std::endl;
00138   }
00139 
00140   // taking what I need: event header
00141   Handle<EcalTBEventHeader> pEventHeader;
00142   const EcalTBEventHeader* evtHeader=0;
00143   try {
00144     iEvent.getByLabel( eventHeaderProducer_ , pEventHeader );
00145     evtHeader = pEventHeader.product(); 
00146   } catch ( std::exception& ex ) {
00147     std::cout << "Error! can't get the event header " <<std::endl;
00148     std::cerr << "Error! can't get the event header " <<std::endl;
00149   }
00150    
00151   // checking everything is there and fine
00152   if ( (!EBRecHits) || (EBRecHits->size() == 0)){ noHits++;  return; }
00153   if (!recTDC)    { noTdc++;     return; }        
00154   if (!recHodo)   { noHodo++;    return; }                  
00155   if (!evtHeader) { noHeader++;  return; }                
00156   tot_events_ok++;
00157 
00158 
00159 
00160   // ---------------------------------------------------------------------
00161   // info on the event
00162   int run   = -999;
00163   int tbm   = -999;
00164   int event = evtHeader->eventNumber();
00165 
00166   // ---------------------------------------------------------------------
00167   // xtal-in-beam  
00168   int nomXtalInBeam  = -999;
00169   int nextXtalInBeam = -999;
00170 
00171   EBDetId xtalInBeamId(1,xtalInBeam, EBDetId::SMCRYSTALMODE); 
00172   if (xtalInBeamId==EBDetId(0)){ return; }
00173   int mySupCry = xtalInBeamId.ic();
00174   int mySupEta = xtalInBeamId.ieta();
00175   int mySupPhi = xtalInBeamId.iphi();
00176 
00177   
00178   // ---------------------------------------------------------------------
00179   // hodoscope information
00180   double x  = recHodo->posX();
00181   double y  = recHodo->posY();
00182   double sx = recHodo->slopeX();
00183   double sy = recHodo->slopeY();
00184   double qx = recHodo->qualX();
00185   double qy = recHodo->qualY();
00186     
00187 
00188   // ---------------------------------------------------------------------
00189   // tdc information
00190   double tdcOffset = recTDC->offset();
00191   
00192 
00193   // ---------------------------------------------------------------------
00194   // Find EBDetId in a 7x7 Matrix
00195   EBDetId Xtals7x7[49];
00196   double energy[49];
00197   int crystal[49];
00198   int allMatrix = 1; 
00199   for (unsigned int icry=0; icry<49; icry++){
00200     unsigned int row    = icry/7;
00201     unsigned int column = icry%7;
00202     try
00203       {
00204         Xtals7x7[icry]=EBDetId(xtalInBeamId.ieta()+column-3, xtalInBeamId.iphi()+row-3, EBDetId::ETAPHIMODE);
00205         
00206         if ( Xtals7x7[icry].ism() == 1){ 
00207           energy[icry] = EBRecHits->find(Xtals7x7[icry])->energy(); 
00208           crystal[icry] = Xtals7x7[icry].ic();
00209         }
00210         if ( Xtals7x7[icry].ism() != 1){ 
00211           energy[icry] = -100.;
00212           crystal[icry] = -100;
00213           allMatrix = 0; 
00214         }
00215       } 
00216     catch (...)
00217       {
00218         // can not construct 7x7 matrix 
00219         energy[icry]  = -100.; 
00220         crystal[icry] = -100;
00221         allMatrix = 0;  
00222       }
00223   }
00224 
00225 
00226   // ---------------------------------------------------------------------
00227   // Looking for the max energy crystal
00228   double maxEne      = -999.;
00229   int maxEneCry      = 9999;
00230   int maxEneInMatrix = -999;
00231   for (int ii=0; ii<49; ii++){ if (energy[ii] > maxEne){ 
00232     maxEne         = energy[ii]; 
00233     maxEneCry      = crystal[ii]; 
00234     maxEneInMatrix = ii;} 
00235   }  
00236   
00237 
00238 
00239   // Position reconstruction - skipped here
00240   double Xcal   = -999.;
00241   double Ycal   = -999.;
00242 
00243   // filling the tree
00244   myTree->fillInfo(run, event, mySupCry, maxEneCry, nomXtalInBeam, nextXtalInBeam, mySupEta, mySupPhi, tbm, x, y, Xcal, Ycal, sx, sy, qx, qy, tdcOffset, allMatrix, energy, crystal);
00245   myTree->store();
00246 }
00247 
00248 #include "FWCore/PluginManager/interface/ModuleDef.h"
00249 #include "FWCore/Framework/interface/MakerMacros.h"
00250 
00251 //define this as a plug-in
00252 
00253 DEFINE_FWK_MODULE(TreeProducerCalibSimul);