TreeMod/src/ObjectCleaningMod.cc

// $Id: ObjectCleaningMod.cc,v 1.6 2008/09/30 16:38:20 sixie Exp $

#include "MitAna/TreeMod/interface/ObjectCleaningMod.h"
#include <TH1D.h>
#include <TH2D.h>
#include "MitAna/DataTree/interface/Names.h"
#include "MitAna/DataCont/interface/ObjArray.h"
#include "MitAna/Utils/interface/IsolationTools.h"
#include "MitCommon/MathTools/interface/MathUtils.h"
#include "MitAna/Utils/interface/IsolationTools.h"

using namespace mithep;

ClassImp(mithep::ObjectCleaningMod)

//--------------------------------------------------------------------------------------------------
  ObjectCleaningMod::ObjectCleaningMod(const char *name, const char *title) : 
  BaseMod(name,title),
  fPrintDebug(false),
  fMCPartName(Names::gkMCPartBrn),
  fMuonName(Names::gkMuonBrn),
  fElectronName(Names::gkElectronBrn),
  fJetName(Names::gkCaloJetBrn),
  fSCJetName(Names::gkSC5JetBrn),
  fMetName(Names::gkCaloMetBrn),
  fGoodElectronsName("GoodElectrons"),        
  fGoodMuonsName("GoodMuons"),        
  fGoodCentralJetsName("GoodCentralJets" ),        
  fGoodForwardJetsName("GoodForwardJets" ),  
  fGoodCentralSCJetsName("GoodCentralSCJets" ),        
  fGoodForwardSCJetsName("GoodForwardSCJets" ),
  fGenLeptonsName("GenLeptons"),
  fParticles(0),
  fMuons(0),
  fElectrons(0)
{
  // Constructor.
}

//--------------------------------------------------------------------------------------------------
void ObjectCleaningMod::Begin()
{
  // Run startup code on the client machine. For this module, we dont do
  // anything here.
}

