CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes
TopologyWorker Class Reference

#include <TopQuarkAnalysis/TopTools/interface/TopologyWorker.h>

Public Member Functions

void clear (void)
 
double get_aplanarity ()
 
double get_centrality ()
 
double get_et0 ()
 
double get_et56 ()
 
double get_h10 ()
 
double get_h20 ()
 
double get_h30 ()
 
double get_h40 ()
 
double get_h50 ()
 
double get_h60 ()
 
double get_ht ()
 
double get_ht3 ()
 
double get_njetW ()
 
double get_sphericity ()
 
double get_sqrts ()
 
int getFast ()
 
double getThMomPower ()
 
TVector3 majorAxis ()
 
TVector3 minorAxis ()
 
double oblateness ()
 
void planes_sphe (double &pnorm, double &p2, double &p3)
 
void planes_sphe_wei (double &pnorm, double &p2, double &p3)
 
void planes_thrust (double &pnorm, double &p2, double &p3)
 
void setFast (int nf)
 
void setPartList (TObjArray *e1, TObjArray *e2)
 
void setThMomPower (double tp)
 
void setVerbose (bool loud)
 
void sumangles (float &sdeta, float &sdr)
 
TVector3 thrust ()
 
TVector3 thrustAxis ()
 
 TopologyWorker ()
 
 TopologyWorker (bool boost)
 
virtual ~TopologyWorker ()
 

Private Member Functions

double CalcEta (int i)
 
double CalcEta2 (int i)
 
void CalcHTstuff ()
 
double CalcPt (int i)
 
double CalcPt2 (int i)
 
void CalcSqrts ()
 
void CalcWmul ()
 
void fowo ()
 
void getetaphi (double px, double py, double pz, double &eta, double &phi)
 
int iPow (int man, int exp)
 
void ludbrb (TMatrix *mom, double the, double phi, double bx, double by, double bz)
 
void sanda ()
 
double sign (double a, double b)
 
double ulAngle (double x, double y)
 

Private Attributes

double m_apl
 
bool m_boost
 
double m_centrality
 
TMatrix m_dAxes
 
double m_dConv
 
double m_dDeltaThPower
 
double m_dOblateness
 
double m_dSphMomPower
 
double m_dThrust [4]
 
double m_et0
 
double m_et56
 
bool m_fowo_called
 
double m_h10
 
double m_h20
 
double m_h30
 
double m_h40
 
double m_h50
 
double m_h60
 
double m_ht
 
double m_ht3
 
int m_iFast
 
int m_iGood
 
TVector3 m_MajorAxis
 
TVector3 m_MinorAxis
 
TMatrix m_mom
 
TMatrix m_mom2
 
double m_njetsweighed
 
int m_np
 
int m_np2
 
TRandom m_random
 
bool m_sanda_called
 
double m_sph
 
double m_sqrts
 
bool m_sumangles_called
 
TVector3 m_Thrust
 
TVector3 m_ThrustAxis
 
bool m_verbose
 

Static Private Attributes

static int m_maxpart = 1000
 

Detailed Description

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>=""> This class contains the topological methods as used in D0 (all hadronic) analyses.

Definition at line 26 of file TopologyWorker.h.

Constructor & Destructor Documentation

TopologyWorker::TopologyWorker ( )
inline

Definition at line 29 of file TopologyWorker.h.

29 {;}
TopologyWorker::TopologyWorker ( bool  boost)

Definition at line 7 of file TopologyWorker.cc.

References m_apl, m_boost, m_dAxes, m_et0, m_et56, m_fowo_called, m_h10, m_h20, m_h30, m_h40, m_ht, m_ht3, m_maxpart, m_mom, m_mom2, m_njetsweighed, m_np, m_np2, m_sanda_called, m_sph, m_sqrts, m_sumangles_called, and m_verbose.

7  :
8  m_dSphMomPower(2.0),m_dDeltaThPower(0),
9  m_iFast(4),m_dConv(0.0001),m_iGood(2)
10 {
11  m_dAxes.ResizeTo(4,4);
12  m_mom.ResizeTo(m_maxpart,6);
13  m_mom2.ResizeTo(m_maxpart,6);
14  m_np=-1;
15  m_np2=-1;
16  m_sanda_called=false;
17  m_fowo_called=false;
18  m_sumangles_called=false;
19  m_verbose=false;
20  m_boost=boost;
21  m_sph=-1;
22  m_apl=-1;
23  m_h10=-1;
24  m_h20=-1;
25  m_h30=-1;
26  m_h40=-1;
27  m_sqrts=0;
28  m_ht=0;
29  m_ht3=0;
30  m_et56=0;
31  m_njetsweighed=0;
32  m_et0=0;
33 }
TopologyWorker::m_dSphMomPower
double m_dSphMomPower
Definition: TopologyWorker.h:88
TopologyWorker::m_maxpart
static int m_maxpart
Definition: TopologyWorker.h:146
TopologyWorker::m_dDeltaThPower
double m_dDeltaThPower
Definition: TopologyWorker.h:91
virtual TopologyWorker::~TopologyWorker ( )
inlinevirtual

Definition at line 31 of file TopologyWorker.h.

31 {;}

Member Function Documentation

double TopologyWorker::CalcEta ( int  i)
inlineprivate

Definition at line 153 of file TopologyWorker.h.

References eta(), getetaphi(), m_mom, and phi.

153 {double eta, phi;getetaphi(m_mom(i,1),m_mom(i,2),m_mom(i,3),eta,phi); return eta;}
i
int i
Definition: DBlmapReader.cc:9
eta
T eta() const
Definition: Basic3DVectorLD.h:168
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.cc:1408
phi
Definition: DDAxes.h:10
double TopologyWorker::CalcEta2 ( int  i)
inlineprivate

Definition at line 154 of file TopologyWorker.h.

References eta(), getetaphi(), m_mom2, and phi.

154 {double eta, phi; getetaphi(m_mom2(i,1),m_mom2(i,2),m_mom2(i,3),eta,phi); return eta;}
i
int i
Definition: DBlmapReader.cc:9
eta
T eta() const
Definition: Basic3DVectorLD.h:168
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.cc:1408
phi
Definition: DDAxes.h:10
void TopologyWorker::CalcHTstuff ( )
private

Definition at line 1508 of file TopologyWorker.cc.

References CalcPt(), CalcPt2(), h, i, j, m_centrality, m_et0, m_et56, m_ht, m_ht3, m_mom, m_np, m_np2, and mathSSE::sqrt().

Referenced by setPartList().

1508  {
1509  m_ht=0;
1510  m_ht3=0;
1511  m_et56=0;
1512  m_et0=0;
1513  double ptlead=0;
1514  double h=0;
1515  for(int i=0; i< m_np; i++){
1516  //cout << i << "/" << m_np << ":" << CalcPt(i) << endl;
1517  m_ht+=CalcPt(i);
1518  h+=m_mom(i,4);
1519  if(i>1)
1520  m_ht3+=CalcPt(i);
1521  if(i==5)
1522  m_et56=sqrt(CalcPt(i)*CalcPt(i-1));
1523  }
1524 
1525  for(int j=0; j< m_np2; j++){
1526  //cout << j << "/" << m_np2 << ":" << CalcPt2(j) << endl;
1527  if(ptlead<CalcPt2(j))
1528  ptlead=CalcPt2(j);
1529 
1530  }
1531  if(m_ht>0.0001){
1532  m_et0=ptlead/m_ht;
1533  //cout << "calculating ETO" << m_et0 << "=" << ptlead << endl;
1534  }
1535  if(h>0.00001)
1536  m_centrality=m_ht/h;
1537 }
i
int i
Definition: DBlmapReader.cc:9
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
TopologyWorker::CalcPt
double CalcPt(int i)
Definition: TopologyWorker.h:151
j
int j
Definition: DBlmapReader.cc:9
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
TopologyWorker::CalcPt2
double CalcPt2(int i)
Definition: TopologyWorker.h:152
double TopologyWorker::CalcPt ( int  i)
inlineprivate

Definition at line 151 of file TopologyWorker.h.

References m_mom, funct::pow(), and mathSSE::sqrt().

Referenced by CalcHTstuff(), and CalcWmul().

151 { return sqrt(pow(m_mom(i,1),2)+pow(m_mom(i,2),2));}
i
int i
Definition: DBlmapReader.cc:9
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
double TopologyWorker::CalcPt2 ( int  i)
inlineprivate

Definition at line 152 of file TopologyWorker.h.

References m_mom2, funct::pow(), and mathSSE::sqrt().

Referenced by CalcHTstuff().

152 { return sqrt(pow(m_mom2(i,1),2)+pow(m_mom2(i,2),2));}
i
int i
Definition: DBlmapReader.cc:9
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void TopologyWorker::CalcSqrts ( )
private

Definition at line 1492 of file TopologyWorker.cc.

References relval_parameters_module::energy, event(), i, m_mom, m_np, m_sqrts, funct::pow(), and mathSSE::sqrt().

Referenced by setPartList().

1492  {
1493  TLorentzVector event(0,0,0,0);
1494  TLorentzVector worker(0,0,0,0);
1495 
1496  for(int i=0; i< m_np; i++){
1497  double energy=m_mom(i,4);
1498  if(m_mom(i,4)<0.00001)
1499  energy=sqrt(pow(m_mom(i,1),2)+pow(m_mom(i,2),2)+pow(m_mom(i,3),2));
1500  // assume massless particle if only TVector3s are provided...
1501  worker.SetXYZT(m_mom(i,1),m_mom(i,2),m_mom(i,3),energy);
1502  event+=worker;
1503  }
1504  m_sqrts=event.M();
1505 }
i
int i
Definition: DBlmapReader.cc:9
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
event
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void TopologyWorker::CalcWmul ( )
private

Definition at line 1467 of file TopologyWorker.cc.

References CalcPt(), m_njetsweighed, m_np, and query::result.

Referenced by setPartList().

1467  {
1468 
1469  Int_t njets = m_np;
1470  double result=0;
1471  for(Int_t ijet=0; ijet<njets-1; ijet++){
1472  double emin=55;
1473  double emax=55;
1474  if(CalcPt(ijet)<55)
1475  emax=CalcPt(ijet);
1476  if(CalcPt(ijet+1)<55)
1477  emin=CalcPt(ijet+1);
1478  result+=0.5 * (emax*emax-emin*emin)*(ijet+1);
1479  }
1480  double elo=15;
1481  if(CalcPt(njets-1)>elo){
1482  elo=CalcPt(njets-1);
1483  }
1484 
1485  result+=0.5 * (elo*elo-(15*15))*(njets);
1486  result/=((55*55)-100)/2.0;
1487  m_njetsweighed=result;
1488 }
TopologyWorker::CalcPt
double CalcPt(int i)
Definition: TopologyWorker.h:151
query.result
tuple result
Definition: query.py:137
void TopologyWorker::clear ( void  )
inline

Definition at line 33 of file TopologyWorker.h.

References m_np, and m_np2.

Referenced by python.Vispa.Views.BoxDecayView.BoxDecayView::closeEvent(), python.Vispa.Views.LineDecayView.LineDecayView::setDataObjects(), and python.Vispa.Views.BoxDecayView.BoxDecayView::updateContent().

33 {m_np=0;m_np2=0;return;}
void TopologyWorker::fowo ( )
private

Definition at line 1270 of file TopologyWorker.cc.

