- Timestamp:
- Nov 17, 2008, 11:07:45 PM (16 years ago)
- Location:
- trunk/Utilities/Fastjet/plugins/CDFCones
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Utilities/Fastjet/plugins/CDFCones/interface/CalTower.hh
r11 r35 1 1 #ifndef _CAL_TOWER_HH_ 2 2 #define _CAL_TOWER_HH_ 3 4 // Modified by X. Rouby for integration in Delphes 3 5 4 6 #include <cmath> … … 8 10 #endif 9 11 10 const double TOWER_THETA[23] = { 3.000, 5.700, 8.400, 11.100, 13.800, 16.500, 19.200, 21.900, 24.600, 27.300, 30.000, 33.524, 11 36.822, 40.261, 43.614, 47.436, 51.790, 56.735, 62.310, 68.516, 75.297, 82.526, 90.000 }; 12 const double pi = acos(-1); 13 14 // CDF data : 23 towers. step=2.7° at the beginning, and after, it changes 15 //const unsigned int ntower = 23; 16 //const double TOWER_THETA[ntower] = { 3.000, 5.700, 8.400, 11.100, 13.800, 16.500, 19.200, 21.900, 24.600, 27.300, 30.000,// step=2.7° 17 // 33.524, 36.822, 40.261, 43.614, 47.436, 51.790, 56.735, 62.310, 68.516, 75.297, 82.526, 90.000 }; 18 19 // CMS HB: 0.087 x 0.087 = 0.087 x 5° 20 // CMS HF: 0.175 x 0.175 = 0.175 x 10° 21 // CMS data 22 const unsigned int ntower = 40; 23 //const double TOWER_ETA[ntower] = { 24 // 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131, 1.218, 1.305, 1.392, 1.479, 1.566, 25 // 1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.868, 2.950, 3.125, 3.300, 3.475, 3.650, 3.825, 4.000, 4.175, 26 // 4.350, 4.525, 4.700} 27 const double TOWER_THETA[ntower] = { 28 1.04221324479 , 1.24151823936 , 1.4789296564 , 1.76172836867 , 2.09858297042 , 2.49981191976 , 2.97769357253 , 3.54683172517 , 29 4.22458468545 , 5.03156565222 , 5.99222070913 , 6.50323511476 , 8.08260215014 , 9.38520669791 , 11.2029938559 , 13.0015703243 , 30 14.7731659079 , 16.5174067026 , 18.2267592981 , 19.9103474053 , 21.6789816306 , 23.5965023504 , 25.6731077732 , 27.9190341722 , 31 30.3443463769 , 32.9586699802 , 35.7708584993 , 38.7885912317 , 42.017901925 , 45.4626450521 , 49.1239157061 , 52.9994507062 , 32 57.0830516071 , 61.3640831827 , 65.8271109013 , 70.4517445074 , 75.2127486134 , 80.0804628398 , 85.02154355, 90 33 }; 12 34 13 35 class CalTower 14 36 { 15 37 public: 16 17 38 double Et,eta,phi; 18 39 int iEta,iPhi; 19 40 20 41 CalTower(): Et(0), eta(0), phi(0), iEta(-1), iPhi(-1) {} 21 CalTower( double Et0, double eta0,double phi0): Et(Et0), eta(eta0), phi(phi0)42 CalTower(const double Et0, const double eta0, const double phi0): Et(Et0), eta(eta0), phi(phi0) 22 43 { 23 if(fabs(eta) < -log(tan(TOWER_THETA[0]*M_PI/180/2))){24 if(eta <= 0){ 25 for( int i = 0; i < 22; i++)26 if(eta < -log(tan((180 - TOWER_THETA[i + 1])*M_PI/180/2))){44 if(fabs(eta) < pseudorapidity(TOWER_THETA[0])) { // |eta| < 3.64 45 if(eta <= 0){ // central particle, but negative eta 46 for(unsigned int i = 0; i < ntower-1; i++) 47 if(eta < pseudorapidity(180 - TOWER_THETA[i + 1])){ // eta < -3 ; -2.61 ; ... ; 0 27 48 iEta = 4 + i; 28 49 break; 29 50 } 30 } 31 else{ 32 for( int i = 0; i < 22; i++)33 if(-eta < -log(tan((180 - TOWER_THETA[i + 1])*M_PI/180/2))){51 } // central particle but negative eta 52 else{ // central particle but positive eta 53 for(unsigned int i = 0; i < ntower-1; i++) 54 if(-eta < pseudorapidity(180 - TOWER_THETA[i + 1])){ // eta > 3 ; 2.61 ; ... ; 0 34 55 iEta = 47 - i; 35 56 break; 36 57 } 37 } 38 if( iEta >= 8 && iEta < 14 || iEta >= 38 && iEta < 44)39 iPhi = int(phi/ 2/M_PI*48)%48;58 } // central particle but positive eta 59 if( (iEta >= 8 && iEta < 14) || (iEta >= 38 && iEta < 44) ) 60 iPhi = int(phi/(2*pi)*48)%48; 40 61 else 41 iPhi = int(phi/ 2/M_PI*24)%24;62 iPhi = int(phi/(2*pi)*24)%24; 42 63 } 43 else{ 64 else{ // particle beyond detector reach : 44 65 iEta = -1; 45 66 iPhi = -1; 46 67 } 47 68 } 48 CalTower( double Et0, double eta0, double phi0, int iEta0,int iPhi0): Et(Et0), eta(eta0), phi(phi0), iEta(iEta0), iPhi(iPhi0) {}69 CalTower(const double Et0, const double eta0, const double phi0, const int iEta0, const int iPhi0): Et(Et0), eta(eta0), phi(phi0), iEta(iEta0), iPhi(iPhi0) {} 49 70 CalTower(const CalTower& c): Et(c.Et), eta(c.eta), phi(c.phi), iEta(c.iEta), iPhi(c.iPhi) {} 50 bool isEqual( CalTowerc)71 bool isEqual(const CalTower& c) 51 72 { 52 73 return Et == c.Et && eta == c.eta && phi == c.phi && iEta == c.iEta && iPhi == c.iPhi; 53 74 } 75 76 // converts theta angle [degrees] into pseudorapidity 77 double pseudorapidity(const double theta) { 78 return -log(tan( (theta/2.) * (pi/180.) )); 79 } 80 54 81 }; 55 82 83 /***** 84 e.g. 85 iEta==4 si -3.642 < eta < -3.000 86 iEta==47 si 3.000 < eta < 3.642 87 iEta==-1 si eta < -3.642 ou si eta > 3.642 88 89 i theta pseudorapidity iEta 90 -------------------------------------- 91 0 3.0 -3.00008274261 4 92 1 5.7 -2.61134904017 5 93 2 8.4 -2.33129451223 6 94 3 11.1 -2.1118548488 7 95 4 13.8 -1.93106912741 8 96 5 16.5 -1.77704423284 9 97 6 19.2 -1.64260787639 10 98 7 21.9 -1.52309848418 11 99 8 24.6 -1.41531379676 12 100 9 27.3 -1.31695789692 13 101 10 30.0 -1.19999982651 14 102 11 33.524 -1.10000635195 15 103 12 36.822 -1.00361552815 16 104 13 40.261 -0.916007799535 17 105 14 43.614 -0.822472442947 18 106 15 47.436 -0.72264587494 19 107 16 51.79 -0.616250691646 20 108 17 56.735 -0.503273260393 21 109 18 62.31 -0.384075299436 22 110 19 68.516 -0.259479460739 23 111 20 75.297 -0.130817437415 24 112 21 82.526 1.11022302463e-16 25 113 114 115 i theta pseudorapidity iEta 116 -------------------------------------- 117 0 3.0 -3.00008274261 47 118 1 5.7 -2.61134904017 46 119 2 8.4 -2.33129451223 45 120 3 11.1 -2.1118548488 44 121 4 13.8 -1.93106912741 43 122 5 16.5 -1.77704423284 42 123 6 19.2 -1.64260787639 41 124 7 21.9 -1.52309848418 40 125 8 24.6 -1.41531379676 39 126 9 27.3 -1.31695789692 38 127 10 30.0 -1.19999982651 37 128 11 33.524 -1.10000635195 36 129 12 36.822 -1.00361552815 35 130 13 40.261 -0.916007799535 34 131 14 43.614 -0.822472442947 33 132 15 47.436 -0.72264587494 32 133 16 51.79 -0.616250691646 31 134 17 56.735 -0.503273260393 30 135 18 62.31 -0.384075299436 29 136 19 68.516 -0.259479460739 28 137 20 75.297 -0.130817437415 27 138 21 82.526 1.11022302463e-16 26 139 *****/ 140 56 141 #endif 142 -
trunk/Utilities/Fastjet/plugins/CDFCones/interface/Centroid.hh
r11 r35 11 11 { 12 12 public: 13 14 13 double Et,eta,phi; 15 14 … … 17 16 Centroid(double centroidEt, double centroidEta, double centroidPhi): Et(centroidEt), eta(centroidEta), phi(centroidPhi) {} 18 17 Centroid(const Centroid& c): Et(c.Et), eta(c.eta), phi(c.phi) {} 19 void add( Centroidc)18 void add(const Centroid& c) 20 19 { 21 20 double newEt = Et + c.Et; … … 33 32 Et = newEt; 34 33 } 35 void subtract( Centroidc)34 void subtract(const Centroid& c) 36 35 { 37 36 double newEt = Et - c.Et; … … 49 48 Et = newEt; 50 49 } 51 bool isEqual( Centroidc)50 bool isEqual(const Centroid& c) 52 51 { 53 52 return Et == c.Et && eta == c.eta && phi == c.phi; -
trunk/Utilities/Fastjet/plugins/CDFCones/interface/Cluster.hh
r11 r35 30 30 pt_tilde = 0.0; 31 31 } 32 void addTower( PhysicsTowerp)32 void addTower(const PhysicsTower& p) 33 33 { 34 34 towerList.push_back(p); … … 37 37 pt_tilde += p.fourVector.pt(); 38 38 } 39 void removeTower( PhysicsTowerp)39 void removeTower(const PhysicsTower& p) 40 40 { 41 41 for(std::vector<PhysicsTower>::iterator towerIter = towerList.begin(); towerIter != towerList.end(); towerIter++) -
trunk/Utilities/Fastjet/plugins/CDFCones/interface/LorentzVector.hh
r11 r35 30 30 return r; 31 31 } 32 void add( LorentzVectorv)32 void add(const LorentzVector& v) 33 33 { 34 34 px += v.px; … … 37 37 E += v.E; 38 38 } 39 void subtract( LorentzVectorv)39 void subtract(const LorentzVector& v) 40 40 { 41 41 px -= v.px; … … 44 44 E -= v.E; 45 45 } 46 bool isEqual( LorentzVectorv)46 bool isEqual(const LorentzVector& v) 47 47 { 48 48 return px == v.px && py == v.py && pz == v.pz && E == v.E; -
trunk/Utilities/Fastjet/plugins/CDFCones/interface/MidPointAlgorithm.hh
r11 r35 15 15 double _coneRadius; 16 16 double _coneAreaFraction; 17 int _maxPairSize;17 unsigned int _maxPairSize; 18 18 int _maxIterations; 19 19 double _overlapThreshold; … … 44 44 std::vector<Cluster>& stableCones, bool reduceConeSize); 45 45 void addClustersToPairs(std::vector<int>& testPair, std::vector< std::vector<int> >& pairs, 46 std::vector< std::vector<bool> >& distanceOK, int maxClustersInPair);46 std::vector< std::vector<bool> >& distanceOK, unsigned int maxClustersInPair); 47 47 void splitAndMerge(std::vector<Cluster>& stableCones, std::vector<Cluster>& jets); 48 48 void run(std::vector<PhysicsTower>& particles, std::vector<Cluster>& jets); -
trunk/Utilities/Fastjet/plugins/CDFCones/interface/PhysicsTower.hh
r11 r35 8 8 { 9 9 public: 10 11 10 LorentzVector fourVector; 12 11 CalTower calTower; 12 /// addition by GPS (2008-08-15) for tracking within fastjet 13 int fjindex; 13 14 14 15 PhysicsTower(): fourVector(LorentzVector()), calTower(CalTower()), fjindex(-1) {} 15 PhysicsTower( LorentzVector v, CalTowerc): fourVector(v), calTower(c), fjindex(-1) {}16 PhysicsTower(const LorentzVector& v, const CalTower& c): fourVector(v), calTower(c), fjindex(-1) {} 16 17 PhysicsTower(const PhysicsTower& p): fourVector(p.fourVector), calTower(p.calTower), fjindex(p.fjindex) {} 17 PhysicsTower( CalTowerc):18 PhysicsTower(const CalTower& c): 18 19 fourVector(LorentzVector(c.Et*cos(c.phi),c.Et*sin(c.phi),c.Et*sinh(c.eta),c.Et*cosh(c.eta))), calTower(c), fjindex(-1) {} 19 PhysicsTower( LorentzVectorv): fourVector(v), calTower(CalTower(v.Et(),v.eta(),v.phi())), fjindex(-1) {}20 PhysicsTower(const LorentzVector& v): fourVector(v), calTower(CalTower(v.Et(),v.eta(),v.phi())), fjindex(-1) {} 20 21 double Et() const {return calTower.Et;} 21 22 double eta() const {return calTower.eta;} … … 23 24 int iEta() const {return calTower.iEta;} 24 25 int iPhi() const {return calTower.iPhi;} 25 bool isEqual( PhysicsTowerp)26 bool isEqual(const PhysicsTower& p) 26 27 { 27 28 return fourVector.isEqual(p.fourVector) && calTower.isEqual(p.calTower); 28 29 } 29 /// addition by GPS (2008-08-15) for tracking within fastjet 30 int fjindex; 30 31 31 }; 32 32 -
trunk/Utilities/Fastjet/plugins/CDFCones/src/MidPointAlgorithm.cc
r11 r35 18 18 std::vector< std::vector<bool> > distanceOK; 19 19 distanceOK.resize(stableCones.size() - 1); 20 for( int nCluster1 = 1; nCluster1 < stableCones.size(); nCluster1++){20 for(unsigned int nCluster1 = 1; nCluster1 < stableCones.size(); nCluster1++){ 21 21 distanceOK[nCluster1 - 1].resize(nCluster1); 22 22 double cluster1Rapidity = stableCones[nCluster1].fourVector.y(); 23 23 double cluster1Phi = stableCones[nCluster1].fourVector.phi(); 24 for( int nCluster2 = 0; nCluster2 < nCluster1; nCluster2++){24 for(unsigned int nCluster2 = 0; nCluster2 < nCluster1; nCluster2++){ 25 25 double cluster2Rapidity = stableCones[nCluster2].fourVector.y(); 26 26 double cluster2Phi = stableCones[nCluster2].fourVector.phi(); … … 37 37 std::vector< std::vector<int> > pairs(0); 38 38 std::vector<int> testPair(0); 39 int maxClustersInPair = _maxPairSize;39 unsigned int maxClustersInPair = _maxPairSize; 40 40 if(!maxClustersInPair) 41 41 maxClustersInPair = stableCones.size(); … … 44 44 // Loop over all combinations. Calculate MidPoint. Make midPointClusters. 45 45 bool reduceConeSize = false; 46 for( int iPair = 0; iPair < pairs.size(); iPair++){46 for(unsigned int iPair = 0; iPair < pairs.size(); iPair++){ 47 47 // Calculate rapidity, phi and pT of MidPoint. 48 48 LorentzVector midPoint(0,0,0,0); 49 for( int iPairMember = 0; iPairMember < pairs[iPair].size(); iPairMember++)49 for(unsigned int iPairMember = 0; iPairMember < pairs[iPair].size(); iPairMember++) 50 50 midPoint.add(stableCones[pairs[iPair][iPairMember]].fourVector); 51 51 iterateCone(midPoint.y(),midPoint.phi(),midPoint.pt(),towers,stableCones,reduceConeSize); … … 116 116 117 117 void MidPointAlgorithm::addClustersToPairs(std::vector<int>& testPair, std::vector< std::vector<int> >& pairs, 118 std::vector< std::vector<bool> >& distanceOK, int maxClustersInPair)118 std::vector< std::vector<bool> >& distanceOK, unsigned int maxClustersInPair) 119 119 { 120 120 // Recursively adds clusters to pairs, triplets, ... whose mid-points are then calculated. … … 124 124 if(testPair.size()) 125 125 nextClusterStart = testPair.back() + 1; 126 for( int nextCluster = nextClusterStart; nextCluster <= distanceOK.size(); nextCluster++){126 for(unsigned int nextCluster = nextClusterStart; nextCluster <= distanceOK.size(); nextCluster++){ 127 127 // Is new SeedCone less than 2*_coneRadius apart from all clusters in testPair? 128 128 bool addCluster = true; 129 for( int iCluster = 0; iCluster < testPair.size() && addCluster; iCluster++)129 for(unsigned int iCluster = 0; iCluster < testPair.size() && addCluster; iCluster++) 130 130 if(!distanceOK[nextCluster - 1][testPair[iCluster]]) 131 131 addCluster = false;
Note:
See TracChangeset
for help on using the changeset viewer.