Mods/src/MVASystematicsMod.cc

// $Id: MVASystematicsMod.cc,v 1.1 2011/12/13 21:13:23 bendavid Exp $

#include <TMath.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TNtuple.h>
#include <TRandom3.h>
#include "MitCommon/MathTools/interface/MathUtils.h"
#include "MitAna/DataUtil/interface/Debug.h"
#include "MitAna/DataTree/interface/Names.h"
#include "MitAna/DataTree/interface/TrackCol.h"
#include "MitAna/DataTree/interface/VertexCol.h"
#include "MitAna/DataTree/interface/PhotonCol.h"
#include "MitPhysics/Mods/interface/MVASystematicsMod.h"
#include "MitAna/DataTree/interface/BeamSpotCol.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"

#include "MitPhysics/Utils/interface/IsolationTools.h"
#include "MitPhysics/Utils/interface/PhotonTools.h"

#include <iostream>
#include <sstream>

using namespace mithep;

ClassImp(mithep::MVASystematicsMod)

//--------------------------------------------------------------------------------------------------
MVASystematicsMod::MVASystematicsMod(const char *name, const char *title) :
  BaseMod             (name,title),

  // define all the Branches to load
  fMCParticleName    (Names::gkMCPartBrn),
  fPVName            (Names::gkPVBeamSpotBrn),
  fEBSCName            (Names::gkBarrelSuperClusterBrn),
  fEESCName            (Names::gkEndcapSuperClusterBrn),  


  // ----------------------------------------
  // collections....
  fMCParticles       (0),
  fPV       (0),
  fEBSC       (0),
  fEESC       (0),
  // -------------------------------------------
  
  fMCR9ScaleEB(1.0),
  fMCR9ScaleEE(1.0),
  
  fIsData            (false),

  fTupleName("hMVAtuple")
  
{
  // Constructor.
}

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

//--------------------------------------------------------------------------------------------------
void MVASystematicsMod::Process()
{  
  IncNEventsProcessed();
  if( !fIsData )
    LoadBranch(fMCParticleName);

  LoadBranch(fPVName);
  LoadBranch(fEBSCName);
  LoadBranch(fEESCName);
  
  const Vertex *vtx = 0;
  if (fPV->GetEntries()>0) vtx = fPV->At(0);
  
  const MCParticle *h = 0;
  const MCParticle *p1 = 0;
  const MCParticle *p2 = 0;
  
  const SuperCluster *sc1 = 0;
  const SuperCluster *sc2 = 0;
  
  Float_t _pth    = -100.;
  Float_t _y      = -100.;
  Float_t _genmass = -100.;
  if( !fIsData ) h = FindHiggsPtAndY(_pth, _y, _genmass);  

  if (h && h->NDaughters()>=2) {
    p1 = h->Daughter(0);
    p2 = h->Daughter(1);
  }
  
  Float_t pt1 = -100;
  Float_t eta1 = -100;
  Float_t phi1 = -100;

  Float_t pt2 = -100;
  Float_t eta2 = -100;
  Float_t phi2 = -100;  
  
  Float_t scet1 = -100;
  Float_t sceta1 = -100;
  Float_t scphi1 = -100;
  Float_t scr91 = -100;
  bool iseb1 = kFALSE;

  Float_t scet2 = -100;
  Float_t sceta2 = -100;
  Float_t scphi2 = -100;  
  Float_t scr92 = -100;  
  bool iseb2 = kFALSE;
  
  if (p1) {
    pt1 = p1->Pt();
    eta1 = p1->Eta();
    phi1 = p1->Phi();
    sc1 = MatchSC(p1,iseb1);
  }

  if (p2) {
    pt2 = p2->Pt();
    eta2 = p2->Eta();
    phi2 = p2->Phi();
    sc2 = MatchSC(p2,iseb2);
  }
  
  if (sc1) {
    double r9scale;
    if   (iseb1) r9scale = fMCR9ScaleEB;
    else r9scale = fMCR9ScaleEE;
      
    // compute the probe momentum wr to the chosen Vtx
    FourVectorM scmom;  
    ThreeVectorC scp;
    if (vtx) scp = sc1->Point() - vtx->Position();
    else scp = sc1->Point();
    scp = scp/scp.R();
    scmom.SetXYZT(sc1->Energy()*scp.X(), sc1->Energy()*scp.Y(), sc1->Energy()*scp.Z(), sc1->Energy());
      
    scet1 = scmom.Pt();
    sceta1 = sc1->Eta();
    scphi1 = sc1->Phi();
    scr91 = r9scale*sc1->R9();
    
  }


  if (sc2) {
    double r9scale;
    if   (iseb2) r9scale = fMCR9ScaleEB;
    else r9scale = fMCR9ScaleEE;
      
    // compute the probe momentum wr to the chosen Vtx
    FourVectorM scmom;  
    ThreeVectorC scp;
    if (vtx) scp = sc2->Point() - vtx->Position();
    else scp = sc2->Point();
    scp = scp/scp.R();
    scmom.SetXYZT(sc2->Energy()*scp.X(), sc2->Energy()*scp.Y(), sc2->Energy()*scp.Z(), sc2->Energy());
      
    scet2 = scmom.Pt();
    sceta2 = sc2->Eta();
    scphi2 = sc2->Phi();
    scr92 = r9scale*sc2->R9();
    
  }

  Float_t fill[] = {_pth, _y, _genmass,pt1,eta1,phi1,pt2,eta2,phi2,scet1,sceta1,scphi1,scr91,scet2,sceta2,scphi2,scr92};

  hMVAtuple->Fill(fill);


  return;
}

