#include <GammaJetAnalysis.h>

Inheritance diagram for cms::GammaJetAnalysis:

Public Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &)
virtual void	beginJob ()
virtual void	endJob ()
	GammaJetAnalysis (const edm::ParameterSet &)
	~GammaJetAnalysis ()
Private Attributes
bool	allowMissingInputs_
int	Code [4000]
double	CutOnEgammaEnergy_
int	EcalClusDet [20]
std::string	ecalInput_
float	ecut [3][3]
int	event
std::string	fOutputFileName
std::string	gammaClus_
float	GammaIsoEcal [9][20]
float	GammaIsoHcal [9][20]
float	GammaRecoEt [20]
float	GammaRecoEta [20]
float	GammaRecoPhi [20]
const CaloGeometry *	geo
std::string	hbheInput_
float	HcalDet [8000]
float	HcalRecoEt [8000]
float	HcalRecoEta [8000]
float	HcalRecoPhi [8000]
std::string	hfInput_
std::string	hoInput_
TFile *	hOutputFile
std::string	jetCalo_
float	JetGenEt [10]
float	JetGenEta [10]
float	JetGenPhi [10]
float	JetGenType [10]
float	JetRecoEt [10]
float	JetRecoEta [10]
float	JetRecoPhi [10]
float	JetRecoType [10]
int	Mother1 [4000]
std::string	myName
std::ofstream *	myout_ecal
std::ofstream *	myout_hcal
std::ofstream *	myout_jet
std::ofstream *	myout_part
std::ofstream *	myout_photon
TTree *	myTree
std::string	nameProd_
int	NumGenJets
int	NumPart
int	NumRecoGamma
int	NumRecoHcal
int	NumRecoJets
int	NumRecoTrack
float	parte [4000]
float	partm [4000]
float	partpx [4000]
float	partpy [4000]
float	partpz [4000]
float	partvt [4000]
float	partvx [4000]
float	partvy [4000]
float	partvz [4000]
float	risol [3]
int	run
int	Status [4000]
float	TrackRecoEt [10]
float	TrackRecoEta [10]
float	TrackRecoPhi [10]
std::string	Tracks_
bool	useMC

Detailed Description

Definition at line 50 of file GammaJetAnalysis.h.

Constructor & Destructor Documentation

cms::GammaJetAnalysis::GammaJetAnalysis ( const edm::ParameterSet & iConfig ) [explicit]

Definition at line 44 of file GammaJetAnalysis.cc.

References edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), and PatBasicFWLiteJetAnalyzer_Selector_cfg::useMC.

{
  // get names of modules, producing object collections

  nameProd_ = iConfig.getUntrackedParameter<std::string>("nameProd");
  jetCalo_ = iConfig.getUntrackedParameter<std::string>("jetCalo","GammaJetJetBackToBackCollection");
  gammaClus_ = iConfig.getUntrackedParameter<std::string>("gammaClus","GammaJetGammaBackToBackCollection");
  ecalInput_=iConfig.getUntrackedParameter<std::string>("ecalInput","GammaJetEcalRecHitCollection");
  hbheInput_ = iConfig.getUntrackedParameter<std::string>("hbheInput");
  hoInput_ = iConfig.getUntrackedParameter<std::string>("hoInput");
  hfInput_ = iConfig.getUntrackedParameter<std::string>("hfInput");
  Tracks_ = iConfig.getUntrackedParameter<std::string>("Tracks","GammaJetTracksCollection");
  CutOnEgammaEnergy_  = iConfig.getParameter<double>("CutOnEgammaEnergy");

  myName = iConfig.getParameter<std::string> ("textout");
  useMC = iConfig.getParameter<bool>("useMCInfo"); 
  allowMissingInputs_=iConfig.getUntrackedParameter<bool>("AllowMissingInputs",false);
  // get name of output file with histogramms
  fOutputFileName = iConfig.getUntrackedParameter<string>("HistOutFile"); 
  risol[0] = 0.5;
  risol[1] = 0.7;
  risol[2] = 1.0;
  
  ecut[0][0] = 0.09;
  ecut[0][1] = 0.18;
  ecut[0][2] = 0.27;
  
  ecut[1][0] = 0.45;
  ecut[1][1] = 0.9;
  ecut[1][2] = 1.35;
  
  ecut[2][0] = 0.5;
  ecut[2][1] = 1.;
  ecut[2][2] = 1.5;
}