//--------------------------------------------------------------------------------------------------
void ObjectCleaningMod::Process()
{
  // Process entries of the tree. For this module, we just load the branches and
  //output debug info or not   

  fNEventsProcessed++;
 
  if (fNEventsProcessed % 1000 == 0 || fPrintDebug) {
    time_t systime;
    systime = time(NULL);

    cerr << endl << "ObjectCleaningMod : Process Event " << fNEventsProcessed << "  Time: " << ctime(&systime) << endl;  
  }
  
 //Muons
  ObjArray<Muon> *GoodMuons = new ObjArray<Muon>; 
  LoadBranch(fMuonName);

  for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
    Muon *mu = fMuons->At(i);
    fAllMuonPtHist->Fill(mu->Pt());
    fAllMuonEtaHist->Fill(mu->Eta());    

    Double_t MuonClass = -1;    
    //Find Class of Muons: Global(0), Standalone(1), Tracker(2)
    //Double_t MuonClass = -1;
    
    if (mu->GlobalTrk())      
      MuonClass = 0;
    else if (mu->StandaloneTrk())

      MuonClass = 1;
    else if (mu->TrackerTrk())
      MuonClass = 2;

    //Define the ID Cuts
    const int nCuts = 4;
    double cutValue[nCuts] = {0.1, 3.0, 3.0, 1.5 };
    bool passCut[nCuts] = {false, false, false, false};
    
    double muonD0 = -0.05;

    muonD0 = mu->BestTrk()->D0();
    if(muonD0 < cutValue[0] &&  MuonClass == 0 ) 
      passCut[0] = true;
    if(mu->IsoR03SumPt() < cutValue[1]) passCut[1] = true;
    if(mu->IsoR03EmEt() + 
       mu->IsoR03HadEt() < cutValue[2]) passCut[2] = true;    
    if(mu->Pt() > 5)
      passCut[3] = true;   

    // if(mu->Trk()->NHits() > 6)
//       passCut[4] = true;

//     if(mu->Trk()->Chi2() / mu->Trk()->Ndof() < 5 )
//       passCut[5] = true;

    // Final decision
    bool allCuts = true;
    for(int c=0; c<nCuts; c++) {
      allCuts = allCuts & passCut[c];
    }

    //make muon ID selection histograms
    fMuonSelection->Fill(-1);    
    //Fill the rest of the muon selection histograms
    for (int k=0;k<3;k++) {
      bool pass = true;
      for (int p=0;p<=k;p++)
        pass = (pass && passCut[p]);
      
      if (pass) {
        fMuonSelection->Fill(k);
      }
    } 
    //Fill histogram for the Good Muons
    if ( allCuts
         && abs(mu->Eta()) < 2.5
      ) {     
      fGoodMuonPtHist->Fill(mu->Pt());
      fGoodMuonEtaHist->Fill(mu->Eta());       
      GoodMuons->Add(mu);
    }                    
  }
  

  //Get Electrons
  LoadBranch(fElectronName);

  //we have to use a vector first and then fill the ObjArray with the vector 
  //contents because ObJArray does not allow removal of duplicates.
  vector<Electron*> GoodElectronsVector;
  vector<Electron*> GoodLooseElectronsVector;
  vector<Electron*> GoodTightElectronsVector;
  vector<Electron*> GoodLikelihoodIDElectronsVector;

  for (UInt_t i=0; i<fElectrons->GetEntries(); ++i) {    
    Electron *e = fElectrons->At(i);   
    fAllElectronPtHist->Fill(e->Pt());
    fAllElectronEtaHist->Fill(e->Eta());               
    
    //from RecoEgamma/ElectronIdentification/src/CurBasedElectronID.cc  : CMSSW 2_1_0
    Int_t InEndcapOrBarrel = (fabs(e->Eta()) < 1.479)? 0 : 1;     
    double sigmaEtaEta = e->CovEtaEta();
    //have to correct sigma_etaeta dependance on eta in the endcap
    if (InEndcapOrBarrel == 1) {
      sigmaEtaEta = sigmaEtaEta - 0.02*(fabs(e->Eta()) - 2.3);
      e->SetCovEtaEta(e->CovEtaEta() - 0.02*(fabs(e->Eta()) - 2.3));
    }

    //Calculate the electron category for determining which cuts to make
    Int_t ElectronCategory;
    Double_t fBrem = (e->PIn() - e->POut())/e->PIn();
    if((fabs(e->Eta())<1.479 && fBrem<0.06) || (fabs(e->Eta())>1.479 && fBrem<0.1)) 
      ElectronCategory=1;
    else if (e->ESuperClusterOverP() < 1.2 && e->ESuperClusterOverP() > 0.8) 
      ElectronCategory=0;
    else 
      ElectronCategory=2;       
    
    // Final decision
    bool allCuts = true;  
    allCuts = e->PassTightID();    

    //Check whether it overlaps with a good muon.
    bool isMuonOverlap = false;
    for (UInt_t j=0; j<GoodMuons->GetEntries();j++) {
      double deltaR = MathUtils::DeltaR(GoodMuons->At(j)->Mom(), e->Mom());     
      if (deltaR < 0.1) {
        isMuonOverlap = true;
        break;   
      }
    }

    //Check whether it overlaps with another electron candidate
    bool isElectronOverlap = false;
    for (UInt_t j=0; j<GoodElectronsVector.size(); j++) {      
      double deltaR = MathUtils::DeltaR(GoodElectronsVector[j]->Mom(), e->Mom());

      if (deltaR < 0.1) {
        isElectronOverlap = true;        
        
        //if there's an overlap and the new one being considered has E/P closer to 1.0 
        //then replace the overlapping one with this new one because it's better.
        //Once we get super cluster info, we should make sure the superclusters match 
        //before declaring it is a duplicate
        //can have one SC matched to many tracks. or one track matched to many SC's. 
        //we should cover all the cases. see SUSYAnalyzer...
        //Si: This will be unnecessary very soon. It is to be removed in reconstruction.
        if (abs(GoodElectronsVector[j]->ESuperClusterOverP() - 1) > 
            abs(e->ESuperClusterOverP() - 1)) {    
          GoodElectronsVector[j] = e;
        }        
        break;
      }
    }

    //These are Good Electrons
    if ( allCuts && !isMuonOverlap && !isElectronOverlap                   
      ) {       
      fGoodElectronPtHist->Fill(e->Pt());
      fGoodElectronEtaHist->Fill(e->Eta());  
      fGoodElectronClassification->Fill(e->Classification());      
      GoodElectronsVector.push_back(fElectrons->At(i));   
    }   

  }     
  //fill the electron ObjArray with the contents of the vector
  //this is necessary because I want to swap out the duplicates, can't be done with ObjArray...
  ObjArray<Electron> *GoodElectrons = new ObjArray<Electron>;
  for (UInt_t j=0; j<GoodElectronsVector.size(); j++) 
    GoodElectrons->Add(GoodElectronsVector[j]);
  
  //Get Jet info
  ObjArray<Jet> *GoodCentralJets = new ObjArray<Jet>; 
  ObjArray<Jet> *GoodForwardJets = new ObjArray<Jet>; 
  LoadBranch(fJetName);
  for (UInt_t i=0; i<fJets->GetEntries(); ++i) {
    Jet *jet = fJets->At(i);
        
    bool isElectronOverlap =  false;
     
    //Check for overlap with an electron
    for (UInt_t j=0; j<GoodElectrons->GetEntries(); j++) {
      double deltaR = MathUtils::DeltaR(GoodElectrons->At(j)->Mom(),jet->Mom());  
      if (deltaR < 0.1) {
        isElectronOverlap = true;                       
        break;           
      }      
    }
    if (isElectronOverlap) continue; //skip this jet if it was an overlap
    
    
    fAllJetPtHist->Fill(jet->Pt());
    fAllJetEtaHist->Fill(jet->Eta());
        
    const int nCuts = 4;
    double cutValue[nCuts] = {1.0, 15., 2.5, 0.2};
    bool passCut[nCuts] = {false, false, false, false};
    
    passCut[0] = (!isElectronOverlap); //This is supposed to check e/ph overlap
    if(jet->Et() > cutValue[1])           passCut[1] = true;
    if(fabs(jet->Eta()) < cutValue[2])    passCut[2] = true;        
    if(jet->Alpha() > cutValue[3] ||
       jet->Et() > 20.) 
      passCut[3] = true; 
    
    // Final decision
    bool passAllCentralJetCuts = true;
    bool passAllForwardJetCuts = true;
    for(int i=0; i<nCuts; i++) {
      passAllCentralJetCuts = passAllCentralJetCuts && passCut[i];
      passAllForwardJetCuts = passAllForwardJetCuts && ((i==2)?(!passCut[i]):passCut[i]);       
    }
    
    //make electron ID selection histogram
    fCentralJetSelection->Fill(-1);
   
    for (int k=0;k<nCuts;k++) {
      bool pass = true;
      for (int p=0;p<=k;p++)
        pass = (pass && passCut[p]);
      
      if (pass) {
        fCentralJetSelection->Fill(k);  
     }
    } 
                
    //Save Good Jets
    if (passAllCentralJetCuts) {
      GoodCentralJets->Add(jet);     
    }
    
    if(passAllForwardJetCuts)
      GoodForwardJets->Add(jet);
    
  } //for all jets
   

