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 <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 32 of file TopologyWorker.h.

Constructor & Destructor Documentation

TopologyWorker::TopologyWorker ( )
inline

Definition at line 35 of file TopologyWorker.h.

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

Definition at line 176 of file TopologyWorker.h.

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.

176  :
177  m_dSphMomPower(2.0),m_dDeltaThPower(0),
178  m_iFast(4),m_dConv(0.0001),m_iGood(2)
179 {
180  m_dAxes.ResizeTo(4,4);
181  m_mom.ResizeTo(m_maxpart,6);
182  m_mom2.ResizeTo(m_maxpart,6);
183  m_np=-1;
184  m_np2=-1;
185  m_sanda_called=false;
186  m_fowo_called=false;
187  m_sumangles_called=false;
188  m_verbose=false;
189  m_boost=boost;
190  m_sph=-1;
191  m_apl=-1;
192  m_h10=-1;
193  m_h20=-1;
194  m_h30=-1;
195  m_h40=-1;
196  m_sqrts=0;
197  m_ht=0;
198  m_ht3=0;
199  m_et56=0;
200  m_njetsweighed=0;
201  m_et0=0;
202 }
TopologyWorker::m_dSphMomPower
double m_dSphMomPower
Definition: TopologyWorker.h:94
TopologyWorker::m_maxpart
static int m_maxpart
Definition: TopologyWorker.h:152
TopologyWorker::m_dDeltaThPower
double m_dDeltaThPower
Definition: TopologyWorker.h:97
virtual TopologyWorker::~TopologyWorker ( )
inlinevirtual

Definition at line 37 of file TopologyWorker.h.

37 {;}

Member Function Documentation

double TopologyWorker::CalcEta ( int  i)
inlineprivate

Definition at line 159 of file TopologyWorker.h.

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