References Exhume::I, m_fowo_called, m_h10, m_h20, m_h30, m_h40, m_h50, m_h60, m_mom, m_np, MultiGaussianStateTransform::N, P, funct::pow(), and mathSSE::sqrt().

Referenced by get_h10(), get_h20(), get_h30(), get_h40(), get_h50(), and get_h60().

1270  {
1271 // 20020830 changed: from p/E to Et/Ettot and include H50 and H60
1272  m_fowo_called=true;
1273  // include fox wolframs
1274  float H10=-1;
1275  float H20=-1;
1276  float H30=-1;
1277  float H40=-1;
1278  float H50=-1;
1279  float H60=-1;
1280  if (1==1) {
1281  float P[1000][6],H0,HD,CTHE;
1282  int N,NP,I,J,I1,I2;
1283  H0=HD=0.;
1284  N=m_np;
1285  NP=0;
1286  for (I=1;I<N+1;I++){
1287  NP++; // start at one
1288  P[NP][1]=m_mom(I-1,1) ;
1289  P[NP][2]=m_mom(I-1,2) ;
1290  P[NP][3]=m_mom(I-1,3) ;
1291  P[NP][4]=m_mom(I-1,4) ;
1292  P[NP][5]=m_mom(I-1,5) ;
1293  }
1294 
1295  N=NP;
1296  NP=0;
1297 
1298  for (I=1;I<N+1;I++) {
1299  NP=NP+1;
1300  for (J=1;J<5;J++) {
1301  P[N+NP][J]=P[I][J];
1302  }
1303 // p/E
1304  P[N+NP][4]=sqrt(pow(P[I][1],2)+pow(P[I][2],2)+pow(P[I][3],2));
1305 // Et/Ettot
1306  P[N+NP][5]=sqrt(pow(P[I][1],2)+pow(P[I][2],2));
1307  H0=H0+P[N+NP][5];
1308  HD=HD+pow(P[N+NP][5],2);
1309  }
1310  H0=H0*H0;
1311 
1312 
1313 
1314  // Very low multiplicities (0 or 1) not considered.
1315  if (NP<2) {
1316  H10=-1.;
1317  H20=-1.;
1318  H30=-1.;
1319  H40=-1.;
1320  H50=-1.;
1321  H60=-1.;
1322  return;
1323  }
1324 
1325  // Calculate H1 - H4.
1326  H10=0.;
1327  H20=0.;
1328  H30=0.;
1329  H40=0.;
1330  H50=0.;
1331  H60=0.;
1332  for (I1=N+1;I1<N+NP+1;I1++) { //130
1333  for (I2=I1+1;I2<N+NP+1;I2++) { // 120
1334  CTHE=(P[I1][1]*P[I2][1]+P[I1][2]*P[I2][2]+P[I1][3]*P[I2][3])/
1335  (P[I1][4]*P[I2][4]);
1336  H10=H10+P[I1][5]*P[I2][5]*CTHE;
1337  double C2=(1.5*CTHE*CTHE-0.5);
1338  H20=H20+P[I1][5]*P[I2][5]*C2;
1339  double C3=(2.5*CTHE*CTHE*CTHE-1.5*CTHE);
1340  H30=H30+P[I1][5]*P[I2][5]*C3;
1341  // use recurrence
1342  double C4=(7*CTHE*C3 - 3*C2)/4.;
1343  double C5=(9*CTHE*C4 - 4*C3)/5.;
1344  double C6=(11*CTHE*C5 - 5*C4)/6.;
1345 // H40=H40+P[I1][5]*P[I2][5]*(4.375*pow(CTHE,4)-3.75*CTHE*CTHE+0.375);
1346 // H50=H50+P[I1][5]*P[I2][5]*
1347 // (63./8.*pow(CTHE,5)-70./8.*CTHE*CTHE*CTHE+15./8.*CTHE);
1348 // H60=H60+P[I1][5]*P[I2][5]*
1349 // (231/16.*pow(CTHE,6)-315./16.*CTHE*CTHE*CTHE*CTHE+105./16.*CTHE*CTHE-5./16.);
1350  H40=H40+P[I1][5]*P[I2][5]*C4;
1351  H50=H50+P[I1][5]*P[I2][5]*C5;
1352  H60=H60+P[I1][5]*P[I2][5]*C6;
1353  } // 120
1354  } // 130
1355 
1356  H10=(HD+2.*H10)/H0;
1357  H20=(HD+2.*H20)/H0;
1358  H30=(HD+2.*H30)/H0;
1359  H40=(HD+2.*H40)/H0;
1360  H50=(HD+2.*H50)/H0;
1361  H60=(HD+2.*H60)/H0;
1362  m_h10=H10;
1363  m_h20=H20;
1364  m_h30=H30;
1365  m_h40=H40;
1366  m_h50=H50;
1367  m_h60=H60;
1368  }
1369 
1370 }
P
#define P
Definition: PFJetBenchmark.cc:33
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
Exhume::I
const std::complex< double > I
Definition: I.h:8
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
double TopologyWorker::get_aplanarity ( )

Definition at line 1403 of file TopologyWorker.cc.

References m_apl, m_sanda_called, and sanda().

1403  {
1404  if (!m_sanda_called) sanda();
1405  return m_apl;
1406 }
double TopologyWorker::get_centrality ( )
inline

Definition at line 72 of file TopologyWorker.h.

References m_centrality.

72 { return m_centrality;}
double TopologyWorker::get_et0 ( )
inline

Definition at line 68 of file TopologyWorker.h.

References m_et0.

68 {return m_et0;}
double TopologyWorker::get_et56 ( )
inline

Definition at line 71 of file TopologyWorker.h.

References m_et56.

71 {return m_et56;}
double TopologyWorker::get_h10 ( )

Definition at line 1372 of file TopologyWorker.cc.

References fowo(), m_fowo_called, and m_h10.

1372  {
1373  if (!m_fowo_called) fowo();
1374  return m_h10;
1375 }
double TopologyWorker::get_h20 ( )

Definition at line 1376 of file TopologyWorker.cc.

References fowo(), m_fowo_called, and m_h20.

1376  {
1377  if (!m_fowo_called) fowo();
1378  return m_h20;
1379 }
double TopologyWorker::get_h30 ( )

Definition at line 1380 of file TopologyWorker.cc.

References fowo(), m_fowo_called, and m_h30.

1380  {
1381  if (!m_fowo_called) fowo();
1382  return m_h30;
1383 }
double TopologyWorker::get_h40 ( )

Definition at line 1384 of file TopologyWorker.cc.

References fowo(), m_fowo_called, and m_h40.

1384  {
1385  if (!m_fowo_called) fowo();
1386  return m_h40;
1387 }
double TopologyWorker::get_h50 ( )

Definition at line 1389 of file TopologyWorker.cc.

References fowo(), m_fowo_called, and m_h50.

1389  {
1390  if (!m_fowo_called) fowo();
1391  return m_h50;
1392 }
double TopologyWorker::get_h60 ( )

Definition at line 1393 of file TopologyWorker.cc.

References fowo(), m_fowo_called, and m_h60.

1393  {
1394  if (!m_fowo_called) fowo();
1395  return m_h60;
1396 }
double TopologyWorker::get_ht ( )
inline

Definition at line 66 of file TopologyWorker.h.

References m_ht.

66 {return m_ht;}
double TopologyWorker::get_ht3 ( )
inline

Definition at line 67 of file TopologyWorker.h.

References m_ht3.

67 {return m_ht3;}
double TopologyWorker::get_njetW ( )
inline

Definition at line 70 of file TopologyWorker.h.

References m_njetsweighed.

70 {return m_njetsweighed;}
double TopologyWorker::get_sphericity ( )

Definition at line 1399 of file TopologyWorker.cc.

References m_sanda_called, m_sph, and sanda().

1399  {
1400  if (!m_sanda_called) sanda();
1401  return m_sph;
1402 }
double TopologyWorker::get_sqrts ( )
inline

Definition at line 69 of file TopologyWorker.h.

References m_sqrts.

69 {return m_sqrts;}
void TopologyWorker::getetaphi ( double  px,
double  py,
double  pz,
double &  eta,
double &  phi 
)
private

Definition at line 1408 of file TopologyWorker.cc.

References funct::log(), L1TEmulatorMonitor_cff::p, pi, mathSSE::sqrt(), and funct::tan().

Referenced by CalcEta(), CalcEta2(), planes_sphe(), planes_sphe_wei(), and sumangles().

1408  {
1409 
1410  double pi=3.1415927;
1411 
1412  double p=sqrt(px*px+py*py+pz*pz);
1413  // get eta and phi
1414  double th=pi/2.;
1415  if (p!=0) {
1416  th=acos(pz/p); // Theta
1417  }
1418  float thg=th;
1419  if (th<=0) {
1420  thg = pi + th;
1421  }
1422  eta=-9.;
1423  if (tan( thg/2.0 )>0.000001) {
1424  eta = -log( tan( thg/2.0 ) );
1425  }
1426  phi = atan2(py,px);
1427  if(phi<=0) phi += 2.0*pi;
1428  return;
1429 }
eta
T eta() const
Definition: Basic3DVectorLD.h:168
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
funct::log
Log< T >::type log(const T &t)
Definition: Log.h:22
L1TEmulatorMonitor_cff.p
tuple p
Definition: L1TEmulatorMonitor_cff.py:48
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
phi
Definition: DDAxes.h:10
Int_t TopologyWorker::getFast ( )

Definition at line 412 of file TopologyWorker.cc.

References m_iFast.

413 {
414  return m_iFast;
415 }
double TopologyWorker::getThMomPower ( )

Definition at line 398 of file TopologyWorker.cc.

References m_dDeltaThPower.

399 {
400  return 1.0 + m_dDeltaThPower;
401 }
TopologyWorker::m_dDeltaThPower
double m_dDeltaThPower
Definition: TopologyWorker.h:91
Int_t TopologyWorker::iPow ( int  man,
int  exp 
)
private

Definition at line 1455 of file TopologyWorker.cc.

References funct::exp(), and gen::k.

Referenced by setPartList().

1456 {
1457  Int_t ans = 1;
1458  for( Int_t k = 0; k < exp; k++)
1459  {
1460  ans = ans*man;
1461  }
1462  return ans;
1463 }
funct::exp
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
void TopologyWorker::ludbrb ( TMatrix *  mom,
double  the,
double  phi,
double  bx,
double  by,
double  bz 
)
private

Definition at line 483 of file TopologyWorker.cc.

References beta, i, j, gen::k, and runTheMatrix::np.

Referenced by setPartList().

