Morgan/src/TotoAnalyzer.cc

#include "UserCode/Morgan/interface/TotoAnalyzer.h"

#include "FWCore/Framework/interface/TriggerNames.h"
#include "DataFormats/Common/interface/TriggerResults.h"

using namespace std;
using namespace reco;
using namespace edm;


TotoAnalyzer::TotoAnalyzer(const edm::ParameterSet& iConfig)
{
        myConfig = iConfig.getParameter<ParameterSet>("myConfig");
}


TotoAnalyzer::~TotoAnalyzer()
{}


// ------------ method called once each job just before starting event loop  ------------
void TotoAnalyzer::beginJob(const edm::EventSetup&)
{

        // Load Config parameters       
        verbosity = myConfig.getUntrackedParameter<int>("verbosity", 0);
        rootFileName_ = myConfig.getUntrackedParameter<string>("RootFileName","noname.root");
        isCSA07Soup = myConfig.getUntrackedParameter<bool>("isCSA07Soup",false);
        doHLT = myConfig.getUntrackedParameter<bool>("doHLT",false);
        doMC = myConfig.getUntrackedParameter<bool>("doMC",false);
        doSignalMC = myConfig.getUntrackedParameter<bool>("doSignalMC",false);
        doTrack = myConfig.getUntrackedParameter<bool>("doTrack",false);
        doJet = myConfig.getUntrackedParameter<bool>("doJet",false);
        doMuon = myConfig.getUntrackedParameter<bool>("doMuon",false);
        doElectron = myConfig.getUntrackedParameter<bool>("doElectron",false);  
        doPhoton = myConfig.getUntrackedParameter<bool>("doPhoton",false);
        doCluster = myConfig.getUntrackedParameter<bool>("doCluster",false);
        doPhotonIsolation = myConfig.getUntrackedParameter<bool>("doPhotonIsolation",false);
        doPhotonConversion = myConfig.getUntrackedParameter<bool>("doPhotonConversion",false);
        
        nTotEvt_ = 0;
        
        // initialize root output file
        if(verbosity>0) cout << "New RootFile " << rootFileName_.c_str() << " is created" << endl;
        rootFile_ = new TFile(rootFileName_.c_str(), "recreate");
        rootFile_->cd();

        runInfos_ = new TRootRun();
        runTree_ = new TTree("runTree", "Global Run Infos");
        runTree_->Branch ("runInfos", "TRootRun", &runInfos_);

        rootEvent = 0;
        eventTree_ = new TTree("eventTree", "Event Infos");
        eventTree_->Branch ("Event", "TRootEvent", &rootEvent);

        if(doHLT)
        {
                hltAnalyzer_ = new HLTAnalyzer(&myConfig);
                hltAnalyzer_->setVerbosity(verbosity);
        }

        if(doSignalMC)
        {                       
                if(verbosity>0) cout << "Signal generator info will be added to rootuple" << endl;
                rootMCSignalEvent = 0;
                cout << "Create MuMuGamma Branch" << endl;
                eventTree_->Branch ("MuMuGamma", "TRootSignalEvent", &rootMCSignalEvent);
        }

        if(doMC)
        {
                if(verbosity>0) cout << "MCParticles info will be added to rootuple" << endl;
                mcParticles = new TClonesArray("TRootParticle", 1000);
                eventTree_->Branch ("MCParticles", "TClonesArray", &mcParticles);
        }
        
        if(doTrack)
        {
                if(verbosity>0) cout << "Tracks info will be added to rootuple" << endl;
                tracks = new TClonesArray("TRootTrack", 1000);
                eventTree_->Branch ("Tracks", "TClonesArray", &tracks);
        }

        if(doJet)
        {
                if(verbosity>0) cout << "Jets info will be added to rootuple" << endl;
                jets = new TClonesArray("TRootJet", 1000);
                eventTree_->Branch ("Jets", "TClonesArray", &jets);
        }

        if(doMuon)
        {
                if(verbosity>0) cout << "Muons info will be added to rootuple" << endl;
                muons = new TClonesArray("TRootMuon", 1000);
                eventTree_->Branch ("Muons", "TClonesArray", &muons);
        }
        
        if(doElectron)
        {
                if(verbosity>0) cout << "Electrons info will be added to rootuple" << endl;
                electrons = new TClonesArray("TRootElectron", 1000);
                eventTree_->Branch ("Electrons", "TClonesArray", &electrons);
        }

        if(doPhoton)
        {
                if(verbosity>0) cout << "Photons info will be added to rootuple" << endl;
                uncorrectedPhotons = new TClonesArray("TRootPhoton", 1000);
                eventTree_->Branch ("uncorrectedPhotons", "TClonesArray", &uncorrectedPhotons);
                correctedPhotons = new TClonesArray("TRootPhoton", 1000);
                eventTree_->Branch ("correctedPhotons", "TClonesArray", &correctedPhotons);
        }
        
        if(doCluster)
        {
                if(verbosity>0) cout << "ECAL Clusters info will be added to rootuple" << endl;
                clusters = new TClonesArray("TRootCluster", 1000);
                eventTree_->Branch ("BasicClusters", "TClonesArray", &clusters);
                superClusters = new TClonesArray("TRootSuperCluster", 1000);
                eventTree_->Branch ("SuperClusters", "TClonesArray", &superClusters);
        }

        if(doPhotonConversion)
        {
                if(verbosity>0) cout << "Conversion Tracks info will be added to rootuple" << endl;
                conversionTracks = new TClonesArray("TRootParticle", 1000);
                eventTree_->Branch ("ConversionTracks", "TClonesArray", &conversionTracks);
        }
}


// ------------ method called once each job just after ending the event loop  ------------
void TotoAnalyzer::endJob()
{
        hltAnalyzer_ ->copySummary(runInfos_);
        runTree_->Fill();

        std::cout << "Total number of events: " << nTotEvt_ << std::endl;
        std::cout << "Closing rootfile " << rootFile_->GetName() << std::endl;
        rootFile_->Write();
        rootFile_->Close();

        // Trigger Summary Tables
        if(doHLT)
        {       
                cout << "Trigger Summary Tables" << endl;
                hltAnalyzer_->printStats();
        }

        delete mcParticles;
        delete tracks;
        delete jets;
        delete muons;
        delete electrons;       
        delete correctedPhotons;
        delete uncorrectedPhotons;
        delete clusters;        
        delete superClusters;   
        delete conversionTracks;        
}


