Utils/src/VertexTools.cc

// $Id: VertexTools.cc,v 1.3 2011/07/01 13:49:01 maxi Exp $

#include "MitPhysics/Utils/interface/VertexTools.h"
#include "MitPhysics/Utils/interface/ElectronTools.h"
#include "MitAna/DataTree/interface/StableData.h"
#include <TFile.h>
#include <TVector3.h>

ClassImp(mithep::VertexTools)

using namespace mithep;

VertexTools* VertexTools::meobject = NULL;

//--------------------------------------------------------------------------------------------------
VertexTools::VertexTools(const char* str)  
{
  // Constructor.
  relname = str;
  reader = new TMVA::Reader( "!Color:!Silent" );    
  reader->AddVariable( "var1", &tmvar1 );
  reader->AddVariable( "var2", &tmvar2 );
  reader->AddVariable( "var3", &tmvar3 );
  reader->AddVariable( "var4", &tmvar4 );
  reader->AddVariable( "var5", &tmvar5 );
  reader->AddVariable( "var6", &tmvar6 ); 
  reader->BookMVA( "BDTG method",relname + TString("/src/MitPhysics/data/TMVAClassification_BDTG.weights.xml").Data());
  reader->BookMVA( "BDTD method",relname + TString("/src/MitPhysics/data/TMVAClassification_BDTD.weights.xml" ).Data());
  //reader->BookMVA( "CFMlpANN method", "/home/maxi/cms/root/TMVAClassification_CFMlpANN.weights.xml" );
  reader->BookMVA( "MLP method", relname + TString("/src/MitPhysics/data/TMVAClassification_MLP.weights.xml").Data());
  reader->BookMVA( "MLPBFGS method",relname + TString("/src/MitPhysics/data/TMVAClassification_MLPBFGS.weights.xml" ).Data());
}

double VertexTools::NewMass(const Photon* ph1, const Photon* ph2, const BaseVertex* vert)
{
  ThreeVector drv1 = (ThreeVector(ph1->CaloPos()) - vert->Position()).Unit();
  //ThreeVector drv1 = (ThreeVector(ph1->SCluster()->Point()) - vert->Position()).Unit();
  FourVector pho1c(drv1.X()*ph1->E(),drv1.Y()*ph1->E(),drv1.Z()*ph1->E(),ph1->E());
  ThreeVector drv2 = (ThreeVector(ph2->CaloPos()) - vert->Position()).Unit();
  FourVector pho2c(drv2.X()*ph2->E(),drv2.Y()*ph2->E(),drv2.Z()*ph2->E(),ph2->E());

  FourVector diboso = pho1c+pho2c;
  return diboso.M();
}

//--------------------------------------------------------------------------------------------------
VertexZarray VertexTools::ExtractZarray(const VertexCol* vcol, float zmin, float zmax,  const BaseVertex  *fBeamSpot)
{
  VertexZarray zs;
  if(vcol == NULL) return zs;

  for(unsigned vv = 0; vv < vcol->GetEntries(); vv++){    
    const Vertex* vert = vcol->At(vv);
    double zpos = vert->Z();
    if(fBeamSpot != NULL) 
      zpos = zpos - fBeamSpot->Z();
    if(zpos > zmin && zpos > zmin)
      zs.push_back(zpos);
  }
  return zs;
}

VertexZarray VertexTools::ExtractZarray(float zmin, float zmax, float step)
{
  VertexZarray zs;
  for(float zpos = zmin; zpos < zmax+step; zpos = zpos+step){
    zs.push_back(zpos);
  }
  return zs;
}

const Vertex* VertexTools::BestVtx( const PFCandidateCol *fPFJets, const VertexCol *c, 
                                    const BaseVertex  *fBeamSpot, FourVector diboso) {
  
  if (!c || !c->GetEntries()) return NULL;
  
  double bestprob = -100.;
  const Vertex* bestvert = NULL;
  for(unsigned vv = 0; vv < c->GetEntries(); vv++){
    
    const Vertex* vert = c->At(vv);
    double zpos = vert->Z();     
    double prob = Prob(fPFJets, zpos, fBeamSpot, diboso);
    if(prob > bestprob){
      bestprob = prob;
      bestvert = vert;
    }
  }
  return bestvert;
}

