PhotonsWithConversionsAnalyzer.cc

Go to the documentation of this file.

#include <iostream>
//
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
//
#include "RecoEgamma/Examples/plugins/PhotonsWithConversionsAnalyzer.h"
#include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
#include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruth.h"
#include "RecoEgamma/EgammaMCTools/interface/ElectronMCTruth.h"
//
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/Track/interface/SimTrack.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/Vertex/interface/SimVertex.h"
#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
//
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
//
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackExtra.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/EgammaCandidates/interface/Conversion.h"
#include "DataFormats/EgammaCandidates/interface/ConversionFwd.h"
#include "DataFormats/EgammaCandidates/interface/Photon.h"
#include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
//
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "MagneticField/Engine/interface/MagneticField.h"

//
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TTree.h"
#include "TVector3.h"
#include "TProfile.h"
//



using namespace std;


PhotonsWithConversionsAnalyzer::PhotonsWithConversionsAnalyzer( const edm::ParameterSet& pset )
  {

  photonCollectionProducer_ = pset.getParameter<std::string>("phoProducer");
  photonCollection_ = pset.getParameter<std::string>("photonCollection");
  //


}



PhotonsWithConversionsAnalyzer::~PhotonsWithConversionsAnalyzer() {


  delete thePhotonMCTruthFinder_;

}


void PhotonsWithConversionsAnalyzer::beginJob()
{


  nEvt_=0;
  nMCPho_=0;
  nMatched_=0;


  thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();



  edm::Service<TFileService> fs;


  h_ErecoEMC_  = fs->make<TH1F>("deltaE","    delta(reco-mc) energy",100,0.,2.);
  h_deltaPhi_ = fs->make<TH1F>("deltaPhi","  delta(reco-mc) phi",100,-0.1, 0.1);
  h_deltaEta_ = fs->make<TH1F>("deltaEta","  delta(reco-mc) eta",100,-0.05, 0.05);

  h_MCphoE_ = fs->make<TH1F>("MCphoE","MC photon energy",100,0.,100.);
  h_MCphoPhi_ = fs->make<TH1F>("MCphoPhi","MC photon phi",40,-3.14, 3.14);
  h_MCphoEta_ = fs->make<TH1F>("MCphoEta","MC photon eta",40,-3., 3.);


  h_MCConvE_ = fs->make<TH1F>("MCConvE","MC converted photon energy",100,0.,100.);
  h_MCConvPt_ = fs->make<TH1F>("MCConvPt","MC converted photon pt",100,0.,100.);
  h_MCConvEta_ = fs->make<TH1F>("MCConvEta","MC converted photon eta",50, 0., 2.5);

  h_scE_ = fs->make<TH1F>("scE","SC Energy ",100,0., 200.);
  h_scEt_ = fs->make<TH1F>("scEt","SC Et ",100,0., 200.);
  h_scEta_ = fs->make<TH1F>("scEta","SC Eta ",40,-3., 3.);
  h_scPhi_ = fs->make<TH1F>("scPhi","SC Phi ",40, -3.14, 3.14);
  //
  h_phoE_ = fs->make<TH1F>("phoE","Photon Energy ",100,0., 200.);
  h_phoEta_ = fs->make<TH1F>("phoEta","Photon Eta ",40,-3., 3.);
  h_phoPhi_ = fs->make<TH1F>("phoPhi","Photon  Phi ",40,  -3.14, 3.14);

  // Recontructed tracks from converted photon candidates
  h2_tk_nHitsVsR_ = fs->make<TH2F>("tknHitsVsR","Tracks Hits vs R  ", 12,0.,120.,20,0.5, 20.5);
  h2_tk_inPtVsR_  = fs->make<TH2F>("tkInPtvsR","Tracks inner Pt vs R  ", 12,0.,120.,100,0., 100.);


  return ;
}