//--------------------------------------------------------------------------------------------------
void MVASystematicsMod::SlaveBegin()
{
  // Run startup code on the computer (slave) doing the actual analysis. Here,
  // we typically initialize histograms and other analysis objects and request
  // branches or objects created by earlier modules.

  if (!fIsData) {
    ReqBranch(fMCParticleName,        fMCParticles);
  }
  
  ReqBranch(fPVName,        fPV);
  ReqBranch(fEBSCName,      fEBSC);
  ReqBranch(fEESCName,      fEESC);
  

  hMVAtuple = new TNtuple(fTupleName.Data(),fTupleName.Data(),"hpt:hy:hm:pt1:eta1:phi1:pt2:eta2:phi2:scet1:sceta1:scphi1:scr91:scet2:sceta2:scphi2:scr92");

  AddOutput(hMVAtuple);
}

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

//--------------------------------------------------------------------------------------------------
void MVASystematicsMod::Terminate()
{
  // Run finishing code on the client computer. For this module, we dont do
  // anything here.
}

// ----------------------------------------------------------------------------------------
// some helpfer functions....
const MCParticle *MVASystematicsMod::FindHiggsPtAndY(Float_t& pt, Float_t& Y, Float_t& mass) {

  const MCParticle *h = 0;
  
  pt = -999.;
  Y  = -999;
  mass = -999.;

  // loop over all GEN particles and look for status 1 photons
  for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
    const MCParticle* p = fMCParticles->At(i);
    if( p->Is(MCParticle::kH) ) {
      pt=p->Pt();
      Y = p->Rapidity();
      mass = p->Mass();
      h = p;
      break;
    }
  }
  
  return h;
}