double VertexTools::BestVtx( const PFCandidateCol *fPFJets, VertexZarray zcol,
                             const BaseVertex  *fBeamSpot, FourVector diboso){
  
  double bestprob = -100.;
  double bestz = -100.;
  for(unsigned vv = 0; vv < zcol.size(); vv++){
    double zpos = zcol[vv];     
    double prob = Prob(fPFJets, zpos, fBeamSpot, diboso);
    if(prob > bestprob){
      bestprob = prob;
      bestz = zpos;
    }
  }
  return bestz;
  
}

double VertexTools::Prob(const PFCandidateCol *fPFJets, double zpos,  
                         const BaseVertex  *fBeamSpot, FourVector diboso){
  
  double bosophi = diboso.Phi();
  double bosopt = diboso.Pt();
  
  Vertex* vert = new Vertex(0,0,zpos);
  double ZVC = vert->Z()-fBeamSpot->Z();
  VertexTools* vtool = VertexTools::instance("");
  
  Double_t sinsum = 0.;
  Double_t cossum = 0.;
  Double_t sumpt = 0.;
  Double_t ntrks = 0.;
  Double_t ntplus = 0.;
  Double_t ortminus = 0.;
  Double_t ortplus = 0.;
  Double_t bdplus = 0.;
  Double_t bdminus = 0.;
  Double_t zmean = 0.;
  Double_t zmeansq = 0.;
  Double_t ww = 0.;
  
  for(unsigned pfj = 0; pfj < fPFJets->GetEntries(); pfj++){
    const PFCandidate* pfca = fPFJets->At(pfj);
    if(! pfca->HasTrk()) continue;
    if(!(pfca->PFType() == PFCandidate::eX ||
         pfca->PFType() == PFCandidate::eHadron || 
         pfca->PFType() == PFCandidate::eElectron || 
         pfca->PFType() == PFCandidate::eMuon) ) continue;
    const Track *t = pfca->Trk();
    if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 0.2) continue;
    
    if(pfca->Pt()<0.3 || pfca->Pt()>200) continue;

    std::vector<const Track*>::iterator itt;
    itt = find ((vtool->excluded).begin(), (vtool->excluded).end(), t);
    if(itt != (vtool->excluded).end()) continue;
    
    zmean = zmean + pfca->Pt()* t->DzCorrected(*fBeamSpot);
    zmeansq = zmeansq + pfca->Pt()*t->DzCorrected(*fBeamSpot)*t->DzCorrected(*fBeamSpot);
    ntrks++;
    sumpt  = sumpt + pfca->Pt()*pfca->Pt();
    ww = ww + pfca->Pt();
    
    sinsum = sinsum + pfca->Pt()*TMath::Sin(t->Phi());
    cossum = cossum + pfca->Pt()*TMath::Cos(t->Phi());
  }
  if(ntrks < 2 || !(sumpt > 0.) ) return 0;
  
  Double_t phim = TMath::ATan2(sinsum,cossum);
  zmean = zmean/ww;
  zmeansq = sqrt(zmeansq/ww - zmean*zmean);
    
  //--------------------
  Double_t bosoproj = 0.;
  Float_t ymean = 0.;
  Float_t ymsq = 0.;
  for(unsigned pfj = 0; pfj < fPFJets->GetEntries(); pfj++){
    const PFCandidate* pfca = fPFJets->At(pfj);
    if(! pfca->HasTrk()) continue;
    if(!(pfca->PFType() == PFCandidate::eX ||
         pfca->PFType() == PFCandidate::eHadron ||
         pfca->PFType() == PFCandidate::eElectron ||
         pfca->PFType() == PFCandidate::eMuon) ) continue;
    
    const Track *t = pfca->Trk();
    if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 0.2) continue;
    //if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 3*zwidth) continue;
    if(pfca->Pt()<0.3 || pfca->Pt()>200) continue;

    std::vector<const Track*>::iterator itt;
    itt = find ((vtool->excluded).begin(), (vtool->excluded).end(), t);
    if(itt != (vtool->excluded).end()) continue;
    
    //Float_t phid = phim - t->Phi();
    Float_t phid = bosophi+3.14 - t->Phi();
    ymean = ymean + pfca->Pt()*TMath::Sin(phid);
    ymsq = ymsq + pow(pfca->Pt()*TMath::Sin(phid),2);
    bosoproj = bosoproj + pfca->Pt()*TMath::Cos(bosophi-t->Phi());
    if(TMath::Sin(phid) > 0.) ortplus = ortplus + pfca->Pt()*TMath::Sin(phid);
    else                      ortminus = ortminus + pfca->Pt()*fabs(TMath::Sin(phid));
    if(TMath::Cos(phid) > 0.) bdplus = bdplus + pfca->Pt()*TMath::Cos(phid);
    else                      bdminus = bdminus + pfca->Pt()*fabs(TMath::Cos(phid));
    if(TMath::Cos(phid) > 0.) ntplus = ntplus + TMath::Cos(phid);
  }
  
  Float_t phimsq = sqrt(ymsq/ntrks - pow(ymean/ntrks,2))*180/3.14;
  // not used... commented out by Fabian
  //Double_t A2n = (ortplus+ortminus) > 0.01 ? (bdplus+bdminus)/(ortplus+ortminus) : (bdplus+bdminus)/0.01;
  Double_t A1n = bdplus/bosopt;
  Double_t A3n = ntplus > 0 ? bdplus/ntplus : 0.;
  //-------------------
  Double_t angle = 180/3.14*TVector2::Phi_0_2pi(bosophi-phim);
  
  vtool->tmvar1 = ntrks;
  vtool->tmvar2 = sumpt;
  vtool->tmvar3 = A1n;
  vtool->tmvar4 = angle;
  vtool->tmvar5 = phimsq;
  vtool->tmvar6 = A3n;
  
  double p1 = vtool->reader->GetProba ( "BDTD method" );
  double p2 = vtool->reader->GetProba ( "BDTG method" );
  double p3 = vtool->reader->GetProba ( "MLP method");
  double p4 = vtool->reader->GetProba ( "MLPBFGS method");
  
  double retval = p1*p2*p3*p4;
  
  return retval;
}