489 {
490  // Ignore "zeroth" elements in rot,pr,dp.
491  // Trying to use physics-like notation.
492  TMatrix rot(4,4);
493  double pr[4];
494  double dp[5];
495  Int_t np = mom->GetNrows();
496  if ( the*the + phi*phi > 1.0E-20 )
497  {
498  rot(1,1) = TMath::Cos(the)*TMath::Cos(phi);
499  rot(1,2) = -TMath::Sin(phi);
500  rot(1,3) = TMath::Sin(the)*TMath::Cos(phi);
501  rot(2,1) = TMath::Cos(the)*TMath::Sin(phi);
502  rot(2,2) = TMath::Cos(phi);
503  rot(2,3) = TMath::Sin(the)*TMath::Sin(phi);
504  rot(3,1) = -TMath::Sin(the);
505  rot(3,2) = 0.0;
506  rot(3,3) = TMath::Cos(the);
507  for ( Int_t i = 0; i < np; i++ )
508  {
509  for ( Int_t j = 1; j < 4; j++ )
510  {
511  pr[j] = (*mom)(i,j);
512  (*mom)(i,j) = 0;
513  }
514  for ( Int_t jb = 1; jb < 4; jb++)
515  {
516  for ( Int_t k = 1; k < 4; k++)
517  {
518  (*mom)(i,jb) = (*mom)(i,jb) + rot(jb,k)*pr[k];
519  }
520  }
521  }
522  double beta = TMath::Sqrt( bx*bx + by*by + bz*bz );
523  if ( beta*beta > 1.0E-20 )
524  {
525  if ( beta > 0.99999999 )
526  {
527  //send message: boost too large, resetting to <~1.0.
528  bx = bx*(0.99999999/beta);
529  by = by*(0.99999999/beta);
530  bz = bz*(0.99999999/beta);
531  beta = 0.99999999;
532  }
533  double gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
534  for ( Int_t i = 0; i < np; i++ )
535  {
536  for ( Int_t j = 1; j < 5; j++ )
537  {
538  dp[j] = (*mom)(i,j);
539  }
540  double bp = bx*dp[1] + by*dp[2] + bz*dp[3];
541  double gbp = gamma*(gamma*bp/(1.0 + gamma) + dp[4]);
542  (*mom)(i,1) = dp[1] + gbp*bx;
543  (*mom)(i,2) = dp[2] + gbp*by;
544  (*mom)(i,3) = dp[3] + gbp*bz;
545  (*mom)(i,4) = gamma*(dp[4] + bp);
546  }
547  }
548  }
549  return;
550 }
beta
const double beta
Definition: CosmicMuonParameters.h:69
i
int i
Definition: DBlmapReader.cc:9
j
int j
Definition: DBlmapReader.cc:9
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
phi
Definition: DDAxes.h:10
TVector3 TopologyWorker::majorAxis ( )

Definition at line 426 of file TopologyWorker.cc.

References m_dAxes, and m_MajorAxis.

Referenced by planes_thrust().

426  {
427  TVector3 m_MajorAxis(m_dAxes(2,1),m_dAxes(2,2),m_dAxes(2,3));
428  return m_MajorAxis;
429 }
TopologyWorker::m_MajorAxis
TVector3 m_MajorAxis
Definition: TopologyWorker.h:111
TVector3 TopologyWorker::minorAxis ( )

Definition at line 432 of file TopologyWorker.cc.

References m_dAxes, and m_MinorAxis.

Referenced by planes_thrust().

432  {
433  TVector3 m_MinorAxis(m_dAxes(3,1),m_dAxes(3,2),m_dAxes(3,3));
434  return m_MinorAxis;
435 }
TopologyWorker::m_MinorAxis
TVector3 m_MinorAxis
Definition: TopologyWorker.h:112
double TopologyWorker::oblateness ( )

Definition at line 444 of file TopologyWorker.cc.

References m_dOblateness.

444  {
445  return m_dOblateness;
446 }
void TopologyWorker::planes_sphe ( double &  pnorm,
double &  p2,
double &  p3 
)

Definition at line 710 of file TopologyWorker.cc.

References a, b, trackerHits::c, funct::cos(), debug_cff::d1, eta(), funct::exp(), getetaphi(), Exhume::I, gen::k, m_mom, m_np, siStripFEDMonitor_P5_cff::Max, siStripFEDMonitor_P5_cff::Min, MultiGaussianStateTransform::N, P, p3, p4, phi, pi, funct::pow(), PWT, SGN, funct::sin(), and mathSSE::sqrt().