cms::GammaJetAnalysis::~GammaJetAnalysis ( )

Definition at line 81 of file GammaJetAnalysis.cc.

{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}

Member Function Documentation

void cms::GammaJetAnalysis::analyze	(	const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)		`[virtual]`

Implements edm::EDAnalyzer.

Definition at line 161 of file GammaJetAnalysis.cc.

References Reference_intrackfit_cff::barrel, gather_cfg::cout, cond::rpcobgas::detid, ExpressReco_HICollisions_FallBack::e, DetId::Ecal, eta(), PV3DBase< T, PVType, FrameType >::eta(), edm::EventID::event(), edm::EventSetup::get(), edm::Event::getAllProvenance(), edm::Event::getByLabel(), DetId::Hcal, i, edm::EventBase::id(), j, metsig::jet, analyzePatCleaning_cfg::jets, PV3DBase< T, PVType, FrameType >::phi(), phi, pos, edm::ESHandle< T >::product(), edm::Handle< T >::product(), dt_offlineAnalysis_common_cff::reco, edm::EventID::run(), CrabTask::run, and mathSSE::sqrt().

{

   edm::ESHandle<CaloGeometry> pG;
   iSetup.get<CaloGeometryRecord>().get(pG);
   geo = pG.product();


  using namespace edm;
  std::vector<Provenance const*> theProvenance;
  iEvent.getAllProvenance(theProvenance);
//  for( std::vector<Provenance const*>::const_iterator ip = theProvenance.begin();
//                                                      ip != theProvenance.end(); ip++)
//  {
//     cout<<" Print all module/label names "<<(**ip).moduleName()<<" "<<(**ip).moduleLabel()<<
//     " "<<(**ip).productInstanceName()<<endl;
//  }
// Load generator information
// write HEPEVT block into file
   run = iEvent.id().run();
   event = iEvent.id().event();
  (*myout_part)<<"Event "<<iEvent.id().run()<<" "<<iEvent.id().event()<<endl;
  (*myout_jet)<<"Event "<<iEvent.id().run()<<" "<<iEvent.id().event()<<endl;
  (*myout_hcal)<<"Event "<<iEvent.id().run()<<" "<<iEvent.id().event()<<endl;
  (*myout_ecal)<<"Event "<<iEvent.id().run()<<" "<<iEvent.id().event()<<endl;
  (*myout_photon)<<"Event "<<iEvent.id().run()<<" "<<iEvent.id().event()<<endl;
  

    std::vector<edm::InputTag>::const_iterator ic;
    int jettype = 0;
    int jetexist = -100;
    int reco = 1;
    double etlost = -100.1;
    
    NumRecoJets = 0;
    
     try {
       
       edm::Handle<reco::CaloJetCollection> jets;
       iEvent.getByLabel(nameProd_, jetCalo_, jets);
       reco::CaloJetCollection::const_iterator jet = jets->begin ();
       cout<<" Size of Calo jets "<<jets->size()<<endl;
       jettype++;
       
       if(jets->size() > 0 )
       {
         int ij = 0;
         for (; jet != jets->end (); jet++) 
         {
            cout<<" Jet et "<<(*jet).et()<<" "<<(*jet).eta()<<" "<<(*jet).phi()<<endl;
            ij++;
            if(ij<4) (*myout_jet)<<jettype<<" "<<reco<<" "<<ij<<" "<<(*jet).et()<<" "<<(*jet).eta()<<" "<<(*jet).phi()
            <<" "<<iEvent.id().event()<<endl;
            jetexist = ij;
            if( NumRecoJets < 8 )
            {
             JetRecoEt[NumRecoJets] = (*jet).et();
             JetRecoEta[NumRecoJets] = (*jet).eta();
             JetRecoPhi[NumRecoJets] = (*jet).phi();
             JetRecoType[NumRecoJets] = jettype;
             NumRecoJets++;
            }
         }
       }
     } catch (cms::Exception& e) { // can't find it!
       if (!allowMissingInputs_) {
         cout<< " Calojets are missed "<<endl;
         throw e;
        }        
     }
   
     cout<<" We filled CaloJet part "<<jetexist<<endl;
     
     if( jetexist < 0 ) (*myout_jet)<<jetexist<<" "<<reco<<" "<<etlost
                         <<" "<<etlost<<" "<<etlost
                         <<" "<<iEvent.id().event()<<endl;
// Load EcalRecHits
  
    std::vector<edm::InputTag>::const_iterator i;
    vector<CaloRecHit> theRecHits;
      
    try {
      
      edm::Handle<EcalRecHitCollection> ec;
      iEvent.getByLabel(nameProd_, ecalInput_,ec);
      
       for(EcalRecHitCollection::const_iterator recHit = (*ec).begin();
                                                recHit != (*ec).end(); ++recHit)
       {
// EcalBarrel = 1, EcalEndcap = 2

         GlobalPoint pos = geo->getPosition(recHit->detid());
         theRecHits.push_back(*recHit);

         if( (*recHit).energy()> ecut[recHit->detid().subdetId()-1][0] )
                    (*myout_ecal)<<recHit->detid().subdetId()<<" "<<(*recHit).energy()<<" "<<pos.phi()<<" "<<pos.eta()
                                 <<" "<<iEvent.id().event()<<endl;
             
       } 
      
    } catch (cms::Exception& e) { // can't find it!
    if (!allowMissingInputs_) {
      cout<<" Ecal collection is missed "<<endl;
      throw e;
     } 
    }

     cout<<" Fill EcalRecHits "<<endl;
//  cout<<" Start to get hbhe "<<endl;
// Hcal Barrel and endcap for isolation
    try {
      edm::Handle<HBHERecHitCollection> hbhe;
      iEvent.getByLabel(nameProd_,hbheInput_,hbhe);

//      (*myout_hcal)<<(*hbhe).size()<<endl;
  for(HBHERecHitCollection::const_iterator hbheItr = (*hbhe).begin();

      hbheItr != (*hbhe).end(); ++hbheItr)
      {
        DetId id = (hbheItr)->detid();
        GlobalPoint pos = geo->getPosition(hbheItr->detid());
        (*myout_hcal)<<id.subdetId()<<" "<<(*hbheItr).energy()<<" "<<pos.phi()<<
                                      " "<<pos.eta()<<" "<<iEvent.id().event()<<endl;
        theRecHits.push_back(*hbheItr);

      }
    } catch (cms::Exception& e) { // can't find it!
      if (!allowMissingInputs_) {
        cout<<" HBHE collection is missed "<<endl;
        throw e;
      }
    }

 
 for(int i = 0; i<9; i++)
 {
    for(int j= 0; j<10; j++) GammaIsoEcal[i][j] = 0.;
 }

// Load Ecal clusters
 jetexist = -100; 
 int barrel = 1;
 NumRecoGamma = 0;
 
 try {
 int ij = 0;
  // Get island super clusters after energy correction
  Handle<reco::SuperClusterCollection> eclus;
  iEvent.getByLabel(nameProd_,gammaClus_, eclus);
  const reco::SuperClusterCollection* correctedSuperClusters=eclus.product();
  // loop over the super clusters and fill the histogram
  for(reco::SuperClusterCollection::const_iterator aClus = correctedSuperClusters->begin();
                                                           aClus != correctedSuperClusters->end(); aClus++) {
    double vet = aClus->energy()/cosh(aClus->eta());
    cout<<" Supercluster " << ij<<" Et "<< vet <<" energy "<<aClus->energy()<<" eta "<<aClus->eta()<<" Cut "<<CutOnEgammaEnergy_<<endl;

    if(vet>CutOnEgammaEnergy_) {
      ij++;
      float gammaiso_ecal[9] = {0.,0.,0.,0.,0.,0.,0.,0.,0.};
     for(vector<CaloRecHit>::iterator it = theRecHits.begin(); it != theRecHits.end(); it++)
      {
           GlobalPoint pos = geo->getPosition(it->detid());
           double eta = pos.eta();
           double phi = pos.phi();
           double deta = fabs(eta-aClus->eta());
           double dphi = fabs(phi-aClus->phi());
           if(dphi>4.*atan(1.)) dphi = 8.*atan(1.)-dphi;
           double dr = sqrt(deta*deta+dphi*dphi);
           
           double rmin = 0.07;
           if( fabs(aClus->eta()) > 1.47 ) rmin = 0.07*(fabs(aClus->eta())-.47)*1.2;
           if( fabs(aClus->eta()) > 2.2 ) rmin = 0.07*(fabs(aClus->eta())-.47)*1.4;
           
           int itype_ecal = 0;
           double ecutn = 0.;
           for (int i = 0; i<3; i++)
           {
             for (int j = 0; j<3; j++)
             {
             
                if(it->detid().det() == DetId::Ecal ) 
                {
                  if(it->detid().subdetId() == 1) ecutn = ecut[0][j];
                  if(it->detid().subdetId() == 2) ecutn = ecut[1][j];
                  if( dr>rmin && dr<risol[i])
                  {
                   if((*it).energy() > ecutn) gammaiso_ecal[itype_ecal] = gammaiso_ecal[itype_ecal]+(*it).energy()/cosh(eta);
                  } 
                }
                
                if(it->detid().det() == DetId::Hcal ) 
                {
                   ecutn = ecut[2][j];
                   if( dr>rmin && dr<risol[i])
                   {
                     if((*it).energy() > ecutn) 
                     {
                        gammaiso_ecal[itype_ecal] = gammaiso_ecal[itype_ecal]+(*it).energy()/cosh(eta);
                     }
                   }
                } 
                jetexist = ij;
                itype_ecal++;
                
             } // Ecal
           } // cycle on iso radii      
      } // cycle on rechits
      
      
// Fill Tree      
           if( NumRecoGamma < 10 ) 
           {
            for (int ii = 0; ii<9 ; ii++)
            {
             GammaIsoEcal[ii][NumRecoGamma] = gammaiso_ecal[ii]; 
            } 
             EcalClusDet[NumRecoGamma] = 1;
             GammaRecoEt[NumRecoGamma] = vet;
             GammaRecoEta[NumRecoGamma] = aClus->eta();
             GammaRecoPhi[NumRecoGamma] = aClus->phi();
             NumRecoGamma++;
            }
    (*myout_photon)<<ij<<" "<<barrel<<" "<<vet<<" "<<aClus->eta()<<" "<<aClus->phi()<<" "<<iEvent.id().event()<<endl;
    (*myout_photon)<<ij<<" "<<gammaiso_ecal[0]<<" "<<gammaiso_ecal[1] <<" "<<gammaiso_ecal[2]<<" "<<gammaiso_ecal[3]
                   <<" "<<gammaiso_ecal[4]<<" "<<gammaiso_ecal[5]<<" "<<gammaiso_ecal[6]<<" "<<gammaiso_ecal[7]<<" "<<gammaiso_ecal[8]<<endl;
      
       jetexist = ij;
    } //vet  
  } // number of superclusters
  } catch (cms::Exception& e) { // can't find it!
    if (!allowMissingInputs_) {
       cout<<" Ecal barrel clusters are missed "<<endl;
       throw e;
    }   
  }

    cout<<" After iso cuts "<<jetexist<<endl;

    double ecluslost = -100.1;
    if(jetexist<0) (*myout_photon)<<jetexist<<" "<<barrel<<" "<<ecluslost<<" "<<ecluslost
                                  <<" "<<ecluslost<<" "<<iEvent.id().event()<<endl;
  
    cout<<" Event is ready "<<endl;
   
    myTree->Fill();
   
} // analyze method

