CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TopologyWorker.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: TopTools
4 // Class: TopologyWorker
5 //
14 #ifndef __TOPTOOLSTOPOLOGYWORKER__
15 #define __TOPTOOLSTOPOLOGYWORKER__
16 
17 #warning The TopologyWorker class is currently not supported anymore! There might be issues in its implementation.
18 #warning If you are still using it or planned to do so, please contact the admins of the corresponding CMSSW package.
19 #warning You can find their e-mail adresses in: cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/TopQuarkAnalysis/TopTools/.admin/
20 
21 #include "TF1.h"
22 #include "TMath.h"
23 #include "TClass.h"
24 #include "TString.h"
25 #include "TRandom.h"
26 #include "TMatrixD.h"
27 #include "TLorentzVector.h"
28 
29 #include <cmath>
30 #include <iostream>
31 
32 class TopologyWorker
33 {
34 public:
35  TopologyWorker(){;}
36  TopologyWorker(bool boost);
37  virtual ~TopologyWorker(){;}
38 
39  void clear(void){m_np=0;m_np2=0;return;}
40 
41  void setPartList(TObjArray* e1, TObjArray* e2);
42  void setVerbose(bool loud){m_verbose=loud; return;}
43 
44  void setThMomPower(double tp);
45  double getThMomPower();
46  void setFast(int nf);
47  int getFast();
48 
49  TVector3 thrustAxis();
50  TVector3 majorAxis();
51  TVector3 minorAxis();
52 
53  TVector3 thrust();
54  // thrust :: Corresponding thrust, major, and minor value.
55 
56  double oblateness();
57  double get_sphericity();
58  double get_aplanarity();
59  double get_h10();
60  double get_h20();
61  double get_h30();
62  double get_h40();
63  double get_h50();
64  double get_h60();
65 
66 
67  void planes_sphe(double& pnorm,double& p2, double& p3);
68  void planes_sphe_wei(double& pnorm,double& p2, double& p3);
69  void planes_thrust(double& pnorm,double& p2, double& p3);
70  void sumangles(float& sdeta, float& sdr);
71 
72  double get_ht() {return m_ht;}
73  double get_ht3() {return m_ht3;}
74  double get_et0() {return m_et0;}
75  double get_sqrts() {return m_sqrts;}
76  double get_njetW() {return m_njetsweighed;}
77  double get_et56() {return m_et56;}
78  double get_centrality() { return m_centrality;}
79 
80 private:
81  bool m_verbose;
82  void getetaphi(double px, double py, double pz, double& eta, double& phi);
83  double ulAngle(double x, double y);
84  double sign(double a, double b);
85  void ludbrb(TMatrix *mom,
86  double the,
87  double phi,
88  double bx,
89  double by,
90  double bz);
91 
92  int iPow(int man, int exp);
93 
94  double m_dSphMomPower;
95  // PARU(41): Power of momentum dependence in sphericity finder.
96 
97  double m_dDeltaThPower;
98  // PARU(42): Power of momentum dependence in thrust finder.
99 
100  int m_iFast;
101  // MSTU(44): # of initial fastest particles choosen to start search.
102 
103  double m_dConv;
104  // PARU(48): Convergence criteria for axis maximization.
105 
106  int m_iGood;
107  // MSTU(45): # different starting configurations that must
108  // converge before axis is accepted as correct.
109 
110  TMatrix m_dAxes;
111  // data: results
112  // m_dAxes[1] is the Thrust axis.
113  // m_dAxes[2] is the Major axis.
114  // m_dAxes[3] is the Minor axis.
115 
116  TVector3 m_ThrustAxis;
117  TVector3 m_MajorAxis;
118  TVector3 m_MinorAxis;
119  TVector3 m_Thrust;
120 
121  TRandom m_random;
122 
123  TMatrix m_mom;
124  TMatrix m_mom2;
125 
126  double m_dThrust[4];
127  double m_dOblateness;
128  int m_np;
129  int m_np2;
130  bool m_sanda_called;
131  bool m_fowo_called;
132  bool m_boost;
133  bool m_sumangles_called;
134  double m_sph;
135  double m_apl;
136  double m_h10;
137  double m_h20;
138  double m_h30;
139  double m_h40;
140  double m_h50;
141  double m_h60;
142  double m_ht;
143  double m_ht3;
144  double m_et0;
145  double m_sqrts;
146  double m_njetsweighed;
147  double m_et56;
148  double m_centrality;
149 
150  void sanda();
151  void fowo();
152  static int m_maxpart;
153 
154  void CalcWmul();
155  void CalcSqrts();
156  void CalcHTstuff();
157  double CalcPt(int i) { return sqrt(pow(m_mom(i,1),2)+pow(m_mom(i,2),2));}
158  double CalcPt2(int i) { return sqrt(pow(m_mom2(i,1),2)+pow(m_mom2(i,2),2));}
159  double CalcEta(int i) {double eta, phi;getetaphi(m_mom(i,1),m_mom(i,2),m_mom(i,3),eta,phi); return eta;}
160  double CalcEta2(int i) {double eta, phi; getetaphi(m_mom2(i,1),m_mom2(i,2),m_mom2(i,3),eta,phi); return eta;}
161 
162 };
163 
164 class LessThan {
165  public :
166  // retrieve tru info MC stuff
167  bool operator () (const TLorentzVector & tl1, const TLorentzVector &
168  tl2)
169  const {
170  return tl2.Pt() < tl1.Pt();
171  }
172 };
173 
174 Int_t TopologyWorker::m_maxpart = 1000;
175 
176 TopologyWorker::TopologyWorker(bool boost):
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 }
203 //______________________________________________________________
204 
205 
206 // Input the particle 3(4)-vector list
207 // e: 3-vector TVector3 ..(px,py,pz) or
208 // 4-vector TLorentzVector ..(px,py,pz,E)
209 // Even input the TLorentzVector, we don't use Energy
210 //
211 void TopologyWorker::setPartList(TObjArray* e1, TObjArray* e2)
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 }
555 //______________________________________________________________
556 
557 // Setting and getting parameters.
558 
559 void TopologyWorker::setThMomPower(double tp)
560 {
561  // Error if sp not positive.
562  if ( tp > 0. ) m_dDeltaThPower = tp - 1.0;
563  return;
564 }
565 //______________________________________________________________
566 
567 double TopologyWorker::getThMomPower()
568 {
569  return 1.0 + m_dDeltaThPower;
570 }
571 //______________________________________________________________
572 
573 void TopologyWorker::setFast(Int_t nf)
574 {
575  // Error if sp not positive.
576  if ( nf > 3 ) m_iFast = nf;
577  return;
578 }
579 //______________________________________________________________
580 
581 Int_t TopologyWorker::getFast()
582 {
583  return m_iFast;
584 }
585 //______________________________________________________________
586 
587 // Returning results
588 
589 TVector3 TopologyWorker::thrustAxis() {
590  TVector3 m_ThrustAxis(m_dAxes(1,1),m_dAxes(1,2),m_dAxes(1,3));
591  return m_ThrustAxis;
592 }
593 //______________________________________________________________
594 
595 TVector3 TopologyWorker::majorAxis() {
596  TVector3 m_MajorAxis(m_dAxes(2,1),m_dAxes(2,2),m_dAxes(2,3));
597  return m_MajorAxis;
598 }
599 //______________________________________________________________
600 
601 TVector3 TopologyWorker::minorAxis() {
602  TVector3 m_MinorAxis(m_dAxes(3,1),m_dAxes(3,2),m_dAxes(3,3));
603  return m_MinorAxis;
604 }
605 //______________________________________________________________
606 
607 TVector3 TopologyWorker::thrust() {
608  TVector3 m_Thrust(m_dThrust[1],m_dThrust[2],m_dThrust[3]);
609  return m_Thrust;
610 }
611 //______________________________________________________________
612 
613 double TopologyWorker::oblateness() {
614  return m_dOblateness;
615 }
616 //______________________________________________________________
617 
618 // utilities(from Jetset):
619 double TopologyWorker::ulAngle(double x, double y)
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 }
640 //______________________________________________________________
641 
642 double TopologyWorker::sign(double a, double b) {
643  if ( b < 0 ) {
644  return -TMath::Abs(a);
645  }
646  else {
647  return TMath::Abs(a);
648  }
649 }
650 //______________________________________________________________
651 
652 void TopologyWorker::ludbrb(TMatrix* mom,
653  double the,
654  double phi,
655  double bx,
656  double by,
657  double bz)
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 }
720 
721 
722 
723 // APLANARITY and SPHERICITY
724 
725 void TopologyWorker::sanda() {
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
875 
876 
877 
878 
879 void TopologyWorker::planes_sphe(double& pnorm, double& p2, double& p3) {
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
1112 
1113 
1114 void TopologyWorker::planes_sphe_wei(double& pnorm, double& p2, double& p3) {
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
1345 
1346 
1347 
1348 void TopologyWorker::planes_thrust(double& pnorm, double& p2, double& p3) {
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
1443 
1444 
1445 void TopologyWorker::fowo() {
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 }
1546 
1547 double TopologyWorker::get_h10() {
1548  if (!m_fowo_called) fowo();
1549  return m_h10;
1550 }
1551 double TopologyWorker::get_h20() {
1552  if (!m_fowo_called) fowo();
1553  return m_h20;
1554 }
1555 double TopologyWorker::get_h30() {
1556  if (!m_fowo_called) fowo();
1557  return m_h30;
1558 }
1559 double TopologyWorker::get_h40() {
1560  if (!m_fowo_called) fowo();
1561  return m_h40;
1562 }
1563 
1564 double TopologyWorker::get_h50() {
1565  if (!m_fowo_called) fowo();
1566  return m_h50;
1567 }
1568 double TopologyWorker::get_h60() {
1569  if (!m_fowo_called) fowo();
1570  return m_h60;
1571 }
1572 
1573 
1574 double TopologyWorker::get_sphericity() {
1575  if (!m_sanda_called) sanda();
1576  return m_sph;
1577 }
1578 double TopologyWorker::get_aplanarity() {
1579  if (!m_sanda_called) sanda();
1580  return m_apl;
1581 }
1582 
1583 void TopologyWorker::getetaphi(double px, double py, double pz, double& eta, double& phi) {
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 }
1605 
1606 
1607 
1608 void TopologyWorker::sumangles(float& sdeta, float& sdr) {
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 }
1627 
1628 //______________________________________________________________
1629 
1630 Int_t TopologyWorker::iPow(Int_t man, Int_t exp)
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 }
1639 
1640 // added by Freya:
1641 
1642 void TopologyWorker::CalcWmul(){
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 }
1664 
1665 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1666 
1667 void TopologyWorker::CalcSqrts(){
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 }
1681 
1682 //++++++++++++++++++++++++++++++++++++++
1683 void TopologyWorker::CalcHTstuff(){
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 }
1713 
1714 #endif
1715 
1716 
1717 
1718 
TopologyWorker::clear
void clear(void)
Definition: TopologyWorker.h:39
Gflash::Z
const double Z[kNumberCalorimeter]
Definition: GflashNameSpace.h:58
TopologyWorker::ulAngle
double ulAngle(double x, double y)
Definition: TopologyWorker.h:619
beta
const double beta
Definition: CosmicMuonParameters.h:69
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
i
int i
Definition: DBlmapReader.cc:9
TopologyWorker::get_njetW
double get_njetW()
Definition: TopologyWorker.h:76
create_public_lumi_plots.exp
tuple exp
Definition: create_public_lumi_plots.py:1087
cuy.ib
int ib
Definition: cuy.py:660
siStripFEDMonitor_P5_cff.Min
tuple Min
Definition: siStripFEDMonitor_P5_cff.py:55
kp
int kp
Definition: CascadeWrapper.h:14
FWPFMaths::sgn
float sgn(float val)
Definition: FWPFMaths.cc:9
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
TopologyWorker::planes_sphe_wei
void planes_sphe_wei(double &pnorm, double &p2, double &p3)
Definition: TopologyWorker.h:1114
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1102
X
#define X(str)
Definition: MuonsGrabber.cc:48
TopologyWorker::thrustAxis
TVector3 thrustAxis()
Definition: TopologyWorker.h:589
P
#define P
Definition: PFJetBenchmark.cc:35
AnalysisDataFormats_SUSYBSMObjects::pa
susybsm::HSCParticle pa
Definition: classes.h:8
eta
T eta() const
Definition: Basic3DVectorLD.h:179
LessThan::operator()
bool operator()(const TLorentzVector &tl1, const TLorentzVector &tl2) const
Definition: TopologyWorker.h:167
TopologyWorker::m_dSphMomPower
double m_dSphMomPower
Definition: TopologyWorker.h:94
TopologyWorker::m_MinorAxis
TVector3 m_MinorAxis
Definition: TopologyWorker.h:118
TopologyWorker::get_et56
double get_et56()
Definition: TopologyWorker.h:77
siStripFEDMonitor_P5_cff.Max
tuple Max
Definition: siStripFEDMonitor_P5_cff.py:56
TopologyWorker::oblateness
double oblateness()
Definition: TopologyWorker.h:613
create_public_lumi_plots.ax
tuple ax
Definition: create_public_lumi_plots.py:1091
TopologyWorker::minorAxis
TVector3 minorAxis()
Definition: TopologyWorker.h:601
pi
const Double_t pi
Definition: trackSplitPlot.h:35
TopologyWorker::m_dThrust
double m_dThrust[4]
Definition: TopologyWorker.h:126
TopologyWorker::CalcEta2
double CalcEta2(int i)
Definition: TopologyWorker.h:160
TopologyWorker::ludbrb
void ludbrb(TMatrix *mom, double the, double phi, double bx, double by, double bz)
Definition: TopologyWorker.h:652
TopologyWorker::get_sqrts
double get_sqrts()
Definition: TopologyWorker.h:75
np
int np
Definition: AMPTWrapper.h:33
TopologyWorker::CalcEta
double CalcEta(int i)
Definition: TopologyWorker.h:159
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
p4
double p4[4]
Definition: TauolaWrapper.h:92
TopologyWorker::CalcPt
double CalcPt(int i)
Definition: TopologyWorker.h:157
query.result
tuple result
Definition: query.py:137
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
j
int j
Definition: DBlmapReader.cc:9
TopologyWorker::planes_sphe
void planes_sphe(double &pnorm, double &p2, double &p3)
Definition: TopologyWorker.h:879
TopologyWorker::setVerbose
void setVerbose(bool loud)
Definition: TopologyWorker.h:42
Exhume::I
const std::complex< double > I
Definition: I.h:8
TopologyWorker::getThMomPower
double getThMomPower()
Definition: TopologyWorker.h:567
TopologyWorker::setThMomPower
void setThMomPower(double tp)
Definition: TopologyWorker.h:559
PWT
double PWT[12]
Definition: herwig.h:108
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
TopologyWorker::sign
double sign(double a, double b)
Definition: TopologyWorker.h:642
p2
double p2[4]
Definition: TauolaWrapper.h:90
TopologyWorker::planes_thrust
void planes_thrust(double &pnorm, double &p2, double &p3)
Definition: TopologyWorker.h:1348
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:81
TopologyWorker::get_sphericity
double get_sphericity()
Definition: TopologyWorker.h:1574
TopologyWorker::m_maxpart
static int m_maxpart
Definition: TopologyWorker.h:152
reco::e1
Float e1
Definition: deltaR.h:22
tmax
static const double tmax[3]
Definition: CastorTimeSlew.cc:7
np2
int np2
Definition: TauolaWrapper.h:100
TopologyWorker::getetaphi
void getetaphi(double px, double py, double pz, double &eta, double &phi)
Definition: TopologyWorker.h:1583
TopologyWorker::setPartList
void setPartList(TObjArray *e1, TObjArray *e2)
Definition: TopologyWorker.h:211
TopologyWorker::thrust
TVector3 thrust()
Definition: TopologyWorker.h:607
TopologyWorker::m_ThrustAxis
TVector3 m_ThrustAxis
Definition: TopologyWorker.h:116
SGN
#define SGN(A)
Definition: Simplify_begin.h:68
N
#define N
Definition: blowfish.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::m_MajorAxis
TVector3 m_MajorAxis
Definition: TopologyWorker.h:117
TopologyWorker::sumangles
void sumangles(float &sdeta, float &sdr)
Definition: TopologyWorker.h:1608
reco::dp
auto dp
Definition: deltaR.h:24
TopologyWorker::get_ht3
double get_ht3()
Definition: TopologyWorker.h:73
b
double b
Definition: hdecay.h:120
dtNoiseDBValidation_cfg.cerr
tuple cerr
Definition: dtNoiseDBValidation_cfg.py:22
reco::e2
Float e2
Definition: deltaR.h:23
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
TopologyWorker::get_aplanarity
double get_aplanarity()
Definition: TopologyWorker.h:1578
TopologyWorker::setFast
void setFast(int nf)
Definition: TopologyWorker.h:573
a
double a
Definition: hdecay.h:121
TopologyWorker::m_dDeltaThPower
double m_dDeltaThPower
Definition: TopologyWorker.h:97
TopologyWorker::get_et0
double get_et0()
Definition: TopologyWorker.h:74
TopologyWorker::majorAxis
TVector3 majorAxis()
Definition: TopologyWorker.h:595
TopologyWorker::get_centrality
double get_centrality()
Definition: TopologyWorker.h:78
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
x
Definition: DDAxes.h:10
TopologyWorker::~TopologyWorker
virtual ~TopologyWorker()
Definition: TopologyWorker.h:37
TopologyWorker::CalcPt2
double CalcPt2(int i)
Definition: TopologyWorker.h:158
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
TopologyWorker::iPow
int iPow(int man, int exp)
Definition: TopologyWorker.h:1630
p3
double p3[4]
Definition: TauolaWrapper.h:91
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5