double VertexTools::VertexWidth(const Vertex* vert,  const BaseVertex  *fBeamSpot){
  double width = 0.;
  double zmean = 0.;
  double zmeansq = 0.;
  double ww = 0.;
  if(vert == NULL) return width;
  for(unsigned i = 0; i < vert->NTracks(); i++){
    const Track *t = vert->Trk(i);
    if(t->Pt() < 0.3 ) continue;
    if(t->Pt() > 500.) continue;
    ww = ww + t->Pt();
    zmean = zmean + t->Pt()*t->DzCorrected(*fBeamSpot);
    zmeansq = zmeansq + t->Pt()*t->DzCorrected(*fBeamSpot)*t->DzCorrected(*fBeamSpot);
  }
  if( !(ww > 0.) ) return width;
  zmean = zmean/ww;

  width = sqrt(zmeansq/ww - zmean*zmean);
  return width;
  
}

void VertexTools::BanThisTrack(const Track* track){
  VertexTools* vtool = VertexTools::instance("");
  (vtool->excluded).push_back(track);
}
void VertexTools::Reset(){
  VertexTools* vtool = VertexTools::instance("");
  (vtool->excluded).clear();
}
Revision:	1.4
Committed:	Mon Jul 4 13:48:50 2011 UTC (13 years, 10 months ago) by fabstoec
Content type:	text/plain
Branch:	MAIN
Changes since 1.3:	+3 -2 lines
Log Message:	fixed warnings
#	User	Rev	Content
1	fabstoec	1.4	// $Id: VertexTools.cc,v 1.3 2011/07/01 13:49:01 maxi Exp $
2	bendavid	1.1
3			#include "MitPhysics/Utils/interface/VertexTools.h"
4			#include "MitPhysics/Utils/interface/ElectronTools.h"
5			#include "MitAna/DataTree/interface/StableData.h"
6			#include <TFile.h>
7	maxi	1.2	#include <TVector3.h>
8	bendavid	1.1
9			ClassImp(mithep::VertexTools)
10
11			using namespace mithep;
12
13	maxi	1.2	VertexTools* VertexTools::meobject = NULL;
14
15	bendavid	1.1	//--------------------------------------------------------------------------------------------------
16	maxi	1.2	VertexTools::VertexTools(const char* str)
17	bendavid	1.1	{
18			// Constructor.
19	maxi	1.2	relname = str;
20			reader = new TMVA::Reader( "!Color:!Silent" );
21			reader->AddVariable( "var1", &tmvar1 );
22			reader->AddVariable( "var2", &tmvar2 );
23			reader->AddVariable( "var3", &tmvar3 );
24			reader->AddVariable( "var4", &tmvar4 );
25			reader->AddVariable( "var5", &tmvar5 );
26			reader->AddVariable( "var6", &tmvar6 );
27			reader->BookMVA( "BDTG method",relname + TString("/src/MitPhysics/data/TMVAClassification_BDTG.weights.xml").Data());
28			reader->BookMVA( "BDTD method",relname + TString("/src/MitPhysics/data/TMVAClassification_BDTD.weights.xml" ).Data());
29			//reader->BookMVA( "CFMlpANN method", "/home/maxi/cms/root/TMVAClassification_CFMlpANN.weights.xml" );
30			reader->BookMVA( "MLP method", relname + TString("/src/MitPhysics/data/TMVAClassification_MLP.weights.xml").Data());
31			reader->BookMVA( "MLPBFGS method",relname + TString("/src/MitPhysics/data/TMVAClassification_MLPBFGS.weights.xml" ).Data());
32			}
33
34			double VertexTools::NewMass(const Photon* ph1, const Photon* ph2, const BaseVertex* vert)
35			{
36	maxi	1.3	ThreeVector drv1 = (ThreeVector(ph1->CaloPos()) - vert->Position()).Unit();
37			//ThreeVector drv1 = (ThreeVector(ph1->SCluster()->Point()) - vert->Position()).Unit();
38	maxi	1.2	FourVector pho1c(drv1.X()ph1->E(),drv1.Y()ph1->E(),drv1.Z()*ph1->E(),ph1->E());
39	maxi	1.3	ThreeVector drv2 = (ThreeVector(ph2->CaloPos()) - vert->Position()).Unit();
40	maxi	1.2	FourVector pho2c(drv2.X()ph2->E(),drv2.Y()ph2->E(),drv2.Z()*ph2->E(),ph2->E());
41
42			FourVector diboso = pho1c+pho2c;
43			return diboso.M();
44	bendavid	1.1	}
45
46			//--------------------------------------------------------------------------------------------------
47	maxi	1.2	VertexZarray VertexTools::ExtractZarray(const VertexCol* vcol, float zmin, float zmax, const BaseVertex *fBeamSpot)
48			{
49			VertexZarray zs;
50			if(vcol == NULL) return zs;
51
52			for(unsigned vv = 0; vv < vcol->GetEntries(); vv++){
53			const Vertex* vert = vcol->At(vv);
54			double zpos = vert->Z();
55			if(fBeamSpot != NULL)
56			zpos = zpos - fBeamSpot->Z();
57			if(zpos > zmin && zpos > zmin)
58			zs.push_back(zpos);
59			}
60			return zs;
61			}
62
63			VertexZarray VertexTools::ExtractZarray(float zmin, float zmax, float step)
64			{
65			VertexZarray zs;
66			for(float zpos = zmin; zpos < zmax+step; zpos = zpos+step){
67			zs.push_back(zpos);
68			}
69			return zs;
70			}
71
72			const Vertex* VertexTools::BestVtx( const PFCandidateCol fPFJets, const VertexCol c,
73			const BaseVertex *fBeamSpot, FourVector diboso) {
74
75			if (!c \|\| !c->GetEntries()) return NULL;
76
77			double bestprob = -100.;
78			const Vertex* bestvert = NULL;
79			for(unsigned vv = 0; vv < c->GetEntries(); vv++){
80
81			const Vertex* vert = c->At(vv);
82			double zpos = vert->Z();
83			double prob = Prob(fPFJets, zpos, fBeamSpot, diboso);
84			if(prob > bestprob){
85			bestprob = prob;
86			bestvert = vert;
87			}
88			}
89			return bestvert;
90			}
91
92			double VertexTools::BestVtx( const PFCandidateCol *fPFJets, VertexZarray zcol,
93			const BaseVertex *fBeamSpot, FourVector diboso){
94
95			double bestprob = -100.;
96			double bestz = -100.;
97			for(unsigned vv = 0; vv < zcol.size(); vv++){
98			double zpos = zcol[vv];
99			double prob = Prob(fPFJets, zpos, fBeamSpot, diboso);
100			if(prob > bestprob){
101			bestprob = prob;
102			bestz = zpos;
103			}
104			}
105			return bestz;
106
107			}
108	bendavid	1.1
109	maxi	1.2	double VertexTools::Prob(const PFCandidateCol *fPFJets, double zpos,
110			const BaseVertex *fBeamSpot, FourVector diboso){
111
112			double bosophi = diboso.Phi();
113			double bosopt = diboso.Pt();
114
115			Vertex* vert = new Vertex(0,0,zpos);
116			double ZVC = vert->Z()-fBeamSpot->Z();
117	maxi	1.3	VertexTools* vtool = VertexTools::instance("");
118	maxi	1.2
119			Double_t sinsum = 0.;
120			Double_t cossum = 0.;
121			Double_t sumpt = 0.;
122			Double_t ntrks = 0.;
123			Double_t ntplus = 0.;
124			Double_t ortminus = 0.;
125			Double_t ortplus = 0.;
126			Double_t bdplus = 0.;
127			Double_t bdminus = 0.;
128			Double_t zmean = 0.;
129			Double_t zmeansq = 0.;
130			Double_t ww = 0.;
131
132			for(unsigned pfj = 0; pfj < fPFJets->GetEntries(); pfj++){
133			const PFCandidate* pfca = fPFJets->At(pfj);
134			if(! pfca->HasTrk()) continue;
135			if(!(pfca->PFType() == PFCandidate::eX \|\|
136			pfca->PFType() == PFCandidate::eHadron \|\|
137			pfca->PFType() == PFCandidate::eElectron \|\|
138			pfca->PFType() == PFCandidate::eMuon) ) continue;
139			const Track *t = pfca->Trk();
140			if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 0.2) continue;
141
142			if(pfca->Pt()<0.3 \|\| pfca->Pt()>200) continue;
143	maxi	1.3
144			std::vector<const Track*>::iterator itt;
145			itt = find ((vtool->excluded).begin(), (vtool->excluded).end(), t);
146			if(itt != (vtool->excluded).end()) continue;
147	maxi	1.2
148			zmean = zmean + pfca->Pt()* t->DzCorrected(*fBeamSpot);
149			zmeansq = zmeansq + pfca->Pt()t->DzCorrected(fBeamSpot)t->DzCorrected(fBeamSpot);
150			ntrks++;
151			sumpt = sumpt + pfca->Pt()*pfca->Pt();
152			ww = ww + pfca->Pt();
153
154			sinsum = sinsum + pfca->Pt()*TMath::Sin(t->Phi());
155			cossum = cossum + pfca->Pt()*TMath::Cos(t->Phi());
156			}
157			if(ntrks < 2 \|\| !(sumpt > 0.) ) return 0;
158
159			Double_t phim = TMath::ATan2(sinsum,cossum);
160			zmean = zmean/ww;
161			zmeansq = sqrt(zmeansq/ww - zmean*zmean);
162
163			//--------------------
164			Double_t bosoproj = 0.;
165			Float_t ymean = 0.;
166			Float_t ymsq = 0.;
167			for(unsigned pfj = 0; pfj < fPFJets->GetEntries(); pfj++){
168			const PFCandidate* pfca = fPFJets->At(pfj);
169			if(! pfca->HasTrk()) continue;
170			if(!(pfca->PFType() == PFCandidate::eX \|\|
171			pfca->PFType() == PFCandidate::eHadron \|\|
172			pfca->PFType() == PFCandidate::eElectron \|\|
173			pfca->PFType() == PFCandidate::eMuon) ) continue;
174
175			const Track *t = pfca->Trk();
176			if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 0.2) continue;
177			//if(fabs(t->DzCorrected(fBeamSpot) - ZVC ) > 3zwidth) continue;
178			if(pfca->Pt()<0.3 \|\| pfca->Pt()>200) continue;
179	maxi	1.3
180			std::vector<const Track*>::iterator itt;
181			itt = find ((vtool->excluded).begin(), (vtool->excluded).end(), t);
182			if(itt != (vtool->excluded).end()) continue;
183	maxi	1.2
184			//Float_t phid = phim - t->Phi();
185			Float_t phid = bosophi+3.14 - t->Phi();
186			ymean = ymean + pfca->Pt()*TMath::Sin(phid);
187			ymsq = ymsq + pow(pfca->Pt()*TMath::Sin(phid),2);
188			bosoproj = bosoproj + pfca->Pt()*TMath::Cos(bosophi-t->Phi());
189			if(TMath::Sin(phid) > 0.) ortplus = ortplus + pfca->Pt()*TMath::Sin(phid);
190			else ortminus = ortminus + pfca->Pt()*fabs(TMath::Sin(phid));
191			if(TMath::Cos(phid) > 0.) bdplus = bdplus + pfca->Pt()*TMath::Cos(phid);
192			else bdminus = bdminus + pfca->Pt()*fabs(TMath::Cos(phid));
193			if(TMath::Cos(phid) > 0.) ntplus = ntplus + TMath::Cos(phid);
194			}
195	bendavid	1.1
196	maxi	1.2	Float_t phimsq = sqrt(ymsq/ntrks - pow(ymean/ntrks,2))*180/3.14;
197	fabstoec	1.4	// not used... commented out by Fabian
198			//Double_t A2n = (ortplus+ortminus) > 0.01 ? (bdplus+bdminus)/(ortplus+ortminus) : (bdplus+bdminus)/0.01;
199	maxi	1.2	Double_t A1n = bdplus/bosopt;
200			Double_t A3n = ntplus > 0 ? bdplus/ntplus : 0.;
201			//-------------------
202			Double_t angle = 180/3.14*TVector2::Phi_0_2pi(bosophi-phim);
203	bendavid	1.1
204	maxi	1.2	vtool->tmvar1 = ntrks;
205			vtool->tmvar2 = sumpt;
206			vtool->tmvar3 = A1n;
207			vtool->tmvar4 = angle;
208			vtool->tmvar5 = phimsq;
209			vtool->tmvar6 = A3n;
210
211			double p1 = vtool->reader->GetProba ( "BDTD method" );
212			double p2 = vtool->reader->GetProba ( "BDTG method" );
213			double p3 = vtool->reader->GetProba ( "MLP method");
214			double p4 = vtool->reader->GetProba ( "MLPBFGS method");
215
216			double retval = p1p2p3*p4;
217
218			return retval;
219			}
220
221			double VertexTools::VertexWidth(const Vertex* vert, const BaseVertex *fBeamSpot){
222			double width = 0.;
223			double zmean = 0.;
224			double zmeansq = 0.;
225			double ww = 0.;
226			if(vert == NULL) return width;
227			for(unsigned i = 0; i < vert->NTracks(); i++){
228			const Track *t = vert->Trk(i);
229			if(t->Pt() < 0.3 ) continue;
230			if(t->Pt() > 500.) continue;
231			ww = ww + t->Pt();
232			zmean = zmean + t->Pt()t->DzCorrected(fBeamSpot);
233			zmeansq = zmeansq + t->Pt()t->DzCorrected(fBeamSpot)t->DzCorrected(fBeamSpot);
234			}
235			if( !(ww > 0.) ) return width;
236			zmean = zmean/ww;
237
238			width = sqrt(zmeansq/ww - zmean*zmean);
239			return width;
240
241			}
242	maxi	1.3
243			void VertexTools::BanThisTrack(const Track* track){
244			VertexTools* vtool = VertexTools::instance("");
245			(vtool->excluded).push_back(track);
246			}
247			void VertexTools::Reset(){
248			VertexTools* vtool = VertexTools::instance("");
249			(vtool->excluded).clear();
250			}