float PhotonsWithConversionsAnalyzer::etaTransformation(  float EtaParticle , float Zvertex)  {

//---Definitions
    const float PI    = 3.1415927;
    //UNUSED const float TWOPI = 2.0*PI;

//---Definitions for ECAL
    const float R_ECAL           = 136.5;
    const float Z_Endcap         = 328.0;
    const float etaBarrelEndcap  = 1.479;

//---ETA correction

    float Theta = 0.0  ;
        float ZEcal = R_ECAL*sinh(EtaParticle)+Zvertex;

    if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
    if(Theta<0.0) Theta = Theta+PI ;
    float ETA = - log(tan(0.5*Theta));

    if( fabs(ETA) > etaBarrelEndcap )
      {
       float Zend = Z_Endcap ;
       if(EtaParticle<0.0 )  Zend = -Zend ;
       float Zlen = Zend - Zvertex ;
       float RR = Zlen/sinh(EtaParticle);
       Theta = atan(RR/Zend);
       if(Theta<0.0) Theta = Theta+PI ;
       ETA = - log(tan(0.5*Theta));
      }
//---Return the result
        return ETA;
//---end
}






void PhotonsWithConversionsAnalyzer::analyze( const edm::Event& e, const edm::EventSetup& )
{


  using namespace edm;
  const float etaPhiDistance=0.01;
  // Fiducial region
  //UNUSED const float TRK_BARL =0.9;
  //UNUSED const float BARL = 1.4442; // DAQ TDR p.290
  //UNUSED const float END_LO = 1.566;
  //UNUSED const float END_HI = 2.5;
 // Electron mass
  //UNUSED const Float_t mElec= 0.000511;


  nEvt_++;
  LogInfo("ConvertedPhotonAnalyzer") << "ConvertedPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
  //  LogDebug("ConvertedPhotonAnalyzer") << "ConvertedPhotonAnalyzer Analyzing event number: "  << e.id() << " Global Counter " << nEvt_ <<"\n";
  std::cout << "ConvertedPhotonAnalyzer Analyzing event number: "  << e.id() << " Global Counter " << nEvt_ <<"\n";


  Handle<reco::PhotonCollection> photonHandle;
  e.getByLabel(photonCollectionProducer_, photonCollection_ , photonHandle);
  const reco::PhotonCollection photonCollection = *(photonHandle.product());
  std::cout  << "ConvertedPhotonAnalyzer  Photons with conversions collection size " << photonCollection.size() << "\n";


  std::cout  << " ConvertedPhotonAnalyzer Looking for MC truth " << "\n";

  //get simtrack info
  std::vector<SimTrack> theSimTracks;
  std::vector<SimVertex> theSimVertices;

  edm::Handle<SimTrackContainer> SimTk;
  edm::Handle<SimVertexContainer> SimVtx;
  e.getByLabel("g4SimHits",SimTk);
  e.getByLabel("g4SimHits",SimVtx);

  theSimTracks.insert(theSimTracks.end(),SimTk->begin(),SimTk->end());
  theSimVertices.insert(theSimVertices.end(),SimVtx->begin(),SimVtx->end());

  std::vector<PhotonMCTruth> mcPhotons=thePhotonMCTruthFinder_->find (theSimTracks,  theSimVertices);
  std::cout << " ConvertedPhotonAnalyzer mcPhotons size " <<  mcPhotons.size() << std::endl;



  // Loop over simulated photons
  //UNUSED int iDet=0;
  //UNUSED int iRadius=-1;
  //UNUSED int indPho=0;

  for ( std::vector<PhotonMCTruth>::const_iterator mcPho=mcPhotons.begin(); mcPho !=mcPhotons.end(); mcPho++) {
    float mcPhi= (*mcPho).fourMomentum().phi();
    float mcEta= (*mcPho).fourMomentum().pseudoRapidity();
    mcEta = etaTransformation(mcEta, (*mcPho).primaryVertex().z() );

    if ( (*mcPho).fourMomentum().et() < 20 ) continue;
    //    if ( ! (  fabs(mcEta) <= BARL || ( fabs(mcEta) >= END_LO && fabs(mcEta) <=END_HI ) ) ) {
    //     continue;
    //} // all ecal fiducial region



    h_MCphoE_->Fill(  (*mcPho).fourMomentum().e());
    h_MCphoEta_->Fill( (*mcPho).fourMomentum().eta());
    h_MCphoPhi_->Fill(  (*mcPho).fourMomentum().phi());



    // keep only visible conversions
    if (  (*mcPho).isAConversion() == 0 ) continue;


    nMCPho_++;

    h_MCConvEta_ ->Fill ( fabs( (*mcPho).fourMomentum().pseudoRapidity()) - 0.001);


    bool REJECTED;


    //std::cout   << " ConvertedPhotonAnalyzer  Starting loop over photon candidates " << "\n";
    for( reco::PhotonCollection::const_iterator  iPho = photonCollection.begin(); iPho != photonCollection.end(); iPho++) {
      REJECTED=false;

      //      std::cout  << " ConvertedPhotonAnalyzer Reco SC energy " << (*iPho).superCluster()->energy() <<  "\n";

      float phiClu=(*iPho).superCluster()->phi();
      float etaClu=(*iPho).superCluster()->eta();
      float deltaPhi = phiClu-mcPhi;
      float deltaEta = etaClu-mcEta;


      if ( deltaPhi > Geom::pi()  ) deltaPhi -= Geom::twoPi();
      if ( deltaPhi < -Geom::pi() ) deltaPhi += Geom::twoPi();
      deltaPhi=pow(deltaPhi,2);
      deltaEta=pow(deltaEta,2);
      float delta =  deltaPhi+deltaEta ;
      if (  delta >= etaPhiDistance  )  REJECTED=true;


      //      if ( ! (  fabs(etaClu) <= BARL || ( fabs(etaClu) >= END_LO && fabs(etaClu) <=END_HI ) ) ) REJECTED=true;

      if ( REJECTED ) continue;
      std::cout << " MATCHED " << std::endl;
      nMatched_++;


      //      std::cout << " ConvertedPhotonAnalyzer Matching candidate " << std::endl;

      // std::cout << " ConvertedPhotonAnalyzer Photons isAconversion " << (*mcPho).isAConversion() << " mcMatchingPhoton energy " <<  (*mcPho).fourMomentum().e()  << " ConvertedPhotonAnalyzer conversion vertex R " <<  (*mcPho).vertex().perp() << " Z " <<  (*mcPho).vertex().z() <<  std::endl;


      h_ErecoEMC_->Fill(   (*iPho).superCluster()->energy()/(*mcPho).fourMomentum().e());
      h_deltaPhi_-> Fill ( (*iPho).superCluster()->position().phi()- mcPhi);
      h_deltaEta_-> Fill ( (*iPho).superCluster()->position().eta()- mcEta);

      h_scE_->Fill( (*iPho).superCluster()->energy() );
      h_scEt_->Fill( (*iPho).superCluster()->energy()/cosh( (*iPho).superCluster()->position().eta()) );
      h_scEta_->Fill( (*iPho).superCluster()->position().eta() );
      h_scPhi_->Fill( (*iPho).superCluster()->position().phi() );


      h_phoE_->Fill( (*iPho).energy() );
      h_phoEta_->Fill( (*iPho).eta() );
      h_phoPhi_->Fill( (*iPho).phi() );

      if ( !(*iPho).hasConversionTracks() ) continue;
      //   std::cout << " This photons has " << (*iPho).conversions().size() << " conversions candidates " << std::endl;
      reco::ConversionRefVector conversions = (*iPho).conversions();
      //std::vector<reco::ConversionRef> conversions = (*iPho).conversions();


      for (unsigned int i=0; i<conversions.size(); i++) {
    //std::cout << " Conversion candidate Energy " << (*iPho).energy() << " number of tracks " << conversions[i]->nTracks() << std::endl;
        std::vector< edm::RefToBase<reco::Track> > tracks = conversions[i]->tracks();

    for (unsigned int i=0; i<tracks.size(); i++) {

      //      std::cout  << " ConvertedPhotonAnalyzer Reco Track charge " <<  tracks[i]->charge() << "  Num of RecHits " << tracks[i]->recHitsSize() << " inner momentum " <<  sqrt ( tracks[i]->innerMomentum().Mag2() )  <<  "\n";


      h2_tk_nHitsVsR_ -> Fill (  (*mcPho).vertex().perp(), tracks[i]->recHitsSize()   );
      h2_tk_inPtVsR_->Fill  (  (*mcPho).vertex().perp(),  sqrt( tracks[i]->innerMomentum().Mag2() ) );
    }

      }


    }

  }


}




void PhotonsWithConversionsAnalyzer::endJob()
{




  //   fOutputFile_->Write() ;
  // fOutputFile_->Close() ;

   edm::LogInfo("ConvertedPhotonAnalyzer") << "Analyzed " << nEvt_  << "\n";
   // std::cout  << "::endJob Analyzed " << nEvt_ << " events " << " with total " << nPho_ << " Photons " << "\n";
   std::cout  << "ConvertedPhotonAnalyzer::endJob Analyzed " << nEvt_ << " events " << "\n";

   return ;
}