710  {
711  float SPH=-1;
712  float APL=-1;
713 // C...Calculate matrix to be diagonalized.
714  float P[1000][6];
715  double SM[4][4],SV[4][4];
716  double PA,PWT,PS,SQ,SR,SP,SMAX,SGN;
717  int NP;
718  int J, J1, J2, I, N, JA, JB, J3, JC, JB1, JB2;
719  JA=JB=JC=0;
720  double RL;
721  float rlu,rlu1;
722  //
723  // --- get the input form GTRACK arrays
724  //
725  N=m_np;
726  NP=0;
727  for (I=1;I<N+1;I++){
728  NP++; // start at one
729  P[NP][1]=m_mom(I-1,1) ;
730  P[NP][2]=m_mom(I-1,2) ;
731  P[NP][3]=m_mom(I-1,3) ;
732  P[NP][4]=m_mom(I-1,4) ;
733  P[NP][5]=0;
734  }
735  //
736  //---
737  //
738  N=NP;
739 
740  for (J1=1;J1<4;J1++) {
741  for (J2=J1;J2<4;J2++) {
742  SM[J1][J2]=0.;
743  }
744  }
745  PS=0.;
746  for (I=1;I<N+1;I++) { // 140
747  PA=sqrt(pow(P[I][1],2)+pow(P[I][2],2)+pow(P[I][3],2));
748  double eta,phi;
749  getetaphi(P[I][1],P[I][2],P[I][3],eta,phi);
750  PWT=exp(-std::fabs(eta));
751  PWT=1.;
752  for (J1=1;J1<4;J1++) { // 130
753  for (J2=J1;J2<4;J2++) { // 120
754  SM[J1][J2]=SM[J1][J2]+PWT*P[I][J1]*P[I][J2];
755  }
756  } // 130
757  PS=PS+PWT*PA*PA;
758  } //140
759 // C...Very low multiplicities (0 or 1) not considered.
760  if(NP<2) {
761  SPH=-1.;
762  APL=-1.;
763  return;
764  }
765  for (J1=1;J1<4;J1++) { // 160
766  for (J2=J1;J2<4;J2++) { // 150
767  SM[J1][J2]=SM[J1][J2]/PS;
768  }
769  } // 160
770 // C...Find eigenvalues to matrix (third degree equation).
771  SQ=(SM[1][1]*SM[2][2]+SM[1][1]*SM[3][3]+SM[2][2]*SM[3][3]
772  -pow(SM[1][2],2)
773  -pow(SM[1][3],2)-pow(SM[2][3],2))/3.-1./9.;
774  SR=-0.5*(SQ+1./9.+SM[1][1]*pow(SM[2][3],2)+SM[2][2]*pow(SM[1][3],2)+SM[3][3]*
775  pow(SM[1][2],2)-SM[1][1]*SM[2][2]*SM[3][3])+SM[1][2]*SM[1][3]*SM[2][3]+1./27.;
776 
777  SP=TMath::Cos(TMath::ACos(TMath::Max(TMath::Min(SR/TMath::Sqrt(-pow(SQ,3)),1.),-1.))/3.);
778 
779  P[N+1][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Max(2.*SP,TMath::Sqrt(3.*(1.-SP*SP))-SP);
780  P[N+3][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Min(2.*SP,-TMath::Sqrt(3.*(1.-SP*SP))-SP);
781  P[N+2][4]=1.-P[N+1][4]-P[N+3][4];
782  if (P[N+2][4]> 1E-5) {
783 
784 // C...Find first and last eigenvector by solving equation system.
785  for (I=1;I<4;I=I+2) { // 240
786  for (J1=1;J1<4;J1++) { // 180
787  SV[J1][J1]=SM[J1][J1]-P[N+I][4];
788  for (J2=J1+1;J2<4;J2++) { // 170
789  SV[J1][J2]=SM[J1][J2];
790  SV[J2][J1]=SM[J1][J2];
791  }
792  } // 180
793  SMAX=0.;
794  for (J1=1;J1<4;J1++) { // 200
795  for (J2=1;J2<4;J2++) { // 190
796  if(std::fabs(SV[J1][J2])>SMAX) { // 190
797  JA=J1;
798  JB=J2;
799  SMAX=std::fabs(SV[J1][J2]);
800  }
801  } // 190
802  } // 200
803  SMAX=0.;
804  for (J3=JA+1;J3<JA+3;J3++) { // 220
805  J1=J3-3*((J3-1)/3);
806  RL=SV[J1][JB]/SV[JA][JB];
807  for (J2=1;J2<4;J2++) { // 210
808  SV[J1][J2]=SV[J1][J2]-RL*SV[JA][J2];
809  if (std::fabs(SV[J1][J2])>SMAX) { // GOTO 210
810  JC=J1;
811  SMAX=std::fabs(SV[J1][J2]);
812  }
813  } // 210
814  } // 220
815  JB1=JB+1-3*(JB/3);
816  JB2=JB+2-3*((JB+1)/3);
817  P[N+I][JB1]=-SV[JC][JB2];
818  P[N+I][JB2]=SV[JC][JB1];
819  P[N+I][JB]=-(SV[JA][JB1]*P[N+I][JB1]+SV[JA][JB2]*P[N+I][JB2])/
820  SV[JA][JB];
821  PA=TMath::Sqrt(pow(P[N+I][1],2)+pow(P[N+I][2],2)+pow(P[N+I][3],2));
822 // make a random number
823  float pa=P[N-1][I];
824  rlu=std::fabs(pa)-std::fabs(int(pa)*1.);
825  rlu1=std::fabs(pa*pa)-std::fabs(int(pa*pa)*1.);
826  SGN=pow((-1.),1.*int(rlu+0.5));
827  for (J=1;J<4;J++) { // 230
828  P[N+I][J]=SGN*P[N+I][J]/PA;
829  } // 230
830  } // 240
831 
832 // C...Middle axis orthogonal to other two. Fill other codes.
833  SGN=pow((-1.),1.*int(rlu1+0.5));
834  P[N+2][1]=SGN*(P[N+1][2]*P[N+3][3]-P[N+1][3]*P[N+3][2]);
835  P[N+2][2]=SGN*(P[N+1][3]*P[N+3][1]-P[N+1][1]*P[N+3][3]);
836  P[N+2][3]=SGN*(P[N+1][1]*P[N+3][2]-P[N+1][2]*P[N+3][1]);
837 
838 // C...Calculate sphericity and aplanarity. Select storing option.
839  SPH=1.5*(P[N+2][4]+P[N+3][4]);
840  APL=1.5*P[N+3][4];
841 
842  } // check 1
843 
844  // so assume now we have Sphericity axis, which one give the minimal Pts
845  double etstot[4];
846  double eltot[4];
847  double sum23=0;
848  double sum22=0;
849  double sum33=0;
850  double pina[4];
851  double ax[4], ay[4], az[4];
852  for (int ia=1;ia<4;ia++) {
853  etstot[ia]=0;
854  eltot[ia]=0;
855  pina[ia]=0;
856  ax[ia]=P[N+ia][1];
857  ay[ia]=P[N+ia][2];
858  az[ia]=P[N+ia][3];
859  ax[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
860  ay[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
861  az[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
862  }
863 
864 
865  for (int k =0 ; k<m_np ; k++) {
866  // double eta,phi;
867  // getetaphi(m_mom(k,1),m_mom(k,2),m_mom(k,3),eta,phi);
868  double W=1.0;
869  for (int ia=1;ia<4;ia++) {
870  double e=sqrt(m_mom(k,1)*m_mom(k,1) +
871  m_mom(k,2)*m_mom(k,2) +
872  m_mom(k,3)*m_mom(k,3));
873  double el=ax[ia]*m_mom(k,1) + ay[ia]*m_mom(k,2) + az[ia]*m_mom(k,3);
874  pina[ia]=el;
875  double ets=(e*e-el*el);
876  etstot[ia]+=ets*W;
877  eltot[ia]+=el*el;
878  }
879  double a2=pina[2];
880  double a3=pina[3];
881  // double h=0.4;
882  //a2=pina[2]*cos(h)+pina[3]*sin(h);
883  //a3=pina[3]*cos(h)-pina[2]*sin(h);
884  sum22+=a2*a2*W;
885  sum23+=a2*a3*W;
886  sum33+=a3*a3*W;
887  }
888 
889 
890 
891  double pi=3.1415927;
892  double phi=pi/2.0;
893  double phip=pi/2.0;
894  double a=sum23;
895  double c=-a;
896  double b=sum22-sum33;
897  double disc=b*b-4*a*c;
898  // cout << " disc " << disc << endl;
899  if (disc>=0) {
900  double x1=(sqrt(disc)-b)/2/a;
901  double x2=(-sqrt(disc)-b)/2/a;
902  phi=atan(x1);
903  phip=atan(x2);
904  if (phi<0) phi=2.*pi+phi;
905  if (phip<0) phip=2.*pi+phip;
906  }
907  double p21=sum22*cos(phi)*cos(phi)+sum33*sin(phi)*sin(phi)+2*sum23*cos(phi)*sin(phi);
908  double p31=sum22*sin(phi)*sin(phi)+sum33*cos(phi)*cos(phi)-2*sum23*cos(phi)*sin(phi);
909 
910  double p22=sum22*cos(phip)*cos(phip)+sum33*sin(phip)*sin(phip)+2*sum23*cos(phip)*sin(phip);
911  double p32=sum22*sin(phip)*sin(phip)+sum33*cos(phip)*cos(phip)-2*sum23*cos(phip)*sin(phip);
912 
913 
914  double d1=std::fabs(p31*p31 - p21*p21);
915  double d2=std::fabs(p32*p32 - p22*p22);
916  //cout << " eltot " << eltot[2] << " " << eltot[3] << endl;
917  //cout << " phi " << phi << " " << phip << endl;
918  //cout << " d " << d1 << " " << d2 << endl;
919  p2=p21;
920  p3=p31;
921  if (d2>d1) {
922  p2=p22;
923  p3=p32;
924  }
925  pnorm=sqrt(eltot[1]);
926  if (p2>p3) {
927  p3=sqrt(p3);
928  p2=sqrt(p2);
929  }else {
930  double p4=p3;
931  p3=sqrt(p2);
932  p2=sqrt(p4);
933  }
934  //cout << " sum2 sum3 " << sqrt(sum22) << " " << sqrt(sum33) << endl;
935  //cout << " p2 p3 " << p2 << " " << p3 << endl;
936  //double els=sqrt(eltot[2]*eltot[2]+eltot[3]*eltot[3]);
937  // cout << " ets els " << (ettot[1]) << " " << els << endl;
938  return;
939 } // end planes_sphe
siStripFEDMonitor_P5_cff.Min
tuple Min
Definition: siStripFEDMonitor_P5_cff.py:55
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
P
#define P
Definition: PFJetBenchmark.cc:33
funct::exp
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
eta
T eta() const
Definition: Basic3DVectorLD.h:168
debug_cff.d1
tuple d1
Definition: debug_cff.py:7
siStripFEDMonitor_P5_cff.Max
tuple Max
Definition: siStripFEDMonitor_P5_cff.py:56
PWT
double PWT[12]
Definition: herwig.h:108
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
p4
double p4[4]
Definition: TauolaWrapper.h:92
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Exhume::I
const std::complex< double > I
Definition: I.h:8
p2
double p2[4]
Definition: TauolaWrapper.h:90
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.cc:1408
SGN
#define SGN(A)
Definition: Simplify_begin.h:68
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
p3
double p3[4]
Definition: TauolaWrapper.h:91
phi
Definition: DDAxes.h:10
void TopologyWorker::planes_sphe_wei ( double &  pnorm,
double &  p2,
double &  p3 
)

Definition at line 942 of file TopologyWorker.cc.

References a, b, trackerHits::c, funct::cos(), debug_cff::d1, eta(), funct::exp(), getetaphi(), Exhume::I, gen::k, m_mom, m_np, siStripFEDMonitor_P5_cff::Max, siStripFEDMonitor_P5_cff::Min, MultiGaussianStateTransform::N, P, p3, p4, phi, pi, funct::pow(), PWT, SGN, funct::sin(), and mathSSE::sqrt().

942  {
943  float SPH=-1;
944  float APL=-1;
945 // C...Calculate matrix to be diagonalized.
946  float P[1000][6];
947  double SM[4][4],SV[4][4];
948  double PA,PWT,PS,SQ,SR,SP,SMAX,SGN;
949  int NP;
950  int J, J1, J2, I, N, JA, JB, J3, JC, JB1, JB2;
951  JA=JB=JC=0;
952  double RL;
953  float rlu,rlu1;
954  //
955  // --- get the input form GTRACK arrays
956  //
957  N=m_np;
958  NP=0;
959  for (I=1;I<N+1;I++){
960  NP++; // start at one
961  P[NP][1]=m_mom(I-1,1) ;
962  P[NP][2]=m_mom(I-1,2) ;
963  P[NP][3]=m_mom(I-1,3) ;
964  P[NP][4]=m_mom(I-1,4) ;
965  P[NP][5]=0;
966  }
967  //
968  //---
969  //
970  N=NP;
971 
972  for (J1=1;J1<4;J1++) {
973  for (J2=J1;J2<4;J2++) {
974  SM[J1][J2]=0.;
975  }
976  }
977  PS=0.;
978  for (I=1;I<N+1;I++) { // 140
979  PA=sqrt(pow(P[I][1],2)+pow(P[I][2],2)+pow(P[I][3],2));
980  // double eta,phi;
981  // getetaphi(P[I][1],P[I][2],P[I][3],eta,phi);
982  // PWT=exp(-std::fabs(eta));
983  PWT=1.;
984  for (J1=1;J1<4;J1++) { // 130
985  for (J2=J1;J2<4;J2++) { // 120
986  SM[J1][J2]=SM[J1][J2]+PWT*P[I][J1]*P[I][J2];
987  }
988  } // 130
989  PS=PS+PWT*PA*PA;
990  } //140
991 // C...Very low multiplicities (0 or 1) not considered.
992  if(NP<2) {
993  SPH=-1.;
994  APL=-1.;
995  return;
996  }
997  for (J1=1;J1<4;J1++) { // 160
998  for (J2=J1;J2<4;J2++) { // 150
999  SM[J1][J2]=SM[J1][J2]/PS;
1000  }
1001  } // 160
1002 // C...Find eigenvalues to matrix (third degree equation).
1003  SQ=(SM[1][1]*SM[2][2]+SM[1][1]*SM[3][3]+SM[2][2]*SM[3][3]
1004  -pow(SM[1][2],2)
1005  -pow(SM[1][3],2)-pow(SM[2][3],2))/3.-1./9.;
1006  SR=-0.5*(SQ+1./9.+SM[1][1]*pow(SM[2][3],2)+SM[2][2]*pow(SM[1][3],2)+SM[3][3]*
1007  pow(SM[1][2],2)-SM[1][1]*SM[2][2]*SM[3][3])+SM[1][2]*SM[1][3]*SM[2][3]+1./27.;
1008 
1009  SP=TMath::Cos(TMath::ACos(TMath::Max(TMath::Min(SR/TMath::Sqrt(-pow(SQ,3)),1.),-1.))/3.);
1010 
1011  P[N+1][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Max(2.*SP,TMath::Sqrt(3.*(1.-SP*SP))-SP);
1012  P[N+3][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Min(2.*SP,-TMath::Sqrt(3.*(1.-SP*SP))-SP);
1013  P[N+2][4]=1.-P[N+1][4]-P[N+3][4];
1014  if (P[N+2][4]> 1E-5) {
1015 
1016 // C...Find first and last eigenvector by solving equation system.
1017  for (I=1;I<4;I=I+2) { // 240
1018  for (J1=1;J1<4;J1++) { // 180
1019  SV[J1][J1]=SM[J1][J1]-P[N+I][4];
1020  for (J2=J1+1;J2<4;J2++) { // 170
1021  SV[J1][J2]=SM[J1][J2];
1022  SV[J2][J1]=SM[J1][J2];
1023  }
1024  } // 180
1025  SMAX=0.;
1026  for (J1=1;J1<4;J1++) { // 200
1027  for (J2=1;J2<4;J2++) { // 190
1028  if(std::fabs(SV[J1][J2])>SMAX) { // 190
1029  JA=J1;
1030  JB=J2;
1031  SMAX=std::fabs(SV[J1][J2]);
1032  }
1033  } // 190
1034  } // 200
1035  SMAX=0.;
1036  for (J3=JA+1;J3<JA+3;J3++) { // 220
1037  J1=J3-3*((J3-1)/3);
1038  RL=SV[J1][JB]/SV[JA][JB];
1039  for (J2=1;J2<4;J2++) { // 210
1040  SV[J1][J2]=SV[J1][J2]-RL*SV[JA][J2];
1041  if (std::fabs(SV[J1][J2])>SMAX) { // GOTO 210
1042  JC=J1;
1043  SMAX=std::fabs(SV[J1][J2]);
1044  }
1045  } // 210
1046  } // 220
1047  JB1=JB+1-3*(JB/3);
1048  JB2=JB+2-3*((JB+1)/3);
1049  P[N+I][JB1]=-SV[JC][JB2];
1050  P[N+I][JB2]=SV[JC][JB1];
1051  P[N+I][JB]=-(SV[JA][JB1]*P[N+I][JB1]+SV[JA][JB2]*P[N+I][JB2])/
1052  SV[JA][JB];
1053  PA=TMath::Sqrt(pow(P[N+I][1],2)+pow(P[N+I][2],2)+pow(P[N+I][3],2));
1054 // make a random number
1055  float pa=P[N-1][I];
1056  rlu=std::fabs(pa)-std::fabs(int(pa)*1.);
1057  rlu1=std::fabs(pa*pa)-std::fabs(int(pa*pa)*1.);
1058  SGN=pow((-1.),1.*int(rlu+0.5));
1059  for (J=1;J<4;J++) { // 230
1060  P[N+I][J]=SGN*P[N+I][J]/PA;
1061  } // 230
1062  } // 240
1063 
1064 // C...Middle axis orthogonal to other two. Fill other codes.
1065  SGN=pow((-1.),1.*int(rlu1+0.5));
1066  P[N+2][1]=SGN*(P[N+1][2]*P[N+3][3]-P[N+1][3]*P[N+3][2]);
1067  P[N+2][2]=SGN*(P[N+1][3]*P[N+3][1]-P[N+1][1]*P[N+3][3]);
1068  P[N+2][3]=SGN*(P[N+1][1]*P[N+3][2]-P[N+1][2]*P[N+3][1]);
1069 
1070 // C...Calculate sphericity and aplanarity. Select storing option.
1071  SPH=1.5*(P[N+2][4]+P[N+3][4]);
1072  APL=1.5*P[N+3][4];
1073 
1074  } // check 1
1075 
1076  // so assume now we have Sphericity axis, which one give the minimal Pts
1077  double etstot[4];
1078  double eltot[4];
1079  double sum23=0;
1080  double sum22=0;
1081  double sum33=0;
1082  double pina[4];
1083  double ax[4], ay[4], az[4];
1084  for (int ia=1;ia<4;ia++) {
1085  etstot[ia]=0;
1086  eltot[ia]=0;
1087  pina[ia]=0;
1088  ax[ia]=P[N+ia][1];
1089  ay[ia]=P[N+ia][2];
1090  az[ia]=P[N+ia][3];
1091  ax[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1092  ay[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1093  az[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1094  }
1095 
1096  for (int k =0 ; k<m_np ; k++) {
1097  double eta,phi;
1098  getetaphi(m_mom(k,1),m_mom(k,2),m_mom(k,3),eta,phi);
1099  double W=exp(-std::fabs(eta*1.0));
1100  for (int ia=1;ia<4;ia++) {
1101  double e=sqrt(m_mom(k,1)*m_mom(k,1) +
1102  m_mom(k,2)*m_mom(k,2) +
1103  m_mom(k,3)*m_mom(k,3));
1104  double el=ax[ia]*m_mom(k,1) + ay[ia]*m_mom(k,2) + az[ia]*m_mom(k,3);
1105  pina[ia]=el;
1106  double ets=(e*e-el*el);
1107  etstot[ia]+=ets*W;
1108  eltot[ia]+=el*el*W;
1109  }
1110  double a2=pina[2];
1111  double a3=pina[3];
1112  // double h=0.4;
1113  //a2=pina[2]*cos(h)+pina[3]*sin(h);
1114  //a3=pina[3]*cos(h)-pina[2]*sin(h);
1115  sum22+=a2*a2*W;
1116  sum23+=a2*a3*W;
1117  sum33+=a3*a3*W;
1118  }
1119 
1120 
1121  double pi=3.1415927;
1122  double phi=pi/2.0;
1123  double phip=pi/2.0;
1124  double a=sum23;
1125  double c=-a;
1126  double b=sum22-sum33;
1127  double disc=b*b-4*a*c;
1128  // cout << " disc " << disc << endl;
1129  if (disc>=0) {
1130  double x1=(sqrt(disc)-b)/2/a;
1131  double x2=(-sqrt(disc)-b)/2/a;
1132  phi=atan(x1);
1133  phip=atan(x2);
1134  if (phi<0) phi=2.*pi+phi;
1135  if (phip<0) phip=2.*pi+phip;
1136  }
1137  double p21=sum22*cos(phi)*cos(phi)+sum33*sin(phi)*sin(phi)+2*sum23*cos(phi)*sin(phi);
1138  double p31=sum22*sin(phi)*sin(phi)+sum33*cos(phi)*cos(phi)-2*sum23*cos(phi)*sin(phi);
1139 
1140  double p22=sum22*cos(phip)*cos(phip)+sum33*sin(phip)*sin(phip)+2*sum23*cos(phip)*sin(phip);
1141  double p32=sum22*sin(phip)*sin(phip)+sum33*cos(phip)*cos(phip)-2*sum23*cos(phip)*sin(phip);
1142 
1143 
1144  double d1=std::fabs(p31*p31 - p21*p21);
1145  double d2=std::fabs(p32*p32 - p22*p22);
1146  //cout << " eltot " << eltot[2] << " " << eltot[3] << endl;
1147  //cout << " phi " << phi << " " << phip << endl;
1148  //cout << " d " << d1 << " " << d2 << endl;
1149  p2=p21;
1150  p3=p31;
1151  if (d2>d1) {
1152  p2=p22;
1153  p3=p32;
1154  }
1155  pnorm=sqrt(eltot[1]);
1156  if (p2>p3) {
1157  p3=sqrt(p3);
1158  p2=sqrt(p2);
1159  }else {
1160  double p4=p3;
1161  p3=sqrt(p2);
1162  p2=sqrt(p4);
1163  }
1164  //cout << " sum2 sum3 " << sqrt(sum22) << " " << sqrt(sum33) << endl;
1165  //cout << " p2 p3 " << p2 << " " << p3 << endl;
1166  //double els=sqrt(eltot[2]*eltot[2]+eltot[3]*eltot[3]);
1167  // cout << " ets els " << (ettot[1]) << " " << els << endl;
1168  return;
1169 } // end planes_sphe
siStripFEDMonitor_P5_cff.Min
tuple Min
Definition: siStripFEDMonitor_P5_cff.py:55
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
P
#define P
Definition: PFJetBenchmark.cc:33
funct::exp
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
eta
T eta() const
Definition: Basic3DVectorLD.h:168
debug_cff.d1
tuple d1
Definition: debug_cff.py:7
siStripFEDMonitor_P5_cff.Max
tuple Max
Definition: siStripFEDMonitor_P5_cff.py:56
PWT
double PWT[12]
Definition: herwig.h:108
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
p4
double p4[4]
Definition: TauolaWrapper.h:92
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Exhume::I
const std::complex< double > I
Definition: I.h:8
p2
double p2[4]
Definition: TauolaWrapper.h:90
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.cc:1408
SGN
#define SGN(A)
Definition: Simplify_begin.h:68
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
p3
double p3[4]
Definition: TauolaWrapper.h:91
phi
Definition: DDAxes.h:10
void TopologyWorker::planes_thrust ( double &  pnorm,
double &  p2,
double &  p3 
)

Definition at line 1173 of file TopologyWorker.cc.

References a, b, trackerHits::c, funct::cos(), debug_cff::d1, gen::k, m_mom, m_np, majorAxis(), minorAxis(), p3, p4, phi, pi, funct::sin(), mathSSE::sqrt(), and thrustAxis().

1173  {
1174  TVector3 thrustaxis=thrustAxis();
1175  TVector3 majoraxis=majorAxis();
1176  TVector3 minoraxis=minorAxis();
1177  // so assume now we have Sphericity axis, which one give the minimal Pts
1178  double etstot[4];
1179  double eltot[4];
1180  double sum23=0;
1181  double sum22=0;
1182  double sum33=0;
1183  double pina[4];
1184  double ax[4], ay[4], az[4];
1185  ax[1]=thrustaxis(0); ay[1]=thrustaxis(1); az[1]=thrustaxis(2);
1186  ax[2]=minoraxis(0); ay[2]=minoraxis(1); az[2]=minoraxis(2);
1187  ax[3]=majoraxis(0); ay[3]=majoraxis(1); az[3]=majoraxis(2);
1188  for (int ia=1;ia<4;ia++) {
1189  etstot[ia]=0;
1190  eltot[ia]=0;
1191  pina[ia]=0;
1192  ax[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1193  ay[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1194  az[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1195  }
1196 
1197  for (int k =0 ; k<m_np ; k++) {
1198  for (int ia=1;ia<4;ia++) {
1199  double e=sqrt(m_mom(k,1)*m_mom(k,1) +
1200  m_mom(k,2)*m_mom(k,2) +
1201  m_mom(k,3)*m_mom(k,3));
1202  double el=ax[ia]*m_mom(k,1) + ay[ia]*m_mom(k,2) + az[ia]*m_mom(k,3);
1203  pina[ia]=el;
1204  double ets=(e*e-el*el);
1205  etstot[ia]+=ets;
1206  eltot[ia]+=std::fabs(el);
1207  }
1208  double a2=pina[2];
1209  double a3=pina[3];
1210  // double h=0.4;
1211  //a2=pina[2]*cos(h)+pina[3]*sin(h);
1212  //a3=pina[3]*cos(h)-pina[2]*sin(h);
1213  sum22+=a2*a2;
1214  sum23+=a2*a3;
1215  sum33+=a3*a3;
1216  }
1217 
1218 
1219  double pi=3.1415927;
1220  double phi=pi/2.0;
1221  double phip=pi/2.0;
1222  double a=sum23;
1223  double c=-a;
1224  double b=sum22-sum33;
1225  double disc=b*b-4*a*c;
1226  // cout << " disc " << disc << endl;
1227  if (disc>=0) {
1228  double x1=(sqrt(disc)-b)/2/a;
1229  double x2=(-sqrt(disc)-b)/2/a;
1230  phi=atan(x1);
1231  phip=atan(x2);
1232  if (phi<0) phi=2.*pi+phi;
1233  if (phip<0) phip=2.*pi+phip;
1234  }
1235  double p21=sum22*cos(phi)*cos(phi)+sum33*sin(phi)*sin(phi)+2*sum23*cos(phi)*sin(phi);
1236  double p31=sum22*sin(phi)*sin(phi)+sum33*cos(phi)*cos(phi)-2*sum23*cos(phi)*sin(phi);
1237 
1238  double p22=sum22*cos(phip)*cos(phip)+sum33*sin(phip)*sin(phip)+2*sum23*cos(phip)*sin(phip);
1239  double p32=sum22*sin(phip)*sin(phip)+sum33*cos(phip)*cos(phip)-2*sum23*cos(phip)*sin(phip);
1240 
1241 
1242  double d1=std::fabs(p31*p31 - p21*p21);
1243  double d2=std::fabs(p32*p32 - p22*p22);
1244  //cout << " eltot " << eltot[2] << " " << eltot[3] << endl;
1245  //cout << " phi " << phi << " " << phip << endl;
1246  //cout << " d " << d1 << " " << d2 << endl;
1247  p2=p21;
1248  p3=p31;
1249  if (d2>d1) {
1250  p2=p22;
1251  p3=p32;
1252  }
1253  pnorm=sqrt(etstot[1]);
1254  if (p2>p3) {
1255  p3=sqrt(p3);
1256  p2=sqrt(p2);
1257  }else {
1258  double p4=p3;
1259  p3=sqrt(p2);
1260  p2=sqrt(p4);
1261  }
1262  //cout << " sum2 sum3 " << sqrt(sum22) << " " << sqrt(sum33) << endl;
1263  //cout << " p2 p3 " << p2 << " " << p3 << endl;
1264  //double els=sqrt(eltot[2]*eltot[2]+eltot[3]*eltot[3]);
1265  // cout << " ets els " << (ettot[1]) << " " << els << endl;
1266  return;
1267 } // end planes_thru
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
TopologyWorker::thrustAxis
TVector3 thrustAxis()
Definition: TopologyWorker.cc:420
debug_cff.d1
tuple d1
Definition: debug_cff.py:7
TopologyWorker::minorAxis
TVector3 minorAxis()
Definition: TopologyWorker.cc:432
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
p4
double p4[4]
Definition: TauolaWrapper.h:92
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
p2
double p2[4]
Definition: TauolaWrapper.h:90
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
TopologyWorker::majorAxis
TVector3 majorAxis()
Definition: TopologyWorker.cc:426
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
p3
double p3[4]
Definition: TauolaWrapper.h:91
phi
Definition: DDAxes.h:10
void TopologyWorker::sanda ( )
private

Definition at line 556 of file TopologyWorker.cc.

References Exhume::I, m_apl, m_mom, m_np, m_sanda_called, m_sph, siStripFEDMonitor_P5_cff::Max, siStripFEDMonitor_P5_cff::Min, MultiGaussianStateTransform::N, P, funct::pow(), PWT, SGN, and mathSSE::sqrt().

Referenced by get_aplanarity(), and get_sphericity().

556  {
557  float SPH=-1;
558  float APL=-1;
559  m_sanda_called=true;
560  //=======================================================================
561  // By M.Vreeswijk, (core was fortran, stolen from somewhere)
562  // Purpose: to perform sphericity tensor analysis to give sphericity
563  // and aplanarity.
564  //
565  // Needs: Array (standard from root-tuples): GTRACK_px, py, pz
566  // The number of tracks in these arrays: GTRACK_ijet
567  // In addition: Array GTRACK_ijet contains a jet number 'ijet'
568  // (if you wish to change this, simply change code)
569  //
570  // Uses: TVector3 and TLorentzVector classes
571  //
572  // Output: Sphericity SPH and Aplanarity APL
573  //=======================================================================
574 // C...Calculate matrix to be diagonalized.
575  float P[1000][6];
576  double SM[4][4],SV[4][4];
577  double PA,PWT,PS,SQ,SR,SP,SMAX,SGN;
578  int NP;
579  int J, J1, J2, I, N, JA, JB, J3, JC, JB1, JB2;
580  JA=JB=JC=0;
581  double RL;
582  float rlu,rlu1;
583  //
584  // --- get the input form GTRACK arrays
585  //
586  N=m_np;
587  NP=0;
588  for (I=1;I<N+1;I++){
589  NP++; // start at one
590  P[NP][1]=m_mom(I-1,1) ;
591  P[NP][2]=m_mom(I-1,2) ;
592  P[NP][3]=m_mom(I-1,3) ;
593  P[NP][4]=m_mom(I-1,4) ;
594  P[NP][5]=0;
595  }
596  //
597  //---
598  //
599  N=NP;
600 
601  for (J1=1;J1<4;J1++) {
602  for (J2=J1;J2<4;J2++) {
603  SM[J1][J2]=0.;
604  }
605  }
606  PS=0.;
607  for (I=1;I<N+1;I++) { // 140
608  PA=sqrt(pow(P[I][1],2)+pow(P[I][2],2)+pow(P[I][3],2));
609  PWT=1.;
610  for (J1=1;J1<4;J1++) { // 130
611  for (J2=J1;J2<4;J2++) { // 120
612  SM[J1][J2]=SM[J1][J2]+PWT*P[I][J1]*P[I][J2];
613  }
614  } // 130
615  PS=PS+PWT*PA*PA;
616  } //140
617 // C...Very low multiplicities (0 or 1) not considered.
618  if(NP<2) {
619  SPH=-1.;
620  APL=-1.;
621  return;
622  }
623  for (J1=1;J1<4;J1++) { // 160
624  for (J2=J1;J2<4;J2++) { // 150
625  SM[J1][J2]=SM[J1][J2]/PS;
626  }
627  } // 160
628 // C...Find eigenvalues to matrix (third degree equation).
629  SQ=(SM[1][1]*SM[2][2]+SM[1][1]*SM[3][3]+SM[2][2]*SM[3][3]
630  -pow(SM[1][2],2)
631  -pow(SM[1][3],2)-pow(SM[2][3],2))/3.-1./9.;
632  SR=-0.5*(SQ+1./9.+SM[1][1]*pow(SM[2][3],2)+SM[2][2]*pow(SM[1][3],2)+SM[3][3]*
633  pow(SM[1][2],2)-SM[1][1]*SM[2][2]*SM[3][3])+SM[1][2]*SM[1][3]*SM[2][3]+1./27.;
634 
635  SP=TMath::Cos(TMath::ACos(TMath::Max(TMath::Min(SR/TMath::Sqrt(-pow(SQ,3)),1.),-1.))/3.);
636 
637  P[N+1][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Max(2.*SP,TMath::Sqrt(3.*(1.-SP*SP))-SP);
638  P[N+3][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Min(2.*SP,-TMath::Sqrt(3.*(1.-SP*SP))-SP);
639  P[N+2][4]=1.-P[N+1][4]-P[N+3][4];
640  if (P[N+2][4]> 1E-5) {
641 
642 // C...Find first and last eigenvector by solving equation system.
643  for (I=1;I<4;I=I+2) { // 240
644  for (J1=1;J1<4;J1++) { // 180
645  SV[J1][J1]=SM[J1][J1]-P[N+I][4];
646  for (J2=J1+1;J2<4;J2++) { // 170
647  SV[J1][J2]=SM[J1][J2];
648  SV[J2][J1]=SM[J1][J2];
649  }
650  } // 180
651  SMAX=0.;
652  for (J1=1;J1<4;J1++) { // 200
653  for (J2=1;J2<4;J2++) { // 190
654  if(std::fabs(SV[J1][J2])>SMAX) { // 190
655  JA=J1;
656  JB=J2;
657  SMAX=std::fabs(SV[J1][J2]);
658  }
659  } // 190
660  } // 200
661  SMAX=0.;
662  for (J3=JA+1;J3<JA+3;J3++) { // 220
663  J1=J3-3*((J3-1)/3);
664  RL=SV[J1][JB]/SV[JA][JB];
665  for (J2=1;J2<4;J2++) { // 210
666  SV[J1][J2]=SV[J1][J2]-RL*SV[JA][J2];
667  if (std::fabs(SV[J1][J2])>SMAX) { // GOTO 210
668  JC=J1;
669  SMAX=std::fabs(SV[J1][J2]);
670  }
671  } // 210
672  } // 220
673  JB1=JB+1-3*(JB/3);
674  JB2=JB+2-3*((JB+1)/3);
675  P[N+I][JB1]=-SV[JC][JB2];
676  P[N+I][JB2]=SV[JC][JB1];
677  P[N+I][JB]=-(SV[JA][JB1]*P[N+I][JB1]+SV[JA][JB2]*P[N+I][JB2])/
678  SV[JA][JB];
679  PA=TMath::Sqrt(pow(P[N+I][1],2)+pow(P[N+I][2],2)+pow(P[N+I][3],2));
680 // make a random number
681  float pa=P[N-1][I];
682  rlu=std::fabs(pa)-std::fabs(int(pa)*1.);
683  rlu1=std::fabs(pa*pa)-std::fabs(int(pa*pa)*1.);
684  SGN=pow((-1.),1.*int(rlu+0.5));
685  for (J=1;J<4;J++) { // 230
686  P[N+I][J]=SGN*P[N+I][J]/PA;
687  } // 230
688  } // 240
689 
690 // C...Middle axis orthogonal to other two. Fill other codes.
691  SGN=pow((-1.),1.*int(rlu1+0.5));
692  P[N+2][1]=SGN*(P[N+1][2]*P[N+3][3]-P[N+1][3]*P[N+3][2]);
693  P[N+2][2]=SGN*(P[N+1][3]*P[N+3][1]-P[N+1][1]*P[N+3][3]);
694  P[N+2][3]=SGN*(P[N+1][1]*P[N+3][2]-P[N+1][2]*P[N+3][1]);
695 
696 // C...Calculate sphericity and aplanarity. Select storing option.
697  SPH=1.5*(P[N+2][4]+P[N+3][4]);
698  APL=1.5*P[N+3][4];
699 
700  } // check 1
701 
702  m_sph=SPH;
703  m_apl=APL;
704  return;
705 } // end sanda
siStripFEDMonitor_P5_cff.Min
tuple Min
Definition: siStripFEDMonitor_P5_cff.py:55
P
#define P
Definition: PFJetBenchmark.cc:33
siStripFEDMonitor_P5_cff.Max
tuple Max
Definition: siStripFEDMonitor_P5_cff.py:56
PWT
double PWT[12]
Definition: herwig.h:108
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
Exhume::I
const std::complex< double > I
Definition: I.h:8
SGN
#define SGN(A)
Definition: Simplify_begin.h:68
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void TopologyWorker::setFast ( int  nf)

Definition at line 404 of file TopologyWorker.cc.

References m_iFast.

405 {
406  // Error if sp not positive.
407  if ( nf > 3 ) m_iFast = nf;
408  return;
409 }
void TopologyWorker::setPartList ( TObjArray *  e1,
TObjArray *  e2 
)

Definition at line 42 of file TopologyWorker.cc.

References CalcHTstuff(), CalcSqrts(), CalcWmul(), benchmark_cfg::cerr, gather_cfg::cout, i, iPow(), j, gen::k, ludbrb(), m_boost, m_dAxes, m_dConv, m_dDeltaThPower, m_dOblateness, m_dThrust, m_fowo_called, m_iFast, m_iGood, m_maxpart, m_mom, m_mom2, m_np, m_np2, m_random, m_sanda_called, m_verbose, siStripFEDMonitor_P5_cff::Min, n, runTheMatrix::np, np2, connectstrParser::o, L1TEmulatorMonitor_cff::p, phi, FWPFMaths::sgn(), sign(), cond::rpcobtemp::temp, tmax, ulAngle(), v, X, and Gflash::Z.

43 {
44  //To make this look like normal physics notation the
45  //zeroth element of each array, mom[i][0], will be ignored
46  //and operations will be on elements 1,2,3...
47  TMatrix mom(m_maxpart,6);
48  TMatrix mom2(m_maxpart,6);
49  double tmax = 0;
50  double phi = 0.;
51  double the = 0.;
52  double sgn;
53  TMatrix fast(m_iFast + 1,6);
54  TMatrix work(11,6);
55  double tdi[4] = {0.,0.,0.,0.};
56  double tds;
57  double tpr[4] = {0.,0.,0.,0.};
58  double thp;
59  double thps;
60  double pxtot=0;
61  double pytot=0;
62  double pztot=0;
63  double etot=0;
64 
65  TMatrix temp(3,5);
66  Int_t np = 0;
67  Int_t numElements = e1->GetEntries();
68  Int_t numElements2 = e2->GetEntries();
69 
70  // trying to sort...
71 
72 
73 
74  m_np=0;
75  for(Int_t elem=0;elem<numElements;elem++) {
76  if(m_verbose){
77  std::cout << "looping over array, element " << elem << std::endl;
78  }
79  TObject* o = (TObject*) e1->At(elem);
80  if(m_verbose){
81  std::cerr << "TopologyWorker:SetPartList(): adding jet " << elem << "." << std::endl;
82  }
83  if (np >= m_maxpart) {
84  printf("Too many particles input to TopologyWorker");
85  return;
86  }
87  if(m_verbose){
88  std::cout << "getting name of object..." << std::endl;
89  }
90  TString nam(o->IsA()->GetName());
91  if(m_verbose){
92  std::cout << "TopologyWorker::setPartList(): object is of type " << nam << std::endl;
93  }
94  if (nam.Contains("TVector3")) {
95  TVector3 v(((TVector3 *) o)->X(),
96  ((TVector3 *) o)->Y(),
97  ((TVector3 *) o)->Z());
98  mom(np,1) = v.X();
99  mom(np,2) = v.Y();
100  mom(np,3) = v.Z();
101  mom(np,4) = v.Mag();
102  }
103  else if (nam.Contains("TLorentzVector")) {
104  TVector3 v(((TLorentzVector *) o)->X(),
105  ((TLorentzVector *) o)->Y(),
106  ((TLorentzVector *) o)->Z());
107  mom(np,1) = v.X();
108  mom(np,2) = v.Y();
109  mom(np,3) = v.Z();
110  mom(np,4) = ((TLorentzVector *) o)->T();
111  }
112  else {
113  printf("TopologyWorker::setEvent input is not a TVector3 or a TLorentzVector\n");
114  continue;
115  }
116 
117 
118  if ( TMath::Abs( m_dDeltaThPower ) <= 0.001 ) {
119  mom(np,5) = 1.0;
120  }
121  else {
122  mom(np,5) = TMath::Power(mom(np,4),m_dDeltaThPower);
123  }
124  tmax = tmax + mom(np,4)*mom(np,5);
125  pxtot+=mom(np,1);
126  pytot+=mom(np,2);
127  pztot+=mom(np,3);
128  etot+=mom(np,4);
129  np++;
130  m_np=np;
131  }
132  Int_t np2=0;
133  // second jet array.... only use some values here.
134  for(Int_t elem=0;elem<numElements2;elem++) {
135  //cout << elem << endl;
136  TObject* o = (TObject*) e2->At(elem);
137  if (np2 >= m_maxpart) {
138  printf("Too many particles input to TopologyWorker");
139  return;
140  }
141 
142  TString nam(o->IsA()->GetName());
143  if (nam.Contains("TVector3")) {
144  TVector3 v(((TVector3 *) o)->X(),
145  ((TVector3 *) o)->Y(),
146  ((TVector3 *) o)->Z());
147  mom2(np2,1) = v.X();
148  mom2(np2,2) = v.Y();
149  mom2(np2,3) = v.Z();
150  mom2(np2,4) = v.Mag();
151  }
152  else if (nam.Contains("TLorentzVector")) {
153  TVector3 v(((TLorentzVector *) o)->X(),
154  ((TLorentzVector *) o)->Y(),
155  ((TLorentzVector *) o)->Z());
156  mom2(np2,1) = v.X();
157  mom2(np2,2) = v.Y();
158  mom2(np2,3) = v.Z();
159  mom2(np2,4) = ((TLorentzVector *) o)->T();
160  // cout << "mom2: " << mom2(np2,1) << ", " << mom2(np2,2)<<", " << mom2(np2,3)<<", " << mom2(np2,4)<< endl;
161  }
162  else {
163  printf("TopologyWorker::setEvent Second vector input is not a TVector3 or a TLorentzVector\n");
164  continue;
165  }
166  np2++;
167  m_np2=np2;
168  }
169  m_mom2=mom2;
170 
171  if (m_boost && m_np>1) {
172  printf("TopologyWorker::setEvent Only boosting first vector so watch out when you do this!!!");
173  TVector3 booze(-pxtot/etot,-pytot/etot,-pztot/etot);
174  TLorentzVector l1;
175  for (int k=0; k<m_np ; k++) {
176  l1.SetPx(mom(k,1));
177  l1.SetPy(mom(k,2));
178  l1.SetPz(mom(k,3));
179  l1.SetE(mom(k,4));
180  l1.Boost(booze);
181  mom(k,1)=l1.Px();
182  mom(k,2)=l1.Py();
183  mom(k,3)=l1.Pz();
184  mom(k,4)=l1.E();
185  }
186  }
187 
188  m_sanda_called=false;
189  m_fowo_called=false;
190  for (int ip=0;ip<m_maxpart;ip++) {
191  for (int id=0;id<6;id++) {
192  m_mom(ip,id)=mom(ip,id);
193  }
194  }
195 
196  if ( np < 2 ) {
197  m_dThrust[1] = -1.0;
198  m_dOblateness = -1.0;
199  return;
200  }
201  // for pass = 1: find thrust axis.
202  // for pass = 2: find major axis.
203  for ( Int_t pass=1; pass < 3; pass++) {
204  if ( pass == 2 ) {
205  phi = ulAngle(m_dAxes(1,1), m_dAxes(1,2));
206  ludbrb( &mom, 0, -phi, 0., 0., 0. );
207  for ( Int_t i = 0; i < 3; i++ ) {
208  for ( Int_t j = 1; j < 4; j++ ) {
209  temp(i,j) = m_dAxes(i+1,j);
210  }
211  temp(i,4) = 0;
212  }
213  ludbrb(&temp,0.,-phi,0.,0.,0.);
214  for ( Int_t ib = 0; ib < 3; ib++ ) {
215  for ( Int_t j = 1; j < 4; j++ ) {
216  m_dAxes(ib+1,j) = temp(ib,j);
217  }
218  }
219  the = ulAngle( m_dAxes(1,3), m_dAxes(1,1) );
220  ludbrb( &mom, -the, 0., 0., 0., 0. );
221  for ( Int_t ic = 0; ic < 3; ic++ ) {
222  for ( Int_t j = 1; j < 4; j++ ) {
223  temp(ic,j) = m_dAxes(ic+1,j);
224  }
225  temp(ic,4) = 0;
226  }
227  ludbrb(&temp,-the,0.,0.,0.,0.);
228  for ( Int_t id = 0; id < 3; id++ ) {
229  for ( Int_t j = 1; j < 4; j++ ) {
230  m_dAxes(id+1,j) = temp(id,j);
231  }
232  }
233  }
234  for ( Int_t ifas = 0; ifas < m_iFast + 1 ; ifas++ ) {
235  fast(ifas,4) = 0.;
236  }
237  // Find the m_iFast highest momentum particles and
238  // put the highest in fast[0], next in fast[1],....fast[m_iFast-1].
239  // fast[m_iFast] is just a workspace.
240  for ( Int_t i = 0; i < np; i++ ) {
241  if ( pass == 2 ) {
242  mom(i,4) = TMath::Sqrt( mom(i,1)*mom(i,1)
243  + mom(i,2)*mom(i,2) );
244  }
245  for ( Int_t ifas = m_iFast - 1; ifas > -1; ifas-- ) {
246  if ( mom(i,4) > fast(ifas,4) ) {
247  for ( Int_t j = 1; j < 6; j++ ) {
248  fast(ifas+1,j) = fast(ifas,j);
249  if ( ifas == 0 ) fast(ifas,j) = mom(i,j);
250  }
251  }
252  else {
253  for ( Int_t j = 1; j < 6; j++ ) {
254  fast(ifas+1,j) = mom(i,j);
255  }
256  break;
257  }
258  }
259  }
260  // Find axis with highest thrust (case 1)/ highest major (case 2).
261  for ( Int_t ie = 0; ie < work.GetNrows(); ie++ ) {
262  work(ie,4) = 0.;
263  }
264  Int_t p = TMath::Min( m_iFast, np ) - 1;
265  // Don't trust Math.pow to give right answer always.
266  // Want nc = 2**p.
267  Int_t nc = iPow(2,p);
268  for ( Int_t n = 0; n < nc; n++ ) {
269  for ( Int_t j = 1; j < 4; j++ ) {
270  tdi[j] = 0.;
271  }
272  for ( Int_t i = 0; i < TMath::Min(m_iFast,n); i++ ) {
273  sgn = fast(i,5);
274  if (iPow(2,(i+1))*((n+iPow(2,i))/iPow(2,(i+1))) >= i+1){
275  sgn = -sgn;
276  }
277  for ( Int_t j = 1; j < 5-pass; j++ ) {
278  tdi[j] = tdi[j] + sgn*fast(i,j);
279  }
280  }
281  tds = tdi[1]*tdi[1] + tdi[2]*tdi[2] + tdi[3]*tdi[3];
282  for ( Int_t iw = TMath::Min(n,9); iw > -1; iw--) {
283  if( tds > work(iw,4) ) {
284  for ( Int_t j = 1; j < 5; j++ ) {
285  work(iw+1,j) = work(iw,j);
286  if ( iw == 0 ) {
287  if ( j < 4 ) {
288  work(iw,j) = tdi[j];
289  }
290  else {
291  work(iw,j) = tds;
292  }
293  }
294  }
295  }
296  else {
297  for ( Int_t j = 1; j < 4; j++ ) {
298  work(iw+1,j) = tdi[j];
299  }
300  work(iw+1,4) = tds;
301  }
302  }
303  }
304  // Iterate direction of axis until stable maximum.
305  m_dThrust[pass] = 0;
306  thp = -99999.;
307  Int_t nagree = 0;
308  for ( Int_t iw = 0;
309  iw < TMath::Min(nc,10) && nagree < m_iGood; iw++ ){
310  thp = 0.;
311  thps = -99999.;
312  while ( thp > thps + m_dConv ) {
313  thps = thp;
314  for ( Int_t j = 1; j < 4; j++ ) {
315  if ( thp <= 1E-10 ) {
316  tdi[j] = work(iw,j);
317  }
318  else {
319  tdi[j] = tpr[j];
320  tpr[j] = 0;
321  }
322  }
323  for ( Int_t i = 0; i < np; i++ ) {
324  sgn = sign(mom(i,5),
325  tdi[1]*mom(i,1) +
326  tdi[2]*mom(i,2) +
327  tdi[3]*mom(i,3));
328  for ( Int_t j = 1; j < 5 - pass; j++ ){
329  tpr[j] = tpr[j] + sgn*mom(i,j);
330  }
331  }
332  thp = TMath::Sqrt(tpr[1]*tpr[1]
333  + tpr[2]*tpr[2]
334  + tpr[3]*tpr[3])/tmax;
335  }
336  // Save good axis. Try new initial axis until enough
337  // tries agree.
338  if ( thp < m_dThrust[pass] - m_dConv ) {
339  break;
340  }
341  if ( thp > m_dThrust[pass] + m_dConv ) {
342  nagree = 0;
343  sgn = iPow( -1, (Int_t)TMath::Nint(m_random.Rndm()) );
344  for ( Int_t j = 1; j < 4; j++ ) {
345  m_dAxes(pass,j) = sgn*tpr[j]/(tmax*thp);
346  }
347  m_dThrust[pass] = thp;
348  }
349  nagree = nagree + 1;
350  }
351  }
352  // Find minor axis and value by orthogonality.
353  sgn = iPow( -1, (Int_t)TMath::Nint(m_random.Rndm()));
354  m_dAxes(3,1) = -sgn*m_dAxes(2,2);
355  m_dAxes(3,2) = sgn*m_dAxes(2,1);
356  m_dAxes(3,3) = 0.;
357  thp = 0.;
358  for ( Int_t i = 0; i < np; i++ ) {
359  thp += mom(i,5)*TMath::Abs(m_dAxes(3,1)*mom(i,1) +
360  m_dAxes(3,2)*mom(i,2));
361  }
362  m_dThrust[3] = thp/tmax;
363  // Rotate back to original coordinate system.
364  for ( Int_t i6 = 0; i6 < 3; i6++ ) {
365  for ( Int_t j = 1; j < 4; j++ ) {
366  temp(i6,j) = m_dAxes(i6+1,j);
367  }
368  temp(i6,4) = 0;
369  }
370  ludbrb(&temp,the,phi,0.,0.,0.);
371  for ( Int_t i7 = 0; i7 < 3; i7++ ) {
372  for ( Int_t j = 1; j < 4; j++ ) {
373  m_dAxes(i7+1,j) = temp(i7,j);
374  }
375  }
376  m_dOblateness = m_dThrust[2] - m_dThrust[3];
377 
378  // more stuff:
379 
380  // calculate weighed jet multiplicity();
381  CalcWmul();
382  CalcHTstuff();
383  CalcSqrts();
384 
385 }
Gflash::Z
const double Z[kNumberCalorimeter]
Definition: GflashNameSpace.h:57
TopologyWorker::ulAngle
double ulAngle(double x, double y)
Definition: TopologyWorker.cc:450
i
int i
Definition: DBlmapReader.cc:9
siStripFEDMonitor_P5_cff.Min
tuple Min
Definition: siStripFEDMonitor_P5_cff.py:55
FWPFMaths::sgn
float sgn(float val)
Definition: FWPFMaths.cc:9
X
#define X(str)
Definition: MuonsGrabber.cc:49
TopologyWorker::m_dThrust
double m_dThrust[4]
Definition: TopologyWorker.h:120
TopologyWorker::ludbrb
void ludbrb(TMatrix *mom, double the, double phi, double bx, double by, double bz)
Definition: TopologyWorker.cc:483
j
int j
Definition: DBlmapReader.cc:9
TopologyWorker::sign
double sign(double a, double b)
Definition: TopologyWorker.cc:473
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
TopologyWorker::m_maxpart
static int m_maxpart
Definition: TopologyWorker.h:146
tmax
static const double tmax[3]
Definition: CastorTimeSlew.cc:7
np2
int np2
Definition: TauolaWrapper.h:100
benchmark_cfg.cerr
tuple cerr
Definition: benchmark_cfg.py:96
L1TEmulatorMonitor_cff.p
tuple p
Definition: L1TEmulatorMonitor_cff.py:48
TopologyWorker::m_dDeltaThPower
double m_dDeltaThPower
Definition: TopologyWorker.h:91
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
v
mathSSE::Vec4< T > v
Definition: newBasic3DVector.h:343
TopologyWorker::iPow
int iPow(int man, int exp)
Definition: TopologyWorker.cc:1455
phi
Definition: DDAxes.h:10
void TopologyWorker::setThMomPower ( double  tp)

Definition at line 390 of file TopologyWorker.cc.

References m_dDeltaThPower.

391 {
392  // Error if sp not positive.
393  if ( tp > 0. ) m_dDeltaThPower = tp - 1.0;
394  return;
395 }
TopologyWorker::m_dDeltaThPower
double m_dDeltaThPower
Definition: TopologyWorker.h:91
void TopologyWorker::setVerbose ( bool  loud)
inline

Definition at line 36 of file TopologyWorker.h.

References m_verbose.

36 {m_verbose=loud; return;}
double TopologyWorker::sign ( double  a,
double  b 
)
private

Definition at line 473 of file TopologyWorker.cc.

Referenced by setPartList(), and ulAngle().

473  {
474  if ( b < 0 ) {
475  return -TMath::Abs(a);
476  }
477  else {
478  return TMath::Abs(a);
479  }
480 }
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
void TopologyWorker::sumangles ( float &  sdeta,
float &  sdr 
)

Definition at line 1433 of file TopologyWorker.cc.

References getetaphi(), gen::k, m_mom, m_np, m_sumangles_called, and mathSSE::sqrt().

1433  {
1434  double eta1,eta2,phi1,phi2,deta,dphi,dr;
1435  m_sumangles_called=true;
1436  sdeta=0;
1437  sdr=0;
1438  for (int k=0;k<m_np;k++){
1439  for (int kp=k;kp<m_np;kp++){
1440  getetaphi(m_mom(k,1) , m_mom(k,2), m_mom(k,3), eta1,phi1);
1441  getetaphi(m_mom(kp,1) , m_mom(kp,2), m_mom(kp,3), eta2,phi2);
1442  dphi=std::fabs(phi1-phi2);
1443  if (dphi>3.1415) dphi=2*3.1415-dphi;
1444  deta=std::fabs(eta1-eta2);
1445  dr=sqrt(dphi*dphi+deta*deta);
1446  sdeta+=deta;
1447  sdr+=dr;
1448  }
1449  }
1450  return;
1451 }
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
gen::k
int k[5][pyjets_maxn]
Definition: Pythia6Hadronizer.cc:86
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.cc:1408
TVector3 TopologyWorker::thrust ( )

Definition at line 438 of file TopologyWorker.cc.

References m_dThrust, and m_Thrust.

438  {
439  TVector3 m_Thrust(m_dThrust[1],m_dThrust[2],m_dThrust[3]);
440  return m_Thrust;
441 }
TopologyWorker::m_dThrust
double m_dThrust[4]
Definition: TopologyWorker.h:120
TVector3 TopologyWorker::thrustAxis ( )

Definition at line 420 of file TopologyWorker.cc.

References m_dAxes, and m_ThrustAxis.

Referenced by planes_thrust().

420  {
421  TVector3 m_ThrustAxis(m_dAxes(1,1),m_dAxes(1,2),m_dAxes(1,3));
422  return m_ThrustAxis;
423 }
TopologyWorker::m_ThrustAxis
TVector3 m_ThrustAxis
Definition: TopologyWorker.h:110
double TopologyWorker::ulAngle ( double  x,
double  y 
)
private

Definition at line 450 of file TopologyWorker.cc.

References Pi, csvReporter::r, and sign().

Referenced by setPartList().

451 {
452  double ulangl = 0;
453  double r = TMath::Sqrt(x*x + y*y);
454  if ( r < 1.0E-20 ) {
455  return ulangl;
456  }
457  if ( TMath::Abs(x)/r < 0.8 ) {
458  ulangl = sign(TMath::ACos(x/r),y);
459  }
460  else {
461  ulangl = TMath::ASin(y/r);
462  if ( x < 0. && ulangl >= 0. ) {
463  ulangl = TMath::Pi() - ulangl;
464  }
465  else if ( x < 0. ) {
466  ulangl = -TMath::Pi() - ulangl;
467  }
468  }
469  return ulangl;
470 }
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
TopologyWorker::sign
double sign(double a, double b)
Definition: TopologyWorker.cc:473
x
Definition: DDAxes.h:10

Member Data Documentation

double TopologyWorker::m_apl
private

Definition at line 129 of file TopologyWorker.h.

Referenced by get_aplanarity(), sanda(), and TopologyWorker().

bool TopologyWorker::m_boost
private

Definition at line 126 of file TopologyWorker.h.

Referenced by setPartList(), and TopologyWorker().

double TopologyWorker::m_centrality
private

Definition at line 142 of file TopologyWorker.h.

Referenced by CalcHTstuff(), and get_centrality().

TMatrix TopologyWorker::m_dAxes
private

Definition at line 104 of file TopologyWorker.h.

Referenced by majorAxis(), minorAxis(), setPartList(), thrustAxis(), and TopologyWorker().

double TopologyWorker::m_dConv
private

Definition at line 97 of file TopologyWorker.h.

Referenced by setPartList().

double TopologyWorker::m_dDeltaThPower
private

Definition at line 91 of file TopologyWorker.h.

Referenced by getThMomPower(), setPartList(), and setThMomPower().

double TopologyWorker::m_dOblateness
private

Definition at line 121 of file TopologyWorker.h.

Referenced by oblateness(), and setPartList().

double TopologyWorker::m_dSphMomPower
private

Definition at line 88 of file TopologyWorker.h.

double TopologyWorker::m_dThrust[4]
private

Definition at line 120 of file TopologyWorker.h.

Referenced by setPartList(), and thrust().

double TopologyWorker::m_et0
private

Definition at line 138 of file TopologyWorker.h.

Referenced by CalcHTstuff(), get_et0(), and TopologyWorker().

double TopologyWorker::m_et56
private

Definition at line 141 of file TopologyWorker.h.

Referenced by CalcHTstuff(), get_et56(), and TopologyWorker().

bool TopologyWorker::m_fowo_called
private

Definition at line 125 of file TopologyWorker.h.

Referenced by fowo(), get_h10(), get_h20(), get_h30(), get_h40(), get_h50(), get_h60(), setPartList(), and TopologyWorker().

double TopologyWorker::m_h10
private

Definition at line 130 of file TopologyWorker.h.

Referenced by fowo(), get_h10(), and TopologyWorker().

double TopologyWorker::m_h20
private

Definition at line 131 of file TopologyWorker.h.

Referenced by fowo(), get_h20(), and TopologyWorker().

double TopologyWorker::m_h30
private

Definition at line 132 of file TopologyWorker.h.

Referenced by fowo(), get_h30(), and TopologyWorker().

double TopologyWorker::m_h40
private

Definition at line 133 of file TopologyWorker.h.

Referenced by fowo(), get_h40(), and TopologyWorker().

double TopologyWorker::m_h50
private

Definition at line 134 of file TopologyWorker.h.

Referenced by fowo(), and get_h50().

double TopologyWorker::m_h60
private

Definition at line 135 of file TopologyWorker.h.

Referenced by fowo(), and get_h60().

double TopologyWorker::m_ht
private

Definition at line 136 of file TopologyWorker.h.

Referenced by CalcHTstuff(), get_ht(), and TopologyWorker().

double TopologyWorker::m_ht3
private

Definition at line 137 of file TopologyWorker.h.

Referenced by CalcHTstuff(), get_ht3(), and TopologyWorker().

int TopologyWorker::m_iFast
private

Definition at line 94 of file TopologyWorker.h.

Referenced by getFast(), setFast(), and setPartList().

int TopologyWorker::m_iGood
private

Definition at line 100 of file TopologyWorker.h.

Referenced by setPartList().

TVector3 TopologyWorker::m_MajorAxis
private

Definition at line 111 of file TopologyWorker.h.

Referenced by majorAxis().

Int_t TopologyWorker::m_maxpart = 1000
staticprivate

Definition at line 146 of file TopologyWorker.h.

Referenced by setPartList(), and TopologyWorker().

TVector3 TopologyWorker::m_MinorAxis
private

Definition at line 112 of file TopologyWorker.h.

Referenced by minorAxis().

TMatrix TopologyWorker::m_mom
private

Definition at line 117 of file TopologyWorker.h.

Referenced by CalcEta(), CalcHTstuff(), CalcPt(), CalcSqrts(), fowo(), planes_sphe(), planes_sphe_wei(), planes_thrust(), sanda(), setPartList(), sumangles(), and TopologyWorker().

TMatrix TopologyWorker::m_mom2
private

Definition at line 118 of file TopologyWorker.h.

Referenced by CalcEta2(), CalcPt2(), setPartList(), and TopologyWorker().

double TopologyWorker::m_njetsweighed
private

Definition at line 140 of file TopologyWorker.h.

Referenced by CalcWmul(), get_njetW(), and TopologyWorker().

int TopologyWorker::m_np
private

Definition at line 122 of file TopologyWorker.h.

Referenced by CalcHTstuff(), CalcSqrts(), CalcWmul(), clear(), fowo(), planes_sphe(), planes_sphe_wei(), planes_thrust(), sanda(), setPartList(), sumangles(), and TopologyWorker().

int TopologyWorker::m_np2
private

Definition at line 123 of file TopologyWorker.h.

Referenced by CalcHTstuff(), clear(), setPartList(), and TopologyWorker().

TRandom TopologyWorker::m_random
private

Definition at line 115 of file TopologyWorker.h.

Referenced by setPartList().

bool TopologyWorker::m_sanda_called
private

Definition at line 124 of file TopologyWorker.h.

Referenced by get_aplanarity(), get_sphericity(), sanda(), setPartList(), and TopologyWorker().

double TopologyWorker::m_sph
private

Definition at line 128 of file TopologyWorker.h.

Referenced by get_sphericity(), sanda(), and TopologyWorker().

double TopologyWorker::m_sqrts
private

Definition at line 139 of file TopologyWorker.h.

Referenced by CalcSqrts(), get_sqrts(), and TopologyWorker().

bool TopologyWorker::m_sumangles_called
private

Definition at line 127 of file TopologyWorker.h.

Referenced by sumangles(), and TopologyWorker().

TVector3 TopologyWorker::m_Thrust
private

Definition at line 113 of file TopologyWorker.h.

Referenced by thrust().

TVector3 TopologyWorker::m_ThrustAxis
private

Definition at line 110 of file TopologyWorker.h.

Referenced by thrustAxis().

bool TopologyWorker::m_verbose
private

Definition at line 75 of file TopologyWorker.h.

Referenced by setPartList(), setVerbose(), and TopologyWorker().

Generated for CMSSW Reference Manual by doxygen 1.8.5