// ------------ method called to for each event  ------------
void TotoAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{

        rootFile_->cd();
        nTotEvt_++;
        if( (verbosity>1) || (verbosity>0 && nTotEvt_%10==0 && nTotEvt_<=100)  || (verbosity>0 && nTotEvt_%100==0 && nTotEvt_>100) )
                cout << endl << endl << "####### TotoAnalyzer - Event " << nTotEvt_ << " #######" << endl;

        
        // Global Event Infos
        rootEvent = new TRootEvent();
        rootEvent->setNb(nTotEvt_);

        
        // CSA07 Process Id and Event Weight
        if(isCSA07Soup)
        {
                edm::Handle< double> weightHandle;
                iEvent.getByLabel ("csaweightproducer","weight", weightHandle);
                double weight = * weightHandle;
                int processId = csa07::csa07ProcessId(iEvent, 1000., "csaweightproducer");
                char * pname= csa07::csa07ProcessName(processId);
                if(verbosity>1) cout << "CSA07 Soup - id=" << processId << " - "<< pname << " - weight=" << weight << endl;
                rootEvent->setCsa07id(processId);
                rootEvent->setCsa07weight(weight);
                rootEvent->setCsa07process(pname);      
        }

        
        // Trigger
        if(doHLT)
        {       
                cout << "Get TriggerResults..." << endl;
                if (nTotEvt_==1) hltAnalyzer_->init(iEvent, rootEvent);
                hltAnalyzer_->process(iEvent, rootEvent);
        }


        // MC Info
        if(doMC)
        {
                rootMCSignalEvent = new TRootSignalEvent();
                if(verbosity>1) cout << "Analysing MC info..." << endl;
                myMCAnalyzer = new MCAnalyzer(&myConfig);
                //myMCAnalyzer->SetVerbosity(1);
                myMCAnalyzer->Process(iEvent, iSetup, rootEvent, rootMCSignalEvent, mcParticles);
                delete myMCAnalyzer;
        }
        
        
        // Get Primary Vertices
        if(verbosity>1) cout << "Analysing primary vertices collection..." << endl;
        edm::Handle< reco::VertexCollection > recoVertex;
        iEvent.getByLabel("offlinePrimaryVerticesFromCTFTracks", recoVertex);
        if(verbosity>1) std::cout << "Number of primary vertices in the event: " << recoVertex->size() << std::endl;
        for (unsigned int j=0; j<recoVertex->size(); j++)
        {
                if(verbosity>1) cout << "offlinePrimaryVerticesFromCTFTracks  vx=" << (*recoVertex)[j].x() << " vy=" << (*recoVertex)[j].y() << " vz=" << (*recoVertex)[j].z() << endl;
                rootEvent->addPrimaryVertex( TVector3( (*recoVertex)[j].x(), (*recoVertex)[j].y(), (*recoVertex)[j].z() ) );
        }


        // Tracks
        if(doTrack)
        {
                if(verbosity>1) cout << "Analysing tracks collection..." << endl;
                edm::Handle<reco::TrackCollection> recoTracks;
                iEvent.getByLabel("ctfWithMaterialTracks", recoTracks);
                if(verbosity>1) std::cout << "Number of tracks in the event: " << recoTracks->size() << std::endl;

                for (unsigned int j=0; j<recoTracks->size(); j++)
                {
                        const reco::HitPattern& hit = (*recoTracks)[j].hitPattern();
                        
                        //TRootTrack localTrack( (*recoTracks)[j].px(), (*recoTracks)[j].py(), (*recoTracks)[j].pz(), (*recoTracks)[j].p(), (*recoTracks)[j].vx(), (*recoTracks)[j].vy(), (*recoTracks)[j].vz(), 0, (*recoTracks)[j].charge() );

                        TRootTrack localTrack( 
                                        (*recoTracks)[j].px()
                                        ,(*recoTracks)[j].py()
                                        ,(*recoTracks)[j].pz()
                                        ,(*recoTracks)[j].p()
                                        ,(*recoTracks)[j].vx()
                                        ,(*recoTracks)[j].vy()
                                        ,(*recoTracks)[j].vz()
                                        ,0
                                        ,(*recoTracks)[j].charge()
                                        ,hit.numberOfValidPixelHits()
                                        ,hit.numberOfValidTrackerHits()
                                        ,(*recoTracks)[j].chi2()
                                        ,(*recoTracks)[j].d0()
                                        ,(*recoTracks)[j].d0Error()
                                        ,(*recoTracks)[j].dz()
                                        ,(*recoTracks)[j].dzError()
                        );

                        new( (*tracks)[j] ) TRootTrack(localTrack); 
                        if(verbosity>2) cout << "Track - "<< localTrack << endl;

                }               
        }
        
        
        // Calo Jet
        if(doJet)
        {
                if(verbosity>1) cout << "Analysing jets collection..." << endl;
                edm::Handle< reco::CaloJetCollection > recoJets;
                iEvent.getByLabel("iterativeCone5CaloJets", recoJets);
                if(verbosity>1) std::cout << "Number of jets in the event: " << recoJets->size() << std::endl;
                for (unsigned int j=0; j<recoJets->size(); j++)
                {
                        if(verbosity>2) cout << "jet pt="<< (*recoJets)[j].pt() << endl;
                        new( (*jets)[j] ) TRootJet( (*recoJets)[j].px(), (*recoJets)[j].py(), (*recoJets)[j].pz(), (*recoJets)[j].energy() );
                }
        }

        
        // Muons
        if(doMuon)
        {
                Float_t sintheta = 0.;
                if(verbosity>1) cout << "Analysing muons collection..." << endl;
                edm::Handle < reco::MuonCollection > recoMuons;
                iEvent.getByLabel ("muons", recoMuons);
                if(verbosity>1) std::cout << "Number of muons in the event: " << recoMuons->size() << std::endl;
                for (unsigned int j=0; j<recoMuons->size(); j++)
                {
                        TRootMuon localMuon( (*recoMuons)[j].px(), (*recoMuons)[j].py(), (*recoMuons)[j].pz(), (*recoMuons)[j].energy(), (*recoMuons)[j].vx(), (*recoMuons)[j].vy(), (*recoMuons)[j].vz(), 13, (*recoMuons)[j].charge() ); 
                        sintheta = sin( localMuon.Theta() );
                        localMuon.setCaloEnergy( (*recoMuons)[j].getCalEnergy().em * sintheta, (*recoMuons)[j].getCalEnergy().emS9 * sintheta,
                                                            (*recoMuons)[j].getCalEnergy().had * sintheta, (*recoMuons)[j].getCalEnergy().hadS9 * sintheta,
                                                            (*recoMuons)[j].getCalEnergy().ho * sintheta, (*recoMuons)[j].getCalEnergy().hoS9 * sintheta, (*recoMuons)[j].getCaloCompatibility() );
                        localMuon.setIsoR03( (*recoMuons)[j].getIsolationR03().emEt, (*recoMuons)[j].getIsolationR03().hadEt, (*recoMuons)[j].getIsolationR03().hoEt,
                                                        (*recoMuons)[j].getIsolationR03().sumPt, (*recoMuons)[j].getIsolationR03().nTracks, (*recoMuons)[j].getIsolationR03().nJets);
                        localMuon.setIsoR05( (*recoMuons)[j].getIsolationR05().emEt, (*recoMuons)[j].getIsolationR05().hadEt, (*recoMuons)[j].getIsolationR05().hoEt,
                                                        (*recoMuons)[j].getIsolationR05().sumPt, (*recoMuons)[j].getIsolationR05().nTracks, (*recoMuons)[j].getIsolationR05().nJets);
                        new( (*muons)[j] ) TRootMuon(localMuon);
                        if(verbosity>2) cout << localMuon << endl;
                }
        }
                        
                        
        // Electrons
        if(doElectron)
        {
                if(verbosity>1) cout << "Analysing electrons collection..." << endl;
                edm::Handle< reco::PixelMatchGsfElectronCollection > recoElectrons;
                iEvent.getByLabel("pixelMatchGsfElectrons", recoElectrons);
                if(verbosity>1) std::cout << "Number of electrons in the event: " << recoElectrons->size() << std::endl;
                for (unsigned int j=0; j<recoElectrons->size(); j++)
                {
                        float elec_e = ((*recoElectrons)[j].superCluster())->energy();
                        float elec_px = (*recoElectrons)[j].px() * elec_e / (*recoElectrons)[j].p();
                        float elec_py = (*recoElectrons)[j].py() * elec_e / (*recoElectrons)[j].p();
                        float elec_pz = (*recoElectrons)[j].pz() * elec_e / (*recoElectrons)[j].p();
                        TRootElectron localElectron( elec_px, elec_py, elec_pz, elec_e, (*recoElectrons)[j].vx(), (*recoElectrons)[j].vy(), (*recoElectrons)[j].vz(), 11, (*recoElectrons)[j].charge() );
                        new( (*electrons)[j] ) TRootElectron(localElectron);
                        if(verbosity>2) cout << localElectron << endl;
                }
        }
                        
                
        // Photons
        if(doPhoton)
        {
                if(verbosity>1) cout << "Analysing photons collection..." << endl;
                myPhotonAnalyzer = new PhotonAnalyzer();
                myPhotonAnalyzer->SetVerbosity(verbosity);
                myPhotonAnalyzer->Process(iEvent, rootEvent, uncorrectedPhotons, "photons");
                myPhotonAnalyzer->Process(iEvent, rootEvent, correctedPhotons, "correctedPhotons");
                delete myPhotonAnalyzer;
        }
                
                
        // Ecal Clusters
        if(doCluster)
        {
                if(verbosity>1) cout << "Analysing BasicClusters collection..." << endl;
                myClusterAnalyzer = new ClusterAnalyzer();
                myClusterAnalyzer->SetVerbosity(verbosity);
                myClusterAnalyzer->Process(iEvent, rootEvent, clusters, "islandBasicClusters","islandBarrelBasicClusters", "islandBarrelShape", 110);
                myClusterAnalyzer->Process(iEvent, rootEvent, clusters, "islandBasicClusters", "islandEndcapBasicClusters", "islandEndcapShape", 120);
                myClusterAnalyzer->Process(iEvent, rootEvent, clusters, "hybridSuperClusters", "", "", 210);
                delete myClusterAnalyzer;
                
                if(verbosity>1) cout << "Analysing SuperClusters collection..." << endl;
                mySClusterAnalyzer = new SuperClusterAnalyzer();
                mySClusterAnalyzer->SetVerbosity(verbosity);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "islandSuperClusters","islandBarrelSuperClusters",110);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "islandSuperClusters","islandEndcapSuperClusters",120);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedIslandEndcapSuperClusters","",121);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedEndcapSuperClustersWithPreshower","",122);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "hybridSuperClusters","",210);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedHybridSuperClusters","",211);
                delete mySClusterAnalyzer;
        }
        
        /*
        cout << "Nb of BC110=" << rootEvent->nBasicClusters(110) << endl;
        cout << "Nb of BC120"= << rootEvent->nBasicClusters(120) << endl;
        cout << "Nb of BC210=" << rootEvent->nBasicClusters(210) << endl; );
        if(verbosity>2) cout << "electron et="<< (*recoElectrons)[j].pt() <
        cout << "Nb of SC110=" << rootEvent->nSuperClusters(110) << endl;
        cout << "Nb of SC120=" << rootEvent->nSuperClusters(120) << endl;
        cout << "Nb of SC121=" << rootEvent->nSuperClusters(121) << endl;Misc
        cout << "Nb of SC122=" << rootEvent->nSuperClusters(122) << endl;
        cout << "Nb of SC210=" << rootEvent->nSuperClusters(210) << endl;
        cout << "Nb of SC211=" << rootEvent->nSuperClusters(211) << endl;
        */
        
        
        if(doCluster)
        {
                ClusterAssociator* myClusterAssociator = new ClusterAssociator();
                myClusterAssociator->SetVerbosity(verbosity);
                myClusterAssociator->Process(superClusters, clusters);
                //myClusterAssociator->PrintSuperClusters(superClusters, clusters,0);  // 0 to print all types SC
                delete myClusterAssociator;
        }
        
        
        if(doPhoton && doCluster)
        {
                PhotonAssociator* myPhotonAssociator = new PhotonAssociator();
                myPhotonAssociator->SetVerbosity(verbosity);
                myPhotonAssociator->AssociateSuperCluster(correctedPhotons, superClusters);
                //myPhotonAssociator->AssociateSuperCluster(uncorrectedPhotons, superClusters);
                if(doPhotonConversion) myPhotonAssociator->AssociateConversion(iEvent, correctedPhotons, superClusters, conversionTracks);
                //if(doPhotonConversion) myPhotonAssociator->AssociateConversion(iEvent, uncorrectedPhotons, superClusters, conversionTracks);
                //myPhotonAssociator->PrintPhotons(correctedPhotons,superClusters,0);  // 0 to print all types associated SC
                myPhotonAssociator->FullPrintPhotons(correctedPhotons,superClusters,clusters,0);  // 0 to print all types associated SC
                delete myPhotonAssociator;
        }
        
        
        if(doPhoton && doPhotonIsolation) //FIXME - Test HCAL rechits availability
        {
                if(verbosity>1) cout << "PhotonIsolator......" << endl;
                
                // Get HB/HE rechits
                edm::Handle<HBHERecHitCollection> hbheHandle;
                iEvent.getByLabel("hbhereco", hbheHandle);
                const HBHERecHitCollection* hbheHits = hbheHandle.product();
  
                // Get HO rechits
                edm::Handle<HORecHitCollection> hoHandle;
                iEvent.getByLabel("horeco", hoHandle);
                const HORecHitCollection* hoHits = hoHandle.product();

                // Get HF rechits
                edm::Handle<HFRecHitCollection> hfHandle;
                iEvent.getByLabel("hfreco", hfHandle);
                const HFRecHitCollection* hfHits = hfHandle.product();

                // Get Calo Geometry
                edm::ESHandle<CaloGeometry> geomHandle;
                (iSetup.get<IdealGeometryRecord>()).get(geomHandle);
                const CaloGeometry* caloGeom = geomHandle.product();

                Float_t ecalIslandIsolation;
                Float_t ecalDoubleConeIsolation;
                Float_t hcalRecHitIsolation;
                Float_t isoTracks;
                Int_t isoNTracks;
                TRootPhoton* localPhoton;

                PhotonIsolator* myPhotonIsolator = new PhotonIsolator(&myConfig);
                for (int iphoton=0; iphoton<correctedPhotons->GetEntriesFast(); iphoton++)
                {
                        localPhoton = (TRootPhoton*)correctedPhotons->At(iphoton);

                        ecalIslandIsolation=-1.;
                        if (doCluster) ecalIslandIsolation = myPhotonIsolator->EcalIslandIsolation(localPhoton, superClusters, clusters);

                        ecalDoubleConeIsolation=-1.;
                        if (doCluster) ecalDoubleConeIsolation = myPhotonIsolator->EcalDoubleConeIsolation(localPhoton, superClusters, clusters);

                        hcalRecHitIsolation=-1.;
                        hcalRecHitIsolation = myPhotonIsolator->HcalRecHitIsolation(localPhoton, hbheHits, hoHits, hfHits, caloGeom);

                        isoTracks=-1.;
                        isoNTracks=-1;
                        if(doTrack) isoTracks = myPhotonIsolator->TrackerIsolation(localPhoton, tracks, isoNTracks);

                        localPhoton->setIsolation(ecalIslandIsolation, ecalDoubleConeIsolation, hcalRecHitIsolation, isoTracks, isoNTracks);
                        if(verbosity>2) cout << "PhotonIsolator for [" << iphoton << "] " << *localPhoton << endl;
                }
                delete myPhotonIsolator;
        }
        
        
        if(verbosity>2) cout << "Filling rootuple..." << endl;
        eventTree_->Fill();
        if(verbosity>2) cout << "Deleting objects..." << endl;
        delete rootEvent;
        if(doSignalMC) delete rootMCSignalEvent;
        if(doMC) (*mcParticles).Delete();
        if(doTrack) (*tracks).Delete();
        if(doJet) (*jets).Delete();
        if(doMuon) (*muons).Delete();
        if(doElectron) (*electrons).Delete();
        if(doPhoton) (*uncorrectedPhotons).Delete();
        if(doPhoton) (*correctedPhotons).Delete();
        if(doCluster) (*clusters).Delete();
        if(doCluster) (*superClusters).Delete();
        if(doPhotonConversion) (*conversionTracks).Delete();
        
}