void cms::GammaJetAnalysis::beginJob ( void ) [virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 89 of file GammaJetAnalysis.cc.

References gather_cfg::cout, event(), and CrabTask::run.

{
   hOutputFile   = new TFile( fOutputFileName.c_str(), "RECREATE" ) ; 
   myTree = new TTree("GammaJet","GammaJet Tree");
   myTree->Branch("run",  &run, "run/I");
   myTree->Branch("event",  &event, "event/I");
   
   NumRecoJets = 0;
   NumGenJets = 0;
   NumRecoGamma = 0;
   NumRecoTrack = 0;
   NumPart = 0;
// Jet block
   myTree->Branch("NumRecoJets", &NumRecoJets, "NumRecoJets/I");   
   myTree->Branch("JetRecoEt",  JetRecoEt, "JetRecoEt[10]/F");
   myTree->Branch("JetRecoEta",  JetRecoEta, "JetRecoEta[10]/F");
   myTree->Branch("JetRecoPhi",  JetRecoPhi, "JetRecoPhi[10]/F");
   myTree->Branch("JetRecoType",  JetRecoType, "JetRecoType[10]/F");
   
// Gamma block for ECAL isolated gammas
   myTree->Branch("NumRecoGamma", &NumRecoGamma, "NumRecoGamma/I");
   myTree->Branch("EcalClusDet", &EcalClusDet, "EcalClusDet[20]/I");
   myTree->Branch("GammaRecoEt",  GammaRecoEt, "GammaRecoEt[20]/F");
   myTree->Branch("GammaRecoEta",  GammaRecoEta, "GammaRecoEta[20]/F");
   myTree->Branch("GammaRecoPhi",  GammaRecoPhi, "GammaRecoPhi[20]/F");
   myTree->Branch("GammaIsoEcal",  GammaIsoEcal, "GammaIsoEcal[9][20]/F");

// Tracks block
   myTree->Branch("NumRecoTrack", &NumRecoTrack, "NumRecoTrack/I");
   myTree->Branch("TrackRecoEt",  TrackRecoEt, "TrackRecoEt[200]/F");
   myTree->Branch("TrackRecoEta",  TrackRecoEta, "TrackRecoEta[200]/F");
   myTree->Branch("TrackRecoPhi",  TrackRecoPhi, "TrackRecoPhi[200]/F");

// end of tree declaration

//   edm::ESHandle<CaloGeometry> pG;
//   iSetup.get<CaloGeometryRecord>().get(pG);
//   geo = pG.product();

  myout_part = new ofstream((myName+"_part.dat").c_str()); 
  if(!myout_part) cout << " Output file not open!!! "<<endl;
  myout_hcal = new ofstream((myName+"_hcal.dat").c_str()); 
  if(!myout_hcal) cout << " Output file not open!!! "<<endl;
  myout_ecal = new ofstream((myName+"_ecal.dat").c_str()); 
  if(!myout_ecal) cout << " Output file not open!!! "<<endl;
  
  myout_jet = new ofstream((myName+"_jet.dat").c_str()); 
  if(!myout_jet) cout << " Output file not open!!! "<<endl;
  myout_photon = new ofstream((myName+"_photon.dat").c_str()); 
  if(!myout_photon) cout << " Output file not open!!! "<<endl;
   
   
}

void cms::GammaJetAnalysis::endJob ( void ) [virtual]

Reimplemented from edm::EDAnalyzer.

Definition at line 143 of file GammaJetAnalysis.cc.

References gather_cfg::cout.

{

   cout << "===== Start writing user histograms =====" << endl;
   hOutputFile->SetCompressionLevel(2);
   hOutputFile->cd();
   myTree->Write();
   hOutputFile->Close() ;
   cout << "===== End writing user histograms =======" << endl; 
}