//Get Siscone 5 Jet info
  ObjArray<Jet> *GoodCentralSCJets = new ObjArray<Jet>; 
  ObjArray<Jet> *GoodForwardSCJets = new ObjArray<Jet>; 
  LoadBranch(fSCJetName);
  for (UInt_t i=0; i<fSCJets->GetEntries(); ++i) {
    Jet *jet = fSCJets->At(i);
        
    bool isElectronOverlap =  false;
     
    //Check for overlap with an electron
    for (UInt_t j=0; j<GoodElectrons->GetEntries(); j++) {
      double deltaR = MathUtils::DeltaR(GoodElectrons->At(j)->Mom(),jet->Mom());  
      if (deltaR < 0.1) {
        isElectronOverlap = true;                       
        break;           
      }      
    }
    if (isElectronOverlap) continue; //skip this jet if it was an overlap
            
    const int nCuts = 4;
    double cutValue[nCuts] = {1.0, 15., 2.5, 0.2};
    bool passCut[nCuts] = {false, false, false, false};
    
    passCut[0] = (!isElectronOverlap); //This is supposed to check e/ph overlap
    if(jet->Et() > cutValue[1])           passCut[1] = true;
    if(fabs(jet->Eta()) < cutValue[2])    passCut[2] = true;        
    if(jet->Alpha() > cutValue[3] ||
       jet->Et() > 20.) 
      passCut[3] = true; 
    
    // Final decision
    bool passAllCentralJetCuts = true;
    bool passAllForwardJetCuts = true;
    for(int i=0; i<nCuts; i++) {
      passAllCentralJetCuts = passAllCentralJetCuts && passCut[i];
      passAllForwardJetCuts = passAllForwardJetCuts && ((i==2)?(!passCut[i]):passCut[i]);       
    }
        
    for (int k=0;k<nCuts;k++) {
      bool pass = true;
      for (int p=0;p<=k;p++)
        pass = (pass && passCut[p]);
            
    }       
    
    //Save Good SCJets
    if (passAllCentralJetCuts) {
      GoodCentralSCJets->Add(jet);      
    }
    
    if(passAllForwardJetCuts)
      GoodForwardSCJets->Add(jet);
    
  } //for all jets    

   //Get MET
  LoadBranch(fMetName);

  Met *met = fMet->At(0); //define the met
  fMetPtHist->Fill(met->Pt());
  fMetPhiHist->Fill(met->Phi());
      
   //Final Summary Debug Output   
  if ( fPrintDebug ) {
    cerr << "Event Dump: " << fNEventsProcessed << endl;
    
    //print out event content to text
    cerr << "Electrons" << endl;
    for (UInt_t i = 0; i < GoodElectrons->GetEntries(); i++) {
      cerr << i << " " << GoodElectrons->At(i)->Pt() << " " << GoodElectrons->At(i)->Eta() 
           << " " << GoodElectrons->At(i)->Phi() << " " 
           << GoodElectrons->At(i)->ESuperClusterOverP() << endl;    
    }   
    
    cerr << "Muons" << endl;
    for (UInt_t i = 0; i < GoodMuons->GetEntries(); i++) {
      cerr << i << " " << GoodMuons->At(i)->Pt() << " " << GoodMuons->At(i)->Eta() 
           << " " << GoodMuons->At(i)->Phi() << endl;    
    }
    
    cerr << "Central Jets" << endl;
    for (UInt_t i = 0; i < GoodCentralJets->GetEntries(); i++) {
      cerr << i << " " << GoodCentralJets->At(i)->Pt() << " " 
           << GoodCentralJets->At(i)->Eta() << " " << GoodCentralJets->At(i)->Phi() << endl;    
    }
    
    cerr << "MET" << endl;        
    cerr << met->Pt() << " " << met->Eta() << " " 
         << met->Phi() << endl;                         
  }   
  
  //Save Objects for Other Modules to use
  AddObjThisEvt(GoodElectrons, fGoodElectronsName.Data());
  AddObjThisEvt(GoodMuons, fGoodMuonsName.Data());
  AddObjThisEvt(GoodCentralJets, fGoodCentralJetsName.Data());
  AddObjThisEvt(GoodForwardJets, fGoodForwardJetsName.Data());
  AddObjThisEvt(GoodCentralSCJets, fGoodCentralSCJetsName.Data());
  AddObjThisEvt(GoodForwardSCJets, fGoodForwardSCJetsName.Data());
  
}