//define this as a plug-in
//DEFINE_ANOTHER_FWK_MODULE(TotoAnalyzer);
Revision:	1.6
Committed:	Mon Jun 23 13:31:11 2008 UTC (16 years, 10 months ago) by mlethuil
Content type:	text/plain
Branch:	MAIN
Changes since 1.5:	+16 -0 lines
Log Message:	CSA07 process id and event weight
#	User	Rev	Content
1	mlethuil	1.1	#include "UserCode/Morgan/interface/TotoAnalyzer.h"
2
3	mlethuil	1.5	#include "FWCore/Framework/interface/TriggerNames.h"
4			#include "DataFormats/Common/interface/TriggerResults.h"
5
6	mlethuil	1.1	using namespace std;
7			using namespace reco;
8			using namespace edm;
9
10
11			TotoAnalyzer::TotoAnalyzer(const edm::ParameterSet& iConfig)
12			{
13			myConfig = iConfig.getParameter<ParameterSet>("myConfig");
14			}
15
16
17			TotoAnalyzer::~TotoAnalyzer()
18			{}
19
20
21
22			// ------------ method called once each job just before starting event loop ------------
23			void TotoAnalyzer::beginJob(const edm::EventSetup&)
24			{
25
26			// Load Config parameters
27			verbosity = myConfig.getUntrackedParameter<int>("verbosity", 0);
28			rootFileName_ = myConfig.getUntrackedParameter<string>("RootFileName","noname.root");
29	mlethuil	1.6	isCSA07Soup = myConfig.getUntrackedParameter<bool>("isCSA07Soup",false);
30	mlethuil	1.5	doHLT = myConfig.getUntrackedParameter<bool>("doHLT",false);
31	mlethuil	1.1	doMC = myConfig.getUntrackedParameter<bool>("doMC",false);
32			doSignalMC = myConfig.getUntrackedParameter<bool>("doSignalMC",false);
33			doTrack = myConfig.getUntrackedParameter<bool>("doTrack",false);
34			doJet = myConfig.getUntrackedParameter<bool>("doJet",false);
35			doMuon = myConfig.getUntrackedParameter<bool>("doMuon",false);
36			doElectron = myConfig.getUntrackedParameter<bool>("doElectron",false);
37			doPhoton = myConfig.getUntrackedParameter<bool>("doPhoton",false);
38			doCluster = myConfig.getUntrackedParameter<bool>("doCluster",false);
39	mlethuil	1.4	doPhotonIsolation = myConfig.getUntrackedParameter<bool>("doPhotonIsolation",false);
40	mlethuil	1.3	doPhotonConversion = myConfig.getUntrackedParameter<bool>("doPhotonConversion",false);
41	mlethuil	1.1
42			nTotEvt_ = 0;
43
44			// initialize root output file
45			if(verbosity>0) cout << "New RootFile " << rootFileName_.c_str() << " is created" << endl;
46			rootFile_ = new TFile(rootFileName_.c_str(), "recreate");
47			rootFile_->cd();
48	mlethuil	1.5
49			runInfos_ = new TRootRun();
50			runTree_ = new TTree("runTree", "Global Run Infos");
51			runTree_->Branch ("runInfos", "TRootRun", &runInfos_);
52
53	mlethuil	1.1	rootEvent = 0;
54	mlethuil	1.5	eventTree_ = new TTree("eventTree", "Event Infos");
55			eventTree_->Branch ("Event", "TRootEvent", &rootEvent);
56
57			if(doHLT)
58			{
59			hltAnalyzer_ = new HLTAnalyzer(&myConfig);
60			hltAnalyzer_->setVerbosity(verbosity);
61			}
62
63	mlethuil	1.1	if(doSignalMC)
64			{
65			if(verbosity>0) cout << "Signal generator info will be added to rootuple" << endl;
66			rootMCSignalEvent = 0;
67			cout << "Create MuMuGamma Branch" << endl;
68	mlethuil	1.5	eventTree_->Branch ("MuMuGamma", "TRootSignalEvent", &rootMCSignalEvent);
69	mlethuil	1.1	}
70
71			if(doMC)
72			{
73			if(verbosity>0) cout << "MCParticles info will be added to rootuple" << endl;
74			mcParticles = new TClonesArray("TRootParticle", 1000);
75	mlethuil	1.5	eventTree_->Branch ("MCParticles", "TClonesArray", &mcParticles);
76	mlethuil	1.1	}
77
78			if(doTrack)
79			{
80			if(verbosity>0) cout << "Tracks info will be added to rootuple" << endl;
81	mlethuil	1.4	tracks = new TClonesArray("TRootTrack", 1000);
82	mlethuil	1.5	eventTree_->Branch ("Tracks", "TClonesArray", &tracks);
83	mlethuil	1.1	}
84
85			if(doJet)
86			{
87			if(verbosity>0) cout << "Jets info will be added to rootuple" << endl;
88			jets = new TClonesArray("TRootJet", 1000);
89	mlethuil	1.5	eventTree_->Branch ("Jets", "TClonesArray", &jets);
90	mlethuil	1.1	}
91
92			if(doMuon)
93			{
94			if(verbosity>0) cout << "Muons info will be added to rootuple" << endl;
95			muons = new TClonesArray("TRootMuon", 1000);
96	mlethuil	1.5	eventTree_->Branch ("Muons", "TClonesArray", &muons);
97	mlethuil	1.1	}
98
99			if(doElectron)
100			{
101			if(verbosity>0) cout << "Electrons info will be added to rootuple" << endl;
102			electrons = new TClonesArray("TRootElectron", 1000);
103	mlethuil	1.5	eventTree_->Branch ("Electrons", "TClonesArray", &electrons);
104	mlethuil	1.1	}
105
106			if(doPhoton)
107			{
108			if(verbosity>0) cout << "Photons info will be added to rootuple" << endl;
109			uncorrectedPhotons = new TClonesArray("TRootPhoton", 1000);
110	mlethuil	1.5	eventTree_->Branch ("uncorrectedPhotons", "TClonesArray", &uncorrectedPhotons);
111	mlethuil	1.1	correctedPhotons = new TClonesArray("TRootPhoton", 1000);
112	mlethuil	1.5	eventTree_->Branch ("correctedPhotons", "TClonesArray", &correctedPhotons);
113	mlethuil	1.1	}
114
115			if(doCluster)
116			{
117			if(verbosity>0) cout << "ECAL Clusters info will be added to rootuple" << endl;
118			clusters = new TClonesArray("TRootCluster", 1000);
119	mlethuil	1.5	eventTree_->Branch ("BasicClusters", "TClonesArray", &clusters);
120	mlethuil	1.1	superClusters = new TClonesArray("TRootSuperCluster", 1000);
121	mlethuil	1.5	eventTree_->Branch ("SuperClusters", "TClonesArray", &superClusters);
122	mlethuil	1.1	}
123	mlethuil	1.3
124			if(doPhotonConversion)
125			{
126			if(verbosity>0) cout << "Conversion Tracks info will be added to rootuple" << endl;
127			conversionTracks = new TClonesArray("TRootParticle", 1000);
128	mlethuil	1.5	eventTree_->Branch ("ConversionTracks", "TClonesArray", &conversionTracks);
129	mlethuil	1.3	}
130	mlethuil	1.1	}
131
132
133			// ------------ method called once each job just after ending the event loop ------------
134			void TotoAnalyzer::endJob()
135			{
136	mlethuil	1.5	hltAnalyzer_ ->copySummary(runInfos_);
137			runTree_->Fill();
138
139	mlethuil	1.1	std::cout << "Total number of events: " << nTotEvt_ << std::endl;
140			std::cout << "Closing rootfile " << rootFile_->GetName() << std::endl;
141			rootFile_->Write();
142			rootFile_->Close();
143	mlethuil	1.5
144			// Trigger Summary Tables
145			if(doHLT)
146			{
147			cout << "Trigger Summary Tables" << endl;
148			hltAnalyzer_->printStats();
149			}
150
151	mlethuil	1.1	delete mcParticles;
152			delete tracks;
153			delete jets;
154			delete muons;
155			delete electrons;
156			delete correctedPhotons;
157			delete uncorrectedPhotons;
158			delete clusters;
159			delete superClusters;
160	mlethuil	1.3	delete conversionTracks;
161	mlethuil	1.1	}
162
163
164
165			// ------------ method called to for each event ------------
166			void TotoAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
167			{
168
169			rootFile_->cd();
170			nTotEvt_++;
171			if( (verbosity>1) \|\| (verbosity>0 && nTotEvt_%10==0 && nTotEvt_<=100) \|\| (verbosity>0 && nTotEvt_%100==0 && nTotEvt_>100) )
172			cout << endl << endl << "####### TotoAnalyzer - Event " << nTotEvt_ << " #######" << endl;
173
174
175			// Global Event Infos
176			rootEvent = new TRootEvent();
177			rootEvent->setNb(nTotEvt_);
178
179	mlethuil	1.6
180			// CSA07 Process Id and Event Weight
181			if(isCSA07Soup)
182			{
183			edm::Handle< double> weightHandle;
184			iEvent.getByLabel ("csaweightproducer","weight", weightHandle);
185			double weight = * weightHandle;
186			int processId = csa07::csa07ProcessId(iEvent, 1000., "csaweightproducer");
187			char * pname= csa07::csa07ProcessName(processId);
188			if(verbosity>1) cout << "CSA07 Soup - id=" << processId << " - "<< pname << " - weight=" << weight << endl;
189			rootEvent->setCsa07id(processId);
190			rootEvent->setCsa07weight(weight);
191			rootEvent->setCsa07process(pname);
192			}
193	mlethuil	1.5
194	mlethuil	1.6
195	mlethuil	1.5	// Trigger
196			if(doHLT)
197			{
198			cout << "Get TriggerResults..." << endl;
199			if (nTotEvt_==1) hltAnalyzer_->init(iEvent, rootEvent);
200			hltAnalyzer_->process(iEvent, rootEvent);
201			}
202
203
204	mlethuil	1.1	// MC Info
205			if(doMC)
206			{
207			rootMCSignalEvent = new TRootSignalEvent();
208			if(verbosity>1) cout << "Analysing MC info..." << endl;
209			myMCAnalyzer = new MCAnalyzer(&myConfig);
210			//myMCAnalyzer->SetVerbosity(1);
211			myMCAnalyzer->Process(iEvent, iSetup, rootEvent, rootMCSignalEvent, mcParticles);
212			delete myMCAnalyzer;
213			}
214
215
216			// Get Primary Vertices
217			if(verbosity>1) cout << "Analysing primary vertices collection..." << endl;
218			edm::Handle< reco::VertexCollection > recoVertex;
219			iEvent.getByLabel("offlinePrimaryVerticesFromCTFTracks", recoVertex);
220			if(verbosity>1) std::cout << "Number of primary vertices in the event: " << recoVertex->size() << std::endl;
221			for (unsigned int j=0; j<recoVertex->size(); j++)
222			{
223			if(verbosity>1) cout << "offlinePrimaryVerticesFromCTFTracks vx=" << (recoVertex)[j].x() << " vy=" << (recoVertex)[j].y() << " vz=" << (*recoVertex)[j].z() << endl;
224			rootEvent->addPrimaryVertex( TVector3( (recoVertex)[j].x(), (recoVertex)[j].y(), (*recoVertex)[j].z() ) );
225			}
226
227
228			// Tracks
229			if(doTrack)
230			{
231			if(verbosity>1) cout << "Analysing tracks collection..." << endl;
232			edm::Handle<reco::TrackCollection> recoTracks;
233			iEvent.getByLabel("ctfWithMaterialTracks", recoTracks);
234			if(verbosity>1) std::cout << "Number of tracks in the event: " << recoTracks->size() << std::endl;
235	mlethuil	1.4
236	mlethuil	1.1	for (unsigned int j=0; j<recoTracks->size(); j++)
237			{
238	mlethuil	1.4	const reco::HitPattern& hit = (*recoTracks)[j].hitPattern();
239
240			//TRootTrack localTrack( (recoTracks)[j].px(), (recoTracks)[j].py(), (recoTracks)[j].pz(), (recoTracks)[j].p(), (recoTracks)[j].vx(), (recoTracks)[j].vy(), (recoTracks)[j].vz(), 0, (recoTracks)[j].charge() );
241
242			TRootTrack localTrack(
243			(*recoTracks)[j].px()
244			,(*recoTracks)[j].py()
245			,(*recoTracks)[j].pz()
246			,(*recoTracks)[j].p()
247			,(*recoTracks)[j].vx()
248			,(*recoTracks)[j].vy()
249			,(*recoTracks)[j].vz()
250			,0
251			,(*recoTracks)[j].charge()
252			,hit.numberOfValidPixelHits()
253			,hit.numberOfValidTrackerHits()
254			,(*recoTracks)[j].chi2()
255			,(*recoTracks)[j].d0()
256			,(*recoTracks)[j].d0Error()
257			,(*recoTracks)[j].dz()
258			,(*recoTracks)[j].dzError()
259			);
260
261			new( (*tracks)[j] ) TRootTrack(localTrack);
262	mlethuil	1.1	if(verbosity>2) cout << "Track - "<< localTrack << endl;
263	mlethuil	1.4
264	mlethuil	1.1	}
265			}
266
267
268			// Calo Jet
269			if(doJet)
270			{
271			if(verbosity>1) cout << "Analysing jets collection..." << endl;
272			edm::Handle< reco::CaloJetCollection > recoJets;
273			iEvent.getByLabel("iterativeCone5CaloJets", recoJets);
274			if(verbosity>1) std::cout << "Number of jets in the event: " << recoJets->size() << std::endl;
275			for (unsigned int j=0; j<recoJets->size(); j++)
276			{
277			if(verbosity>2) cout << "jet pt="<< (*recoJets)[j].pt() << endl;
278			new( (jets)[j] ) TRootJet( (recoJets)[j].px(), (recoJets)[j].py(), (recoJets)[j].pz(), (*recoJets)[j].energy() );
279			}
280			}
281
282
283			// Muons
284			if(doMuon)
285			{
286			Float_t sintheta = 0.;
287			if(verbosity>1) cout << "Analysing muons collection..." << endl;
288			edm::Handle < reco::MuonCollection > recoMuons;
289			iEvent.getByLabel ("muons", recoMuons);
290			if(verbosity>1) std::cout << "Number of muons in the event: " << recoMuons->size() << std::endl;
291			for (unsigned int j=0; j<recoMuons->size(); j++)
292			{
293			TRootMuon localMuon( (recoMuons)[j].px(), (recoMuons)[j].py(), (recoMuons)[j].pz(), (recoMuons)[j].energy(), (recoMuons)[j].vx(), (recoMuons)[j].vy(), (recoMuons)[j].vz(), 13, (recoMuons)[j].charge() );
294			sintheta = sin( localMuon.Theta() );
295			localMuon.setCaloEnergy( (recoMuons)[j].getCalEnergy().em sintheta, (recoMuons)[j].getCalEnergy().emS9 sintheta,
296			(recoMuons)[j].getCalEnergy().had sintheta, (recoMuons)[j].getCalEnergy().hadS9 sintheta,
297			(recoMuons)[j].getCalEnergy().ho sintheta, (recoMuons)[j].getCalEnergy().hoS9 sintheta, (*recoMuons)[j].getCaloCompatibility() );
298			localMuon.setIsoR03( (recoMuons)[j].getIsolationR03().emEt, (recoMuons)[j].getIsolationR03().hadEt, (*recoMuons)[j].getIsolationR03().hoEt,
299			(recoMuons)[j].getIsolationR03().sumPt, (recoMuons)[j].getIsolationR03().nTracks, (*recoMuons)[j].getIsolationR03().nJets);
300			localMuon.setIsoR05( (recoMuons)[j].getIsolationR05().emEt, (recoMuons)[j].getIsolationR05().hadEt, (*recoMuons)[j].getIsolationR05().hoEt,
301			(recoMuons)[j].getIsolationR05().sumPt, (recoMuons)[j].getIsolationR05().nTracks, (*recoMuons)[j].getIsolationR05().nJets);
302			new( (*muons)[j] ) TRootMuon(localMuon);
303			if(verbosity>2) cout << localMuon << endl;
304			}
305			}
306
307
308			// Electrons
309			if(doElectron)
310			{
311			if(verbosity>1) cout << "Analysing electrons collection..." << endl;
312			edm::Handle< reco::PixelMatchGsfElectronCollection > recoElectrons;
313			iEvent.getByLabel("pixelMatchGsfElectrons", recoElectrons);
314			if(verbosity>1) std::cout << "Number of electrons in the event: " << recoElectrons->size() << std::endl;
315			for (unsigned int j=0; j<recoElectrons->size(); j++)
316			{
317			float elec_e = ((*recoElectrons)[j].superCluster())->energy();
318			float elec_px = (recoElectrons)[j].px() elec_e / (*recoElectrons)[j].p();
319			float elec_py = (recoElectrons)[j].py() elec_e / (*recoElectrons)[j].p();
320			float elec_pz = (recoElectrons)[j].pz() elec_e / (*recoElectrons)[j].p();
321			TRootElectron localElectron( elec_px, elec_py, elec_pz, elec_e, (recoElectrons)[j].vx(), (recoElectrons)[j].vy(), (recoElectrons)[j].vz(), 11, (recoElectrons)[j].charge() );
322			new( (*electrons)[j] ) TRootElectron(localElectron);
323			if(verbosity>2) cout << localElectron << endl;
324			}
325			}
326
327
328			// Photons
329			if(doPhoton)
330			{
331			if(verbosity>1) cout << "Analysing photons collection..." << endl;
332			myPhotonAnalyzer = new PhotonAnalyzer();
333			myPhotonAnalyzer->SetVerbosity(verbosity);
334			myPhotonAnalyzer->Process(iEvent, rootEvent, uncorrectedPhotons, "photons");
335			myPhotonAnalyzer->Process(iEvent, rootEvent, correctedPhotons, "correctedPhotons");
336			delete myPhotonAnalyzer;
337			}
338
339
340			// Ecal Clusters
341			if(doCluster)
342			{
343			if(verbosity>1) cout << "Analysing BasicClusters collection..." << endl;
344			myClusterAnalyzer = new ClusterAnalyzer();
345			myClusterAnalyzer->SetVerbosity(verbosity);
346			myClusterAnalyzer->Process(iEvent, rootEvent, clusters, "islandBasicClusters","islandBarrelBasicClusters", "islandBarrelShape", 110);
347			myClusterAnalyzer->Process(iEvent, rootEvent, clusters, "islandBasicClusters", "islandEndcapBasicClusters", "islandEndcapShape", 120);
348			myClusterAnalyzer->Process(iEvent, rootEvent, clusters, "hybridSuperClusters", "", "", 210);
349			delete myClusterAnalyzer;
350
351			if(verbosity>1) cout << "Analysing SuperClusters collection..." << endl;
352			mySClusterAnalyzer = new SuperClusterAnalyzer();
353			mySClusterAnalyzer->SetVerbosity(verbosity);
354			mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "islandSuperClusters","islandBarrelSuperClusters",110);
355			mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "islandSuperClusters","islandEndcapSuperClusters",120);
356			mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedIslandEndcapSuperClusters","",121);
357			mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedEndcapSuperClustersWithPreshower","",122);
358			mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "hybridSuperClusters","",210);
359			mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedHybridSuperClusters","",211);
360			delete mySClusterAnalyzer;
361			}
362
363			/*
364			cout << "Nb of BC110=" << rootEvent->nBasicClusters(110) << endl;
365			cout << "Nb of BC120"= << rootEvent->nBasicClusters(120) << endl;
366			cout << "Nb of BC210=" << rootEvent->nBasicClusters(210) << endl; );
367			if(verbosity>2) cout << "electron et="<< (*recoElectrons)[j].pt() <
368			cout << "Nb of SC110=" << rootEvent->nSuperClusters(110) << endl;
369			cout << "Nb of SC120=" << rootEvent->nSuperClusters(120) << endl;
370			cout << "Nb of SC121=" << rootEvent->nSuperClusters(121) << endl;Misc
371			cout << "Nb of SC122=" << rootEvent->nSuperClusters(122) << endl;
372			cout << "Nb of SC210=" << rootEvent->nSuperClusters(210) << endl;
373			cout << "Nb of SC211=" << rootEvent->nSuperClusters(211) << endl;
374			*/
375
376
377
378			if(doCluster)
379			{
380			ClusterAssociator* myClusterAssociator = new ClusterAssociator();
381			myClusterAssociator->SetVerbosity(verbosity);
382			myClusterAssociator->Process(superClusters, clusters);
383			//myClusterAssociator->PrintSuperClusters(superClusters, clusters,0); // 0 to print all types SC
384			delete myClusterAssociator;
385			}
386
387
388			if(doPhoton && doCluster)
389			{
390			PhotonAssociator* myPhotonAssociator = new PhotonAssociator();
391			myPhotonAssociator->SetVerbosity(verbosity);
392	mlethuil	1.2	myPhotonAssociator->AssociateSuperCluster(correctedPhotons, superClusters);
393			//myPhotonAssociator->AssociateSuperCluster(uncorrectedPhotons, superClusters);
394	mlethuil	1.3	if(doPhotonConversion) myPhotonAssociator->AssociateConversion(iEvent, correctedPhotons, superClusters, conversionTracks);
395			//if(doPhotonConversion) myPhotonAssociator->AssociateConversion(iEvent, uncorrectedPhotons, superClusters, conversionTracks);
396	mlethuil	1.1	//myPhotonAssociator->PrintPhotons(correctedPhotons,superClusters,0); // 0 to print all types associated SC
397			myPhotonAssociator->FullPrintPhotons(correctedPhotons,superClusters,clusters,0); // 0 to print all types associated SC
398			delete myPhotonAssociator;
399			}
400
401
402	mlethuil	1.4	if(doPhoton && doPhotonIsolation) //FIXME - Test HCAL rechits availability
403	mlethuil	1.1	{
404			if(verbosity>1) cout << "PhotonIsolator......" << endl;
405
406			// Get HB/HE rechits
407			edm::Handle<HBHERecHitCollection> hbheHandle;
408			iEvent.getByLabel("hbhereco", hbheHandle);
409			const HBHERecHitCollection* hbheHits = hbheHandle.product();
410
411			// Get HO rechits
412			edm::Handle<HORecHitCollection> hoHandle;
413			iEvent.getByLabel("horeco", hoHandle);
414			const HORecHitCollection* hoHits = hoHandle.product();
415
416			// Get HF rechits
417			edm::Handle<HFRecHitCollection> hfHandle;
418			iEvent.getByLabel("hfreco", hfHandle);
419			const HFRecHitCollection* hfHits = hfHandle.product();
420
421			// Get Calo Geometry
422			edm::ESHandle<CaloGeometry> geomHandle;
423			(iSetup.get<IdealGeometryRecord>()).get(geomHandle);
424			const CaloGeometry* caloGeom = geomHandle.product();
425
426			Float_t ecalIslandIsolation;
427			Float_t ecalDoubleConeIsolation;
428			Float_t hcalRecHitIsolation;
429	mlethuil	1.4	Float_t isoTracks;
430			Int_t isoNTracks;
431			TRootPhoton* localPhoton;
432
433	mlethuil	1.1	PhotonIsolator* myPhotonIsolator = new PhotonIsolator(&myConfig);
434			for (int iphoton=0; iphoton<correctedPhotons->GetEntriesFast(); iphoton++)
435			{
436	mlethuil	1.4	localPhoton = (TRootPhoton*)correctedPhotons->At(iphoton);
437
438	mlethuil	1.1	ecalIslandIsolation=-1.;
439	mlethuil	1.4	if (doCluster) ecalIslandIsolation = myPhotonIsolator->EcalIslandIsolation(localPhoton, superClusters, clusters);
440
441	mlethuil	1.1	ecalDoubleConeIsolation=-1.;
442	mlethuil	1.4	if (doCluster) ecalDoubleConeIsolation = myPhotonIsolator->EcalDoubleConeIsolation(localPhoton, superClusters, clusters);
443
444			hcalRecHitIsolation=-1.;
445	mlethuil	1.1	hcalRecHitIsolation = myPhotonIsolator->HcalRecHitIsolation(localPhoton, hbheHits, hoHits, hfHits, caloGeom);
446	mlethuil	1.4
447			isoTracks=-1.;
448			isoNTracks=-1;
449			if(doTrack) isoTracks = myPhotonIsolator->TrackerIsolation(localPhoton, tracks, isoNTracks);
450
451			localPhoton->setIsolation(ecalIslandIsolation, ecalDoubleConeIsolation, hcalRecHitIsolation, isoTracks, isoNTracks);
452			if(verbosity>2) cout << "PhotonIsolator for [" << iphoton << "] " << *localPhoton << endl;
453	mlethuil	1.1	}
454			delete myPhotonIsolator;
455			}
456
457
458
459			if(verbosity>2) cout << "Filling rootuple..." << endl;
460	mlethuil	1.5	eventTree_->Fill();
461	mlethuil	1.1	if(verbosity>2) cout << "Deleting objects..." << endl;
462			delete rootEvent;
463			if(doSignalMC) delete rootMCSignalEvent;
464			if(doMC) (*mcParticles).Delete();
465			if(doTrack) (*tracks).Delete();
466			if(doJet) (*jets).Delete();
467			if(doMuon) (*muons).Delete();
468			if(doElectron) (*electrons).Delete();
469			if(doPhoton) (*uncorrectedPhotons).Delete();
470			if(doPhoton) (*correctedPhotons).Delete();
471			if(doCluster) (*clusters).Delete();
472			if(doCluster) (*superClusters).Delete();
473	mlethuil	1.3	if(doPhotonConversion) (*conversionTracks).Delete();
474	mlethuil	1.1
475			}
476
477
478
479			//define this as a plug-in
480			//DEFINE_ANOTHER_FWK_MODULE(TotoAnalyzer);