Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    GenPurposeSkimmerData
00004 // Class:      GenPurposeSkimmerData
00005 // 
00040 //
00041 // Original Author:  Nikolaos Rompotis
00042 //         Created:  Thu Oct 16 17:11:55 CEST 2008
00043 // $Id: GenPurposeSkimmerData.cc,v 1.4 2011/10/13 14:39:54 eulisse Exp $
00044 //
00045 //
00046 
00047 #include "ElectroWeakAnalysis/WENu/interface/GenPurposeSkimmerData.h"
00048 
00049 //
00050 //
00051 //
00052 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00053 #include "DataFormats/TrackReco/interface/Track.h"
00054 #include "DataFormats/PatCandidates/interface/Muon.h"
00055 #include "FWCore/Common/interface/TriggerNames.h"
00056 //
00057 //
00058 GenPurposeSkimmerData::GenPurposeSkimmerData(const edm::ParameterSet& ps)
00059 
00060 {
00061 //
00062 //   I N P U T      P A R A M E T E R S
00063 //
00064   // output file name
00065   outputFile_ = ps.getUntrackedParameter<std::string>("outputfile");
00066   //
00067   // Electron Collection
00068   ElectronCollection_=ps.getUntrackedParameter<edm::InputTag>("ElectronCollection");
00069   //
00070   // MC:
00071   //MCCollection_ = ps.getUntrackedParameter<edm::InputTag>("MCCollection");
00072   //MCMatch_Deta_ = ps.getUntrackedParameter<double>("MCMatch_Deta",0.1);
00073   //MCMatch_Dphi_ = ps.getUntrackedParameter<double>("MCMatch_Dphi",0.35);
00074   //
00075   // MET Collections:
00076   MetCollectionTag_ = ps.getUntrackedParameter<edm::InputTag>("MetCollectionTag");
00077   mcMetCollectionTag_ = ps.getUntrackedParameter<edm::InputTag>("mcMetCollectionTag");
00078   t1MetCollectionTag_ = ps.getUntrackedParameter<edm::InputTag>("t1MetCollectionTag");
00079   pfMetCollectionTag_ = ps.getUntrackedParameter<edm::InputTag>("pfMetCollectionTag");
00080   tcMetCollectionTag_ = ps.getUntrackedParameter<edm::InputTag>("tcMetCollectionTag");
00081   //  genMetCollectionTag_ = ps.getUntrackedParameter<edm::InputTag>("genMetCollectionTag");
00082   //
00083   // HLT parameters:
00084   // allow info for 2 paths and 2 filters
00085   // ---------------------------------------------------------------------------
00086   HLTCollectionE29_= ps.getUntrackedParameter<edm::InputTag>("HLTCollectionE29");
00087   HLTCollectionE31_= ps.getUntrackedParameter<edm::InputTag>("HLTCollectionE31");
00088   HLTTriggerResultsE29_ = ps.getUntrackedParameter<edm::InputTag>("HLTTriggerResultsE29");
00089   HLTTriggerResultsE31_ = ps.getUntrackedParameter<edm::InputTag>("HLTTriggerResultsE31");
00090   //HLTPath_ = ps.getUntrackedParameter<std::string>("HLTPath","HLT_Ele15_LW_L1R");
00091   //HLTFilterType_ =ps.getUntrackedParameter<edm::InputTag>("HLTFilterType");
00092   //
00093   // matching HLT objects to electrons
00094   ProbeHLTObjMaxDR= ps.getUntrackedParameter<double>("ProbeHLTObjMaxDR",0.2);
00095   //
00096   // ----------------------------------------------------------------------------
00097   //
00098   // detector geometry
00099   //
00100   BarrelMaxEta = ps.getUntrackedParameter<double>("BarrelMaxEta");
00101   EndcapMinEta = ps.getUntrackedParameter<double>("EndcapMinEta");
00102   EndcapMaxEta = ps.getUntrackedParameter<double>("EndcapMaxEta");
00103   // 
00104   ctfTracksTag_ = ps.getUntrackedParameter<edm::InputTag>("ctfTracksTag");
00105   corHybridsc_  = ps.getUntrackedParameter<edm::InputTag>("corHybridsc");
00106   multi5x5sc_   = ps.getUntrackedParameter<edm::InputTag>("multi5x5sc");
00107 
00108 }
00109 
00110 
00111 GenPurposeSkimmerData::~GenPurposeSkimmerData()
00112 {
00113  
00114    // do anything here that needs to be done at desctruction time
00115    // (e.g. close files, deallocate resources etc.)
00116 
00117 }
00118 
00119 
00120 //
00121 // member functions
00122 //
00123 
00124 // ------------ method called to for each event  ------------
00125 void
00126 GenPurposeSkimmerData::analyze(const edm::Event& evt, const edm::EventSetup& es)
00127 {
00128   // MC Collection ------------------------------------------------
00129   
00130   //  edm::Handle<reco::GenParticleCollection> pGenPart;
00131   //  evt.getByLabel(MCCollection_, pGenPart);
00132   //  if ( not  pGenPart.isValid() ) {
00133   //    std::cout <<"Error! Can't get "<<MCCollection_.label() << std::endl;
00134   //    return;
00135   //  }
00136   
00137   //  const reco::GenParticleCollection *McCand = pGenPart.product();
00138   
00139   // GsF Electron Collection ---------------------------------------
00140   edm::Handle<pat::ElectronCollection> pElectrons;
00141 
00142   try{
00143     evt.getByLabel(ElectronCollection_, pElectrons);
00144   }
00145   catch (cms::Exception)
00146     {
00147       edm::LogError("")<< "Error! Can't get ElectronCollection by label. ";
00148     }
00149   // ***********************************************************************
00150   // check which trigger has accepted the event ****************************
00151   // ***********************************************************************
00152   //
00153   // path allocation: first 10 paths belong to the low lum menu, the rest
00154   // in the high lum one
00155   //
00156   // Low Luminosity Menu (8e29)
00157   //
00158   /*
00159   edm::Handle<edm::TriggerResults> HLTResultsE29;
00160   evt.getByLabel(HLTTriggerResultsE29_, HLTResultsE29);
00161   if (not HLTResultsE29.isValid()) {
00162     std::cout << "HLT Results with label: " << HLTTriggerResultsE29_ 
00163               << " not found" << std::endl;
00164     return;
00165   }
00166   //
00167   edm::Handle<trigger::TriggerEvent> pHLTe29;
00168   evt.getByLabel(HLTCollectionE29_, pHLTe29);
00169   if (not pHLTe29.isValid()) {
00170     std::cout << "HLT Results with label: " << HLTCollectionE29_
00171               << " not found" << std::endl;
00172     return;
00173   }
00174   //
00175   int sum = 0;
00176   //
00177   for (int iT=0; iT<10; ++iT) {
00178     event_HLTPath[iT] = 0;
00179     numberOfHLTFilterObjects[iT] =0;
00180     //
00181     const edm::TriggerNames & triggerNames = evt.triggerNames(*HLTResultsE29);
00182     unsigned int trigger_size = HLTResultsE29->size();
00183     unsigned int trigger_position = triggerNames.triggerIndex(HLTPath_[iT]);
00184     if (trigger_position < trigger_size ) 
00185       event_HLTPath[iT] = (int) HLTResultsE29->accept(trigger_position);
00186     //
00187     numberOfHLTFilterObjects[iT] = 0;
00188     // check explicitly that the filter is there
00189     const int nF(pHLTe29->sizeFilters());
00190     const int filterInd = pHLTe29->filterIndex(HLTFilterType_[iT]);
00191     if (nF != filterInd) {
00192       const trigger::Vids& VIDS (pHLTe29->filterIds(filterInd));
00193       const trigger::Keys& KEYS(pHLTe29->filterKeys(filterInd));
00194       const int nI(VIDS.size());
00195       const int nK(KEYS.size());
00196       numberOfHLTFilterObjects[iT] = (nI>nK)? nI:nK;
00197     }
00198     //if (iT==2) // HLT_Ele15_LW_L1R only this trigger is required
00199       sum += numberOfHLTFilterObjects[iT];
00200   }
00201   //
00202   // High Luminosity Menu (1e31) DISABLED - only low lumi level
00203   //
00204   edm::Handle<edm::TriggerResults> HLTResultsE31;
00205   evt.getByLabel(HLTTriggerResultsE31_, HLTResultsE31);
00206   if (not HLTResultsE31.isValid()) {
00207       std::cout << "HLT Results with label: " << HLTTriggerResultsE31_ 
00208             << " not found" << std::endl;
00209     return;
00210   }
00212   edm::Handle<trigger::TriggerEvent> pHLTe31;
00213   evt.getByLabel(HLTCollectionE31_, pHLTe31);
00214   if (not pHLTe31.isValid()) {
00215     std::cout << "HLT Results with label: " << HLTCollectionE31_
00216               << " not found" << std::endl;
00217     return;
00218   }
00220   for (int iT=10; iT<25; ++iT) {
00221     event_HLTPath[iT] = 0;
00222     numberOfHLTFilterObjects[iT] =0;
00223     //
00224     const edm::TriggerNames & triggerNames = evt.triggerNames(*HLTResultsE31);
00225     unsigned int trigger_size = HLTResultsE31->size();
00226     unsigned int trigger_position = triggerNames.triggerIndex(HLTPath_[iT]);
00227     if (trigger_position < trigger_size ) 
00228       event_HLTPath[iT] = (int) HLTResultsE31->accept(trigger_position);
00229     //
00230     numberOfHLTFilterObjects[iT] = 0;
00231     // check explicitly that the filter is there
00232     const int nF(pHLTe31->sizeFilters());
00233     const int filterInd = pHLTe31->filterIndex(HLTFilterType_[iT]);
00234     if (nF != filterInd) {
00235       const trigger::Vids& VIDS (pHLTe31->filterIds(filterInd));
00236       const trigger::Keys& KEYS(pHLTe31->filterKeys(filterInd));
00237       const int nI(VIDS.size());
00238       const int nK(KEYS.size());
00239       numberOfHLTFilterObjects[iT] = (nI>nK)? nI:nK;
00240     }
00241     // not needed
00242     sum += numberOfHLTFilterObjects[iT];
00243   }
00244   if (sum == 0) { 
00245     //std::cout << "No trigger found in this event..." << std::endl;
00246     return;
00247   }
00248   */
00249   //std::cout << "HLT objects: #" << sum << std::endl;
00250   // print out the triggers that exist in this event
00251     // comment this out if you want to see the names of the existing triggers
00252   edm::Handle<trigger::TriggerEvent> pHLTe29;
00253   evt.getByLabel(HLTCollectionE29_, pHLTe29);
00254   if (not pHLTe29.isValid()){
00255     std::cout << "Error!!! HLT is missing!" << std::endl;
00256     return;
00257   } /*
00258   else {
00259     // check explicitly that the filter is there
00260     const int nF(pHLTe29->sizeFilters());
00261     for (int filterInd=0; filterInd< nF; ++filterInd) {
00262       const trigger::Vids& VIDS (pHLTe29->filterIds(filterInd));
00263       const trigger::Keys& KEYS(pHLTe29->filterKeys(filterInd));
00264       const int nI(VIDS.size());
00265       const int nK(KEYS.size());
00266       int   nObjects = (nI>nK)? nI:nK;     
00267       const edm::InputTag filterTag = pHLTe29->filterTag(filterInd);
00268       std::cout << "Found filter with name " << filterTag
00269                 << " and #objects: #" << nObjects << std::endl;
00270     }
00271   }
00272     */
00273   // *********************************************************************
00274   // MET Collections:
00275   //
00276   edm::Handle<reco::CaloMETCollection> caloMET;
00277   evt.getByLabel(MetCollectionTag_, caloMET);  
00278   //
00279   edm::Handle<pat::METCollection> t1MET;
00280   evt.getByLabel(t1MetCollectionTag_, t1MET);
00281   //
00282   edm::Handle<pat::METCollection> mcMET;
00283   evt.getByLabel(mcMetCollectionTag_, mcMET);
00284   //
00285   edm::Handle<reco::METCollection> tcMET;
00286   evt.getByLabel(tcMetCollectionTag_, tcMET);
00287   //
00288   edm::Handle<reco::PFMETCollection> pfMET;
00289   evt.getByLabel(pfMetCollectionTag_, pfMET);
00290   //
00291   //  edm::Handle<reco::GenMETCollection> genMET;
00292   //  evt.getByLabel(genMetCollectionTag_, genMET);
00293   //
00294   // initialize the MET variables ........................................
00295   event_MET     = -99.;   event_MET_phi = -99.;    event_MET_sig = -99.;
00296   event_mcMET     = -99.;   event_mcMET_phi = -99.;    event_mcMET_sig = -99.;
00297   event_tcMET   = -99.;   event_tcMET_phi = -99.;  event_tcMET_sig = -99.;
00298   event_pfMET   = -99.;   event_pfMET_phi = -99.;  event_pfMET_sig = -99.;
00299   event_t1MET   = -99.;   event_t1MET_phi = -99.;  event_t1MET_sig = -99.;
00300   //
00301   //event_genMET  = -99.;   event_genMET_phi= -99.;  event_genMET_sig = -99.;
00302   //
00303   // get the values, if they are available
00304   if ( caloMET.isValid() ) {
00305     const reco::CaloMETRef MET(caloMET, 0);
00306     event_MET = MET->et();  event_MET_phi = MET->phi();
00307     event_MET_sig = MET->mEtSig();
00308   }
00309   else {
00310     std::cout << "caloMET not valid: input Tag: " << MetCollectionTag_
00311               << std::endl;
00312   }
00313   if ( tcMET.isValid() ) {
00314     const reco::METRef MET(tcMET, 0);
00315     event_tcMET = MET->et();  event_tcMET_phi = MET->phi();
00316     event_tcMET_sig = MET->mEtSig();
00317   }
00318   if ( pfMET.isValid() ) {
00319     const reco::PFMETRef MET(pfMET, 0);
00320     event_pfMET = MET->et();  event_pfMET_phi = MET->phi();
00321     event_pfMET_sig = MET->mEtSig();
00322   }
00323   if ( t1MET.isValid() ) {
00324     const pat::METRef MET(t1MET, 0);
00325     event_t1MET = MET->et();  event_t1MET_phi = MET->phi();
00326     event_t1MET_sig = MET->mEtSig();
00327   }
00328   if ( mcMET.isValid() ) {
00329     const pat::METRef MET(mcMET, 0);
00330     event_mcMET = MET->et();  event_mcMET_phi = MET->phi();
00331     event_mcMET_sig = MET->mEtSig();
00332   }
00333 
00334   //  if ( genMET.isValid() ) {
00335   //    const reco::GenMETRef MET(genMET, 0);
00336   //    event_genMET = MET->et();  event_genMET_phi = MET->phi();
00337   //    event_genMET_sig = MET->mEtSig();
00338   //  }
00339 
00340   //  std::cout << "t1MET: " << event_t1MET  << " twikiT1MET: " 
00341   //        << event_twikiT1MET  << ", calo="<<event_MET  << std::endl;
00342   //
00343   // some supercluster collections ...........................................
00344   // correcyedHybridSuperClusters
00345   //InputTag corHybridsc("correctedHybridSuperClusters","",InputTagEnding_);
00346   edm::Handle<reco::SuperClusterCollection> SC1;
00347   evt.getByLabel(corHybridsc_,SC1);
00348   const reco::SuperClusterCollection *sc1 = SC1.product();
00349   // multi5x5SuperClustersWithPreshower
00350   //edm::InputTag multi5x5sc("multi5x5SuperClustersWithPreshower",
00351   //                       "", InputTagEnding_);
00352   edm::Handle<reco::SuperClusterCollection> SC2;
00353   evt.getByLabel(multi5x5sc_,SC2);
00354   const reco::SuperClusterCollection *sc2 = SC2.product();
00355   //
00356   const int n1 =  sc1->size();
00357   const int n2 =  sc2->size();
00358   //std::cout << "SC found: hybrid: " << n1 << ", multi5x5: " 
00359   //        << n2 << std::endl;
00360   // keep details of the 5 highest ET superclusters
00361   for (int i=0; i<5; ++i) {
00362     sc_hybrid_et[i] = -9999.;
00363     sc_hybrid_eta[i] = -9999.;
00364     sc_hybrid_phi[i] = -9999.;
00365     //
00366     sc_multi5x5_et[i] = -9999.;
00367     sc_multi5x5_eta[i] = -9999.;
00368     sc_multi5x5_phi[i] = -9999.;
00369     //
00370   }
00371   // sort the energies of the first sc
00372   std::vector<double> ETsc1;
00373   std::vector<reco::SuperCluster>::const_iterator sc;
00374   for (sc = sc1->begin(); sc !=  sc1->end(); ++sc) {
00375     reco::SuperCluster mySc = *sc;
00376     double scEt = mySc.energy()/(cosh(mySc.eta()));
00377     ETsc1.push_back(scEt);
00378 
00379   }
00380   int *sorted1 = new int[n1];
00381   double *et1 = new double[n1];
00382   for (int i=0; i<n1; ++i) {
00383     et1[i] = ETsc1[i];
00384   }
00385   // array sorted now has the indices of the highest ET electrons
00386   TMath::Sort(n1, et1, sorted1, true);
00387   // .........................................................................
00388   std::vector<double> ETsc2;
00389   for (sc = sc2->begin(); sc !=  sc2->end(); ++sc) {
00390     reco::SuperCluster mySc = *sc;
00391     double scEt = mySc.energy()/(cosh(mySc.eta()));
00392     ETsc2.push_back(scEt);
00393 
00394   }
00395   int *sorted2 = new int[n2];
00396   double *et2 = new double[n2];
00397   for (int i=0; i<n2; ++i) {
00398     et2[i] = ETsc2[i];
00399   }
00400   // array sorted now has the indices of the highest ET electrons
00401   TMath::Sort(n2, et2, sorted2, true);
00402   //
00403   //
00404   for( int probeSc = 0; probeSc < n1; ++probeSc)
00405     {
00406       //std::cout<<"sorted["<< probeIt<< "]=" << sorted[probeIt] << std::endl;
00407       // break if you have more than the appropriate number of electrons
00408       if (probeSc >= 5) break;
00409       //
00410       int sc_index = sorted1[probeSc];
00411       std::vector<reco::SuperCluster>::const_iterator
00412         Rprobe = sc1->begin() + sc_index;
00413       //
00414       reco::SuperCluster sc0 = *Rprobe;
00415       // now keep the relevant stuff:
00416       sc_hybrid_et[probeSc] =  sc0.energy()/(cosh(sc0.eta()));
00417       sc_hybrid_eta[probeSc] = sc0.eta();
00418       sc_hybrid_phi[probeSc] = sc0.phi();
00419     }
00420   // .........................................................................
00421   for( int probeSc = 0; probeSc < n2; ++probeSc)
00422     {
00423       //std::cout<<"sorted["<< probeIt<< "]=" << sorted[probeIt] << std::endl;
00424       // break if you have more than the appropriate number of electrons
00425       if (probeSc >= 5) break;
00426       //
00427       int sc_index = sorted2[probeSc];
00428       std::vector<reco::SuperCluster>::const_iterator
00429         Rprobe = sc2->begin() + sc_index;
00430       //
00431       reco::SuperCluster sc0 = *Rprobe;
00432       // now keep the relevant stuff:
00433       sc_multi5x5_et[probeSc] =  sc0.energy()/(cosh(sc0.eta()));
00434       sc_multi5x5_eta[probeSc] = sc0.eta();
00435       sc_multi5x5_phi[probeSc] = sc0.phi();
00436     }
00437   delete [] sorted1;  delete [] sorted2;
00438   delete [] et1;     delete [] et2;
00440   //  edm::InputTag ctfTracksTag("generalTracks", "", InputTagEnding_);
00441   edm::Handle<reco::TrackCollection> ctfTracks;
00442   evt.getByLabel(ctfTracksTag_, ctfTracks);
00443   const reco::TrackCollection *ctf = ctfTracks.product();
00444   reco::TrackCollection::const_iterator tr;
00445   const int ntracks =  ctf->size();
00446   //
00447   // get the beam spot for the parameter of the track
00448   edm::Handle<reco::BeamSpot> pBeamSpot;
00449   evt.getByLabel("offlineBeamSpot", pBeamSpot);
00450   const reco::BeamSpot *bspot = pBeamSpot.product();
00451   const math::XYZPoint bspotPosition = bspot->position();
00452   //
00453   for (int i=0; i<20; ++i) {
00454     ctf_track_pt[i] = -9999.;
00455     ctf_track_eta[i] = -9999.;
00456     ctf_track_phi[i] = -9999.;
00457     ctf_track_vx[i] = -9999.; ctf_track_vy[i]=-9999.; ctf_track_vz[i] =-9999.;
00458     ctf_track_tip[i] = -9999.;    ctf_track_tip_bs[i] = -9999.;
00459   }
00460   //
00461   std::vector<double> ETtrack;
00462   for (tr = ctf->begin(); tr !=  ctf->end(); ++tr) {
00463     reco::Track mySc = *tr;
00464     double scEt = mySc.pt();
00465     ETtrack.push_back(scEt);
00466   }
00467   int *sortedTr = new int[ntracks];
00468   double *etTr = new double[ntracks];
00469   for (int i=0; i<ntracks; ++i) {
00470     etTr[i] = ETtrack[i];
00471   }
00472   // array sorted now has the indices of the highest ET electrons
00473   TMath::Sort(ntracks, etTr, sortedTr, true);
00474   //
00475   for( int probeSc = 0; probeSc < ntracks; ++probeSc)
00476     {
00477       //std::cout<<"sorted["<< probeIt<< "]=" << sorted[probeIt] << std::endl;
00478       // break if you have more than the appropriate number of electrons
00479       if (probeSc >= 20) break;
00480       //
00481       int sc_index = sortedTr[probeSc];
00482       std::vector<reco::Track>::const_iterator
00483         Rprobe = ctf->begin() + sc_index;
00484       //
00485       reco::Track sc0 = *Rprobe;
00486       // now keep the relevant stuff:
00487       ctf_track_pt[probeSc] =  sc0.pt();
00488       ctf_track_eta[probeSc] = sc0.eta();
00489       ctf_track_phi[probeSc] = sc0.phi();
00490       ctf_track_vx[probeSc] = sc0.vx();
00491       ctf_track_vy[probeSc] = sc0.vy();
00492       ctf_track_vz[probeSc] = sc0.vz();
00493       ctf_track_tip[probeSc] = -sc0.dxy();
00494       ctf_track_tip_bs[probeSc] = -sc0.dxy(bspotPosition);
00495     }
00496   delete [] sortedTr; delete [] etTr;
00497   //
00498   // keep 4 of the selectedLayer1Muons for reference
00499   edm::Handle<pat::MuonCollection> pMuons;
00500   evt.getByLabel("selectedLayer1Muons", pMuons);
00501   const pat::MuonCollection *pmuon = pMuons.product();
00502   pat::MuonCollection::const_iterator muon;
00503   const int nmuons =  pMuons->size();
00504   //
00505   for (int i=0; i<4; ++i) {
00506     muon_pt[i] = -9999.;
00507     muon_eta[i] = -9999.;
00508     muon_phi[i] = -9999.;
00509     muon_vx[i] = -9999.; muon_vy[i] = -9999.; muon_vz[i] = -9999.;
00510     muon_tip[i] = -9999.;    muon_tip_bs[i] = -9999.;
00511   }
00512   //
00513   std::vector<double> ETmuons;
00514   for (muon = pmuon->begin(); muon !=  pmuon->end(); ++muon) {
00515     pat::Muon mySc = *muon;
00516     double scEt = mySc.track()->pt();
00517     ETmuons.push_back(scEt);
00518   }
00519   int *sortedMu = new int[nmuons];
00520   double *etMu = new double[nmuons];
00521   for (int i=0; i<nmuons; ++i) {
00522     etMu[i] = ETmuons[i];
00523   }
00524   // array sorted now has the indices of the highest ET electrons
00525   TMath::Sort(nmuons, etMu, sortedMu, true);
00526   //
00527   for( int probeSc = 0; probeSc < nmuons; ++probeSc)
00528     {
00529       //std::cout<<"sorted["<< probeIt<< "]=" << sorted[probeIt] << std::endl;
00530       // break if you have more than the appropriate number of electrons
00531       if (probeSc >= 4) break;
00532       //
00533       int sc_index = sortedMu[probeSc];
00534       std::vector<pat::Muon>::const_iterator
00535         Rprobe = pmuon->begin() + sc_index;
00536       //
00537       pat::Muon sc0 = *Rprobe;
00538       // now keep the relevant stuff:
00539       muon_pt[probeSc] =  sc0.track()->pt();
00540       muon_eta[probeSc] = sc0.track()->eta();
00541       muon_phi[probeSc] = sc0.track()->phi();
00542       muon_vx[probeSc] = sc0.track()->vx();
00543       muon_vy[probeSc] = sc0.track()->vy();
00544       muon_vz[probeSc] = sc0.track()->vz();
00545       muon_tip[probeSc] = -sc0.track()->dxy();
00546       muon_tip_bs[probeSc] = -sc0.track()->dxy(bspotPosition);
00547     }
00548   delete [] sortedMu; delete [] etMu;
00549   //
00550   if (n1+n2+ntracks == 0) {
00551     std::cout << "Return: no sc in this event" << std::endl;
00552     return;
00553   }
00554   // /////////////////////////////////////////////////////////////////////////
00555   // electron details
00557   const int MAX_PROBES = 4;
00558   for(int i =0; i < MAX_PROBES; i++){
00559     probe_ele_eta_for_tree[i] = -99.0;
00560     probe_ele_et_for_tree[i] = -99.0;
00561     probe_ele_phi_for_tree[i] = -99.0;
00562     probe_ele_Xvertex_for_tree[i] = -99.0;
00563     probe_ele_Yvertex_for_tree[i] = -99.0;
00564     probe_ele_Zvertex_for_tree[i] = -99.0;
00565     probe_ele_tip[i] = -999.;    
00566 
00567     probe_sc_eta_for_tree[i] = -99.0;
00568     probe_sc_et_for_tree[i] = -99.0;
00569     probe_sc_phi_for_tree[i] = -99.0;
00570     
00571     probe_charge_for_tree[i] = -99;
00572     probe_sc_pass_fiducial_cut[i] = 0;
00573     probe_classification_index_for_tree[i]=-99; 
00574     //
00575     // probe isolation values ............
00576     probe_isolation_value[i] = 999.0;
00577     probe_iso_user[i] = 999.0;
00578     probe_ecal_isolation_value[i] = 999;
00579     probe_ecal_iso_user[i] = 999;
00580     probe_hcal_isolation_value[i] = 999;
00581     probe_hcal_iso_user[i] = 999;
00582 
00583     probe_ele_hoe[i]  = 999.;
00584     probe_ele_shh[i]  = 999.;
00585     probe_ele_sihih[i] = 999.;
00586     probe_ele_dhi[i]  = 999.;
00587     probe_ele_dfi[i]  = 999.;
00588     probe_ele_eop[i]  = 999.;
00589     probe_ele_pin[i]  = 999.;
00590     probe_ele_pout[i] = 999.;
00591     probe_ele_e5x5[i] = 999.;
00592     probe_ele_e2x5[i] = 999.;
00593     probe_ele_e1x5[i] = 999.;
00594 
00595     //
00596     //
00597     //for (int j=0; j<25; ++j) {
00598     //  probe_pass_trigger_cut[i][j]=0;
00599     //}
00600     //probe_hlt_matched_dr[i]=0;
00601     //probe_mc_matched[i] = 0;
00602     //probe_mc_matched_deta[i] = 999.;
00603     //probe_mc_matched_dphi[i] = 999.;
00604     //probe_mc_matched_denergy[i] = 999.;
00605     //probe_mc_matched_mother[i] = 999;
00606     //
00607     //
00608   }
00609   const pat::ElectronCollection *electrons= pElectrons.product();
00610   
00611 
00612   elec_number_in_event = electrons->size();
00613   //std::cout << "In this event " << elec_number_in_event << 
00614   //  " electrons were found" << std::endl;
00615   //  if (elec_number_in_event == 0) return;
00616  
00617   std::vector<pat::ElectronRef> UniqueElectrons;
00618   // edm::LogInfo("") << "Starting loop over electrons.";
00619   int index =0;
00620   //***********************************************************************
00621   // NEW METHOD by D WARDROPE implemented 26.05.08 ************************
00622   //************* DUPLICATE ******  REMOVAL *******************************
00623   // 02.06.08: due to a bug in the hybrid algorithm that affects detid ****
00624   //           we change detid matching to superCluster ref matching ******
00625   for(pat::ElectronCollection::const_iterator 
00626         elec = electrons->begin(); elec != electrons->end();++elec) {
00627     const pat::ElectronRef  electronRef(pElectrons, index);
00628     //Remove duplicate electrons which share a supercluster
00629     pat::ElectronCollection::const_iterator BestDuplicate = elec;
00630     int index2 = 0;
00631     for(pat::ElectronCollection::const_iterator
00632           elec2 = electrons->begin();
00633         elec2 != electrons->end(); ++elec2)
00634       {
00635         if(elec != elec2)
00636           {
00637             if( elec->superCluster() == elec2->superCluster())
00638               {
00639                 if(fabs(BestDuplicate->eSuperClusterOverP()-1.)
00640                    >= fabs(elec2->eSuperClusterOverP()-1.))
00641                   {
00642                     BestDuplicate = elec2;
00643                   }
00644               }
00645           }
00646         ++index2;
00647       }
00648     if(BestDuplicate == elec) UniqueElectrons.push_back(electronRef);
00649     ++index;
00650   }
00651   //
00652   // debugging: store electrons after duplicate removal
00653   elec_1_duplicate_removal = UniqueElectrons.size();
00654   //std::cout << "In this event there are " << elec_1_duplicate_removal 
00655   //        << " electrons" << std::endl;
00656   //
00657   //
00658   // duplicate removal is done now:
00659   //           the electron collection is in UniqueElectrons
00660   //
00661   // run over probes - now probe electrons and store
00662   //
00663   // the electron collection is now 
00664   // vector<reco::PixelMatchGsfElectronRef>   UniqueElectrons
00665   std::vector<double> ETs;
00666   std::vector<pat::ElectronRef>::const_iterator  elec;
00667   for (elec = UniqueElectrons.begin(); elec !=  UniqueElectrons.end(); ++elec) {
00668     pat::ElectronRef probeEle;
00669     probeEle = *elec;
00670     double probeEt = probeEle->caloEnergy()/(cosh(probeEle->caloPosition().eta()));
00671     ETs.push_back(probeEt);
00672 
00673   }
00674   int *sorted = new int[elec_1_duplicate_removal];
00675   double *et = new double[elec_1_duplicate_removal];
00676   //std::cout << "Elecs: " << elec_1_duplicate_removal << std::endl;
00677   for (int i=0; i<elec_1_duplicate_removal; ++i) {
00678     et[i] = ETs[i];
00679     //std::cout << "et["<< i << "]=" << et[i] << std::endl;
00680   }
00681   // array sorted now has the indices of the highest ET electrons
00682   TMath::Sort(elec_1_duplicate_removal, et, sorted, true);
00683   //
00684   //
00685   for( int probeIt = 0; probeIt < elec_1_duplicate_removal; ++probeIt)
00686     {
00687       //std::cout<<"sorted["<< probeIt<< "]=" << sorted[probeIt] << std::endl;
00688       // break if you have more than the appropriate number of electrons
00689       if (probeIt >= MAX_PROBES) break;
00690       //
00691       int elec_index = sorted[probeIt];
00692       std::vector<pat::ElectronRef>::const_iterator
00693         Rprobe = UniqueElectrons.begin() + elec_index;
00694       //
00695       pat::ElectronRef probeEle;
00696       probeEle = *Rprobe;
00697       double probeEt = probeEle->caloEnergy()/(cosh(probeEle->caloPosition().eta()));
00698       probe_sc_eta_for_tree[probeIt] = probeEle->caloPosition().eta();
00699       probe_sc_phi_for_tree[probeIt] = probeEle->caloPosition().phi();
00700       probe_sc_et_for_tree[probeIt] = probeEt;
00701       // fiducial cut ...............................
00702       if(fabs(probeEle->caloPosition().eta()) < BarrelMaxEta || 
00703          (fabs(probeEle->caloPosition().eta()) > EndcapMinEta && 
00704           fabs(probeEle->caloPosition().eta()) < EndcapMaxEta)){
00705         probe_sc_pass_fiducial_cut[probeIt] = 1;
00706       }
00707       //
00708       probe_charge_for_tree[probeIt] = probeEle->charge();
00709       probe_ele_eta_for_tree[probeIt] = probeEle->eta();
00710       probe_ele_et_for_tree[probeIt] = probeEle->et();
00711       probe_ele_phi_for_tree[probeIt] =probeEle->phi();
00712       probe_ele_Xvertex_for_tree[probeIt] =probeEle->vx();
00713       probe_ele_Yvertex_for_tree[probeIt] =probeEle->vy();
00714       probe_ele_Zvertex_for_tree[probeIt] =probeEle->vz();
00715       probe_classification_index_for_tree[probeIt] = 
00716         probeEle->classification();
00717       double ProbeTIP = probeEle->gsfTrack()->d0();
00718       probe_ele_tip[probeIt] = ProbeTIP;
00719       // isolation ..................................
00720       // these are the default values: trk 03, ecal, hcal 04
00721       // I know that there is a more direct way, but in this way it
00722       // is clearer what you get each time :P
00723       probe_isolation_value[probeIt] = probeEle->dr03IsolationVariables().tkSumPt;
00724       probe_ecal_isolation_value[probeIt] = probeEle->dr04IsolationVariables().ecalRecHitSumEt;
00725       probe_hcal_isolation_value[probeIt] = 
00726         probeEle->dr04IsolationVariables().hcalDepth1TowerSumEt + 
00727         probeEle->dr04IsolationVariables().hcalDepth2TowerSumEt;
00728       // one extra isos:
00729       probe_iso_user[probeIt] = probeEle->dr04IsolationVariables().tkSumPt;
00730       probe_ecal_iso_user[probeIt] = probeEle->dr03IsolationVariables().ecalRecHitSumEt;
00731       probe_hcal_iso_user[probeIt] = 
00732         probeEle->dr03IsolationVariables().hcalDepth1TowerSumEt + 
00733         probeEle->dr03IsolationVariables().hcalDepth2TowerSumEt;
00734       // ele id variables
00735       double hOverE = probeEle->hadronicOverEm();
00736       double deltaPhiIn = probeEle->deltaPhiSuperClusterTrackAtVtx();
00737       double deltaEtaIn = probeEle->deltaEtaSuperClusterTrackAtVtx();
00738       double eOverP = probeEle->eSuperClusterOverP();
00739       double pin  = probeEle->trackMomentumAtVtx().R(); 
00740       double pout = probeEle->trackMomentumOut().R(); 
00741       double sigmaee = probeEle->scSigmaEtaEta();
00742       double sigma_IetaIeta = probeEle->scSigmaIEtaIEta();
00743       // correct if in endcaps
00744       if( fabs (probeEle->caloPosition().eta()) > 1.479 )  {
00745         sigmaee = sigmaee - 0.02*(fabs(probeEle->caloPosition().eta()) -2.3);
00746       }
00747       //
00748       //double e5x5, e2x5Right, e2x5Left, e2x5Top, e2x5Bottom, e1x5;
00749       double e5x5, e2x5, e1x5;
00750       e5x5 = probeEle->scE5x5();
00751       e1x5 = probeEle->scE1x5();
00752       e2x5 = probeEle->scE2x5Max();
00753       //
00754       // electron ID variables
00755       probe_ele_hoe[probeIt] = hOverE;
00756       probe_ele_shh[probeIt] = sigmaee;
00757       probe_ele_sihih[probeIt] = sigma_IetaIeta;
00758       probe_ele_dfi[probeIt] = deltaPhiIn;
00759       probe_ele_dhi[probeIt] = deltaEtaIn;
00760       probe_ele_eop[probeIt] = eOverP;
00761       probe_ele_pin[probeIt] = pin;
00762       probe_ele_pout[probeIt] = pout;
00763       probe_ele_e5x5[probeIt] = e5x5;
00764       probe_ele_e2x5[probeIt] = e2x5;
00765       probe_ele_e1x5[probeIt] = e1x5;
00766  
00767       //
00768       // HLT filter ------------------------------------------------------
00769       //
00770       //
00771       // low luminosity filters
00772       /*************************************************************
00773       for (int filterNum=0; filterNum<10; ++filterNum) {
00774         int trigger_int_probe = 0;
00775         
00776         //double hlt_matched_dr   = -1.;
00777         const int nF(pHLTe29->sizeFilters());
00778         //
00779         // default (tag) trigger filter
00780         //
00781         // find how many relevant
00782         const int iF = pHLTe29->filterIndex(HLTFilterType_[filterNum]);
00783         // loop over these objects to see whether they match
00784         const trigger::TriggerObjectCollection& TOC(pHLTe29->getObjects());
00785         if (nF != iF) {
00786           // find how many objects there are
00787           const trigger::Keys& KEYS(pHLTe29->filterKeys(iF));
00788           const int nK(KEYS.size());
00789           for (int iTrig = 0;iTrig <nK; ++iTrig ) {
00790             const trigger::TriggerObject& TO(TOC[KEYS[iTrig]]);
00791             //std::cout << "--> filter: "<< HLTFilterType_[filterNum]  <<" TO id: " << TO.id() << std::endl;
00792             // this is better to be left out: HLT matching is with an HLT object
00793             // and we don't care what this object is
00794             //if (abs(TO.id())==11 ) { // demand it to be an electron
00795             double dr_ele_HLT = 
00796               reco::deltaR(probeEle->eta(), probeEle->phi(), TO.eta(), TO.phi());
00797             if (fabs(dr_ele_HLT) < ProbeHLTObjMaxDR) {++trigger_int_probe;
00798             //hlt_matched_dr = dr_ele_HLT;
00799             }
00800             //}
00801           }
00802         }
00803         //
00804         if(trigger_int_probe>0) probe_pass_trigger_cut[probeIt][filterNum] = 1;
00805         //probe_hlt_matched_dr[probeIt] = hlt_matched_dr;
00806       }
00807       // high lumi filters
00808       for (int filterNum=10; filterNum<25; ++filterNum) {
00809         int trigger_int_probe = 0;
00810         
00811         //double hlt_matched_dr   = -1.;
00812         const int nF(pHLTe31->sizeFilters());
00813         //
00814         // default (tag) trigger filter
00815         //
00816         // find how many relevant
00817         const int iF = pHLTe31->filterIndex(HLTFilterType_[filterNum]);
00818         // loop over these objects to see whether they match
00819         const trigger::TriggerObjectCollection& TOC(pHLTe31->getObjects());
00820         if (nF != iF) {
00821           // find how many objects there are
00822           const trigger::Keys& KEYS(pHLTe31->filterKeys(iF));
00823           const int nK(KEYS.size());
00824           for (int iTrig = 0;iTrig <nK; ++iTrig ) {
00825             const trigger::TriggerObject& TO(TOC[KEYS[iTrig]]);
00826             //if (abs(TO.id())==11 ) { // demand it to be an electron
00827             double dr_ele_HLT = 
00828               reco::deltaR(probeEle->eta(), probeEle->phi(), TO.eta(), TO.phi());
00829             if (fabs(dr_ele_HLT) < ProbeHLTObjMaxDR) {++trigger_int_probe;
00830             //hlt_matched_dr = dr_ele_HLT;
00831             }
00832           }
00833         }
00834       
00835         //
00836         if(trigger_int_probe>0) probe_pass_trigger_cut[probeIt][filterNum] = 1;
00837         //probe_hlt_matched_dr[probeIt] = hlt_matched_dr;
00838       }
00839       ******************************************/
00840       // ------------------------------------------------------------------
00841       //
00842       // MC Matching ......................................................
00843       // check whether these electrons are matched to a MC electron
00844       /*
00845       int mc_index = 0;
00846       int matched = 0; int mother_id = 999;
00847       double deta_matched = 999.;      double dphi_matched = 999.;
00848       double denergy_matched = 999.;
00849       for(reco::GenParticleCollection::const_iterator   McParticle = 
00850             McCand->begin(); McParticle != McCand->end();  ++McParticle)
00851         {
00852           // check only for electrons
00853           if(abs(McParticle->pdgId())==11 && McParticle->status()==1) {
00854             mc_index++;
00855             // check whether it matches a gsf electron
00856             double deta = McParticle->eta() - probeEle->eta();
00857             double dphi = McParticle->phi() - probeEle->phi();
00858             if ( fabs(deta) < MCMatch_Deta_  && fabs(dphi) < MCMatch_Dphi_){
00859               ++matched;
00860               deta_matched = deta; dphi_matched = dphi;
00861               denergy_matched = McParticle->energy() - probeEle->caloEnergy();
00862               // find the mother of the MC electron
00863               const reco::Candidate *mum;
00864               bool mother_finder = true;
00865               if (abs(McParticle->mother()->pdgId()) != 11)
00866                 mum = McParticle->mother();
00867               else if (abs(McParticle->mother()->mother()->pdgId())!= 11)
00868                 mum = McParticle->mother()->mother();
00869               else {
00870                 edm::LogInfo("info") << "Going too far to find the mum";
00871                 mother_finder = false;
00872               }          
00873               if (mother_finder) {
00874                 mother_id = mum->pdgId();
00875               }
00876             }
00877           }
00878         }
00879       probe_mc_matched[probeIt] = matched;
00880       probe_mc_matched_deta[probeIt] = deta_matched;
00881       probe_mc_matched_dphi[probeIt] = dphi_matched;
00882       probe_mc_matched_denergy[probeIt] = denergy_matched;
00883       probe_mc_matched_mother[probeIt] = mother_id;
00884       */
00885     }
00886   
00887   probe_tree->Fill();
00888   ++ tree_fills_;
00889   delete []  sorted;
00890   delete []  et;
00891 }
00892 
00893 
00894 // ------------ method called once each job just before starting event loop  --
00895 void 
00896 GenPurposeSkimmerData::beginJob()
00897 {
00898   //std::cout << "In beginJob()" << std::endl;
00899   TString filename_histo = outputFile_;
00900   histofile = new TFile(filename_histo,"RECREATE");
00901   tree_fills_ = 0;
00902 
00903   probe_tree =  new TTree("probe_tree","Tree to store probe variables");
00904 
00905   //probe_tree->Branch("probe_ele_eta",probe_ele_eta_for_tree,"probe_ele_eta[4]/D");
00906   //probe_tree->Branch("probe_ele_phi",probe_ele_phi_for_tree,"probe_ele_phi[4]/D");
00907   //probe_tree->Branch("probe_ele_et",probe_ele_et_for_tree,"probe_ele_et[4]/D");
00908   probe_tree->Branch("probe_ele_tip",probe_ele_tip,"probe_ele_tip[4]/D");
00909   probe_tree->Branch("probe_ele_vertex_x",probe_ele_Xvertex_for_tree,
00910                      "probe_ele_vertex_x[4]/D");
00911   probe_tree->Branch("probe_ele_vertex_y",probe_ele_Yvertex_for_tree,
00912                      "probe_ele_vertex_y[4]/D");
00913   probe_tree->Branch("probe_ele_vertex_z",probe_ele_Zvertex_for_tree,
00914                      "probe_ele_vertex_z[4]/D");
00915   probe_tree->Branch("probe_sc_eta",probe_sc_eta_for_tree,"probe_sc_eta[4]/D");
00916   probe_tree->Branch("probe_sc_phi",probe_sc_phi_for_tree,"probe_sc_phi[4]/D");
00917   probe_tree->Branch("probe_sc_et",probe_sc_et_for_tree,"probe_sc_et[4]/D");
00918 
00919   // trigger related variables
00920   //probe_tree->Branch("probe_trigger_cut",probe_pass_trigger_cut,"probe_trigger_cut[4][25]/I");
00921   //probe_tree->Branch("probe_hlt_matched_dr", probe_hlt_matched_dr,"probe_hlt_matched_dr[4]/D");
00922   // mc matching to electrons
00923   //  probe_tree->Branch("probe_mc_matched",probe_mc_matched,"probe_mc_matched[4]/I");
00924   //probe_tree->Branch("probe_mc_matched_deta",probe_mc_matched_deta,
00925   //         "probe_mc_matched_deta[4]/D");
00926   //probe_tree->Branch("probe_mc_matched_dphi",probe_mc_matched_dphi,
00927   //                 "probe_mc_matched_dphi[4]/D");
00928   //probe_tree->Branch("probe_mc_matched_denergy",probe_mc_matched_denergy,
00929   //                 "probe_mc_matched_denergy[4]/D");
00930   //probe_tree->Branch("probe_mc_matched_mother",probe_mc_matched_mother,
00931   //                 "probe_mc_matched_mother[4]/I");
00932   //
00933   probe_tree->Branch("probe_charge",probe_charge_for_tree,"probe_charge[4]/I");
00934   //probe_tree->Branch("probe_sc_fiducial_cut",probe_sc_pass_fiducial_cut,
00935   //                 "probe_sc_fiducial_cut[4]/I");
00936 
00937 
00938 
00939   //probe_tree->Branch("probe_classification",
00940   //        probe_classification_index_for_tree,"probe_classification[4]/I");
00941   //
00942   // Isolation related variables ........................................
00943   //
00944   probe_tree->Branch("probe_isolation_value",probe_isolation_value, "probe_isolation_value[4]/D");
00945   probe_tree->Branch("probe_ecal_isolation_value",probe_ecal_isolation_value, "probe_ecal_isolation_value[4]/D");
00946   probe_tree->Branch("probe_hcal_isolation_value",probe_hcal_isolation_value,"probe_hcal_isolation_value[4]/D");
00947   //
00948   probe_tree->Branch("probe_iso_user",     probe_iso_user,      "probe_iso_user[4]/D");
00949   probe_tree->Branch("probe_ecal_iso_user",probe_ecal_iso_user, "probe_ecal_iso_user[4]/D");
00950   probe_tree->Branch("probe_hcal_iso_user",probe_hcal_iso_user, "probe_hcal_iso_user[4]/D");
00951 
00952   //......................................................................
00953   // Electron ID Related variables .......................................
00954   probe_tree->Branch("probe_ele_hoe",probe_ele_hoe, "probe_ele_hoe[4]/D");
00955   //probe_tree->Branch("probe_ele_shh",probe_ele_shh, "probe_ele_shh[4]/D");
00956   probe_tree->Branch("probe_ele_sihih",probe_ele_sihih,"probe_ele_sihih[4]/D");
00957   probe_tree->Branch("probe_ele_dfi",probe_ele_dfi, "probe_ele_dfi[4]/D");
00958   probe_tree->Branch("probe_ele_dhi",probe_ele_dhi, "probe_ele_dhi[4]/D");
00959   probe_tree->Branch("probe_ele_eop",probe_ele_eop, "probe_ele_eop[4]/D");
00960   probe_tree->Branch("probe_ele_pin",probe_ele_pin, "probe_ele_pin[4]/D");
00961   probe_tree->Branch("probe_ele_pout",probe_ele_pout, "probe_ele_pout[4]/D");
00962   // probe_tree->Branch("probe_ele_e5x5",probe_ele_e5x5, "probe_ele_e5x5[4]/D");
00963   //probe_tree->Branch("probe_ele_e2x5",probe_ele_e2x5, "probe_ele_e2x5[4]/D");
00964   //probe_tree->Branch("probe_ele_e1x5",probe_ele_e1x5, "probe_ele_e1x5[4]/D");
00965 
00966   //.......................................................................
00967   //
00968   // each entry for each trigger path
00969   //probe_tree->Branch("event_HLTPath",event_HLTPath,"event_HLTPath[25]/I");
00970   //probe_tree->Branch("numberOfHLTFilterObjects", numberOfHLTFilterObjects,
00971   //                 "numberOfHLTFilterObjects[25]/I");
00972   //
00973   // debugging info:
00974   //probe_tree->Branch("elec_number_in_event",&elec_number_in_event,"elec_number_in_event/I");
00975   probe_tree->Branch("elec_1_duplicate_removal",&elec_1_duplicate_removal,"elec_1_duplicate_removal/I");
00976   //
00977 
00978   // Missing ET in the event
00979   probe_tree->Branch("event_MET",&event_MET,"event_MET/D");
00980   probe_tree->Branch("event_MET_phi",&event_MET_phi,"event_MET_phi/D");
00981   //  probe_tree->Branch("event_MET_sig",&event_MET_sig,"event_MET_sig/D");
00982   probe_tree->Branch("event_mcMET",&event_mcMET,"event_mcMET/D");
00983   probe_tree->Branch("event_mcMET_phi",&event_mcMET_phi,"event_mcMET_phi/D");
00984   //
00985   probe_tree->Branch("event_tcMET",&event_tcMET,"event_tcMET/D");
00986   probe_tree->Branch("event_tcMET_phi",&event_tcMET_phi,"event_tcMET_phi/D");
00987   //  probe_tree->Branch("event_tcMET_sig",&event_tcMET_sig,"event_tcMET_sig/D");
00988 
00989   probe_tree->Branch("event_pfMET",&event_pfMET,"event_pfMET/D");
00990   probe_tree->Branch("event_pfMET_phi",&event_pfMET_phi,"event_pfMET_phi/D");
00991   //  probe_tree->Branch("event_pfMET_sig",&event_pfMET_sig,"event_pfMET_sig/D");
00992 
00993   //  probe_tree->Branch("event_genMET",&event_genMET,"event_genMET/D");
00994   //  probe_tree->Branch("event_genMET_phi",&event_genMET_phi, "event_genMET_phi/D");
00995   //  probe_tree->Branch("event_genMET_sig",&event_genMET_sig, "event_genMET_sig/D");
00996   //..... type 1 corrected MET
00997   probe_tree->Branch("event_t1MET", &event_t1MET, "event_t1MET/D");
00998   probe_tree->Branch("event_t1MET_phi", &event_t1MET_phi,"event_t1MET_phi/D");
00999   //probe_tree->Branch("event_t1MET_sig",&event_t1MET_sig,"event_t1MET_sig/D");
01000 
01001   //
01002   // some sc related variables
01003   probe_tree->Branch("sc_hybrid_et", sc_hybrid_et, "sc_hybrid_et[5]/D");
01004   probe_tree->Branch("sc_hybrid_eta", sc_hybrid_eta, "sc_hybrid_eta[5]/D");
01005   probe_tree->Branch("sc_hybrid_phi", sc_hybrid_phi, "sc_hybrid_phi[5]/D");
01006   //
01007   probe_tree->Branch("sc_multi5x5_et",sc_multi5x5_et, "sc_multi5x5_et[5]/D");
01008   probe_tree->Branch("sc_multi5x5_eta",sc_multi5x5_eta,"sc_multi5x5_eta[5]/D");
01009   probe_tree->Branch("sc_multi5x5_phi",sc_multi5x5_phi,"sc_multi5x5_phi[5]/D");
01010   // /////////////////////////////////////////////////////////////////////////
01011   // general tracks in the event: keep 20 tracks
01012   probe_tree->Branch("ctf_track_pt",  ctf_track_pt,  "ctf_track_pt[20]/D");
01013   probe_tree->Branch("ctf_track_eta", ctf_track_eta, "ctf_track_eta[20]/D");
01014   probe_tree->Branch("ctf_track_phi", ctf_track_phi, "ctf_track_phi[20]/D");
01015   probe_tree->Branch("ctf_track_vx", ctf_track_vx, "ctf_track_vx[20]/D");
01016   probe_tree->Branch("ctf_track_vy", ctf_track_vy, "ctf_track_vy[20]/D");
01017   probe_tree->Branch("ctf_track_vz", ctf_track_vz, "ctf_track_vz[20]/D");
01018   probe_tree->Branch("ctf_track_tip", ctf_track_tip, "ctf_track_tip[20]/D");
01019   probe_tree->Branch("ctf_track_tip_bs", ctf_track_tip_bs, 
01020                      "ctf_track_tip_bs[20]/D");
01021   //
01022   probe_tree->Branch("muon_pt",  muon_pt,  "muon_pt[4]/D");
01023   probe_tree->Branch("muon_eta", muon_eta, "muon_eta[4]/D");
01024   probe_tree->Branch("muon_phi", muon_phi, "muon_phi[4]/D");
01025   probe_tree->Branch("muon_vx", muon_vx, "muon_vx[4]/D");
01026   probe_tree->Branch("muon_vy", muon_vy, "muon_vy[4]/D");
01027   probe_tree->Branch("muon_vz", muon_vz, "muon_vz[4]/D");
01028   probe_tree->Branch("muon_tip", muon_tip, "muon_tip[4]/D");
01029   probe_tree->Branch("muon_tip_bs", muon_tip_bs, "muon_tip_bs[4]/D");
01030 
01031 }
01032 
01033 // ------------ method called once each job just after ending the event loop  -
01034 void 
01035 GenPurposeSkimmerData::endJob() {
01036   //std::cout << "In endJob()" << std::endl;
01037   if (tree_fills_ == 0) {
01038     std::cout << "Empty tree: no output..." << std::endl;
01039     return;
01040   }
01041   //probe_tree->Print();
01042   histofile->Write();
01043   histofile->Close();
01044 
01045 }
01046 
01047 
01048 //define this as a plug-in
01049 DEFINE_FWK_MODULE(GenPurposeSkimmerData);
/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/ElectroWeakAnalysis/WENu/src/GenPurposeSkimmerData.cc