//--------------------------------------------------------------------------------------------------
void ObjectCleaningMod::SlaveBegin()
{
  // Run startup code on the computer (slave) doing the actual analysis. Here,
  // we typically initialize histograms and other analysis objects and request
  // branches. For this module, we request a branch of the MitTree.

  ReqBranch(fMCPartName,                fParticles);
  ReqBranch(fMuonName,                  fMuons);
  ReqBranch(fElectronName,              fElectrons);
  ReqBranch(fJetName,                   fJets);
  ReqBranch(fSCJetName,                 fSCJets);
  ReqBranch(fMetName,                   fMet);


  fAllMuonPtHist                   = new TH1D("hAllMuonPtHist",";p_{t};#",200,0.,200.);
  fAllMuonEtaHist                  = new TH1D("hAllMuonEtaHist",";#eta;#",100,-5.,5.);
  fGoodMuonPtHist                  = new TH1D("hGoodMuonPtHist",";p_{t};#",200,0.,200.);
  fGoodMuonEtaHist                 = new TH1D("hGoodMuonEtaHist",";#eta;#",21,-5.,5.);
  fMuonSelection                   = new TH1D("hMuonSelection", ";MuonSelection;#",4,-1.5,2.5 ) ;  
  AddOutput(fAllMuonPtHist);
  AddOutput(fAllMuonEtaHist);
  AddOutput(fGoodMuonPtHist);
  AddOutput(fGoodMuonEtaHist);
  AddOutput(fMuonSelection);
 
  fAllElectronPtHist                     = new TH1D("hAllElectronPtHist",";p_{t};#",100,0.,200.);  
  fAllElectronEtaHist                    = new TH1D("hAllElectronEtaHist",";#eta;#",100,-5.,5.);
  fGoodElectronPtHist                    = new TH1D("hGoodElectronPtHist",";p_{t};#",25,0.,200.);
  fGoodElectronEtaHist                   = new TH1D("hGoodElectronEtaHist",";#eta;#",21,-5.,5.);
  fGoodElectronClassification            = new TH1D("hGoodElectronClassification",
                                                    ";Good Electron Classification;#",51,0,50 ) ; 
  AddOutput(fAllElectronPtHist);
  AddOutput(fAllElectronEtaHist);
  AddOutput(fGoodElectronPtHist);
  AddOutput(fGoodElectronEtaHist);
  AddOutput(fGoodElectronClassification);

  //Jet Plots
  fAllJetPtHist                   = new TH1D("hAllJetPtHist",";All Jet p_{t};#",100,0.,200.);
  fAllJetEtaHist                  = new TH1D("hAllJetEtaHist",";All Jet #eta;#",160,-8.,8.);
  fCentralJetSelection            = new TH1D("hCentralJetSelection",
                                             ";CentralJetSelection;#",5,-1.5,3.5 ) ; 
  AddOutput(fAllJetPtHist);
  AddOutput(fAllJetEtaHist);
  AddOutput(fCentralJetSelection);

  //MET Plots
  fMetPtHist                = new TH1D("hMetPtHist",";p_{t};#",30,0.,300.);
  fMetPhiHist               = new TH1D("hMetPhiHist",";#phi;#",28,-3.5,3.5);
  AddOutput(fMetPtHist);
  AddOutput(fMetPhiHist);


}

//--------------------------------------------------------------------------------------------------
void ObjectCleaningMod::SlaveTerminate()
{
  // Run finishing code on the computer (slave) that did the analysis. For this
  // module, we dont do anything here.

}

