Morgan/src/TotoAnalyzer.cc

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

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


TotoAnalyzer::TotoAnalyzer(const edm::ParameterSet& iConfig)
{
        myConfig_ = iConfig.getParameter<ParameterSet>("myConfig");
        dataType_ = myConfig_.getUntrackedParameter<string>("dataType","unknown");
        cout << "dataType: " << dataType_ << endl;
        if( dataType_=="RECO" )                         producersNames_ = iConfig.getParameter<ParameterSet>("producersNamesRECO");
        else if( dataType_=="AOD" )                     producersNames_ = iConfig.getParameter<ParameterSet>("producersNamesAOD");
        else if( dataType_=="PATAOD" )  producersNames_ = iConfig.getParameter<ParameterSet>("producersNamesPATAOD");
        else if( dataType_=="PAT" )                     producersNames_ = iConfig.getParameter<ParameterSet>("producersNamesPAT");
        else { cout << "TotoAnalyzer::TotoAnalyzer...   dataType is unknown...  exiting..." << endl; exit(1); }
        
        runInfos_ = 0;
        rootEvent = 0;
        hltAnalyzer_ = 0;
        mcParticles = 0;
        genJets = 0;
        genMETs = 0;
        rootMuMuGammaEvent = 0;
        rootMCTopTop = 0;
        mcPhotons = 0;
        vertices = 0;
        tracks = 0;
        jets = 0;
        muons = 0;
        electrons = 0;
        photons = 0;
        clusters = 0;
        superClusters = 0;
        conversionTracks = 0;
        met = 0;
}


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);
        doJetMC = myConfig_.getUntrackedParameter<bool>("doJetMC",false);
        doMETMC = myConfig_.getUntrackedParameter<bool>("doMETMC",false);
        doPDFInfo = myConfig_.getUntrackedParameter<bool>("doPDFInfo",false);
        doSignalMuMuGamma = myConfig_.getUntrackedParameter<bool>("doSignalMuMuGamma",false);
        doSignalTopTop = myConfig_.getUntrackedParameter<bool>("doSignalTopTop",false);
        doPhotonConversionMC = myConfig_.getUntrackedParameter<bool>("doPhotonConversionMC",false);
        drawMCTree = myConfig_.getUntrackedParameter<bool>("drawMCTree",false);
        doPrimaryVertex = myConfig_.getUntrackedParameter<bool>("doPrimaryVertex",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);
        doPhotonConversion = myConfig_.getUntrackedParameter<bool>("doPhotonConversion",false);
        doPhotonIsolation = myConfig_.getUntrackedParameter<bool>("doPhotonIsolation",false);
        doMET = myConfig_.getUntrackedParameter<bool>("doMET",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(producersNames_);
                hltAnalyzer_->setVerbosity(verbosity);
        }

        if(doMC)
        {
                if(verbosity>0) cout << "MC Particles info will be added to rootuple" << endl;
                mcParticles = new TClonesArray("TRootMCParticle", 1000);
                eventTree_->Branch ("MCParticles", "TClonesArray", &mcParticles);
        }
        
        if(doJetMC)
        {
                if(verbosity>0) cout << "genJets info will be added to rootuple" << endl;
                genJets = new TClonesArray("TRootParticle", 1000);
                eventTree_->Branch ("genJets", "TClonesArray", &genJets);
        }

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

        if(doSignalMuMuGamma)
        {                       
                if(verbosity>0) cout << "MC info for Z -> mu mu gamma  will be added to rootuple" << endl;
                rootMuMuGammaEvent = 0;
                cout << "Create MuMuGamma Branch" << endl;
                eventTree_->Branch ("MuMuGamma", "TRootSignalEvent", &rootMuMuGammaEvent);
        }

        if(doSignalTopTop)
        {
                if(verbosity>0) cout << "MC info for Top Top will be added to rootuple" << endl;
                rootMCTopTop = new TClonesArray("TRootTopTop", 1000);
                eventTree_->Branch ("rootMCTopTop", "TClonesArray", &rootMCTopTop);
        }

        if(doPhotonConversionMC)
        {
                if(verbosity>0) cout << "Converted MC Photons info will be added to rootuple" << endl;
                mcPhotons = new TClonesArray("TRootMCPhoton", 1000);
                eventTree_->Branch ("MCPhotons", "TClonesArray", &mcPhotons);
        }
        
        if(doPrimaryVertex)
        {
                if(verbosity>0) cout << "Vertices info will be added to rootuple" << endl;
                vertices = new TClonesArray("TRootVertex", 1000);
                eventTree_->Branch ("Vertices", "TClonesArray", &vertices);
        }

        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;
                photons = new TClonesArray("TRootPhoton", 1000);
                eventTree_->Branch ("Photons", "TClonesArray", &photons);
        }
        
        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("TRootTrack", 1000);
                eventTree_->Branch ("ConversionTracks", "TClonesArray", &conversionTracks);
                        
                conversionLikelihoodCalculator_ = new ConversionLikelihoodCalculator();
                std::string weightsString =  myConfig_.getUntrackedParameter<string>("conversionLikelihoodWeightsFile","RecoEgamma/EgammaTools/data/TMVAnalysis_Likelihood.weights.txt");
                edm::FileInPath weightsFile(weightsString.c_str() );
                conversionLikelihoodCalculator_->setWeightsFile(weightsFile.fullPath().c_str());
        }

        if(doMET)
        {
                if(verbosity>0) cout << "MET info will be added to rootuple" << endl;
                met = new TClonesArray("TRootMET", 1000);
                eventTree_->Branch ("MET", "TClonesArray", &met);
        }

}