159 {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
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.h:1583
double TopologyWorker::CalcEta2 ( int  i)
inlineprivate

Definition at line 160 of file TopologyWorker.h.

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

160 {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
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.h:1583
void TopologyWorker::CalcHTstuff ( )
private

Definition at line 1683 of file TopologyWorker.h.

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().

1683  {
1684  m_ht=0;
1685  m_ht3=0;
1686  m_et56=0;
1687  m_et0=0;
1688  double ptlead=0;
1689  double h=0;
1690  for(int i=0; i< m_np; i++){
1691  //cout << i << "/" << m_np << ":" << CalcPt(i) << endl;
1692  m_ht+=CalcPt(i);
1693  h+=m_mom(i,4);
1694  if(i>1)
1695  m_ht3+=CalcPt(i);
1696  if(i==5)
1697  m_et56=sqrt(CalcPt(i)*CalcPt(i-1));
1698  }
1699 
1700  for(int j=0; j< m_np2; j++){
1701  //cout << j << "/" << m_np2 << ":" << CalcPt2(j) << endl;
1702  if(ptlead<CalcPt2(j))
1703  ptlead=CalcPt2(j);
1704 
1705  }
1706  if(m_ht>0.0001){
1707  m_et0=ptlead/m_ht;
1708  //cout << "calculating ETO" << m_et0 << "=" << ptlead << endl;
1709  }
1710  if(h>0.00001)
1711  m_centrality=m_ht/h;
1712 }
i
int i
Definition: DBlmapReader.cc:9
h
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
Definition: L1TUtmAlgorithmRcd.h:4
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
TopologyWorker::CalcPt
double CalcPt(int i)
Definition: TopologyWorker.h:157
j
int j
Definition: DBlmapReader.cc:9
TopologyWorker::CalcPt2
double CalcPt2(int i)
Definition: TopologyWorker.h:158
double TopologyWorker::CalcPt ( int  i)
inlineprivate

Definition at line 157 of file TopologyWorker.h.

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

Referenced by CalcHTstuff(), and CalcWmul().

157 { 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:18
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 158 of file TopologyWorker.h.

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

Referenced by CalcHTstuff().

158 { 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:18
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void TopologyWorker::CalcSqrts ( )
private

Definition at line 1667 of file TopologyWorker.h.

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

Referenced by setPartList().

1667  {
1668  TLorentzVector event(0,0,0,0);
1669  TLorentzVector worker(0,0,0,0);
1670 
1671  for(int i=0; i< m_np; i++){
1672  double energy=m_mom(i,4);
1673  if(m_mom(i,4)<0.00001)
1674  energy=sqrt(pow(m_mom(i,1),2)+pow(m_mom(i,2),2)+pow(m_mom(i,3),2));
1675  // assume massless particle if only TVector3s are provided...
1676  worker.SetXYZT(m_mom(i,1),m_mom(i,2),m_mom(i,3),energy);
1677  event+=worker;
1678  }
1679  m_sqrts=event.M();
1680 }
i
int i
Definition: DBlmapReader.cc:9
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
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 1642 of file TopologyWorker.h.

References CalcPt(), m_njetsweighed, m_np, HLT_FULL_cff::njets, and mps_fire::result.

Referenced by setPartList().

1642  {
1643 
1644  Int_t njets = m_np;
1645  double result=0;
1646  for(Int_t ijet=0; ijet<njets-1; ijet++){
1647  double emin=55;
1648  double emax=55;
1649  if(CalcPt(ijet)<55)
1650  emax=CalcPt(ijet);
1651  if(CalcPt(ijet+1)<55)
1652  emin=CalcPt(ijet+1);
1653  result+=0.5 * (emax*emax-emin*emin)*(ijet+1);
1654  }
1655  double elo=15;
1656  if(CalcPt(njets-1)>elo){
1657  elo=CalcPt(njets-1);
1658  }
1659 
1660  result+=0.5 * (elo*elo-(15*15))*(njets);
1661  result/=((55*55)-100)/2.0;
1662  m_njetsweighed=result;
1663 }
HLT_FULL_cff.njets
tuple njets
Definition: HLT_FULL_cff.py:9578
mps_fire.result
tuple result
Definition: mps_fire.py:95
TopologyWorker::CalcPt
double CalcPt(int i)
Definition: TopologyWorker.h:157
void TopologyWorker::clear ( void  )
inline

Definition at line 39 of file TopologyWorker.h.

References m_np, and m_np2.

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

Definition at line 1445 of file TopologyWorker.h.

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

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

1445  {
1446 // 20020830 changed: from p/E to Et/Ettot and include H50 and H60
1447  m_fowo_called=true;
1448  // include fox wolframs
1449  float H10=-1;
1450  float H20=-1;
1451  float H30=-1;
1452  float H40=-1;
1453  float H50=-1;
1454  float H60=-1;
1455  if (1==1) {
1456  float P[1000][6],H0,HD,CTHE;
1457  int N,NP,I,J,I1,I2;
1458  H0=HD=0.;
1459  N=m_np;
1460  NP=0;
1461  for (I=1;I<N+1;I++){
1462  NP++; // start at one
1463  P[NP][1]=m_mom(I-1,1) ;
1464  P[NP][2]=m_mom(I-1,2) ;
1465  P[NP][3]=m_mom(I-1,3) ;
1466  P[NP][4]=m_mom(I-1,4) ;
1467  P[NP][5]=m_mom(I-1,5) ;
1468  }
1469 
1470  N=NP;
1471  NP=0;
1472 
1473  for (I=1;I<N+1;I++) {
1474  NP=NP+1;
1475  for (J=1;J<5;J++) {
1476  P[N+NP][J]=P[I][J];
1477  }
1478 // p/E
1479  P[N+NP][4]=sqrt(pow(P[I][1],2)+pow(P[I][2],2)+pow(P[I][3],2));
1480 // Et/Ettot
1481  P[N+NP][5]=sqrt(pow(P[I][1],2)+pow(P[I][2],2));
1482  H0=H0+P[N+NP][5];
1483  HD=HD+pow(P[N+NP][5],2);
1484  }
1485  H0=H0*H0;
1486 
1487 
1488 
1489  // Very low multiplicities (0 or 1) not considered.
1490  if (NP<2) {
1491  H10=-1.;
1492  H20=-1.;
1493  H30=-1.;
1494  H40=-1.;
1495  H50=-1.;
1496  H60=-1.;
1497  return;
1498  }
1499 
1500  // Calculate H1 - H4.
1501  H10=0.;
1502  H20=0.;
1503  H30=0.;
1504  H40=0.;
1505  H50=0.;
1506  H60=0.;
1507  for (I1=N+1;I1<N+NP+1;I1++) { //130
1508  for (I2=I1+1;I2<N+NP+1;I2++) { // 120
1509  CTHE=(P[I1][1]*P[I2][1]+P[I1][2]*P[I2][2]+P[I1][3]*P[I2][3])/
1510  (P[I1][4]*P[I2][4]);
1511  H10=H10+P[I1][5]*P[I2][5]*CTHE;
1512  double C2=(1.5*CTHE*CTHE-0.5);
1513  H20=H20+P[I1][5]*P[I2][5]*C2;
1514  double C3=(2.5*CTHE*CTHE*CTHE-1.5*CTHE);
1515  H30=H30+P[I1][5]*P[I2][5]*C3;
1516  // use recurrence
1517  double C4=(7*CTHE*C3 - 3*C2)/4.;
1518  double C5=(9*CTHE*C4 - 4*C3)/5.;
1519  double C6=(11*CTHE*C5 - 5*C4)/6.;
1520 // H40=H40+P[I1][5]*P[I2][5]*(4.375*pow(CTHE,4)-3.75*CTHE*CTHE+0.375);
1521 // H50=H50+P[I1][5]*P[I2][5]*
1522 // (63./8.*pow(CTHE,5)-70./8.*CTHE*CTHE*CTHE+15./8.*CTHE);
1523 // H60=H60+P[I1][5]*P[I2][5]*
1524 // (231/16.*pow(CTHE,6)-315./16.*CTHE*CTHE*CTHE*CTHE+105./16.*CTHE*CTHE-5./16.);
1525  H40=H40+P[I1][5]*P[I2][5]*C4;
1526  H50=H50+P[I1][5]*P[I2][5]*C5;
1527  H60=H60+P[I1][5]*P[I2][5]*C6;
1528  } // 120
1529  } // 130
1530 
1531  H10=(HD+2.*H10)/H0;
1532  H20=(HD+2.*H20)/H0;
1533  H30=(HD+2.*H30)/H0;
1534  H40=(HD+2.*H40)/H0;
1535  H50=(HD+2.*H50)/H0;
1536  H60=(HD+2.*H60)/H0;
1537  m_h10=H10;
1538  m_h20=H20;
1539  m_h30=H30;
1540  m_h40=H40;
1541  m_h50=H50;
1542  m_h60=H60;
1543  }
1544 
1545 }
P
#define P
Definition: PFJetBenchmark.cc:37
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
Exhume::I
const std::complex< double > I
Definition: I.h:8
N
#define N
Definition: blowfish.cc:9
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
double TopologyWorker::get_aplanarity ( )

Definition at line 1578 of file TopologyWorker.h.

References m_apl, m_sanda_called, and sanda().

1578  {
1579  if (!m_sanda_called) sanda();
1580  return m_apl;
1581 }
double TopologyWorker::get_centrality ( )
inline

Definition at line 78 of file TopologyWorker.h.

References m_centrality.

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

Definition at line 74 of file TopologyWorker.h.

References m_et0.

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

Definition at line 77 of file TopologyWorker.h.

References m_et56.

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

Definition at line 1547 of file TopologyWorker.h.

References fowo(), m_fowo_called, and m_h10.

1547  {
1548  if (!m_fowo_called) fowo();
1549  return m_h10;
1550 }
double TopologyWorker::get_h20 ( )

Definition at line 1551 of file TopologyWorker.h.

References fowo(), m_fowo_called, and m_h20.

1551  {
1552  if (!m_fowo_called) fowo();
1553  return m_h20;
1554 }
double TopologyWorker::get_h30 ( )

Definition at line 1555 of file TopologyWorker.h.

References fowo(), m_fowo_called, and m_h30.

1555  {
1556  if (!m_fowo_called) fowo();
1557  return m_h30;
1558 }
double TopologyWorker::get_h40 ( )

Definition at line 1559 of file TopologyWorker.h.

References fowo(), m_fowo_called, and m_h40.

1559  {
1560  if (!m_fowo_called) fowo();
1561  return m_h40;
1562 }
double TopologyWorker::get_h50 ( )

Definition at line 1564 of file TopologyWorker.h.

References fowo(), m_fowo_called, and m_h50.

1564  {
1565  if (!m_fowo_called) fowo();
1566  return m_h50;
1567 }
double TopologyWorker::get_h60 ( )

Definition at line 1568 of file TopologyWorker.h.

References fowo(), m_fowo_called, and m_h60.

1568  {
1569  if (!m_fowo_called) fowo();
1570  return m_h60;
1571 }
double TopologyWorker::get_ht ( )
inline

Definition at line 72 of file TopologyWorker.h.

References m_ht.

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

Definition at line 73 of file TopologyWorker.h.

References m_ht3.

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

Definition at line 76 of file TopologyWorker.h.

References m_njetsweighed.

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

Definition at line 1574 of file TopologyWorker.h.

References m_sanda_called, m_sph, and sanda().

1574  {
1575  if (!m_sanda_called) sanda();
1576  return m_sph;
1577 }
double TopologyWorker::get_sqrts ( )
inline

Definition at line 75 of file TopologyWorker.h.

References m_sqrts.

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

Definition at line 1583 of file TopologyWorker.h.

References dqm-mbProfile::log, AlCaHLTBitMon_ParallelJobs::p, pi, mathSSE::sqrt(), and funct::tan().

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

1583  {
1584 
1585  double pi=3.1415927;
1586 
1587  double p=sqrt(px*px+py*py+pz*pz);
1588  // get eta and phi
1589  double th=pi/2.;
1590  if (p!=0) {
1591  th=acos(pz/p); // Theta
1592  }
1593  float thg=th;
1594  if (th<=0) {
1595  thg = pi + th;
1596  }
1597  eta=-9.;
1598  if (tan( thg/2.0 )>0.000001) {
1599  eta = -log( tan( thg/2.0 ) );
1600  }
1601  phi = atan2(py,px);
1602  if(phi<=0) phi += 2.0*pi;
1603  return;
1604 }
pi
const Double_t pi
Definition: trackSplitPlot.h:36
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
Int_t TopologyWorker::getFast ( )

Definition at line 581 of file TopologyWorker.h.

References m_iFast.

582 {
583  return m_iFast;
584 }
double TopologyWorker::getThMomPower ( )

Definition at line 567 of file TopologyWorker.h.

References m_dDeltaThPower.

568 {
569  return 1.0 + m_dDeltaThPower;
570 }
TopologyWorker::m_dDeltaThPower
double m_dDeltaThPower
Definition: TopologyWorker.h:97
Int_t TopologyWorker::iPow ( int  man,
int  exp 
)
private

Definition at line 1630 of file TopologyWorker.h.

References create_public_lumi_plots::exp, and relval_2017::k.

Referenced by setPartList().

1631 {
1632  Int_t ans = 1;
1633  for( Int_t k = 0; k < exp; k++)
1634  {
1635  ans = ans*man;
1636  }
1637  return ans;
1638 }
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1088
relval_2017.k
list k
Definition: relval_2017.py:23
void TopologyWorker::ludbrb ( TMatrix *  mom,
double  the,
double  phi,
double  bx,
double  by,
double  bz 
)
private

Definition at line 652 of file TopologyWorker.h.

References beta, reco::dp, i, j, relval_2017::k, np, and makeMuonMisalignmentScenario::rot.

Referenced by setPartList().

658 {
659  // Ignore "zeroth" elements in rot,pr,dp.
660  // Trying to use physics-like notation.
661  TMatrix rot(4,4);
662  double pr[4];
663  double dp[5];
664  Int_t np = mom->GetNrows();
665  if ( the*the + phi*phi > 1.0E-20 )
666  {
667  rot(1,1) = TMath::Cos(the)*TMath::Cos(phi);
668  rot(1,2) = -TMath::Sin(phi);
669  rot(1,3) = TMath::Sin(the)*TMath::Cos(phi);
670  rot(2,1) = TMath::Cos(the)*TMath::Sin(phi);
671  rot(2,2) = TMath::Cos(phi);
672  rot(2,3) = TMath::Sin(the)*TMath::Sin(phi);
673  rot(3,1) = -TMath::Sin(the);
674  rot(3,2) = 0.0;
675  rot(3,3) = TMath::Cos(the);
676  for ( Int_t i = 0; i < np; i++ )
677  {
678  for ( Int_t j = 1; j < 4; j++ )
679  {
680  pr[j] = (*mom)(i,j);
681  (*mom)(i,j) = 0;
682  }
683  for ( Int_t jb = 1; jb < 4; jb++)
684  {
685  for ( Int_t k = 1; k < 4; k++)
686  {
687  (*mom)(i,jb) = (*mom)(i,jb) + rot(jb,k)*pr[k];
688  }
689  }
690  }
691  double beta = TMath::Sqrt( bx*bx + by*by + bz*bz );
692  if ( beta*beta > 1.0E-20 )
693  {
694  if ( beta > 0.99999999 )
695  {
696  //send message: boost too large, resetting to <~1.0.
697  bx = bx*(0.99999999/beta);
698  by = by*(0.99999999/beta);
699  bz = bz*(0.99999999/beta);
700  beta = 0.99999999;
701  }
702  double gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
703  for ( Int_t i = 0; i < np; i++ )
704  {
705  for ( Int_t j = 1; j < 5; j++ )
706  {
707  dp[j] = (*mom)(i,j);
708  }
709  double bp = bx*dp[1] + by*dp[2] + bz*dp[3];
710  double gbp = gamma*(gamma*bp/(1.0 + gamma) + dp[4]);
711  (*mom)(i,1) = dp[1] + gbp*bx;
712  (*mom)(i,2) = dp[2] + gbp*by;
713  (*mom)(i,3) = dp[3] + gbp*bz;
714  (*mom)(i,4) = gamma*(dp[4] + bp);
715  }
716  }
717  }
718  return;
719 }
beta
const double beta
Definition: CosmicMuonParameters.h:69
i
int i
Definition: DBlmapReader.cc:9
relval_2017.k
list k
Definition: relval_2017.py:23
np
int np
Definition: AMPTWrapper.h:33
j
int j
Definition: DBlmapReader.cc:9
reco::dp
auto dp
Definition: deltaR.h:22
TVector3 TopologyWorker::majorAxis ( )

Definition at line 595 of file TopologyWorker.h.

References m_dAxes, and m_MajorAxis.

Referenced by planes_thrust().

595  {
596  TVector3 m_MajorAxis(m_dAxes(2,1),m_dAxes(2,2),m_dAxes(2,3));
597  return m_MajorAxis;
598 }
TopologyWorker::m_MajorAxis
TVector3 m_MajorAxis
Definition: TopologyWorker.h:117
TVector3 TopologyWorker::minorAxis ( )

Definition at line 601 of file TopologyWorker.h.

References m_dAxes, and m_MinorAxis.

Referenced by planes_thrust().

601  {
602  TVector3 m_MinorAxis(m_dAxes(3,1),m_dAxes(3,2),m_dAxes(3,3));
603  return m_MinorAxis;
604 }
TopologyWorker::m_MinorAxis
TVector3 m_MinorAxis
Definition: TopologyWorker.h:118
double TopologyWorker::oblateness ( )

Definition at line 613 of file TopologyWorker.h.

References m_dOblateness.

613  {
614  return m_dOblateness;
615 }
void TopologyWorker::planes_sphe ( double &  pnorm,
double &  p2,
double &  p3 
)

Definition at line 879 of file TopologyWorker.h.

References a, create_public_lumi_plots::ax, b, EnergyCorrector::c, funct::cos(), alignCSCRings::e, eta, create_public_lumi_plots::exp, getetaphi(), Exhume::I, relval_2017::k, m_mom, m_np, Max(), Min(), N, P, p3, p4, AnalysisDataFormats_SUSYBSMObjects::pa, phi, pi, funct::pow(), PWT, SGN, funct::sin(), and mathSSE::sqrt().

879  {
880  //float SPH=-1; //FIXME: commented out since gcc461 complained that this variable is set but unused
881  //float APL=-1; //FIXME: commented out since gcc461 complained that this variable is set but unused
882 // C...Calculate matrix to be diagonalized.
883  float P[1000][6];
884  double SM[4][4],SV[4][4];
885  double PA,PWT,PS,SQ,SR,SP,SMAX,SGN;
886  int NP;
887  int J, J1, J2, I, N, JA, JB, J3, JC, JB1, JB2;
888  JA=JB=JC=0;
889  double RL;
890  float rlu,rlu1;
891  //
892  // --- get the input form GTRACK arrays
893  //
894  N=m_np;
895  NP=0;
896  for (I=1;I<N+1;I++){
897  NP++; // start at one
898  P[NP][1]=m_mom(I-1,1) ;
899  P[NP][2]=m_mom(I-1,2) ;
900  P[NP][3]=m_mom(I-1,3) ;
901  P[NP][4]=m_mom(I-1,4) ;
902  P[NP][5]=0;
903  }
904  //
905  //---
906  //
907  N=NP;
908 
909  for (J1=1;J1<4;J1++) {
910  for (J2=J1;J2<4;J2++) {
911  SM[J1][J2]=0.;
912  }
913  }
914  PS=0.;
915  for (I=1;I<N+1;I++) { // 140
916  PA=sqrt(pow(P[I][1],2)+pow(P[I][2],2)+pow(P[I][3],2));
917  double eta,phi;
918  getetaphi(P[I][1],P[I][2],P[I][3],eta,phi);
919  PWT=exp(-std::fabs(eta));
920  PWT=1.;
921  for (J1=1;J1<4;J1++) { // 130
922  for (J2=J1;J2<4;J2++) { // 120
923  SM[J1][J2]=SM[J1][J2]+PWT*P[I][J1]*P[I][J2];
924  }
925  } // 130
926  PS=PS+PWT*PA*PA;
927  } //140
928 // C...Very low multiplicities (0 or 1) not considered.
929  if(NP<2) {
930  //SPH=-1.; //FIXME: commented out since gcc461 complained that this variable is set but unused
931  //APL=-1.; //FIXME: commented out since gcc461 complained that this variable is set but unused
932  return;
933  }
934  for (J1=1;J1<4;J1++) { // 160
935  for (J2=J1;J2<4;J2++) { // 150
936  SM[J1][J2]=SM[J1][J2]/PS;
937  }
938  } // 160
939 // C...Find eigenvalues to matrix (third degree equation).
940  SQ=(SM[1][1]*SM[2][2]+SM[1][1]*SM[3][3]+SM[2][2]*SM[3][3]
941  -pow(SM[1][2],2)
942  -pow(SM[1][3],2)-pow(SM[2][3],2))/3.-1./9.;
943  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]*
944  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.;
945 
946  SP=TMath::Cos(TMath::ACos(TMath::Max(TMath::Min(SR/TMath::Sqrt(-pow(SQ,3)),1.),-1.))/3.);
947 
948  P[N+1][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Max(2.*SP,TMath::Sqrt(3.*(1.-SP*SP))-SP);
949  P[N+3][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Min(2.*SP,-TMath::Sqrt(3.*(1.-SP*SP))-SP);
950  P[N+2][4]=1.-P[N+1][4]-P[N+3][4];
951  if (P[N+2][4]> 1E-5) {
952 
953 // C...Find first and last eigenvector by solving equation system.
954  for (I=1;I<4;I=I+2) { // 240
955  for (J1=1;J1<4;J1++) { // 180
956  SV[J1][J1]=SM[J1][J1]-P[N+I][4];
957  for (J2=J1+1;J2<4;J2++) { // 170
958  SV[J1][J2]=SM[J1][J2];
959  SV[J2][J1]=SM[J1][J2];
960  }
961  } // 180
962  SMAX=0.;
963  for (J1=1;J1<4;J1++) { // 200
964  for (J2=1;J2<4;J2++) { // 190
965  if(std::fabs(SV[J1][J2])>SMAX) { // 190
966  JA=J1;
967  JB=J2;
968  SMAX=std::fabs(SV[J1][J2]);
969  }
970  } // 190
971  } // 200
972  SMAX=0.;
973  for (J3=JA+1;J3<JA+3;J3++) { // 220
974  J1=J3-3*((J3-1)/3);
975  RL=SV[J1][JB]/SV[JA][JB];
976  for (J2=1;J2<4;J2++) { // 210
977  SV[J1][J2]=SV[J1][J2]-RL*SV[JA][J2];
978  if (std::fabs(SV[J1][J2])>SMAX) { // GOTO 210
979  JC=J1;
980  SMAX=std::fabs(SV[J1][J2]);
981  }
982  } // 210
983  } // 220
984  JB1=JB+1-3*(JB/3);
985  JB2=JB+2-3*((JB+1)/3);
986  P[N+I][JB1]=-SV[JC][JB2];
987  P[N+I][JB2]=SV[JC][JB1];
988  P[N+I][JB]=-(SV[JA][JB1]*P[N+I][JB1]+SV[JA][JB2]*P[N+I][JB2])/
989  SV[JA][JB];
990  PA=TMath::Sqrt(pow(P[N+I][1],2)+pow(P[N+I][2],2)+pow(P[N+I][3],2));
991 // make a random number
992  float pa=P[N-1][I];
993  rlu=std::fabs(pa)-std::fabs(int(pa)*1.);
994  rlu1=std::fabs(pa*pa)-std::fabs(int(pa*pa)*1.);
995  SGN=pow((-1.),1.*int(rlu+0.5));
996  for (J=1;J<4;J++) { // 230
997  P[N+I][J]=SGN*P[N+I][J]/PA;
998  } // 230
999  } // 240
1000 
1001 // C...Middle axis orthogonal to other two. Fill other codes.
1002  SGN=pow((-1.),1.*int(rlu1+0.5));
1003  P[N+2][1]=SGN*(P[N+1][2]*P[N+3][3]-P[N+1][3]*P[N+3][2]);
1004  P[N+2][2]=SGN*(P[N+1][3]*P[N+3][1]-P[N+1][1]*P[N+3][3]);
1005  P[N+2][3]=SGN*(P[N+1][1]*P[N+3][2]-P[N+1][2]*P[N+3][1]);
1006 
1007 // C...Calculate sphericity and aplanarity. Select storing option.
1008  //SPH=1.5*(P[N+2][4]+P[N+3][4]); //FIXME: commented out since gcc461 complained that this variable is set but unused
1009  //APL=1.5*P[N+3][4]; //FIXME: commented out since gcc461 complained that this variable is set but unused
1010 
1011  } // check 1
1012 
1013  // so assume now we have Sphericity axis, which one give the minimal Pts
1014  double etstot[4];
1015  double eltot[4];
1016  double sum23=0;
1017  double sum22=0;
1018  double sum33=0;
1019  double pina[4];
1020  double ax[4], ay[4], az[4];
1021  for (int ia=1;ia<4;ia++) {
1022  etstot[ia]=0;
1023  eltot[ia]=0;
1024  pina[ia]=0;
1025  ax[ia]=P[N+ia][1];
1026  ay[ia]=P[N+ia][2];
1027  az[ia]=P[N+ia][3];
1028  ax[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1029  ay[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1030  az[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1031  }
1032 
1033 
1034  for (int k =0 ; k<m_np ; k++) {
1035  // double eta,phi;
1036  // getetaphi(m_mom(k,1),m_mom(k,2),m_mom(k,3),eta,phi);
1037  double W=1.0;
1038  for (int ia=1;ia<4;ia++) {
1039  double e=sqrt(m_mom(k,1)*m_mom(k,1) +
1040  m_mom(k,2)*m_mom(k,2) +
1041  m_mom(k,3)*m_mom(k,3));
1042  double el=ax[ia]*m_mom(k,1) + ay[ia]*m_mom(k,2) + az[ia]*m_mom(k,3);
1043  pina[ia]=el;
1044  double ets=(e*e-el*el);
1045  etstot[ia]+=ets*W;
1046  eltot[ia]+=el*el;
1047  }
1048  double a2=pina[2];
1049  double a3=pina[3];
1050  // double h=0.4;
1051  //a2=pina[2]*cos(h)+pina[3]*sin(h);
1052  //a3=pina[3]*cos(h)-pina[2]*sin(h);
1053  sum22+=a2*a2*W;
1054  sum23+=a2*a3*W;
1055  sum33+=a3*a3*W;
1056  }
1057 
1058 
1059 
1060  double pi=3.1415927;
1061  double phi=pi/2.0;
1062  double phip=pi/2.0;
1063  double a=sum23;
1064  double c=-a;
1065  double b=sum22-sum33;
1066  double disc=b*b-4*a*c;
1067  // cout << " disc " << disc << endl;
1068  if (disc>=0) {
1069  double x1=(sqrt(disc)-b)/2/a;
1070  double x2=(-sqrt(disc)-b)/2/a;
1071  phi=atan(x1);
1072  phip=atan(x2);
1073  if (phi<0) phi=2.*pi+phi;
1074  if (phip<0) phip=2.*pi+phip;
1075  }
1076  double p21=sum22*cos(phi)*cos(phi)+sum33*sin(phi)*sin(phi)+2*sum23*cos(phi)*sin(phi);
1077  double p31=sum22*sin(phi)*sin(phi)+sum33*cos(phi)*cos(phi)-2*sum23*cos(phi)*sin(phi);
1078 
1079  double p22=sum22*cos(phip)*cos(phip)+sum33*sin(phip)*sin(phip)+2*sum23*cos(phip)*sin(phip);
1080  double p32=sum22*sin(phip)*sin(phip)+sum33*cos(phip)*cos(phip)-2*sum23*cos(phip)*sin(phip);
1081 
1082 
1083  double d1=std::fabs(p31*p31 - p21*p21);
1084  double d2=std::fabs(p32*p32 - p22*p22);
1085  //cout << " eltot " << eltot[2] << " " << eltot[3] << endl;
1086  //cout << " phi " << phi << " " << phip << endl;
1087  //cout << " d " << d1 << " " << d2 << endl;
1088  p2=p21;
1089  p3=p31;
1090  if (d2>d1) {
1091  p2=p22;
1092  p3=p32;
1093  }
1094  pnorm=sqrt(eltot[1]);
1095  if (p2>p3) {
1096  p3=sqrt(p3);
1097  p2=sqrt(p2);
1098  }else {
1099  double p4=p3;
1100  p3=sqrt(p2);
1101  p2=sqrt(p4);
1102  }
1103  //cout << " sum2 sum3 " << sqrt(sum22) << " " << sqrt(sum33) << endl;
1104  //cout << " p2 p3 " << p2 << " " << p3 << endl;
1105  //double els=sqrt(eltot[2]*eltot[2]+eltot[3]*eltot[3]);
1106  // cout << " ets els " << (ettot[1]) << " " << els << endl;
1107 
1108  //m_sph=SPH; //FIXME: shouldn't the local variables be used to reset the class member variables accordingly?
1109  //m_apl=APL; //FIXME: shouldn't the local variables be used to reset the class member variables accordingly?
1110  return;
1111 } // end planes_sphe
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1088
relval_2017.k
list k
Definition: relval_2017.py:23
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
P
#define P
Definition: PFJetBenchmark.cc:37
AnalysisDataFormats_SUSYBSMObjects::pa
susybsm::HSCParticle pa
Definition: classes.h:8
Min
T Min(T a, T b)
Definition: MathUtil.h:39
create_public_lumi_plots.ax
tuple ax
Definition: create_public_lumi_plots.py:1092
pi
const Double_t pi
Definition: trackSplitPlot.h:36
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
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
PWT
double PWT[12]
Definition: herwig.h:108
p2
double p2[4]
Definition: TauolaWrapper.h:90
Max
T Max(T a, T b)
Definition: MathUtil.h:44
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.h:1583
SGN
#define SGN(A)
Definition: Simplify_begin.h:68
N
#define N
Definition: blowfish.cc:9
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
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
void TopologyWorker::planes_sphe_wei ( double &  pnorm,
double &  p2,
double &  p3 
)

Definition at line 1114 of file TopologyWorker.h.

References a, create_public_lumi_plots::ax, b, EnergyCorrector::c, funct::cos(), alignCSCRings::e, eta, create_public_lumi_plots::exp, getetaphi(), Exhume::I, relval_2017::k, m_mom, m_np, Max(), Min(), N, P, p3, p4, AnalysisDataFormats_SUSYBSMObjects::pa, phi, pi, funct::pow(), PWT, SGN, funct::sin(), and mathSSE::sqrt().

1114  {
1115  //float SPH=-1; //FIXME: commented out since gcc461 complained that this variable is set but unused
1116  //float APL=-1; //FIXME: commented out since gcc461 complained that this variable is set but unused
1117 // C...Calculate matrix to be diagonalized.
1118  float P[1000][6];
1119  double SM[4][4],SV[4][4];
1120  double PA,PWT,PS,SQ,SR,SP,SMAX,SGN;
1121  int NP;
1122  int J, J1, J2, I, N, JA, JB, J3, JC, JB1, JB2;
1123  JA=JB=JC=0;
1124  double RL;
1125  float rlu,rlu1;
1126  //
1127  // --- get the input form GTRACK arrays
1128  //
1129  N=m_np;
1130  NP=0;
1131  for (I=1;I<N+1;I++){
1132  NP++; // start at one
1133  P[NP][1]=m_mom(I-1,1) ;
1134  P[NP][2]=m_mom(I-1,2) ;
1135  P[NP][3]=m_mom(I-1,3) ;
1136  P[NP][4]=m_mom(I-1,4) ;
1137  P[NP][5]=0;
1138  }
1139  //
1140  //---
1141  //
1142  N=NP;
1143 
1144  for (J1=1;J1<4;J1++) {
1145  for (J2=J1;J2<4;J2++) {
1146  SM[J1][J2]=0.;
1147  }
1148  }
1149  PS=0.;
1150  for (I=1;I<N+1;I++) { // 140
1151  PA=sqrt(pow(P[I][1],2)+pow(P[I][2],2)+pow(P[I][3],2));
1152  // double eta,phi;
1153  // getetaphi(P[I][1],P[I][2],P[I][3],eta,phi);
1154  // PWT=exp(-std::fabs(eta));
1155  PWT=1.;
1156  for (J1=1;J1<4;J1++) { // 130
1157  for (J2=J1;J2<4;J2++) { // 120
1158  SM[J1][J2]=SM[J1][J2]+PWT*P[I][J1]*P[I][J2];
1159  }
1160  } // 130
1161  PS=PS+PWT*PA*PA;
1162  } //140
1163 // C...Very low multiplicities (0 or 1) not considered.
1164  if(NP<2) {
1165  //SPH=-1.; //FIXME: commented out since gcc461 complained that this variable is set but unused
1166  //APL=-1.; //FIXME: commented out since gcc461 complained that this variable is set but unused
1167  return;
1168  }
1169  for (J1=1;J1<4;J1++) { // 160
1170  for (J2=J1;J2<4;J2++) { // 150
1171  SM[J1][J2]=SM[J1][J2]/PS;
1172  }
1173  } // 160
1174 // C...Find eigenvalues to matrix (third degree equation).
1175  SQ=(SM[1][1]*SM[2][2]+SM[1][1]*SM[3][3]+SM[2][2]*SM[3][3]
1176  -pow(SM[1][2],2)
1177  -pow(SM[1][3],2)-pow(SM[2][3],2))/3.-1./9.;
1178  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]*
1179  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.;
1180 
1181  SP=TMath::Cos(TMath::ACos(TMath::Max(TMath::Min(SR/TMath::Sqrt(-pow(SQ,3)),1.),-1.))/3.);
1182 
1183  P[N+1][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Max(2.*SP,TMath::Sqrt(3.*(1.-SP*SP))-SP);
1184  P[N+3][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Min(2.*SP,-TMath::Sqrt(3.*(1.-SP*SP))-SP);
1185  P[N+2][4]=1.-P[N+1][4]-P[N+3][4];
1186  if (P[N+2][4]> 1E-5) {
1187 
1188 // C...Find first and last eigenvector by solving equation system.
1189  for (I=1;I<4;I=I+2) { // 240
1190  for (J1=1;J1<4;J1++) { // 180
1191  SV[J1][J1]=SM[J1][J1]-P[N+I][4];
1192  for (J2=J1+1;J2<4;J2++) { // 170
1193  SV[J1][J2]=SM[J1][J2];
1194  SV[J2][J1]=SM[J1][J2];
1195  }
1196  } // 180
1197  SMAX=0.;
1198  for (J1=1;J1<4;J1++) { // 200
1199  for (J2=1;J2<4;J2++) { // 190
1200  if(std::fabs(SV[J1][J2])>SMAX) { // 190
1201  JA=J1;
1202  JB=J2;
1203  SMAX=std::fabs(SV[J1][J2]);
1204  }
1205  } // 190
1206  } // 200
1207  SMAX=0.;
1208  for (J3=JA+1;J3<JA+3;J3++) { // 220
1209  J1=J3-3*((J3-1)/3);
1210  RL=SV[J1][JB]/SV[JA][JB];
1211  for (J2=1;J2<4;J2++) { // 210
1212  SV[J1][J2]=SV[J1][J2]-RL*SV[JA][J2];
1213  if (std::fabs(SV[J1][J2])>SMAX) { // GOTO 210
1214  JC=J1;
1215  SMAX=std::fabs(SV[J1][J2]);
1216  }
1217  } // 210
1218  } // 220
1219  JB1=JB+1-3*(JB/3);
1220  JB2=JB+2-3*((JB+1)/3);
1221  P[N+I][JB1]=-SV[JC][JB2];
1222  P[N+I][JB2]=SV[JC][JB1];
1223  P[N+I][JB]=-(SV[JA][JB1]*P[N+I][JB1]+SV[JA][JB2]*P[N+I][JB2])/
1224  SV[JA][JB];
1225  PA=TMath::Sqrt(pow(P[N+I][1],2)+pow(P[N+I][2],2)+pow(P[N+I][3],2));
1226 // make a random number
1227  float pa=P[N-1][I];
1228  rlu=std::fabs(pa)-std::fabs(int(pa)*1.);
1229  rlu1=std::fabs(pa*pa)-std::fabs(int(pa*pa)*1.);
1230  SGN=pow((-1.),1.*int(rlu+0.5));
1231  for (J=1;J<4;J++) { // 230
1232  P[N+I][J]=SGN*P[N+I][J]/PA;
1233  } // 230
1234  } // 240
1235 
1236 // C...Middle axis orthogonal to other two. Fill other codes.
1237  SGN=pow((-1.),1.*int(rlu1+0.5));
1238  P[N+2][1]=SGN*(P[N+1][2]*P[N+3][3]-P[N+1][3]*P[N+3][2]);
1239  P[N+2][2]=SGN*(P[N+1][3]*P[N+3][1]-P[N+1][1]*P[N+3][3]);
1240  P[N+2][3]=SGN*(P[N+1][1]*P[N+3][2]-P[N+1][2]*P[N+3][1]);
1241 
1242 // C...Calculate sphericity and aplanarity. Select storing option.
1243  //SPH=1.5*(P[N+2][4]+P[N+3][4]); //FIXME: commented out since gcc461 complained that this variable is set but unused
1244  //APL=1.5*P[N+3][4]; //FIXME: commented out since gcc461 complained that this variable is set but unused
1245 
1246  } // check 1
1247 
1248  // so assume now we have Sphericity axis, which one give the minimal Pts
1249  double etstot[4];
1250  double eltot[4];
1251  double sum23=0;
1252  double sum22=0;
1253  double sum33=0;
1254  double pina[4];
1255  double ax[4], ay[4], az[4];
1256  for (int ia=1;ia<4;ia++) {
1257  etstot[ia]=0;
1258  eltot[ia]=0;
1259  pina[ia]=0;
1260  ax[ia]=P[N+ia][1];
1261  ay[ia]=P[N+ia][2];
1262  az[ia]=P[N+ia][3];
1263  ax[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1264  ay[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1265  az[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1266  }
1267 
1268  for (int k =0 ; k<m_np ; k++) {
1269  double eta,phi;
1270  getetaphi(m_mom(k,1),m_mom(k,2),m_mom(k,3),eta,phi);
1271  double W=exp(-std::fabs(eta*1.0));
1272  for (int ia=1;ia<4;ia++) {
1273  double e=sqrt(m_mom(k,1)*m_mom(k,1) +
1274  m_mom(k,2)*m_mom(k,2) +
1275  m_mom(k,3)*m_mom(k,3));
1276  double el=ax[ia]*m_mom(k,1) + ay[ia]*m_mom(k,2) + az[ia]*m_mom(k,3);
1277  pina[ia]=el;
1278  double ets=(e*e-el*el);
1279  etstot[ia]+=ets*W;
1280  eltot[ia]+=el*el*W;
1281  }
1282  double a2=pina[2];
1283  double a3=pina[3];
1284  // double h=0.4;
1285  //a2=pina[2]*cos(h)+pina[3]*sin(h);
1286  //a3=pina[3]*cos(h)-pina[2]*sin(h);
1287  sum22+=a2*a2*W;
1288  sum23+=a2*a3*W;
1289  sum33+=a3*a3*W;
1290  }
1291 
1292 
1293  double pi=3.1415927;
1294  double phi=pi/2.0;
1295  double phip=pi/2.0;
1296  double a=sum23;
1297  double c=-a;
1298  double b=sum22-sum33;
1299  double disc=b*b-4*a*c;
1300  // cout << " disc " << disc << endl;
1301  if (disc>=0) {
1302  double x1=(sqrt(disc)-b)/2/a;
1303  double x2=(-sqrt(disc)-b)/2/a;
1304  phi=atan(x1);
1305  phip=atan(x2);
1306  if (phi<0) phi=2.*pi+phi;
1307  if (phip<0) phip=2.*pi+phip;
1308  }
1309  double p21=sum22*cos(phi)*cos(phi)+sum33*sin(phi)*sin(phi)+2*sum23*cos(phi)*sin(phi);
1310  double p31=sum22*sin(phi)*sin(phi)+sum33*cos(phi)*cos(phi)-2*sum23*cos(phi)*sin(phi);
1311 
1312  double p22=sum22*cos(phip)*cos(phip)+sum33*sin(phip)*sin(phip)+2*sum23*cos(phip)*sin(phip);
1313  double p32=sum22*sin(phip)*sin(phip)+sum33*cos(phip)*cos(phip)-2*sum23*cos(phip)*sin(phip);
1314 
1315 
1316  double d1=std::fabs(p31*p31 - p21*p21);
1317  double d2=std::fabs(p32*p32 - p22*p22);
1318  //cout << " eltot " << eltot[2] << " " << eltot[3] << endl;
1319  //cout << " phi " << phi << " " << phip << endl;
1320  //cout << " d " << d1 << " " << d2 << endl;
1321  p2=p21;
1322  p3=p31;
1323  if (d2>d1) {
1324  p2=p22;
1325  p3=p32;
1326  }
1327  pnorm=sqrt(eltot[1]);
1328  if (p2>p3) {
1329  p3=sqrt(p3);
1330  p2=sqrt(p2);
1331  }else {
1332  double p4=p3;
1333  p3=sqrt(p2);
1334  p2=sqrt(p4);
1335  }
1336  //cout << " sum2 sum3 " << sqrt(sum22) << " " << sqrt(sum33) << endl;
1337  //cout << " p2 p3 " << p2 << " " << p3 << endl;
1338  //double els=sqrt(eltot[2]*eltot[2]+eltot[3]*eltot[3]);
1339  // cout << " ets els " << (ettot[1]) << " " << els << endl;
1340 
1341  //m_sph=SPH; //FIXME: shouldn't the local variables be used to reset the class member variables accordingly?
1342  //m_apl=APL; //FIXME: shouldn't the local variables be used to reset the class member variables accordingly?
1343  return;
1344 } // end planes_sphe
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1088
relval_2017.k
list k
Definition: relval_2017.py:23
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
P
#define P
Definition: PFJetBenchmark.cc:37
AnalysisDataFormats_SUSYBSMObjects::pa
susybsm::HSCParticle pa
Definition: classes.h:8
Min
T Min(T a, T b)
Definition: MathUtil.h:39
create_public_lumi_plots.ax
tuple ax
Definition: create_public_lumi_plots.py:1092
pi
const Double_t pi
Definition: trackSplitPlot.h:36
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
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
PWT
double PWT[12]
Definition: herwig.h:108
p2
double p2[4]
Definition: TauolaWrapper.h:90
Max
T Max(T a, T b)
Definition: MathUtil.h:44
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.h:1583
SGN
#define SGN(A)
Definition: Simplify_begin.h:68
N
#define N
Definition: blowfish.cc:9
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
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
void TopologyWorker::planes_thrust ( double &  pnorm,
double &  p2,
double &  p3 
)

Definition at line 1348 of file TopologyWorker.h.

References a, create_public_lumi_plots::ax, b, EnergyCorrector::c, funct::cos(), alignCSCRings::e, relval_2017::k, m_mom, m_np, majorAxis(), minorAxis(), p3, p4, phi, pi, funct::sin(), mathSSE::sqrt(), and thrustAxis().

1348  {
1349  TVector3 thrustaxis=thrustAxis();
1350  TVector3 majoraxis=majorAxis();
1351  TVector3 minoraxis=minorAxis();
1352  // so assume now we have Sphericity axis, which one give the minimal Pts
1353  double etstot[4];
1354  double eltot[4];
1355  double sum23=0;
1356  double sum22=0;
1357  double sum33=0;
1358  double pina[4];
1359  double ax[4], ay[4], az[4];
1360  ax[1]=thrustaxis(0); ay[1]=thrustaxis(1); az[1]=thrustaxis(2);
1361  ax[2]=minoraxis(0); ay[2]=minoraxis(1); az[2]=minoraxis(2);
1362  ax[3]=majoraxis(0); ay[3]=majoraxis(1); az[3]=majoraxis(2);
1363  for (int ia=1;ia<4;ia++) {
1364  etstot[ia]=0;
1365  eltot[ia]=0;
1366  pina[ia]=0;
1367  ax[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1368  ay[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1369  az[ia]/=sqrt(ax[ia]*ax[ia]+ay[ia]*ay[ia]+az[ia]*az[ia]);
1370  }
1371 
1372  for (int k =0 ; k<m_np ; k++) {
1373  for (int ia=1;ia<4;ia++) {
1374  double e=sqrt(m_mom(k,1)*m_mom(k,1) +
1375  m_mom(k,2)*m_mom(k,2) +
1376  m_mom(k,3)*m_mom(k,3));
1377  double el=ax[ia]*m_mom(k,1) + ay[ia]*m_mom(k,2) + az[ia]*m_mom(k,3);
1378  pina[ia]=el;
1379  double ets=(e*e-el*el);
1380  etstot[ia]+=ets;
1381  eltot[ia]+=std::fabs(el);
1382  }
1383  double a2=pina[2];
1384  double a3=pina[3];
1385  // double h=0.4;
1386  //a2=pina[2]*cos(h)+pina[3]*sin(h);
1387  //a3=pina[3]*cos(h)-pina[2]*sin(h);
1388  sum22+=a2*a2;
1389  sum23+=a2*a3;
1390  sum33+=a3*a3;
1391  }
1392 
1393 
1394  double pi=3.1415927;
1395  double phi=pi/2.0;
1396  double phip=pi/2.0;
1397  double a=sum23;
1398  double c=-a;
1399  double b=sum22-sum33;
1400  double disc=b*b-4*a*c;
1401  // cout << " disc " << disc << endl;
1402  if (disc>=0) {
1403  double x1=(sqrt(disc)-b)/2/a;
1404  double x2=(-sqrt(disc)-b)/2/a;
1405  phi=atan(x1);
1406  phip=atan(x2);
1407  if (phi<0) phi=2.*pi+phi;
1408  if (phip<0) phip=2.*pi+phip;
1409  }
1410  double p21=sum22*cos(phi)*cos(phi)+sum33*sin(phi)*sin(phi)+2*sum23*cos(phi)*sin(phi);
1411  double p31=sum22*sin(phi)*sin(phi)+sum33*cos(phi)*cos(phi)-2*sum23*cos(phi)*sin(phi);
1412 
1413  double p22=sum22*cos(phip)*cos(phip)+sum33*sin(phip)*sin(phip)+2*sum23*cos(phip)*sin(phip);
1414  double p32=sum22*sin(phip)*sin(phip)+sum33*cos(phip)*cos(phip)-2*sum23*cos(phip)*sin(phip);
1415 
1416 
1417  double d1=std::fabs(p31*p31 - p21*p21);
1418  double d2=std::fabs(p32*p32 - p22*p22);
1419  //cout << " eltot " << eltot[2] << " " << eltot[3] << endl;
1420  //cout << " phi " << phi << " " << phip << endl;
1421  //cout << " d " << d1 << " " << d2 << endl;
1422  p2=p21;
1423  p3=p31;
1424  if (d2>d1) {
1425  p2=p22;
1426  p3=p32;
1427  }
1428  pnorm=sqrt(etstot[1]);
1429  if (p2>p3) {
1430  p3=sqrt(p3);
1431  p2=sqrt(p2);
1432  }else {
1433  double p4=p3;
1434  p3=sqrt(p2);
1435  p2=sqrt(p4);
1436  }
1437  //cout << " sum2 sum3 " << sqrt(sum22) << " " << sqrt(sum33) << endl;
1438  //cout << " p2 p3 " << p2 << " " << p3 << endl;
1439  //double els=sqrt(eltot[2]*eltot[2]+eltot[3]*eltot[3]);
1440  // cout << " ets els " << (ettot[1]) << " " << els << endl;
1441  return;
1442 } // end planes_thru
relval_2017.k
list k
Definition: relval_2017.py:23
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
TopologyWorker::thrustAxis
TVector3 thrustAxis()
Definition: TopologyWorker.h:589
create_public_lumi_plots.ax
tuple ax
Definition: create_public_lumi_plots.py:1092
TopologyWorker::minorAxis
TVector3 minorAxis()
Definition: TopologyWorker.h:601
pi
const Double_t pi
Definition: trackSplitPlot.h:36
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
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
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
TopologyWorker::majorAxis
TVector3 majorAxis()
Definition: TopologyWorker.h:595
p3
double p3[4]
Definition: TauolaWrapper.h:91
void TopologyWorker::sanda ( )
private

Definition at line 725 of file TopologyWorker.h.

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

Referenced by get_aplanarity(), and get_sphericity().

725  {
726  float SPH=-1;
727  float APL=-1;
728  m_sanda_called=true;
729  //=======================================================================
730  // By M.Vreeswijk, (core was fortran, stolen from somewhere)
731  // Purpose: to perform sphericity tensor analysis to give sphericity
732  // and aplanarity.
733  //
734  // Needs: Array (standard from root-tuples): GTRACK_px, py, pz
735  // The number of tracks in these arrays: GTRACK_ijet
736  // In addition: Array GTRACK_ijet contains a jet number 'ijet'
737  // (if you wish to change this, simply change code)
738  //
739  // Uses: TVector3 and TLorentzVector classes
740  //
741  // Output: Sphericity SPH and Aplanarity APL
742  //=======================================================================
743 // C...Calculate matrix to be diagonalized.
744  float P[1000][6];
745  double SM[4][4],SV[4][4];
746  double PA,PWT,PS,SQ,SR,SP,SMAX,SGN;
747  int NP;
748  int J, J1, J2, I, N, JA, JB, J3, JC, JB1, JB2;
749  JA=JB=JC=0;
750  double RL;
751  float rlu,rlu1;
752  //
753  // --- get the input form GTRACK arrays
754  //
755  N=m_np;
756  NP=0;
757  for (I=1;I<N+1;I++){
758  NP++; // start at one
759  P[NP][1]=m_mom(I-1,1) ;
760  P[NP][2]=m_mom(I-1,2) ;
761  P[NP][3]=m_mom(I-1,3) ;
762  P[NP][4]=m_mom(I-1,4) ;
763  P[NP][5]=0;
764  }
765  //
766  //---
767  //
768  N=NP;
769 
770  for (J1=1;J1<4;J1++) {
771  for (J2=J1;J2<4;J2++) {
772  SM[J1][J2]=0.;
773  }
774  }
775  PS=0.;
776  for (I=1;I<N+1;I++) { // 140
777  PA=sqrt(pow(P[I][1],2)+pow(P[I][2],2)+pow(P[I][3],2));
778  PWT=1.;
779  for (J1=1;J1<4;J1++) { // 130
780  for (J2=J1;J2<4;J2++) { // 120
781  SM[J1][J2]=SM[J1][J2]+PWT*P[I][J1]*P[I][J2];
782  }
783  } // 130
784  PS=PS+PWT*PA*PA;
785  } //140
786 // C...Very low multiplicities (0 or 1) not considered.
787  if(NP<2) {
788  SPH=-1.;
789  APL=-1.;
790  return;
791  }
792  for (J1=1;J1<4;J1++) { // 160
793  for (J2=J1;J2<4;J2++) { // 150
794  SM[J1][J2]=SM[J1][J2]/PS;
795  }
796  } // 160
797 // C...Find eigenvalues to matrix (third degree equation).
798  SQ=(SM[1][1]*SM[2][2]+SM[1][1]*SM[3][3]+SM[2][2]*SM[3][3]
799  -pow(SM[1][2],2)
800  -pow(SM[1][3],2)-pow(SM[2][3],2))/3.-1./9.;
801  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]*
802  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.;
803 
804  SP=TMath::Cos(TMath::ACos(TMath::Max(TMath::Min(SR/TMath::Sqrt(-pow(SQ,3)),1.),-1.))/3.);
805 
806  P[N+1][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Max(2.*SP,TMath::Sqrt(3.*(1.-SP*SP))-SP);
807  P[N+3][4]=1./3.+TMath::Sqrt(-SQ)*TMath::Min(2.*SP,-TMath::Sqrt(3.*(1.-SP*SP))-SP);
808  P[N+2][4]=1.-P[N+1][4]-P[N+3][4];
809  if (P[N+2][4]> 1E-5) {
810 
811 // C...Find first and last eigenvector by solving equation system.
812  for (I=1;I<4;I=I+2) { // 240
813  for (J1=1;J1<4;J1++) { // 180
814  SV[J1][J1]=SM[J1][J1]-P[N+I][4];
815  for (J2=J1+1;J2<4;J2++) { // 170
816  SV[J1][J2]=SM[J1][J2];
817  SV[J2][J1]=SM[J1][J2];
818  }
819  } // 180
820  SMAX=0.;
821  for (J1=1;J1<4;J1++) { // 200
822  for (J2=1;J2<4;J2++) { // 190
823  if(std::fabs(SV[J1][J2])>SMAX) { // 190
824  JA=J1;
825  JB=J2;
826  SMAX=std::fabs(SV[J1][J2]);
827  }
828  } // 190
829  } // 200
830  SMAX=0.;
831  for (J3=JA+1;J3<JA+3;J3++) { // 220
832  J1=J3-3*((J3-1)/3);
833  RL=SV[J1][JB]/SV[JA][JB];
834  for (J2=1;J2<4;J2++) { // 210
835  SV[J1][J2]=SV[J1][J2]-RL*SV[JA][J2];
836  if (std::fabs(SV[J1][J2])>SMAX) { // GOTO 210
837  JC=J1;
838  SMAX=std::fabs(SV[J1][J2]);
839  }
840  } // 210
841  } // 220
842  JB1=JB+1-3*(JB/3);
843  JB2=JB+2-3*((JB+1)/3);
844  P[N+I][JB1]=-SV[JC][JB2];
845  P[N+I][JB2]=SV[JC][JB1];
846  P[N+I][JB]=-(SV[JA][JB1]*P[N+I][JB1]+SV[JA][JB2]*P[N+I][JB2])/
847  SV[JA][JB];
848  PA=TMath::Sqrt(pow(P[N+I][1],2)+pow(P[N+I][2],2)+pow(P[N+I][3],2));
849 // make a random number
850  float pa=P[N-1][I];
851  rlu=std::fabs(pa)-std::fabs(int(pa)*1.);
852  rlu1=std::fabs(pa*pa)-std::fabs(int(pa*pa)*1.);
853  SGN=pow((-1.),1.*int(rlu+0.5));
854  for (J=1;J<4;J++) { // 230
855  P[N+I][J]=SGN*P[N+I][J]/PA;
856  } // 230
857  } // 240
858 
859 // C...Middle axis orthogonal to other two. Fill other codes.
860  SGN=pow((-1.),1.*int(rlu1+0.5));
861  P[N+2][1]=SGN*(P[N+1][2]*P[N+3][3]-P[N+1][3]*P[N+3][2]);
862  P[N+2][2]=SGN*(P[N+1][3]*P[N+3][1]-P[N+1][1]*P[N+3][3]);
863  P[N+2][3]=SGN*(P[N+1][1]*P[N+3][2]-P[N+1][2]*P[N+3][1]);
864 
865 // C...Calculate sphericity and aplanarity. Select storing option.
866  SPH=1.5*(P[N+2][4]+P[N+3][4]);
867  APL=1.5*P[N+3][4];
868 
869  } // check 1
870 
871  m_sph=SPH;
872  m_apl=APL;
873  return;
874 } // end sanda
P
#define P
Definition: PFJetBenchmark.cc:37
AnalysisDataFormats_SUSYBSMObjects::pa
susybsm::HSCParticle pa
Definition: classes.h:8
Min
T Min(T a, T b)
Definition: MathUtil.h:39
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
Exhume::I
const std::complex< double > I
Definition: I.h:8
PWT
double PWT[12]
Definition: herwig.h:108
Max
T Max(T a, T b)
Definition: MathUtil.h:44
SGN
#define SGN(A)
Definition: Simplify_begin.h:68
N
#define N
Definition: blowfish.cc:9
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 573 of file TopologyWorker.h.

References m_iFast.

574 {
575  // Error if sp not positive.
576  if ( nf > 3 ) m_iFast = nf;
577  return;
578 }
void TopologyWorker::setPartList ( TObjArray *  e1,
TObjArray *  e2 
)

Definition at line 211 of file TopologyWorker.h.

References Abs(), CalcHTstuff(), CalcSqrts(), CalcWmul(), ecal_dqm_sourceclient-live_cfg::cerr, gather_cfg::cout, HTMLExport::elem(), i, cuy::ib, iPow(), j, relval_2017::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, Min(), gen::n, np, np2, connectstrParser::o, AlCaHLTBitMon_ParallelJobs::p, phi, FWPFMaths::sgn(), sign(), groupFilesInBlocks::temp, tmax, ulAngle(), findQualityFiles::v, edmIntegrityCheck::work(), X, and Gflash::Z.

212 {
213  //To make this look like normal physics notation the
214  //zeroth element of each array, mom[i][0], will be ignored
215  //and operations will be on elements 1,2,3...
216  TMatrix mom(m_maxpart,6);
217  TMatrix mom2(m_maxpart,6);
218  double tmax = 0;
219  double phi = 0.;
220  double the = 0.;
221  double sgn;
222  TMatrix fast(m_iFast + 1,6);
223  TMatrix work(11,6);
224  double tdi[4] = {0.,0.,0.,0.};
225  double tds;
226  double tpr[4] = {0.,0.,0.,0.};
227  double thp;
228  double thps;
229  double pxtot=0;
230  double pytot=0;
231  double pztot=0;
232  double etot=0;
233 
234  TMatrix temp(3,5);
235  Int_t np = 0;
236  Int_t numElements = e1->GetEntries();
237  Int_t numElements2 = e2->GetEntries();
238 
239  // trying to sort...
240 
241 
242 
243  m_np=0;
244  for(Int_t elem=0;elem<numElements;elem++) {
245  if(m_verbose){
246  std::cout << "looping over array, element " << elem << std::endl;
247  }
248  TObject* o = (TObject*) e1->At(elem);
249  if(m_verbose){
250  std::cerr << "TopologyWorker:SetPartList(): adding jet " << elem << "." << std::endl;
251  }
252  if (np >= m_maxpart) {
253  printf("Too many particles input to TopologyWorker");
254  return;
255  }
256  if(m_verbose){
257  std::cout << "getting name of object..." << std::endl;
258  }
259  TString nam(o->IsA()->GetName());
260  if(m_verbose){
261  std::cout << "TopologyWorker::setPartList(): object is of type " << nam << std::endl;
262  }
263  if (nam.Contains("TVector3")) {
264  TVector3 v(((TVector3 *) o)->X(),
265  ((TVector3 *) o)->Y(),
266  ((TVector3 *) o)->Z());
267  mom(np,1) = v.X();
268  mom(np,2) = v.Y();
269  mom(np,3) = v.Z();
270  mom(np,4) = v.Mag();
271  }
272  else if (nam.Contains("TLorentzVector")) {
273  TVector3 v(((TLorentzVector *) o)->X(),
274  ((TLorentzVector *) o)->Y(),
275  ((TLorentzVector *) o)->Z());
276  mom(np,1) = v.X();
277  mom(np,2) = v.Y();
278  mom(np,3) = v.Z();
279  mom(np,4) = ((TLorentzVector *) o)->T();
280  }
281  else {
282  printf("TopologyWorker::setEvent input is not a TVector3 or a TLorentzVector\n");
283  continue;
284  }
285 
286 
287  if ( TMath::Abs( m_dDeltaThPower ) <= 0.001 ) {
288  mom(np,5) = 1.0;
289  }
290  else {
291  mom(np,5) = TMath::Power(mom(np,4),m_dDeltaThPower);
292  }
293  tmax = tmax + mom(np,4)*mom(np,5);
294  pxtot+=mom(np,1);
295  pytot+=mom(np,2);
296  pztot+=mom(np,3);
297  etot+=mom(np,4);
298  np++;
299  m_np=np;
300  }
301  Int_t np2=0;
302  // second jet array.... only use some values here.
303  for(Int_t elem=0;elem<numElements2;elem++) {
304  //cout << elem << endl;
305  TObject* o = (TObject*) e2->At(elem);
306  if (np2 >= m_maxpart) {
307  printf("Too many particles input to TopologyWorker");
308  return;
309  }
310 
311  TString nam(o->IsA()->GetName());
312  if (nam.Contains("TVector3")) {
313  TVector3 v(((TVector3 *) o)->X(),
314  ((TVector3 *) o)->Y(),
315  ((TVector3 *) o)->Z());
316  mom2(np2,1) = v.X();
317  mom2(np2,2) = v.Y();
318  mom2(np2,3) = v.Z();
319  mom2(np2,4) = v.Mag();
320  }
321  else if (nam.Contains("TLorentzVector")) {
322  TVector3 v(((TLorentzVector *) o)->X(),
323  ((TLorentzVector *) o)->Y(),
324  ((TLorentzVector *) o)->Z());
325  mom2(np2,1) = v.X();
326  mom2(np2,2) = v.Y();
327  mom2(np2,3) = v.Z();
328  mom2(np2,4) = ((TLorentzVector *) o)->T();
329  // cout << "mom2: " << mom2(np2,1) << ", " << mom2(np2,2)<<", " << mom2(np2,3)<<", " << mom2(np2,4)<< endl;
330  }
331  else {
332  printf("TopologyWorker::setEvent Second vector input is not a TVector3 or a TLorentzVector\n");
333  continue;
334  }
335  np2++;
336  m_np2=np2;
337  }
338  m_mom2=mom2;
339 
340  if (m_boost && m_np>1) {
341  printf("TopologyWorker::setEvent Only boosting first vector so watch out when you do this!!!");
342  TVector3 booze(-pxtot/etot,-pytot/etot,-pztot/etot);
343  TLorentzVector l1;
344  for (int k=0; k<m_np ; k++) {
345  l1.SetPx(mom(k,1));
346  l1.SetPy(mom(k,2));
347  l1.SetPz(mom(k,3));
348  l1.SetE(mom(k,4));
349  l1.Boost(booze);
350  mom(k,1)=l1.Px();
351  mom(k,2)=l1.Py();
352  mom(k,3)=l1.Pz();
353  mom(k,4)=l1.E();
354  }
355  }
356 
357  m_sanda_called=false;
358  m_fowo_called=false;
359  for (int ip=0;ip<m_maxpart;ip++) {
360  for (int id=0;id<6;id++) {
361  m_mom(ip,id)=mom(ip,id);
362  }
363  }
364 
365  if ( np < 2 ) {
366  m_dThrust[1] = -1.0;
367  m_dOblateness = -1.0;
368  return;
369  }
370  // for pass = 1: find thrust axis.
371  // for pass = 2: find major axis.
372  for ( Int_t pass=1; pass < 3; pass++) {
373  if ( pass == 2 ) {
374  phi = ulAngle(m_dAxes(1,1), m_dAxes(1,2));
375  ludbrb( &mom, 0, -phi, 0., 0., 0. );
376  for ( Int_t i = 0; i < 3; i++ ) {
377  for ( Int_t j = 1; j < 4; j++ ) {
378  temp(i,j) = m_dAxes(i+1,j);
379  }
380  temp(i,4) = 0;
381  }
382  ludbrb(&temp,0.,-phi,0.,0.,0.);
383  for ( Int_t ib = 0; ib < 3; ib++ ) {
384  for ( Int_t j = 1; j < 4; j++ ) {
385  m_dAxes(ib+1,j) = temp(ib,j);
386  }
387  }
388  the = ulAngle( m_dAxes(1,3), m_dAxes(1,1) );
389  ludbrb( &mom, -the, 0., 0., 0., 0. );
390  for ( Int_t ic = 0; ic < 3; ic++ ) {
391  for ( Int_t j = 1; j < 4; j++ ) {
392  temp(ic,j) = m_dAxes(ic+1,j);
393  }
394  temp(ic,4) = 0;
395  }
396  ludbrb(&temp,-the,0.,0.,0.,0.);
397  for ( Int_t id = 0; id < 3; id++ ) {
398  for ( Int_t j = 1; j < 4; j++ ) {
399  m_dAxes(id+1,j) = temp(id,j);
400  }
401  }
402  }
403  for ( Int_t ifas = 0; ifas < m_iFast + 1 ; ifas++ ) {
404  fast(ifas,4) = 0.;
405  }
406  // Find the m_iFast highest momentum particles and
407  // put the highest in fast[0], next in fast[1],....fast[m_iFast-1].
408  // fast[m_iFast] is just a workspace.
409  for ( Int_t i = 0; i < np; i++ ) {
410  if ( pass == 2 ) {
411  mom(i,4) = TMath::Sqrt( mom(i,1)*mom(i,1)
412  + mom(i,2)*mom(i,2) );
413  }
414  for ( Int_t ifas = m_iFast - 1; ifas > -1; ifas-- ) {
415  if ( mom(i,4) > fast(ifas,4) ) {
416  for ( Int_t j = 1; j < 6; j++ ) {
417  fast(ifas+1,j) = fast(ifas,j);
418  if ( ifas == 0 ) fast(ifas,j) = mom(i,j);
419  }
420  }
421  else {
422  for ( Int_t j = 1; j < 6; j++ ) {
423  fast(ifas+1,j) = mom(i,j);
424  }
425  break;
426  }
427  }
428  }
429  // Find axis with highest thrust (case 1)/ highest major (case 2).
430  for ( Int_t ie = 0; ie < work.GetNrows(); ie++ ) {
431  work(ie,4) = 0.;
432  }
433  Int_t p = TMath::Min( m_iFast, np ) - 1;
434  // Don't trust Math.pow to give right answer always.
435  // Want nc = 2**p.
436  Int_t nc = iPow(2,p);
437  for ( Int_t n = 0; n < nc; n++ ) {
438  for ( Int_t j = 1; j < 4; j++ ) {
439  tdi[j] = 0.;
440  }
441  for ( Int_t i = 0; i < TMath::Min(m_iFast,n); i++ ) {
442  sgn = fast(i,5);
443  if (iPow(2,(i+1))*((n+iPow(2,i))/iPow(2,(i+1))) >= i+1){
444  sgn = -sgn;
445  }
446  for ( Int_t j = 1; j < 5-pass; j++ ) {
447  tdi[j] = tdi[j] + sgn*fast(i,j);
448  }
449  }
450  tds = tdi[1]*tdi[1] + tdi[2]*tdi[2] + tdi[3]*tdi[3];
451  for ( Int_t iw = TMath::Min(n,9); iw > -1; iw--) {
452  if( tds > work(iw,4) ) {
453  for ( Int_t j = 1; j < 5; j++ ) {
454  work(iw+1,j) = work(iw,j);
455  if ( iw == 0 ) {
456  if ( j < 4 ) {
457  work(iw,j) = tdi[j];
458  }
459  else {
460  work(iw,j) = tds;
461  }
462  }
463  }
464  }
465  else {
466  for ( Int_t j = 1; j < 4; j++ ) {
467  work(iw+1,j) = tdi[j];
468  }
469  work(iw+1,4) = tds;
470  }
471  }
472  }
473  // Iterate direction of axis until stable maximum.
474  m_dThrust[pass] = 0;
475  thp = -99999.;
476  Int_t nagree = 0;
477  for ( Int_t iw = 0;
478  iw < TMath::Min(nc,10) && nagree < m_iGood; iw++ ){
479  thp = 0.;
480  thps = -99999.;
481  while ( thp > thps + m_dConv ) {
482  thps = thp;
483  for ( Int_t j = 1; j < 4; j++ ) {
484  if ( thp <= 1E-10 ) {
485  tdi[j] = work(iw,j);
486  }
487  else {
488  tdi[j] = tpr[j];
489  tpr[j] = 0;
490  }
491  }
492  for ( Int_t i = 0; i < np; i++ ) {
493  sgn = sign(mom(i,5),
494  tdi[1]*mom(i,1) +
495  tdi[2]*mom(i,2) +
496  tdi[3]*mom(i,3));
497  for ( Int_t j = 1; j < 5 - pass; j++ ){
498  tpr[j] = tpr[j] + sgn*mom(i,j);
499  }
500  }
501  thp = TMath::Sqrt(tpr[1]*tpr[1]
502  + tpr[2]*tpr[2]
503  + tpr[3]*tpr[3])/tmax;
504  }
505  // Save good axis. Try new initial axis until enough
506  // tries agree.
507  if ( thp < m_dThrust[pass] - m_dConv ) {
508  break;
509  }
510  if ( thp > m_dThrust[pass] + m_dConv ) {
511  nagree = 0;
512  sgn = iPow( -1, (Int_t)TMath::Nint(m_random.Rndm()) );
513  for ( Int_t j = 1; j < 4; j++ ) {
514  m_dAxes(pass,j) = sgn*tpr[j]/(tmax*thp);
515  }
516  m_dThrust[pass] = thp;
517  }
518  nagree = nagree + 1;
519  }
520  }
521  // Find minor axis and value by orthogonality.
522  sgn = iPow( -1, (Int_t)TMath::Nint(m_random.Rndm()));
523  m_dAxes(3,1) = -sgn*m_dAxes(2,2);
524  m_dAxes(3,2) = sgn*m_dAxes(2,1);
525  m_dAxes(3,3) = 0.;
526  thp = 0.;
527  for ( Int_t i = 0; i < np; i++ ) {
528  thp += mom(i,5)*TMath::Abs(m_dAxes(3,1)*mom(i,1) +
529  m_dAxes(3,2)*mom(i,2));
530  }
531  m_dThrust[3] = thp/tmax;
532  // Rotate back to original coordinate system.
533  for ( Int_t i6 = 0; i6 < 3; i6++ ) {
534  for ( Int_t j = 1; j < 4; j++ ) {
535  temp(i6,j) = m_dAxes(i6+1,j);
536  }
537  temp(i6,4) = 0;
538  }
539  ludbrb(&temp,the,phi,0.,0.,0.);
540  for ( Int_t i7 = 0; i7 < 3; i7++ ) {
541  for ( Int_t j = 1; j < 4; j++ ) {
542  m_dAxes(i7+1,j) = temp(i7,j);
543  }
544  }
545  m_dOblateness = m_dThrust[2] - m_dThrust[3];
546 
547  // more stuff:
548 
549  // calculate weighed jet multiplicity();
550  CalcWmul();
551  CalcHTstuff();
552  CalcSqrts();
553 
554 }
Gflash::Z
const double Z[kNumberCalorimeter]
Definition: GflashNameSpace.h:58
TopologyWorker::ulAngle
double ulAngle(double x, double y)
Definition: TopologyWorker.h:619
i
int i
Definition: DBlmapReader.cc:9
cuy.ib
int ib
Definition: cuy.py:660
FWPFMaths::sgn
float sgn(float val)
Definition: FWPFMaths.cc:9
relval_2017.k
list k
Definition: relval_2017.py:23
X
#define X(str)
Definition: MuonsGrabber.cc:48
Min
T Min(T a, T b)
Definition: MathUtil.h:39
ecal_dqm_sourceclient-live_cfg.cerr
tuple cerr
Definition: ecal_dqm_sourceclient-live_cfg.py:49
TopologyWorker::m_dThrust
double m_dThrust[4]
Definition: TopologyWorker.h:126
TopologyWorker::ludbrb
void ludbrb(TMatrix *mom, double the, double phi, double bx, double by, double bz)
Definition: TopologyWorker.h:652
np
int np
Definition: AMPTWrapper.h:33
Abs
T Abs(T a)
Definition: MathUtil.h:49
j
int j
Definition: DBlmapReader.cc:9
TopologyWorker::sign
double sign(double a, double b)
Definition: TopologyWorker.h:642
TopologyWorker::m_maxpart
static int m_maxpart
Definition: TopologyWorker.h:152
reco::e1
Float e1
Definition: deltaR.h:20
tmax
static const double tmax[3]
Definition: CastorTimeSlew.cc:7
np2
int np2
Definition: TauolaWrapper.h:100
reco::e2
Float e2
Definition: deltaR.h:21
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
TopologyWorker::m_dDeltaThPower
double m_dDeltaThPower
Definition: TopologyWorker.h:97
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:145
TopologyWorker::iPow
int iPow(int man, int exp)
Definition: TopologyWorker.h:1630
void TopologyWorker::setThMomPower ( double  tp)

Definition at line 559 of file TopologyWorker.h.

References m_dDeltaThPower.

560 {
561  // Error if sp not positive.
562  if ( tp > 0. ) m_dDeltaThPower = tp - 1.0;
563  return;
564 }
TopologyWorker::m_dDeltaThPower
double m_dDeltaThPower
Definition: TopologyWorker.h:97
void TopologyWorker::setVerbose ( bool  loud)
inline

Definition at line 42 of file TopologyWorker.h.

References m_verbose.

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

Definition at line 642 of file TopologyWorker.h.

References Abs().

Referenced by setPartList(), and ulAngle().

642  {
643  if ( b < 0 ) {
644  return -TMath::Abs(a);
645  }
646  else {
647  return TMath::Abs(a);
648  }
649 }
Abs
T Abs(T a)
Definition: MathUtil.h:49
b
double b
Definition: hdecay.h:120
a
double a
Definition: hdecay.h:121
void TopologyWorker::sumangles ( float &  sdeta,
float &  sdr 
)

Definition at line 1608 of file TopologyWorker.h.

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

1608  {
1609  double eta1,eta2,phi1,phi2,deta,dphi,dr;
1610  m_sumangles_called=true;
1611  sdeta=0;
1612  sdr=0;
1613  for (int k=0;k<m_np;k++){
1614  for (int kp=k;kp<m_np;kp++){
1615  getetaphi(m_mom(k,1) , m_mom(k,2), m_mom(k,3), eta1,phi1);
1616  getetaphi(m_mom(kp,1) , m_mom(kp,2), m_mom(kp,3), eta2,phi2);
1617  dphi=std::fabs(phi1-phi2);
1618  if (dphi>3.1415) dphi=2*3.1415-dphi;
1619  deta=std::fabs(eta1-eta2);
1620  dr=sqrt(dphi*dphi+deta*deta);
1621  sdeta+=deta;
1622  sdr+=dr;
1623  }
1624  }
1625  return;
1626 }
kp
int kp
Definition: CascadeWrapper.h:14
relval_2017.k
list k
Definition: relval_2017.py:23
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.h:1583
TVector3 TopologyWorker::thrust ( )

Definition at line 607 of file TopologyWorker.h.

References m_dThrust, and m_Thrust.

607  {
608  TVector3 m_Thrust(m_dThrust[1],m_dThrust[2],m_dThrust[3]);
609  return m_Thrust;
610 }
TopologyWorker::m_dThrust
double m_dThrust[4]
Definition: TopologyWorker.h:126
TVector3 TopologyWorker::thrustAxis ( )

Definition at line 589 of file TopologyWorker.h.

References m_dAxes, and m_ThrustAxis.

Referenced by planes_thrust().

589  {
590  TVector3 m_ThrustAxis(m_dAxes(1,1),m_dAxes(1,2),m_dAxes(1,3));
591  return m_ThrustAxis;
592 }
TopologyWorker::m_ThrustAxis
TVector3 m_ThrustAxis
Definition: TopologyWorker.h:116
double TopologyWorker::ulAngle ( double  x,
double  y 
)
private

Definition at line 619 of file TopologyWorker.h.

References Abs(), Pi, alignCSCRings::r, and sign().

Referenced by setPartList().

620 {
621  double ulangl = 0;
622  double r = TMath::Sqrt(x*x + y*y);
623  if ( r < 1.0E-20 ) {
624  return ulangl;
625  }
626  if ( TMath::Abs(x)/r < 0.8 ) {
627  ulangl = sign(TMath::ACos(x/r),y);
628  }
629  else {
630  ulangl = TMath::ASin(y/r);
631  if ( x < 0. && ulangl >= 0. ) {
632  ulangl = TMath::Pi() - ulangl;
633  }
634  else if ( x < 0. ) {
635  ulangl = -TMath::Pi() - ulangl;
636  }
637  }
638  return ulangl;
639 }
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
Abs
T Abs(T a)
Definition: MathUtil.h:49
TopologyWorker::sign
double sign(double a, double b)
Definition: TopologyWorker.h:642

Member Data Documentation

double TopologyWorker::m_apl
private

Definition at line 135 of file TopologyWorker.h.

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

bool TopologyWorker::m_boost
private

Definition at line 132 of file TopologyWorker.h.

Referenced by setPartList(), and TopologyWorker().

double TopologyWorker::m_centrality
private

Definition at line 148 of file TopologyWorker.h.

Referenced by CalcHTstuff(), and get_centrality().

TMatrix TopologyWorker::m_dAxes
private

Definition at line 110 of file TopologyWorker.h.

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

double TopologyWorker::m_dConv
private

Definition at line 103 of file TopologyWorker.h.

Referenced by setPartList().

double TopologyWorker::m_dDeltaThPower
private

Definition at line 97 of file TopologyWorker.h.

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

double TopologyWorker::m_dOblateness
private

Definition at line 127 of file TopologyWorker.h.

Referenced by oblateness(), and setPartList().

double TopologyWorker::m_dSphMomPower
private

Definition at line 94 of file TopologyWorker.h.

double TopologyWorker::m_dThrust[4]
private

Definition at line 126 of file TopologyWorker.h.

Referenced by setPartList(), and thrust().

double TopologyWorker::m_et0
private

Definition at line 144 of file TopologyWorker.h.

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

double TopologyWorker::m_et56
private

Definition at line 147 of file TopologyWorker.h.

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

bool TopologyWorker::m_fowo_called
private

Definition at line 131 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 136 of file TopologyWorker.h.

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

double TopologyWorker::m_h20
private

Definition at line 137 of file TopologyWorker.h.

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

double TopologyWorker::m_h30
private

Definition at line 138 of file TopologyWorker.h.

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

double TopologyWorker::m_h40
private

Definition at line 139 of file TopologyWorker.h.

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

double TopologyWorker::m_h50
private

Definition at line 140 of file TopologyWorker.h.

Referenced by fowo(), and get_h50().

double TopologyWorker::m_h60
private

Definition at line 141 of file TopologyWorker.h.

Referenced by fowo(), and get_h60().

double TopologyWorker::m_ht
private

Definition at line 142 of file TopologyWorker.h.

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

double TopologyWorker::m_ht3
private

Definition at line 143 of file TopologyWorker.h.

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

int TopologyWorker::m_iFast
private

Definition at line 100 of file TopologyWorker.h.

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

int TopologyWorker::m_iGood
private

Definition at line 106 of file TopologyWorker.h.

Referenced by setPartList().

TVector3 TopologyWorker::m_MajorAxis
private

Definition at line 117 of file TopologyWorker.h.

Referenced by majorAxis().

Int_t TopologyWorker::m_maxpart = 1000
staticprivate

Definition at line 152 of file TopologyWorker.h.

Referenced by setPartList(), and TopologyWorker().

TVector3 TopologyWorker::m_MinorAxis
private

Definition at line 118 of file TopologyWorker.h.

Referenced by minorAxis().

TMatrix TopologyWorker::m_mom
private

Definition at line 123 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 124 of file TopologyWorker.h.

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

double TopologyWorker::m_njetsweighed
private

Definition at line 146 of file TopologyWorker.h.

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

int TopologyWorker::m_np
private

Definition at line 128 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 129 of file TopologyWorker.h.

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

TRandom TopologyWorker::m_random
private

Definition at line 121 of file TopologyWorker.h.

Referenced by setPartList().

bool TopologyWorker::m_sanda_called
private

Definition at line 130 of file TopologyWorker.h.

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

double TopologyWorker::m_sph
private

Definition at line 134 of file TopologyWorker.h.

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

double TopologyWorker::m_sqrts
private

Definition at line 145 of file TopologyWorker.h.

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

bool TopologyWorker::m_sumangles_called
private

Definition at line 133 of file TopologyWorker.h.

Referenced by sumangles(), and TopologyWorker().

TVector3 TopologyWorker::m_Thrust
private

Definition at line 119 of file TopologyWorker.h.

Referenced by thrust().

TVector3 TopologyWorker::m_ThrustAxis
private

Definition at line 116 of file TopologyWorker.h.

Referenced by thrustAxis().

bool TopologyWorker::m_verbose
private

Definition at line 81 of file TopologyWorker.h.

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

Generated for CMSSW Reference Manual by doxygen 1.8.5