Utils/src/VertexTools.cc

// $Id: VertexTools.cc,v 1.1 2011/05/16 13:26:48 bendavid 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->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->SCluster()->Point()) - 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();
  
  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;
    
    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;
    
    //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;
  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);
  VertexTools* vtool = VertexTools::instance("");
  
  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;
  
}
Revision:	1.2
Committed:	Thu Jun 16 02:53:02 2011 UTC (13 years, 10 months ago) by maxi
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+212 -6 lines
Log Message:	first version
#	User	Rev	Content
1	maxi	1.2	// $Id: VertexTools.cc,v 1.1 2011/05/16 13:26:48 bendavid 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			ThreeVector drv1 = (ThreeVector(ph1->SCluster()->Point()) - vert->Position()).Unit();
37			FourVector pho1c(drv1.X()ph1->E(),drv1.Y()ph1->E(),drv1.Z()*ph1->E(),ph1->E());
38			ThreeVector drv2 = (ThreeVector(ph2->SCluster()->Point()) - vert->Position()).Unit();
39			FourVector pho2c(drv2.X()ph2->E(),drv2.Y()ph2->E(),drv2.Z()*ph2->E(),ph2->E());
40
41			FourVector diboso = pho1c+pho2c;
42			return diboso.M();
43	bendavid	1.1	}
44
45			//--------------------------------------------------------------------------------------------------
46	maxi	1.2	VertexZarray VertexTools::ExtractZarray(const VertexCol* vcol, float zmin, float zmax, const BaseVertex *fBeamSpot)
47			{
48			VertexZarray zs;
49			if(vcol == NULL) return zs;
50
51			for(unsigned vv = 0; vv < vcol->GetEntries(); vv++){
52			const Vertex* vert = vcol->At(vv);
53			double zpos = vert->Z();
54			if(fBeamSpot != NULL)
55			zpos = zpos - fBeamSpot->Z();
56			if(zpos > zmin && zpos > zmin)
57			zs.push_back(zpos);
58			}
59			return zs;
60			}
61
62			VertexZarray VertexTools::ExtractZarray(float zmin, float zmax, float step)
63			{
64			VertexZarray zs;
65			for(float zpos = zmin; zpos < zmax+step; zpos = zpos+step){
66			zs.push_back(zpos);
67			}
68			return zs;
69			}
70
71			const Vertex* VertexTools::BestVtx( const PFCandidateCol fPFJets, const VertexCol c,
72			const BaseVertex *fBeamSpot, FourVector diboso) {
73
74			if (!c \|\| !c->GetEntries()) return NULL;
75
76			double bestprob = -100.;
77			const Vertex* bestvert = NULL;
78			for(unsigned vv = 0; vv < c->GetEntries(); vv++){
79
80			const Vertex* vert = c->At(vv);
81			double zpos = vert->Z();
82			double prob = Prob(fPFJets, zpos, fBeamSpot, diboso);
83			if(prob > bestprob){
84			bestprob = prob;
85			bestvert = vert;
86			}
87			}
88			return bestvert;
89			}
90
91			double VertexTools::BestVtx( const PFCandidateCol *fPFJets, VertexZarray zcol,
92			const BaseVertex *fBeamSpot, FourVector diboso){
93
94			double bestprob = -100.;
95			double bestz = -100.;
96			for(unsigned vv = 0; vv < zcol.size(); vv++){
97			double zpos = zcol[vv];
98			double prob = Prob(fPFJets, zpos, fBeamSpot, diboso);
99			if(prob > bestprob){
100			bestprob = prob;
101			bestz = zpos;
102			}
103			}
104			return bestz;
105
106			}
107	bendavid	1.1
108	maxi	1.2	double VertexTools::Prob(const PFCandidateCol *fPFJets, double zpos,
109			const BaseVertex *fBeamSpot, FourVector diboso){
110
111			double bosophi = diboso.Phi();
112			double bosopt = diboso.Pt();
113
114			Vertex* vert = new Vertex(0,0,zpos);
115			double ZVC = vert->Z()-fBeamSpot->Z();
116
117			Double_t sinsum = 0.;
118			Double_t cossum = 0.;
119			Double_t sumpt = 0.;
120			Double_t ntrks = 0.;
121			Double_t ntplus = 0.;
122			Double_t ortminus = 0.;
123			Double_t ortplus = 0.;
124			Double_t bdplus = 0.;
125			Double_t bdminus = 0.;
126			Double_t zmean = 0.;
127			Double_t zmeansq = 0.;
128			Double_t ww = 0.;
129
130			for(unsigned pfj = 0; pfj < fPFJets->GetEntries(); pfj++){
131			const PFCandidate* pfca = fPFJets->At(pfj);
132			if(! pfca->HasTrk()) continue;
133			if(!(pfca->PFType() == PFCandidate::eX \|\|
134			pfca->PFType() == PFCandidate::eHadron \|\|
135			pfca->PFType() == PFCandidate::eElectron \|\|
136			pfca->PFType() == PFCandidate::eMuon) ) continue;
137			const Track *t = pfca->Trk();
138			if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 0.2) continue;
139
140			if(pfca->Pt()<0.3 \|\| pfca->Pt()>200) continue;
141
142			zmean = zmean + pfca->Pt()* t->DzCorrected(*fBeamSpot);
143			zmeansq = zmeansq + pfca->Pt()t->DzCorrected(fBeamSpot)t->DzCorrected(fBeamSpot);
144			ntrks++;
145			sumpt = sumpt + pfca->Pt()*pfca->Pt();
146			ww = ww + pfca->Pt();
147
148			sinsum = sinsum + pfca->Pt()*TMath::Sin(t->Phi());
149			cossum = cossum + pfca->Pt()*TMath::Cos(t->Phi());
150			}
151			if(ntrks < 2 \|\| !(sumpt > 0.) ) return 0;
152
153			Double_t phim = TMath::ATan2(sinsum,cossum);
154			zmean = zmean/ww;
155			zmeansq = sqrt(zmeansq/ww - zmean*zmean);
156
157			//--------------------
158			Double_t bosoproj = 0.;
159			Float_t ymean = 0.;
160			Float_t ymsq = 0.;
161			for(unsigned pfj = 0; pfj < fPFJets->GetEntries(); pfj++){
162			const PFCandidate* pfca = fPFJets->At(pfj);
163			if(! pfca->HasTrk()) continue;
164			if(!(pfca->PFType() == PFCandidate::eX \|\|
165			pfca->PFType() == PFCandidate::eHadron \|\|
166			pfca->PFType() == PFCandidate::eElectron \|\|
167			pfca->PFType() == PFCandidate::eMuon) ) continue;
168
169			const Track *t = pfca->Trk();
170			if(fabs(t->DzCorrected(*fBeamSpot) - ZVC ) > 0.2) continue;
171			//if(fabs(t->DzCorrected(fBeamSpot) - ZVC ) > 3zwidth) continue;
172			if(pfca->Pt()<0.3 \|\| pfca->Pt()>200) continue;
173
174			//Float_t phid = phim - t->Phi();
175			Float_t phid = bosophi+3.14 - t->Phi();
176			ymean = ymean + pfca->Pt()*TMath::Sin(phid);
177			ymsq = ymsq + pow(pfca->Pt()*TMath::Sin(phid),2);
178			bosoproj = bosoproj + pfca->Pt()*TMath::Cos(bosophi-t->Phi());
179			if(TMath::Sin(phid) > 0.) ortplus = ortplus + pfca->Pt()*TMath::Sin(phid);
180			else ortminus = ortminus + pfca->Pt()*fabs(TMath::Sin(phid));
181			if(TMath::Cos(phid) > 0.) bdplus = bdplus + pfca->Pt()*TMath::Cos(phid);
182			else bdminus = bdminus + pfca->Pt()*fabs(TMath::Cos(phid));
183			if(TMath::Cos(phid) > 0.) ntplus = ntplus + TMath::Cos(phid);
184			}
185	bendavid	1.1
186	maxi	1.2	Float_t phimsq = sqrt(ymsq/ntrks - pow(ymean/ntrks,2))*180/3.14;
187			Double_t A2n = (ortplus+ortminus) > 0.01 ? (bdplus+bdminus)/(ortplus+ortminus) : (bdplus+bdminus)/0.01;
188			Double_t A1n = bdplus/bosopt;
189			Double_t A3n = ntplus > 0 ? bdplus/ntplus : 0.;
190			//-------------------
191			Double_t angle = 180/3.14*TVector2::Phi_0_2pi(bosophi-phim);
192			VertexTools* vtool = VertexTools::instance("");
193	bendavid	1.1
194	maxi	1.2	vtool->tmvar1 = ntrks;
195			vtool->tmvar2 = sumpt;
196			vtool->tmvar3 = A1n;
197			vtool->tmvar4 = angle;
198			vtool->tmvar5 = phimsq;
199			vtool->tmvar6 = A3n;
200
201			double p1 = vtool->reader->GetProba ( "BDTD method" );
202			double p2 = vtool->reader->GetProba ( "BDTG method" );
203			double p3 = vtool->reader->GetProba ( "MLP method");
204			double p4 = vtool->reader->GetProba ( "MLPBFGS method");
205
206			double retval = p1p2p3*p4;
207
208			return retval;
209			}
210
211			double VertexTools::VertexWidth(const Vertex* vert, const BaseVertex *fBeamSpot){
212			double width = 0.;
213			double zmean = 0.;
214			double zmeansq = 0.;
215			double ww = 0.;
216			if(vert == NULL) return width;
217			for(unsigned i = 0; i < vert->NTracks(); i++){
218			const Track *t = vert->Trk(i);
219			if(t->Pt() < 0.3 ) continue;
220			if(t->Pt() > 500.) continue;
221			ww = ww + t->Pt();
222			zmean = zmean + t->Pt()t->DzCorrected(fBeamSpot);
223			zmeansq = zmeansq + t->Pt()t->DzCorrected(fBeamSpot)t->DzCorrected(fBeamSpot);
224			}
225			if( !(ww > 0.) ) return width;
226			zmean = zmean/ww;
227
228			width = sqrt(zmeansq/ww - zmean*zmean);
229			return width;
230
231			}