// ----------------------------------------------------------------------------------------
// some helpfer functions....
const SuperCluster *MVASystematicsMod::MatchSC(const MCParticle *p, bool &iseb) {

  iseb = kFALSE;
  
  for(UInt_t i=0; i<fEBSC->GetEntries(); ++i) {
    iseb = kTRUE;
    const SuperCluster *sc = fEBSC->At(i);
    if (MathUtils::DeltaR(sc,p)<0.2) return sc;
  }

  for(UInt_t i=0; i<fEESC->GetEntries(); ++i) {
    iseb = kFALSE;
    const SuperCluster *sc = fEESC->At(i);
    if (MathUtils::DeltaR(sc,p)<0.2) return sc;
  }

  return 0;
}
Revision:	1.2
Committed:	Thu Mar 22 15:54:07 2012 UTC (13 years, 1 month ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_029c, Mit_029b, Mit_029a, Mit_028a, Mit_028, Mit_027, Mit_027a, Mit_025e, Mit_025d, HEAD
Changes since 1.1:	+143 -7 lines
Log Message:	final synchronization for moriond mva analysis
#	User	Rev	Content
1	bendavid	1.2	// $Id: MVASystematicsMod.cc,v 1.1 2011/12/13 21:13:23 bendavid Exp $
2	bendavid	1.1
3			#include <TMath.h>
4			#include <TH1D.h>
5			#include <TH2D.h>
6			#include <TNtuple.h>
7			#include <TRandom3.h>
8			#include "MitCommon/MathTools/interface/MathUtils.h"
9			#include "MitAna/DataUtil/interface/Debug.h"
10			#include "MitAna/DataTree/interface/Names.h"
11			#include "MitAna/DataTree/interface/TrackCol.h"
12			#include "MitAna/DataTree/interface/VertexCol.h"
13			#include "MitAna/DataTree/interface/PhotonCol.h"
14			#include "MitPhysics/Mods/interface/MVASystematicsMod.h"
15			#include "MitAna/DataTree/interface/BeamSpotCol.h"
16			#include "FWCore/ParameterSet/interface/FileInPath.h"
17
18			#include "MitPhysics/Utils/interface/IsolationTools.h"
19			#include "MitPhysics/Utils/interface/PhotonTools.h"
20
21			#include <iostream>
22			#include <sstream>
23
24			using namespace mithep;
25
26			ClassImp(mithep::MVASystematicsMod)
27
28			//--------------------------------------------------------------------------------------------------
29			MVASystematicsMod::MVASystematicsMod(const char name, const char title) :
30			BaseMod (name,title),
31
32			// define all the Branches to load
33			fMCParticleName (Names::gkMCPartBrn),
34	bendavid	1.2	fPVName (Names::gkPVBeamSpotBrn),
35			fEBSCName (Names::gkBarrelSuperClusterBrn),
36			fEESCName (Names::gkEndcapSuperClusterBrn),
37
38	bendavid	1.1
39			// ----------------------------------------
40			// collections....
41			fMCParticles (0),
42	bendavid	1.2	fPV (0),
43			fEBSC (0),
44			fEESC (0),
45	bendavid	1.1	// -------------------------------------------
46	bendavid	1.2
47			fMCR9ScaleEB(1.0),
48			fMCR9ScaleEE(1.0),
49
50	bendavid	1.1	fIsData (false),
51
52			fTupleName("hMVAtuple")
53
54			{
55			// Constructor.
56			}
57
58			//--------------------------------------------------------------------------------------------------
59			void MVASystematicsMod::Begin()
60			{
61			// Run startup code on the client machine. For this module, we dont do anything here.
62			}
63
64			//--------------------------------------------------------------------------------------------------
65			void MVASystematicsMod::Process()
66			{
67			IncNEventsProcessed();
68			if( !fIsData )
69			LoadBranch(fMCParticleName);
70
71	bendavid	1.2	LoadBranch(fPVName);
72			LoadBranch(fEBSCName);
73			LoadBranch(fEESCName);
74
75			const Vertex *vtx = 0;
76			if (fPV->GetEntries()>0) vtx = fPV->At(0);
77
78			const MCParticle *h = 0;
79			const MCParticle *p1 = 0;
80			const MCParticle *p2 = 0;
81
82			const SuperCluster *sc1 = 0;
83			const SuperCluster *sc2 = 0;
84
85	bendavid	1.1	Float_t _pth = -100.;
86			Float_t _y = -100.;
87			Float_t _genmass = -100.;
88	bendavid	1.2	if( !fIsData ) h = FindHiggsPtAndY(_pth, _y, _genmass);
89
90			if (h && h->NDaughters()>=2) {
91			p1 = h->Daughter(0);
92			p2 = h->Daughter(1);
93			}
94
95			Float_t pt1 = -100;
96			Float_t eta1 = -100;
97			Float_t phi1 = -100;
98
99			Float_t pt2 = -100;
100			Float_t eta2 = -100;
101			Float_t phi2 = -100;
102
103			Float_t scet1 = -100;
104			Float_t sceta1 = -100;
105			Float_t scphi1 = -100;
106			Float_t scr91 = -100;
107			bool iseb1 = kFALSE;
108
109			Float_t scet2 = -100;
110			Float_t sceta2 = -100;
111			Float_t scphi2 = -100;
112			Float_t scr92 = -100;
113			bool iseb2 = kFALSE;
114
115			if (p1) {
116			pt1 = p1->Pt();
117			eta1 = p1->Eta();
118			phi1 = p1->Phi();
119			sc1 = MatchSC(p1,iseb1);
120			}
121
122			if (p2) {
123			pt2 = p2->Pt();
124			eta2 = p2->Eta();
125			phi2 = p2->Phi();
126			sc2 = MatchSC(p2,iseb2);
127			}
128
129			if (sc1) {
130			double r9scale;
131			if (iseb1) r9scale = fMCR9ScaleEB;
132			else r9scale = fMCR9ScaleEE;
133
134			// compute the probe momentum wr to the chosen Vtx
135			FourVectorM scmom;
136			ThreeVectorC scp;
137			if (vtx) scp = sc1->Point() - vtx->Position();
138			else scp = sc1->Point();
139			scp = scp/scp.R();
140			scmom.SetXYZT(sc1->Energy()scp.X(), sc1->Energy()scp.Y(), sc1->Energy()*scp.Z(), sc1->Energy());
141
142			scet1 = scmom.Pt();
143			sceta1 = sc1->Eta();
144			scphi1 = sc1->Phi();
145			scr91 = r9scale*sc1->R9();
146
147			}
148
149	bendavid	1.1
150	bendavid	1.2	if (sc2) {
151			double r9scale;
152			if (iseb2) r9scale = fMCR9ScaleEB;
153			else r9scale = fMCR9ScaleEE;
154
155			// compute the probe momentum wr to the chosen Vtx
156			FourVectorM scmom;
157			ThreeVectorC scp;
158			if (vtx) scp = sc2->Point() - vtx->Position();
159			else scp = sc2->Point();
160			scp = scp/scp.R();
161			scmom.SetXYZT(sc2->Energy()scp.X(), sc2->Energy()scp.Y(), sc2->Energy()*scp.Z(), sc2->Energy());
162
163			scet2 = scmom.Pt();
164			sceta2 = sc2->Eta();
165			scphi2 = sc2->Phi();
166			scr92 = r9scale*sc2->R9();
167
168			}
169
170			Float_t fill[] = {_pth, _y, _genmass,pt1,eta1,phi1,pt2,eta2,phi2,scet1,sceta1,scphi1,scr91,scet2,sceta2,scphi2,scr92};
171
172			hMVAtuple->Fill(fill);
173	bendavid	1.1
174
175			return;
176			}
177
178			//--------------------------------------------------------------------------------------------------
179			void MVASystematicsMod::SlaveBegin()
180			{
181			// Run startup code on the computer (slave) doing the actual analysis. Here,
182			// we typically initialize histograms and other analysis objects and request
183			// branches or objects created by earlier modules.
184
185			if (!fIsData) {
186			ReqBranch(fMCParticleName, fMCParticles);
187			}
188	bendavid	1.2
189			ReqBranch(fPVName, fPV);
190			ReqBranch(fEBSCName, fEBSC);
191			ReqBranch(fEESCName, fEESC);
192
193	bendavid	1.1
194	bendavid	1.2	hMVAtuple = new TNtuple(fTupleName.Data(),fTupleName.Data(),"hpt:hy:hm:pt1:eta1:phi1:pt2:eta2:phi2:scet1:sceta1:scphi1:scr91:scet2:sceta2:scphi2:scr92");
195	bendavid	1.1
196			AddOutput(hMVAtuple);
197			}
198
199			//--------------------------------------------------------------------------------------------------
200			void MVASystematicsMod::SlaveTerminate()
201			{
202			// Run finishing code on the computer (slave) that did the analysis. For this
203			// module, we dont do anything here.
204			}
205
206			//--------------------------------------------------------------------------------------------------
207			void MVASystematicsMod::Terminate()
208			{
209			// Run finishing code on the client computer. For this module, we dont do
210			// anything here.
211			}
212
213			// ----------------------------------------------------------------------------------------
214			// some helpfer functions....
215	bendavid	1.2	const MCParticle *MVASystematicsMod::FindHiggsPtAndY(Float_t& pt, Float_t& Y, Float_t& mass) {
216	bendavid	1.1
217	bendavid	1.2	const MCParticle *h = 0;
218
219	bendavid	1.1	pt = -999.;
220			Y = -999;
221			mass = -999.;
222
223			// loop over all GEN particles and look for status 1 photons
224			for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
225			const MCParticle* p = fMCParticles->At(i);
226			if( p->Is(MCParticle::kH) ) {
227			pt=p->Pt();
228			Y = p->Rapidity();
229			mass = p->Mass();
230	bendavid	1.2	h = p;
231	bendavid	1.1	break;
232			}
233			}
234
235	bendavid	1.2	return h;
236			}
237
238			// ----------------------------------------------------------------------------------------
239			// some helpfer functions....
240			const SuperCluster MVASystematicsMod::MatchSC(const MCParticle p, bool &iseb) {
241
242			iseb = kFALSE;
243
244			for(UInt_t i=0; i<fEBSC->GetEntries(); ++i) {
245			iseb = kTRUE;
246			const SuperCluster *sc = fEBSC->At(i);
247			if (MathUtils::DeltaR(sc,p)<0.2) return sc;
248			}
249
250			for(UInt_t i=0; i<fEESC->GetEntries(); ++i) {
251			iseb = kFALSE;
252			const SuperCluster *sc = fEESC->At(i);
253			if (MathUtils::DeltaR(sc,p)<0.2) return sc;
254			}
255
256			return 0;
257			}