//--------------------------------------------------------------------------------------------------
void ObjectCleaningMod::Terminate()
{
  // Run finishing code on the client computer. For this module, we dont do
  // anything here.
}
Revision:	1.7
Committed:	Tue Oct 14 05:12:49 2008 UTC (16 years, 7 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.6:	+1 -1 lines
State:	*FILE REMOVED*
Log Message:	Moved to MitPhysics/Mods
#	User	Rev	Content
1	loizides	1.7	// $Id: ObjectCleaningMod.cc,v 1.6 2008/09/30 16:38:20 sixie Exp $
2	sixie	1.1
3			#include "MitAna/TreeMod/interface/ObjectCleaningMod.h"
4			#include <TH1D.h>
5			#include <TH2D.h>
6			#include "MitAna/DataTree/interface/Names.h"
7			#include "MitAna/DataCont/interface/ObjArray.h"
8	sixie	1.4	#include "MitAna/Utils/interface/IsolationTools.h"
9	sixie	1.2	#include "MitCommon/MathTools/interface/MathUtils.h"
10	sixie	1.6	#include "MitAna/Utils/interface/IsolationTools.h"
11	sixie	1.4
12	sixie	1.1	using namespace mithep;
13
14			ClassImp(mithep::ObjectCleaningMod)
15
16			//--------------------------------------------------------------------------------------------------
17	sixie	1.6	ObjectCleaningMod::ObjectCleaningMod(const char name, const char title) :
18	sixie	1.1	BaseMod(name,title),
19	sixie	1.6	fPrintDebug(false),
20	sixie	1.1	fMCPartName(Names::gkMCPartBrn),
21			fMuonName(Names::gkMuonBrn),
22			fElectronName(Names::gkElectronBrn),
23			fJetName(Names::gkCaloJetBrn),
24	sixie	1.6	fSCJetName(Names::gkSC5JetBrn),
25	sixie	1.1	fMetName(Names::gkCaloMetBrn),
26	sixie	1.6	fGoodElectronsName("GoodElectrons"),
27			fGoodMuonsName("GoodMuons"),
28			fGoodCentralJetsName("GoodCentralJets" ),
29			fGoodForwardJetsName("GoodForwardJets" ),
30			fGoodCentralSCJetsName("GoodCentralSCJets" ),
31			fGoodForwardSCJetsName("GoodForwardSCJets" ),
32			fGenLeptonsName("GenLeptons"),
33	sixie	1.1	fParticles(0),
34			fMuons(0),
35	sixie	1.6	fElectrons(0)
36	sixie	1.1	{
37			// Constructor.
38			}
39
40			//--------------------------------------------------------------------------------------------------
41			void ObjectCleaningMod::Begin()
42			{
43			// Run startup code on the client machine. For this module, we dont do
44			// anything here.
45			}
46
47			//--------------------------------------------------------------------------------------------------
48			void ObjectCleaningMod::Process()
49			{
50			// Process entries of the tree. For this module, we just load the branches and
51	sixie	1.6	//output debug info or not
52	sixie	1.1
53			fNEventsProcessed++;
54
55	sixie	1.6	if (fNEventsProcessed % 1000 == 0 \|\| fPrintDebug) {
56			time_t systime;
57			systime = time(NULL);
58	sixie	1.1
59	sixie	1.6	cerr << endl << "ObjectCleaningMod : Process Event " << fNEventsProcessed << " Time: " << ctime(&systime) << endl;
60	sixie	1.1	}
61
62	sixie	1.6	//Muons
63	sixie	1.1	ObjArray<Muon> *GoodMuons = new ObjArray<Muon>;
64			LoadBranch(fMuonName);
65	sixie	1.6
66	sixie	1.1	for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
67			Muon *mu = fMuons->At(i);
68			fAllMuonPtHist->Fill(mu->Pt());
69	sixie	1.6	fAllMuonEtaHist->Fill(mu->Eta());
70
71			Double_t MuonClass = -1;
72			//Find Class of Muons: Global(0), Standalone(1), Tracker(2)
73			//Double_t MuonClass = -1;
74
75			if (mu->GlobalTrk())
76			MuonClass = 0;
77			else if (mu->StandaloneTrk())
78
79			MuonClass = 1;
80			else if (mu->TrackerTrk())
81			MuonClass = 2;
82
83	sixie	1.1	//Define the ID Cuts
84	sixie	1.2	const int nCuts = 4;
85	sixie	1.6	double cutValue[nCuts] = {0.1, 3.0, 3.0, 1.5 };
86			bool passCut[nCuts] = {false, false, false, false};
87	sixie	1.1
88			double muonD0 = -0.05;
89	sixie	1.6
90	sixie	1.1	muonD0 = mu->BestTrk()->D0();
91	sixie	1.6	if(muonD0 < cutValue[0] && MuonClass == 0 )
92	sixie	1.1	passCut[0] = true;
93			if(mu->IsoR03SumPt() < cutValue[1]) passCut[1] = true;
94			if(mu->IsoR03EmEt() +
95	sixie	1.6	mu->IsoR03HadEt() < cutValue[2]) passCut[2] = true;
96			if(mu->Pt() > 5)
97			passCut[3] = true;
98
99			// if(mu->Trk()->NHits() > 6)
100			// passCut[4] = true;
101
102			// if(mu->Trk()->Chi2() / mu->Trk()->Ndof() < 5 )
103			// passCut[5] = true;
104	sixie	1.2
105	sixie	1.1	// Final decision
106			bool allCuts = true;
107			for(int c=0; c<nCuts; c++) {
108			allCuts = allCuts & passCut[c];
109			}
110	sixie	1.6
111	sixie	1.1	//make muon ID selection histograms
112	sixie	1.6	fMuonSelection->Fill(-1);
113	sixie	1.1	//Fill the rest of the muon selection histograms
114			for (int k=0;k<3;k++) {
115			bool pass = true;
116			for (int p=0;p<=k;p++)
117			pass = (pass && passCut[p]);
118
119			if (pass) {
120			fMuonSelection->Fill(k);
121			}
122			}
123			//Fill histogram for the Good Muons
124			if ( allCuts
125	sixie	1.6	&& abs(mu->Eta()) < 2.5
126	sixie	1.1	) {
127			fGoodMuonPtHist->Fill(mu->Pt());
128			fGoodMuonEtaHist->Fill(mu->Eta());
129			GoodMuons->Add(mu);
130			}
131			}
132
133	sixie	1.6
134	sixie	1.1	//Get Electrons
135			LoadBranch(fElectronName);
136	sixie	1.6
137	sixie	1.1	//we have to use a vector first and then fill the ObjArray with the vector
138			//contents because ObJArray does not allow removal of duplicates.
139			vector<Electron*> GoodElectronsVector;
140	sixie	1.6	vector<Electron*> GoodLooseElectronsVector;
141			vector<Electron*> GoodTightElectronsVector;
142			vector<Electron*> GoodLikelihoodIDElectronsVector;
143
144			for (UInt_t i=0; i<fElectrons->GetEntries(); ++i) {
145			Electron *e = fElectrons->At(i);
146	sixie	1.1	fAllElectronPtHist->Fill(e->Pt());
147	sixie	1.6	fAllElectronEtaHist->Fill(e->Eta());
148	sixie	1.1
149	sixie	1.6	//from RecoEgamma/ElectronIdentification/src/CurBasedElectronID.cc : CMSSW 2_1_0
150			Int_t InEndcapOrBarrel = (fabs(e->Eta()) < 1.479)? 0 : 1;
151			double sigmaEtaEta = e->CovEtaEta();
152			//have to correct sigma_etaeta dependance on eta in the endcap
153			if (InEndcapOrBarrel == 1) {
154			sigmaEtaEta = sigmaEtaEta - 0.02*(fabs(e->Eta()) - 2.3);
155			e->SetCovEtaEta(e->CovEtaEta() - 0.02*(fabs(e->Eta()) - 2.3));
156			}
157
158			//Calculate the electron category for determining which cuts to make
159			Int_t ElectronCategory;
160			Double_t fBrem = (e->PIn() - e->POut())/e->PIn();
161			if((fabs(e->Eta())<1.479 && fBrem<0.06) \|\| (fabs(e->Eta())>1.479 && fBrem<0.1))
162			ElectronCategory=1;
163			else if (e->ESuperClusterOverP() < 1.2 && e->ESuperClusterOverP() > 0.8)
164			ElectronCategory=0;
165			else
166			ElectronCategory=2;
167	sixie	1.1
168			// Final decision
169	sixie	1.2	bool allCuts = true;
170	sixie	1.6	allCuts = e->PassTightID();
171
172	sixie	1.1	//Check whether it overlaps with a good muon.
173			bool isMuonOverlap = false;
174			for (UInt_t j=0; j<GoodMuons->GetEntries();j++) {
175	loizides	1.3	double deltaR = MathUtils::DeltaR(GoodMuons->At(j)->Mom(), e->Mom());
176	sixie	1.1	if (deltaR < 0.1) {
177			isMuonOverlap = true;
178			break;
179			}
180			}
181	sixie	1.6
182	sixie	1.1	//Check whether it overlaps with another electron candidate
183			bool isElectronOverlap = false;
184	sixie	1.2	for (UInt_t j=0; j<GoodElectronsVector.size(); j++) {
185	loizides	1.3	double deltaR = MathUtils::DeltaR(GoodElectronsVector[j]->Mom(), e->Mom());
186	sixie	1.2
187	sixie	1.1	if (deltaR < 0.1) {
188			isElectronOverlap = true;
189
190			//if there's an overlap and the new one being considered has E/P closer to 1.0
191			//then replace the overlapping one with this new one because it's better.
192			//Once we get super cluster info, we should make sure the superclusters match
193			//before declaring it is a duplicate
194			//can have one SC matched to many tracks. or one track matched to many SC's.
195			//we should cover all the cases. see SUSYAnalyzer...
196	sixie	1.6	//Si: This will be unnecessary very soon. It is to be removed in reconstruction.
197	sixie	1.1	if (abs(GoodElectronsVector[j]->ESuperClusterOverP() - 1) >
198			abs(e->ESuperClusterOverP() - 1)) {
199			GoodElectronsVector[j] = e;
200			}
201			break;
202			}
203			}
204	sixie	1.2
205	sixie	1.1	//These are Good Electrons
206	sixie	1.6	if ( allCuts && !isMuonOverlap && !isElectronOverlap
207	sixie	1.1	) {
208			fGoodElectronPtHist->Fill(e->Pt());
209			fGoodElectronEtaHist->Fill(e->Eta());
210	sixie	1.6	fGoodElectronClassification->Fill(e->Classification());
211	sixie	1.1	GoodElectronsVector.push_back(fElectrons->At(i));
212	sixie	1.6	}
213
214			}
215	sixie	1.1	//fill the electron ObjArray with the contents of the vector
216	sixie	1.6	//this is necessary because I want to swap out the duplicates, can't be done with ObjArray...
217	sixie	1.1	ObjArray<Electron> *GoodElectrons = new ObjArray<Electron>;
218			for (UInt_t j=0; j<GoodElectronsVector.size(); j++)
219			GoodElectrons->Add(GoodElectronsVector[j]);
220
221			//Get Jet info
222			ObjArray<Jet> *GoodCentralJets = new ObjArray<Jet>;
223			ObjArray<Jet> *GoodForwardJets = new ObjArray<Jet>;
224			LoadBranch(fJetName);
225			for (UInt_t i=0; i<fJets->GetEntries(); ++i) {
226			Jet *jet = fJets->At(i);
227
228	sixie	1.6	bool isElectronOverlap = false;
229	sixie	1.1
230			//Check for overlap with an electron
231			for (UInt_t j=0; j<GoodElectrons->GetEntries(); j++) {
232	loizides	1.3	double deltaR = MathUtils::DeltaR(GoodElectrons->At(j)->Mom(),jet->Mom());
233	sixie	1.1	if (deltaR < 0.1) {
234			isElectronOverlap = true;
235			break;
236			}
237			}
238			if (isElectronOverlap) continue; //skip this jet if it was an overlap
239
240
241			fAllJetPtHist->Fill(jet->Pt());
242			fAllJetEtaHist->Fill(jet->Eta());
243	sixie	1.6
244	sixie	1.1	const int nCuts = 4;
245	sixie	1.6	double cutValue[nCuts] = {1.0, 15., 2.5, 0.2};
246	sixie	1.1	bool passCut[nCuts] = {false, false, false, false};
247
248			passCut[0] = (!isElectronOverlap); //This is supposed to check e/ph overlap
249	sixie	1.6	if(jet->Et() > cutValue[1]) passCut[1] = true;
250			if(fabs(jet->Eta()) < cutValue[2]) passCut[2] = true;
251			if(jet->Alpha() > cutValue[3] \|\|
252			jet->Et() > 20.)
253			passCut[3] = true;
254	sixie	1.1
255			// Final decision
256			bool passAllCentralJetCuts = true;
257			bool passAllForwardJetCuts = true;
258			for(int i=0; i<nCuts; i++) {
259	sixie	1.6	passAllCentralJetCuts = passAllCentralJetCuts && passCut[i];
260			passAllForwardJetCuts = passAllForwardJetCuts && ((i==2)?(!passCut[i]):passCut[i]);
261	sixie	1.1	}
262
263			//make electron ID selection histogram
264			fCentralJetSelection->Fill(-1);
265	sixie	1.6
266	sixie	1.1	for (int k=0;k<nCuts;k++) {
267			bool pass = true;
268			for (int p=0;p<=k;p++)
269			pass = (pass && passCut[p]);
270
271			if (pass) {
272			fCentralJetSelection->Fill(k);
273	sixie	1.6	}
274	sixie	1.1	}
275	sixie	1.6
276			//Save Good Jets
277			if (passAllCentralJetCuts) {
278			GoodCentralJets->Add(jet);
279			}
280	sixie	1.1
281	sixie	1.6	if(passAllForwardJetCuts)
282			GoodForwardJets->Add(jet);
283	sixie	1.1
284	sixie	1.6	} //for all jets
285
286
287			//Get Siscone 5 Jet info
288			ObjArray<Jet> *GoodCentralSCJets = new ObjArray<Jet>;
289			ObjArray<Jet> *GoodForwardSCJets = new ObjArray<Jet>;
290			LoadBranch(fSCJetName);
291			for (UInt_t i=0; i<fSCJets->GetEntries(); ++i) {
292			Jet *jet = fSCJets->At(i);
293
294			bool isElectronOverlap = false;
295
296			//Check for overlap with an electron
297			for (UInt_t j=0; j<GoodElectrons->GetEntries(); j++) {
298			double deltaR = MathUtils::DeltaR(GoodElectrons->At(j)->Mom(),jet->Mom());
299			if (deltaR < 0.1) {
300			isElectronOverlap = true;
301			break;
302			}
303	sixie	1.1	}
304	sixie	1.6	if (isElectronOverlap) continue; //skip this jet if it was an overlap
305
306			const int nCuts = 4;
307			double cutValue[nCuts] = {1.0, 15., 2.5, 0.2};
308			bool passCut[nCuts] = {false, false, false, false};
309
310			passCut[0] = (!isElectronOverlap); //This is supposed to check e/ph overlap
311			if(jet->Et() > cutValue[1]) passCut[1] = true;
312			if(fabs(jet->Eta()) < cutValue[2]) passCut[2] = true;
313			if(jet->Alpha() > cutValue[3] \|\|
314			jet->Et() > 20.)
315			passCut[3] = true;
316	sixie	1.1
317	sixie	1.6	// Final decision
318			bool passAllCentralJetCuts = true;
319			bool passAllForwardJetCuts = true;
320			for(int i=0; i<nCuts; i++) {
321			passAllCentralJetCuts = passAllCentralJetCuts && passCut[i];
322			passAllForwardJetCuts = passAllForwardJetCuts && ((i==2)?(!passCut[i]):passCut[i]);
323			}
324
325			for (int k=0;k<nCuts;k++) {
326			bool pass = true;
327			for (int p=0;p<=k;p++)
328			pass = (pass && passCut[p]);
329
330			}
331	sixie	1.1
332	sixie	1.6	//Save Good SCJets
333	sixie	1.1	if (passAllCentralJetCuts) {
334	sixie	1.6	GoodCentralSCJets->Add(jet);
335	sixie	1.1	}
336
337			if(passAllForwardJetCuts)
338	sixie	1.6	GoodForwardSCJets->Add(jet);
339	sixie	1.1
340	sixie	1.6	} //for all jets
341
342	sixie	1.1	//Get MET
343	sixie	1.6	LoadBranch(fMetName);
344
345			Met *met = fMet->At(0); //define the met
346			fMetPtHist->Fill(met->Pt());
347			fMetPhiHist->Fill(met->Phi());
348	sixie	1.1
349	sixie	1.2	//Final Summary Debug Output
350	sixie	1.6	if ( fPrintDebug ) {
351			cerr << "Event Dump: " << fNEventsProcessed << endl;
352
353			//print out event content to text
354			cerr << "Electrons" << endl;
355			for (UInt_t i = 0; i < GoodElectrons->GetEntries(); i++) {
356			cerr << i << " " << GoodElectrons->At(i)->Pt() << " " << GoodElectrons->At(i)->Eta()
357			<< " " << GoodElectrons->At(i)->Phi() << " "
358			<< GoodElectrons->At(i)->ESuperClusterOverP() << endl;
359			}
360
361			cerr << "Muons" << endl;
362			for (UInt_t i = 0; i < GoodMuons->GetEntries(); i++) {
363			cerr << i << " " << GoodMuons->At(i)->Pt() << " " << GoodMuons->At(i)->Eta()
364			<< " " << GoodMuons->At(i)->Phi() << endl;
365			}
366
367			cerr << "Central Jets" << endl;
368			for (UInt_t i = 0; i < GoodCentralJets->GetEntries(); i++) {
369			cerr << i << " " << GoodCentralJets->At(i)->Pt() << " "
370			<< GoodCentralJets->At(i)->Eta() << " " << GoodCentralJets->At(i)->Phi() << endl;
371			}
372
373			cerr << "MET" << endl;
374			cerr << met->Pt() << " " << met->Eta() << " "
375			<< met->Phi() << endl;
376			}
377
378			//Save Objects for Other Modules to use
379			AddObjThisEvt(GoodElectrons, fGoodElectronsName.Data());
380			AddObjThisEvt(GoodMuons, fGoodMuonsName.Data());
381			AddObjThisEvt(GoodCentralJets, fGoodCentralJetsName.Data());
382			AddObjThisEvt(GoodForwardJets, fGoodForwardJetsName.Data());
383			AddObjThisEvt(GoodCentralSCJets, fGoodCentralSCJetsName.Data());
384			AddObjThisEvt(GoodForwardSCJets, fGoodForwardSCJetsName.Data());
385
386	sixie	1.1	}
387
388
389			//--------------------------------------------------------------------------------------------------
390			void ObjectCleaningMod::SlaveBegin()
391			{
392			// Run startup code on the computer (slave) doing the actual analysis. Here,
393			// we typically initialize histograms and other analysis objects and request
394			// branches. For this module, we request a branch of the MitTree.
395
396			ReqBranch(fMCPartName, fParticles);
397			ReqBranch(fMuonName, fMuons);
398			ReqBranch(fElectronName, fElectrons);
399			ReqBranch(fJetName, fJets);
400	sixie	1.6	ReqBranch(fSCJetName, fSCJets);
401	sixie	1.1	ReqBranch(fMetName, fMet);
402
403
404	sixie	1.2	fAllMuonPtHist = new TH1D("hAllMuonPtHist",";p_{t};#",200,0.,200.);
405			fAllMuonEtaHist = new TH1D("hAllMuonEtaHist",";#eta;#",100,-5.,5.);
406			fGoodMuonPtHist = new TH1D("hGoodMuonPtHist",";p_{t};#",200,0.,200.);
407	sixie	1.1	fGoodMuonEtaHist = new TH1D("hGoodMuonEtaHist",";#eta;#",21,-5.,5.);
408			fMuonSelection = new TH1D("hMuonSelection", ";MuonSelection;#",4,-1.5,2.5 ) ;
409			AddOutput(fAllMuonPtHist);
410			AddOutput(fAllMuonEtaHist);
411			AddOutput(fGoodMuonPtHist);
412			AddOutput(fGoodMuonEtaHist);
413			AddOutput(fMuonSelection);
414	sixie	1.6
415	sixie	1.1	fAllElectronPtHist = new TH1D("hAllElectronPtHist",";p_{t};#",100,0.,200.);
416	sixie	1.2	fAllElectronEtaHist = new TH1D("hAllElectronEtaHist",";#eta;#",100,-5.,5.);
417	sixie	1.1	fGoodElectronPtHist = new TH1D("hGoodElectronPtHist",";p_{t};#",25,0.,200.);
418			fGoodElectronEtaHist = new TH1D("hGoodElectronEtaHist",";#eta;#",21,-5.,5.);
419			fGoodElectronClassification = new TH1D("hGoodElectronClassification",
420			";Good Electron Classification;#",51,0,50 ) ;
421			AddOutput(fAllElectronPtHist);
422			AddOutput(fAllElectronEtaHist);
423			AddOutput(fGoodElectronPtHist);
424			AddOutput(fGoodElectronEtaHist);
425			AddOutput(fGoodElectronClassification);
426
427			//Jet Plots
428			fAllJetPtHist = new TH1D("hAllJetPtHist",";All Jet p_{t};#",100,0.,200.);
429			fAllJetEtaHist = new TH1D("hAllJetEtaHist",";All Jet #eta;#",160,-8.,8.);
430			fCentralJetSelection = new TH1D("hCentralJetSelection",
431			";CentralJetSelection;#",5,-1.5,3.5 ) ;
432			AddOutput(fAllJetPtHist);
433			AddOutput(fAllJetEtaHist);
434			AddOutput(fCentralJetSelection);
435	sixie	1.6
436	sixie	1.1	//MET Plots
437	sixie	1.6	fMetPtHist = new TH1D("hMetPtHist",";p_{t};#",30,0.,300.);
438			fMetPhiHist = new TH1D("hMetPhiHist",";#phi;#",28,-3.5,3.5);
439			AddOutput(fMetPtHist);
440			AddOutput(fMetPhiHist);
441
442	sixie	1.1
443			}
444
445			//--------------------------------------------------------------------------------------------------
446			void ObjectCleaningMod::SlaveTerminate()
447			{
448			// Run finishing code on the computer (slave) that did the analysis. For this
449			// module, we dont do anything here.
450
451			}
452
453			//--------------------------------------------------------------------------------------------------
454			void ObjectCleaningMod::Terminate()
455			{
456			// Run finishing code on the client computer. For this module, we dont do
457			// anything here.
458			}