// ------------ method called once each job just after ending the event loop  ------------
void TotoAnalyzer::endJob()
{
        // Trigger Summary Tables
        if(doHLT)
        {
                cout << "Trigger Summary Tables" << endl;
                hltAnalyzer_ ->copySummary(runInfos_);
                hltAnalyzer_->printStats();
        }

        runTree_->Fill();

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

        if(doPhotonConversion) delete conversionLikelihoodCalculator_;
        if(doHLT) delete hltAnalyzer_;
}


// ------------ 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 - Cumulated Events " << nTotEvt_
                        << " - Run " << iEvent.id().run() 
                        << " - Event " << iEvent.id().event() 
                        << " #######" << endl;

        // Global Event Infos
        rootEvent = new TRootEvent();
        rootEvent->setNb(nTotEvt_);
        rootEvent->setEventId(iEvent.id().event());
        rootEvent->setRunId(iEvent.id().run());

        /*
        // 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
        rootEvent->setGlobalHLT(true);
        if(doHLT)
        {       
                cout << endl << "Get TriggerResults..." << endl;
                if (nTotEvt_==1) hltAnalyzer_->init(iEvent, rootEvent);
                hltAnalyzer_->process(iEvent, rootEvent);
        }
        //if ( ! rootEvent->passGlobalHLT() ) return;


        // MC Info
        if(doMC)
        {
                if(verbosity>1) cout << endl << "Analysing MC info..." << endl;
                MCAnalyzer* myMCAnalyzer = new MCAnalyzer(myConfig_, producersNames_);
                myMCAnalyzer->SetVerbosity(verbosity);
                if (drawMCTree) myMCAnalyzer->DrawMCTree(iEvent, iSetup, myConfig_, producersNames_);
                if ( dataType_=="RECO" && doPDFInfo ) myMCAnalyzer->PDFInfo(iEvent, rootEvent);
                myMCAnalyzer->ProcessMCParticle(iEvent, mcParticles);
                if(doJetMC) myMCAnalyzer->ProcessGenJets(iEvent, genJets);
                if(doMETMC) myMCAnalyzer->ProcessGenMETs(iEvent, genMETs);
                if (doSignalMuMuGamma) 
                {
                        rootMuMuGammaEvent = new TRootSignalEvent();
                        myMCAnalyzer->ProcessMuMuGammaEvent(iEvent, rootMuMuGammaEvent);
                }
                if (doSignalTopTop) myMCAnalyzer->ProcessTopTopEvent(iEvent, rootMCTopTop);
                if (doPhotonConversionMC) myMCAnalyzer->ProcessConvertedPhoton(iEvent, mcPhotons);
                delete myMCAnalyzer;
        }

        
        // Get Primary Vertices
        if(doPrimaryVertex)
        {
                if(verbosity>1) cout << endl << "Analysing primary vertices collection..." << endl;
                VertexAnalyzer* myVertexAnalyzer = new VertexAnalyzer(producersNames_, verbosity);
                myVertexAnalyzer->Process(iEvent, vertices);
                myVertexAnalyzer->SelectPrimary(rootEvent, vertices);
                delete myVertexAnalyzer;
        }
        
        // Tracks
        if(doTrack)
        {
                if(verbosity>1) cout << endl << "Analysing tracks collection..." << endl;
                TrackAnalyzer* myTrackAnalyzer = new TrackAnalyzer(producersNames_, verbosity);
                myTrackAnalyzer->Process(iEvent, tracks);
                delete myTrackAnalyzer;
        }
        

        // Calo Jet
        if(doJet)
        {
                if(verbosity>1) cout << endl << "Analysing jets collection..." << endl;
                JetAnalyzer* myJetAnalyzer = new JetAnalyzer(producersNames_, myConfig_, verbosity);
                myJetAnalyzer->Process(iEvent, jets);
                delete myJetAnalyzer;
        }

        
        // Muons
        if(doMuon)
        {
                if(verbosity>1) cout << endl << "Analysing muons collection..." << endl;
                MuonAnalyzer* myMuonAnalyzer = new MuonAnalyzer(producersNames_, myConfig_, verbosity);
                if(doPrimaryVertex) myMuonAnalyzer->initIPCalculator(iEvent, iSetup, rootEvent);
                myMuonAnalyzer->Process(iEvent, muons);
                delete myMuonAnalyzer;
        }

                        
        // Lazy Tools to calculate Cluster shape variables
        EcalClusterLazyTools* lazyTools = 0;
        if( (dataType_=="RECO" || dataType_=="AOD" || dataType_=="PATAOD") && ( doElectron || doPhoton || doCluster)  )
        {
                if(verbosity>1) cout << endl << "Loading egamma LazyTools..." << endl;
                edm::InputTag reducedBarrelEcalRecHitCollection_ = producersNames_.getParameter<edm::InputTag>("reducedBarrelEcalRecHitCollection");
                edm::InputTag reducedEndcapEcalRecHitCollection_ = producersNames_.getParameter<edm::InputTag>("reducedEndcapEcalRecHitCollection");
                // FIXME - Test availability of reducedEcalRecHits...
                lazyTools = new EcalClusterLazyTools( iEvent, iSetup, reducedBarrelEcalRecHitCollection_, reducedEndcapEcalRecHitCollection_ );
        }

        // Electrons
        if(doElectron)
        {
                if(verbosity>1) cout << endl << "Analysing electrons collection..." << endl;
                ElectronAnalyzer* myElectronAnalyzer = new ElectronAnalyzer(producersNames_, myConfig_, verbosity);
                myElectronAnalyzer->Process(iEvent, electrons, *lazyTools);
                delete myElectronAnalyzer;
        }


        // Photons
        if(doPhoton)
        {
                if(verbosity>1) cout << endl << "Analysing photons collection..." << endl;
                PhotonAnalyzer* myPhotonAnalyzer = new PhotonAnalyzer(producersNames_, myConfig_, verbosity);
                // FIXME - conversionTracks may not be present...
                myPhotonAnalyzer->Process(iEvent, iSetup, rootEvent, photons, conversionTracks, conversionLikelihoodCalculator_, *lazyTools);
                delete myPhotonAnalyzer;
        }

        
        // Ecal Clusters
        if(doCluster)
        {
                if(verbosity>1) cout << endl << "Analysing BasicClusters collection..." << endl;
                ClusterAnalyzer* myClusterAnalyzer = new ClusterAnalyzer(producersNames_, verbosity);
                //myClusterAnalyzer->SetVerbosity(verbosity);
                myClusterAnalyzer->Process(iEvent, rootEvent, *lazyTools, clusters, "hybridSuperClusters","hybridBarrelBasicClusters", 110);
                myClusterAnalyzer->Process(iEvent, rootEvent, *lazyTools, clusters, "multi5x5BasicClusters", "multi5x5EndcapBasicClusters", 320);
                delete myClusterAnalyzer;
                                
                if(verbosity>1) cout << endl << "Analysing SuperClusters collection..." << endl;
                SuperClusterAnalyzer* mySClusterAnalyzer = new SuperClusterAnalyzer(producersNames_, verbosity);
                //mySClusterAnalyzer->SetVerbosity(verbosity);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "hybridSuperClusters","",210);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedHybridSuperClusters","",211);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "multi5x5SuperClusters","multi5x5EndcapSuperClusters",320);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "multi5x5SuperClustersWithPreshower","",323);
                mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedMulti5x5SuperClustersWithPreshower","",322);
                delete mySClusterAnalyzer;
        }
        
        
        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(photons, superClusters);
                delete myPhotonAssociator;
        }
        
        
        if(doPhoton && doPhotonIsolation)
        {
                if(verbosity>1) cout << endl << "Calculating Photon isolation..." << endl;

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

                PhotonIsolator* myPhotonIsolator = new PhotonIsolator(&myConfig_, &producersNames_);
                for (int iphoton=0; iphoton<photons->GetEntriesFast(); iphoton++)
                {
                        localPhoton = (TRootPhoton*)photons->At(iphoton);

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

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

                        hcalRecHitIsolation=-1.;
                        //FIXME - Test HCAL rechits availability - LoadHcalRecHits() returning bool
                        myPhotonIsolator->LoadHcalRecHits(iEvent, iSetup);
                        hcalRecHitIsolation = myPhotonIsolator->HcalRecHitIsolation(localPhoton);

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

                        localPhoton->setIsolation(ecalIslandIsolation, ecalDoubleConeIsolation, hcalRecHitIsolation, isoTracks, isoNTracks);
                        if(verbosity>4) cout << "   After isolation - ["<< setw(3) << iphoton << "] " << *localPhoton << endl;
                }
                delete myPhotonIsolator;
        }


        // Print Photons / SuperClusters / BasicClusters arborescence
        if(doPhoton)
        {
                PhotonAssociator* myPhotonAssociator = new PhotonAssociator();
                if(verbosity>2) myPhotonAssociator->FullPrintPhotons(photons,superClusters,clusters,0);
                delete myPhotonAssociator;
        }


        // MET 
        if(doMET)
        {
                if(verbosity>1) cout << endl << "Analysing Missing Et..." << endl;
                METAnalyzer* myMETAnalyzer = new METAnalyzer(producersNames_, myConfig_, verbosity);
                myMETAnalyzer->Process(iEvent, met);
                delete myMETAnalyzer;
        }


        // Associate recoParticles to mcParticles
        if(doMC)
        {
                MCAssociator* myMCAssociator = new MCAssociator(producersNames_, verbosity);
                myMCAssociator->init(iEvent, mcParticles);
                if(doPhoton) myMCAssociator->matchGenParticlesTo(photons);
                if(doElectron) myMCAssociator->matchGenParticlesTo(electrons);
                if(doMuon) myMCAssociator->matchGenParticlesTo(muons);
                if(doMET) myMCAssociator->matchGenParticlesTo(met);
                if(doJet) myMCAssociator->matchGenParticlesTo(jets);
                if(doJet && doJetMC) myMCAssociator->processGenJets(genJets, jets);
                if(verbosity>4) cout << endl << "Printing recoParticles / mcParticles associations... " << endl;
                if(verbosity>4 && doPhoton) myMCAssociator->printParticleAssociation(photons);
                if(verbosity>4 && doElectron) myMCAssociator->printParticleAssociation(electrons);
                if(verbosity>4 && doMuon) myMCAssociator->printParticleAssociation(muons);
                if(verbosity>4 && doMET) myMCAssociator->printParticleAssociation(met);
                if(verbosity>4 && doJet) myMCAssociator->printJetAssociation(jets);
                delete myMCAssociator;
        }


        if(verbosity>1) cout << endl << "Filling rootuple..." << endl;
        eventTree_->Fill();
        if(verbosity>1) cout << endl << "Deleting objects..." << endl;
        delete rootEvent;
        if(doMC) (*mcParticles).Delete();
        if(doJetMC) (*genJets).Delete();
        if(doMETMC) (*genMETs).Delete();
        if(doSignalMuMuGamma) delete rootMuMuGammaEvent;
        if(doSignalTopTop) (*rootMCTopTop).Delete();
        if(doPhotonConversionMC) (*mcPhotons).Delete();
        if(doPrimaryVertex) (*vertices).Delete();
        if(doTrack) (*tracks).Delete();
        if(doJet) (*jets).Delete();
        if(doMuon) (*muons).Delete();
        if(doElectron) (*electrons).Delete();
        if(doPhoton) (*photons).Delete();
        if(doCluster) (*clusters).Delete();
        if(doCluster) (*superClusters).Delete();
        if(doPhotonConversion) (*conversionTracks).Delete();
        if(doMET) (*met).Delete();
        if(verbosity>0) cout << endl;
}

Revision:	1.29
Committed:	Fri Apr 17 15:17:38 2009 UTC (16 years ago) by lethuill
Content type:	text/plain
Branch:	MAIN
Changes since 1.28:	+1 -0 lines
Log Message:	Add significance of the 3D impact parameter for leptons. Calculation done in new LeptonAnalyzer class (Base for MuonAnalyzer and ElectronAnalyzer)